MONTPEL: A multi-component Penman-Monteith
energy balance model

Rami Albasha, Loic Manceau, Heidi Webber, Michaël Chelle, Bruce Kimball,
Pierre Martre

► To cite this version:

Rami Albasha, Loic Manceau, Heidi Webber, Michaël Chelle, Bruce Kimball, et al. MONTPEL: A multi-component Penman-Monteith energy balance model. Agricultural and Forest Meteorology, 2024, 358, pp.110221. 10.1016/j.agrformet.2024.110221. hal-04737238

HAL Id: hal-04737238

https://hal.inrae.fr/hal-04737238v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# MONTPEL: A multi-component Penman-Monteith energy balance model

Rami Albasha<sup>a,b,*</sup>, Loïc Manceau<sup>a</sup>, Heidi Webber<sup>c,d</sup>, Michaël Chelle<sup>e</sup>, Bruce Kimball<sup>f</sup>, Pierre Martre<sup>a,*</sup>

## ARTICLE INFO

### Keywords:

Atmospheric stability correction
Canopy temperature
Crop model
Energy balance
Wheat
Python package

## ABSTRACT

Mechanistic modelling is gradually replacing empiricism in crop models, focusing on leaf-level physiological processes. This shift necessitates simulating crop surface temperature at infra-canopy sub-daily scales but many crop models still rely on empirical formulations for canopy temperature estimation, typically on a daily basis. We developed MONTPEL, a multi-component Penman-Monteith model that allows simulating the crop energy balance with flexible canopy representations (“BigLeaf” vs. “Layered”, “Lumped” vs. “Sunlit-Shaded”) and accounts for atmospheric stability conditions. We analyzed the model behavior, sensitivity and accuracy, using measurements from four wheat (Triticum aestivum L.) experiments conducted under varying pedoclimatic and water stress conditions. Measurements included hourly energy balance terms (total net radiation, soil heat flux, sensible and latent energy fluxes), hourly temperature of the canopy surface or of leaves at different depths inside the canopy, and sunlit and shaded leaf temperatures around solar noon at different dates. MONTPEL reproduced the measured energy balance terms with a root mean square error (RMSE) between 21 and 87 Wm⁻² and a coefficient of determination (R²) exceeding 0.65. The model’s accuracy in simulating canopy temperature, with RMSE ≤ 2.2 °C and R² ≥ 0.92, remained consistent regardless of measurement scale. Adjusting the aerodynamic resistance for atmospheric stability minimized simulated canopy temperature errors, notably in semi-arid conditions. Crop latent energy flux and temperature were most sensitive to the maximal stomatal conductance (g_s max) parameter. However, using a single g_s max value across the simulated experiments yielded satisfactory results, suggesting a weak sensitivity to the temporal and site-to-site variability of g_s max. Distinguishing sunlit from shaded canopy fractions systematically resulted in lower latent energy fluxes compared to “Lumped” canopy representation results. Analysis identified limitations in the multi-component approach, particularly an unrealistic uniform temperature shift across leaf layers when soil surface temperature changes.

## 1. Introduction

Crop models are becoming essential tools for predicting the impact of climate change on crop yield (Chenu et al., 2017; Guarin and Asseng, 2022; Rezaei et al., 2023). This progress has been achieved thanks to a shift in crop modelling paradigm where mechanistic description of physiological processes is gradually replacing empiricism (Boote et al., 2013; Wang et al., 2019). However, despite the increasing complexity of modelled formalisms, the crop temperature is still empirically-determined or supplanted by air temperature in many crop

models (e.g. APSIM, Keating et al. 2003; APSIM-X, Holzworth et al., 2014; STICS, Brisson et al. 2008, Beaudoin et al., 2023; EPIC, Williams et al., 1989). Such approximation is penalizing further modelling improvements for predicting the effects of water and heat stresses on crop development and growth (Rezaei et al., 2015; Maiorano et al., 2017, Wang et al., 2019; Guarin and Asseng, 2022).

Providing the crop modelling community with a flexible, stand-alone and physically-based energy balance model would contribute toward overcoming the current obstacles to crop modelling developments. Flexibility refers herein to the capability of the crop energy balance

<footnote><sup>a</sup> LEPSE, Univ. Montpellier, INRAE, Institut Agro Montpellier, Montpellier, France</footnote>

<footnote><sup>b</sup> itk society, Clapiers, France</footnote>

<footnote><sup>c</sup> Leibniz Centre for Agricultural Landscape Research, Müncheberg, Germany</footnote>

<footnote><sup>d</sup> Brandenburg University of Technology, Faculty of Environment and Natural Sciences, Cottbus, Germany</footnote>

<footnote><sup>e</sup> ECOSYS, INRAE, Université Paris-Saclay, Palaiseau, France</footnote>

<footnote><sup>f</sup> USDA, Agricultural Research Service, U.S. Arid-Land Agricultural Research Center, Maricopa, AZ, USA</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: rami.albasha@itk.fr (R. Albasha), pierre.martre@inrae.fr (P. Martre).</footnote>
model to adapt to different canopy representations in crop models. Indeed, the (uniform) crop canopy structure may be represented by one leaf layer (i.e. BigLeaf canopies) or as a succession of overlying leaf layers (i.e. Layered canopies). Concerning the canopy fluxes (e.g. latent and sensible heat fluxes), they may be calculated by distinguishing fluxes from sunlit and shaded leaf fractions (i.e. Sunlit-Shaded canopies), or by calculating fluxes from “average” leaves (i.e. Lumped canopies). Combined representations of canopy structure and fluxes schemes yield four possible crop categories: BigLeaf Lumped (e.g. CropSyst, Stöckle et al., 2003), BigLeaf Sunlit-Shaded (e.g. GECROS, Yin and van Laar, 2005), Layered Lumped (e.g. Sirius Quality, Martre et al., 2006) and Layered Sunlit-Shaded (e.g. Leuning and Wang, 1998).

This diversity in canopy representation requires calculating all energy balance processes in a bottom-up approach. That is to say, energy fluxes and resistances must be calculated on the leaf level then analytically scaled upwards, to the leaf layer or the canopy level and for lumped or sunlit and shaded leaf fractions. Raupach (1995), and later Alves et al. (1998) and Furon et al. (2007), thoroughly discussed the theoretical basis for a sound upscaling of surface and aerodynamic resistances through canopies. Raupach (1995) indicated that the most physically-sound procedure to scale elementary resistances (leaf) to higher levels (layer or canopy) is to perform the reciprocal of the sum of the conductances of all leaves, assuming that resistances act in parallel. That is to say, fluxes emanating from leaves are additive. This procedure is generally followed in the energy balance modules of Land Surface Models (e.g. Community Land Model CLM5; Lawrence et al., 2019) but has found few echoes in crop models (with the exception of the GECROS model, Yin and van Laar, 2005), probably due to the complexity of the resulting equations. Indeed, a non-negligible number of crop energy balance models employs theoretically unsound resistance scaling procedures, especially those which distinguish between sunlit and shaded leaf fractions of the canopy (e.g. Irmak et al. 2008; Zhang et al., 2011; Ding et al., 2014; Luo et al., 2018; Cuadra et al., 2021). For instance, stomatal conductance at the canopy-level for sunlit and shaded leaf fractions is commonly obtained by multiplying leaf-level average stomatal conductance by the surface area of the sunlit and shaded fractions, respectively (Irmak et al. 2008). This procedure was experimentally shown to fail in reproducing measured canopy surface conductance under dense canopies (Rochette et al., 1991) since it implicitly assumes no micro-climate gradients through canopy depth which contrasts with measurements (Leuning et al., 1995; Chelle, 2005; Chelle and Cellier, 2009).

There is an increasing body of evidence that accounting for the atmospheric stability conditions in crop energy balance models improves the crop temperature predictions, regardless of the canopy representation (Webber et al., 2017, Webber et al., 2018). The atmospheric stability describes the ability of the atmospheric boundary layer to facilitate vertical motions, which affects the crop aerodynamic resistance. When an air parcel with a certain temperature meets a hotter surface, the temperature of the air parcel increases. Since air parcels show adiabatic behavior (weakly exchange temperature with surrounding parcels), their temperature will be higher than the surrounding air and hence their density will be less than the outside air. Consequently, the air parcel will likely rise upwards and will continue rising since at each height its temperature will be higher (and density lesser) than the outside, resulting in an accelerated upwards motion. Under such conditions, the atmosphere is described as “unstable”. Conversely, when the crop or the ground surface is cooler than the overlying air, air parcels in the direct vicinity of the crop will reduce their temperature and consequently increase their density. Consequently, these parcels will tend to go back downwards whenever any force would push them upwards; the atmosphere is then “stable”. Under the special case where the temperature of air parcels is equal to the outside air, at each height above the canopy or ground surface, the atmosphere is “neutral”.

Correction for stability conditions is generally omitted in crop models (e.g. APSIM, Keating et al. 2003; APSIM-X, Holzworth et al.,

2014; GECROS, Yin and van Laar, 2005; STICS, Brisson et al. 2008, Beaudoin et al., 2023). The few early studies which examined the effect of the atmospheric stability correction on the performance of crop models, reported no significant improvements (e.g. van Zyl and Jager, 1987; Kjelgaard et al., 1996). Such conclusions were probably due to the fact that most crop models simulate with a daily time step whereby stability correction may not matter when averaged over days due to the diurnal cycle of stability conditions. During day hours under clear sky, the boundary layer of the atmosphere usually undergoes unstable conditions, and the atmospheric correction reduces the aerodynamic resistance. In contrast, during the night, the stability correction increases the aerodynamic resistance as the atmospheric condition is often stable. On a daily basis, stability correction may therefore have a low impact on the average daily aerodynamic resistance, and hence on daily crop temperature and transpiration. Shuttleworth (2007) recalled that this low impact is an implicit assumption of the daily time step Penman-Monteith equation, which probably explains why stability correction is found in crop models only when the energy balance simulations of the latter simulate on sub-daily time steps (e.g. Grant et al., 2011; Kimball et al., 2015; Webber et al., 2016; Cuadra et al., 2021).

We developed in this study a multi-component energy balance model (MONTPEL) based on the generic framework proposed by Lhomme et al. (2013). The proposed model allows simulating surface temperature, latent and sensible heat fluxes with different canopy representations in order to facilitate its coupling with crop growth models (i.e. those distinguishing sunlit from shaded leaves or considering only Lumped leaves, considering a single BigLeaf layer or multiple layers). To the best of our knowledge, this is the first example of a flexible energy balance modelling approach. The energy balance models proposed by Leuning et al. (1995) and Wang and Leuning (1998), split energy fluxes between sunlit and shaded leaf fractions but do not consider soil resistances in the overall resistance network. The ecosys model (Grant et al., 2011) distinguishes sunlit from shaded leaf fractions but it only runs on BigLeaf basis. Kimball et al. (2015) proposed a similar model that includes stability correction functions but runs only on BigLeaf canopies. The energy balance model developed by Webber et al. (2016) and implemented in the Lintul5 crop growth model is a BigLeaf Lumped (no Sunlit-Shaded leaf distinction) model. The GECROS model (Yin and van Laar, 2005) is a Sunlit-Shaded based model that follows Raupach (1995) scaling procedure but does not consider correction for stability conditions nor distinguish between different layers in simulation of processes. In addition, GECROS has no ‘link’ between crop components whereby the temperature of one component influences that of other components. Finally, the recent energy balance developments in the DSSAT-CSM-CROPGRO model (Cuadra et al., 2021) does not allow for multi-layer simulations and the scaling procedure of leaf-to-canopy resistances do not follow that proposed by Raupach (1995) as discussed above.

We evaluated the performance of MONTPEL against a set of four independent field datasets collected in three contrasting pedo-climatic environments. This evaluation focused only on the energy balance terms and the crop temperature, taking the crop variables (e.g. crop height, leaf area index, etc.) as inputs. Such independent evaluation should help to ensure robust performance as no error compensation between the different simulated processes can take place. This independent evaluation allowed analyzing the effect of canopy representation (BigLeaf vs Layered and Lumped vs Sunlit-Shaded) on predicted crop temperature, latent and sensible heat fluxes. Finally, the added value of including the effect of aerodynamic resistance correction for atmospheric stability conditions on the model performance and theoretical issues related to the degree of complexity that is needed for an accurate prediction of canopy temperature are discussed.
## 2. Materials and methods

### 2.1. Model description

MONTPEL relies on a discretized form of the Penman-Monteith model for multi-component canopies proposed by Lhomme et al. (2013). The crop parts considered in the surface energy balance (i.e. canopy and soil) are separated into components that are distinguishable by their physical properties (e.g. vegetation) or state variables (e.g. absorbed irradiance). A component can represent a Lumped leaf layer, a sunlit or shaded leaf layer fraction, the whole canopy, the sunlit or shaded fraction of the whole canopy, or the soil surface, depending on the aims of a particular analysis. Lhomme et al. (2013)'s formulation ensures that the energy balance of one component is linked to the energy balance of all the other components. This means that the temperature of one crop component influences the temperature of all the other crop components.

The multi-component approach assumes no resistance to vertical diffusion within the canopy. Resistances within the canopy are the bulk

surface resistances ($r_{s, i}$, h m⁻¹) and the bulk boundary layer resistances ($r_{a, i}$, h m⁻¹), both of which being in series while their equivalent resistance (i.e. $r_{s, i} + r_{a, i}$) is in parallel (Fig. 1).

The soil component has its own formula for surface resistance ($r_{s, soil}$, h m⁻¹) and bulk aerodynamic resistance ($r_{a, soil}$, h m⁻¹), which is the sum of the boundary layer resistance and the aerodynamic resistance between the substrate and the source height ($z_m$, m) determined from the integral of the reciprocal of eddy diffusivity following Choudhury and Monteith (1988). The total sensible and latent energy fluxes of the crop (respectively $\lambda E^n$ and $H^n$) are each encountered by the aerodynamic resistance ($r_{a, 0}$, h m⁻¹), which is assumed to be the same for heat and vapor transfer.

For the sake of brevity, we only describe in this section the upscaling procedure from leaf to layer (or canopy) processes and the atmospheric stability correction equations. We provide a detailed description of all model equations in Appendix A and model variables, parameters, and constants (including their nominal values) in Supplementary Tables S1 to S3.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** (a) Schema of the network of sensible and latent energy flux resistances for a crop separated into *n* leaf and soil components (adapted from Lhomme et al., 2013). $\lambda E^n$ (W m⁻² ground) is the total latent energy flux of the crop, $\lambda E_i$ and $\lambda E_{soil}$ (W m⁻² ground) are the latent energy fluxes of the canopy component *i* and the soil component, respectively, $H^n$ (W m⁻² ground) is the total sensible energy flux of the crop, $H_i$ and $H_{soil}$ (W m⁻² ground) are the sensible energy fluxes of the canopy component *i* and the soil component, respectively, $e_a$ (kPa) is vapor pressure at reference height ($z_r$, m), $e_m$ and $e_m^*$ (kPa) are respectively the actual and saturated vapor pressure at the canopy source height ($z_m$, m), $r_{a, 0}$ (h m⁻¹) is the aerodynamic resistance, $r_{a, i}$ (h m⁻¹) is the boundary layer resistance of the canopy component *i*, $r_{a, soil}$ (h m⁻¹) is the aerodynamic resistance between the ground surface and the canopy source height, $r_{s, i}$ and $r_{s, soil}$ (h m⁻¹) are the surface resistance of the canopy component *i* and the soil component, respectively, $T_a$ (K) is air temperature at reference height, $T_m$ (K) is the aerodynamic temperature at the canopy source height, $T_{s, i}$ and $T_{s, soil}$ (K) are the surface temperature of the canopy component *i* and the soil component, respectively, (b) possible representation of the canopy as BigLeaf Lumped, BigLeaf Sunlit-Shaded, Layered Lumped and Layered Sunlit-Shaded; the resistances network in (a) applies to each representation in (b).</figcaption></figure>
## 2.1.1. Leaf surface and boundary layer resistances upscaling

For both leaf surface and boundary layer resistances, the component-level resistance ($r_i$, m m⁻¹) is derived from the reciprocal of the corresponding leaf-level conductance (g, m h⁻¹) by assuming that fluxes emanating from leaves are additive (Raupach, 1995):

$$ \mathbf{r}_i = \left( \int_{L_t}^{L_u} \int_{i}^{i} \mathbf{g} \, \varphi \, d\mathbf{L} \right)^{-1} \qquad (1) $$

where $L_u$, i and $L_t$, i ($m^2_{leaf}$, $m^2_{ground}$) are the cumulative downwards leaf area indices at the top and the bottom of the canopy component i, and $\varphi$ (−) is the leaf surface area fraction relative to entire leaf area (1 for Lumped leaves, $\le 1$ for sunlit and shaded leaf fractions).

If simulations are performed for a BigLeaf canopy, $L_u$, i and $L_t$, i are set to 0 and $L_t$, respectively, where $L_t$ ($m^2_{leaf}$, $m^2_{ground}$) is the total leaf area index of the canopy. For a Lumped canopy representation, $\varphi$ is set to 1, otherwise if the canopy is Sunlit-Shaded the value of $\varphi$ changes between 0 and 1 distinctly for the sunlit and shaded leaf fractions, depending on the sun inclination (see Eqs. A12 and A53 for details).

In order to calculate the bulk surface resistance ($r_s$, i, h⁻¹ m), g in Eq. 1 is set to the stomatal conductance at the cumulative leaf area index L ($g_s$, L, m h⁻¹) which is described assuming that leaf surface temperature within each leaf layer is constant (Leuning et al., 1995; Kelliher et al., 1995; Wang and Leuning, 1998). That is to say, solar shortwave irradiance is assumed to be the main driver of the variation of surface conductance inside the canopy component i (see Eqs. Eq. A51 to Eq. A58 for details).

The bulk boundary layer conductance is composed of both the free ($g_b$, free, L, m h⁻¹) and forced ($g_b$, forced, L, m h⁻¹) convection conductances. $g_b$, free, L is a function of leaf-to-air temperature depression (Monteith and Unsworth, 2013) and is hence constant within each canopy component i (since leaves within a canopy component are assumed to have equal temperature). The temperature depression is defined as the difference in temperature between the leaf (or the canopy) surface and the ambient air at the measurement height. $g_b$, forced, L is a function, among other factors, of wind speed in the direct vicinity of the leaf (Jones, 1992) which changes with depth inside the canopy, expressed in terms of the cumulative leaf area index L (see Eqs. Eq. A39 to Eq. A48 for details).

## 2.1.2. Aerodynamic resistance and correction for stability conditions

The sensible and latent energy fluxes between the canopy source and reference heights are largely determined by forced convection ($r_{a,0}$, forced, h m⁻¹) under moderately stable to unstable conditions and by free convection ($r_{a,0}$, free, h m⁻¹) under strongly stable conditions (Monteith and Unsworth, 2013). $r_{a,0}$, forced and $r_{a,0}$, free are assumed to act in parallel and their relative contribution to the total aerodynamic resistance ($r_{a,0}$, h m⁻¹) is weighted by an empirical factor that is a function of the stability condition expressed by the Richardson number ($R_i$, −) and the Monin-Obukhov length ($\xi_{MO}$, m). See equations Eq. A28 to Eq. A38 for details.

## 2.1.3. Model resolution scheme

The surface temperature of individual components, $T_s$, i controls all energy balance processes at a given time step. Therefore, solving the crop energy balance implies determining $T_s$, i.

We developed an iterative solution to solve the energy balance

<figure><img src="image_4.png" /><figcaption>Fig. 2. Flowcharts of the canopy energy balance model MONTPEL. (a) Overall calculation procedure. All state variables are first initialized assuming neutral atmospheric conditions (b). The model is then run once to solve the energy balance (c). If stability correction is to be performed, the stability-related state variables are updated (d) and the aerodynamic resistance is then updated (a) and the transient energy balance is run with the new values of the state variables (c). This procedure is repeated until convergence is reached (a), otherwise the aerodynamic resistance is forced by assuming neutral stability conditions (e) and the energy balance is calculated for a last time with the forced aerodynamic resistance value. Refer to Supplementary Tables S2 and S3 for variables description.</figcaption></figure>
equations (Eq. A1 to Eq. A59), in which $T_s$, $i$ is the state variable. This solution relies on three levels of iteration loops. The first (inner) level (Fig. 2c) determines the energy balance values calculated with a fixed aerodynamic resistance. This is the only required loop if stability correction is not required. The second (outer) level (Fig. 2a) solves energy balance with varying aerodynamic resistance values that are adjusted for atmospheric stability. The third (intermediate) level (Fig. 2d) is required for stabilizing intermediate stability correction variables (correction functions and friction velocity).

We provided the detailed description of the resolution scheme in the appendix a.3. Resolution scheme.

## 2.2. Analysis of the model behavior

Our analysis of the behavior of MONTPEL considered three aspects. First, the appropriate implementation of the upscaling procedure, that is to say, that the model must simulate the same canopy latent and sensible heat fluxes regardless of the number of components and the thickness of the leaf layer. Second, the simulated latent and sensible heat fluxes, in addition to surface temperature, must respond to environmental stimuli as “expected” (e.g. an increasing soil water deficit must increase both surface temperature and sensible energy fluxes and decrease latent energy fluxes). Finally, this analysis examined how canopy representation (Lumped vs. Sunlit-Shaded) impacts the simulated fluxes.

Model simulations were conducted for a sunny day using the four possible canopy representations (BigLeaf Lumped, Layered Lumped, BigLeaf Sunlit-Shaded, and Layered Sunlit-Shaded). $L_t$ was set at $4 m^2_{leaf} m^{-2}_{ground}$. The canopy was divided into four layers with the same leaf area index of $1 m^2_{leaf} m^{-2}_{ground}$. The same canopy representations and weather data were then used to run the model under well-watered (WW) and water-deficit (WD) conditions.

### 2.2.1. Sensitivity analysis and error source classification

We performed a sensitivity analysis of the model by using the Extended Fourier Amplitude Sensitivity Test (eFAST, Satelli et al., 1999).

We conducted the sensitivity analysis over four canopy representations (BigLeaf Lumped, BigLeaf Sunlit-Shaded, Layered Lumped and Layered Sunlit-Shaded), four combinations of high (H) and low (L) incident solar radiation and air vapor pressure deficit conditions (HH: 775 W $m^{-2}_{ground}$, 3.54 kPa; LH: 416 W $m^{-2}_{ground}$, 1.78 kPa; HL: 914 W $m^{-2}_{ground}$, 0.57 kPa; LL: 416 W $m^{-2}_{ground}$, 0.18 kPa), and three soil water conditions, respectively well-watered ($\psi_{soil} = 0$ MP), mild water deficit ($\psi_{soil} = -0.3$ MPa) and severe water deficit ($\psi_{soil} = -3$ MPa). We varied all model parameters by 20 % around their nominal values (see Supplementary Table S1 for parameter values).

We examined the effect of model parameter values on canopy total latent energy flux and surface temperature of the whole canopy (BigLeaf representation) and of its upper and lower leaf layers (Layered representation) for average leaves (Lumped representation) or for sunlit and shaded leaves (Sunlit-Shaded representation).

We further analyzed the error of simulated vs. measured canopy temperature. First, we used a linear regression against key weather inputs and model state variables. Second, we used a non-parametric supervised learning method to construct a regression tree using the tree module provided in the `sklearn` package in Python™ (Pedregosa et al., 2011). The temperature measurements used in this analysis come from the experiments detailed in the following section (2.2.2).

### 2.2.2. Field experiments for model evaluation

The model accuracy was evaluated by comparing model outputs (crop net radiation, latent and sensible heat flux, soil heat flux and canopy surface temperature) against four wheat (Triticum aestivum L.) experiment datasets collected on deep loam or shallow sandy soils, and

semi-arid or temperate climatic conditions. Datasets included well-watered, water-deficit, and well fertilized treatments. Experiment details are summarized in Supplementary Table S4 and each experiment is briefly described below.

Two among the four experiments were conducted at Maricopa, AZ, USA. This site has a semi-arid climate and a clay loam soil for the top $m$ with more sandy soil further down. The first experiment (Kimball et al. 2017) was a Free-Air CO$_2$ Enrichment experiment (FACE), hereafter named “Maricopa FACE”. It was conducted during the 1992–1993, 1993–1994, 1995–1996 and 1996–1997 growing seasons with the spring wheat cultivar Yecora Rojo. The experiments included WW and WD conditions in 1993 and 1994, Nitrogen (N) sufficient and N deficit conditions in 1995–1996 and 1996–1997. The second experiment (Wall et al., 2011; Ottman et al., 2012; Martre et al. 2018) was a free-air heated experiment, known as the Hot Serial Cereal (HSC) experiment, hereafter named “Maricopa HSC”. In this experiment, the spring wheat cultivar Yecora Rojo was sown every about six weeks from March 2007 through January 2009 under WW and N sufficient conditions. In this study, we used the data from six sowing dates of the unheated reference treatments, which were only equipped with dummy heaters.

The third dataset, hereafter named “Grignon” (Bernard et al. 2022), comes from an experiment conducted at Thiverval-Grignon, France, during the 2011–2012 growing season. The experimental site has a temperate climate and deep (about 2 m) silt loam soil. The experiment was performed with the winter wheat cultivar Trémie sown at high (250 seeds $m^{-2}$) and low (180 seeds $m^{-2}$) densities. The plots were rainfed but did not experience any significant water deficit. All plots were grown under sufficient N conditions.

Finally, the fourth dataset, referred to as “Branschweig FACE” (Manderscheid et al. 2020) comes from a FACE experiment conducted at Braunschweig, Germany with the winter wheat cultivar Batis, during the 2013–2014 and 2014–2015 growing seasons. The site has a temperate climate and a shallow (about 0.5 m) loamy sand soil. All plots were grown under WW and N sufficient conditions.

These four datasets provide variables that allow comparing model outputs to measurements at both the crop and leaf-layer scale. The Maricopa FACE dataset includes hourly crop energy balance variables (total net radiation, soil heat flux, latent, and sensible heat fluxes), hourly canopy surface temperature from stationary mast-mounted infrared thermometers (IRTs), as well as midday sunlit and shaded leaf temperature from hand-held portable IRTs. The Maricopa HSC dataset include essentially hourly crop temperature measurements. The Grignon dataset includes measurements of leaf temperature at several depths across the canopy, and finally the Braunschweig dataset include additional measurements on canopy surface temperature.

The evaluation of the model focused on the impact of ambient weather conditions and soil water status on the crop energy balance. Therefore, the elevated CO$_2$, elevated heat, and N deficit treatments data were not used. In addition, in order to avoid over-estimating temperatures due to the hotspot effect (Kuusk, 1991), we excluded the measurements taken at solar noon with stationary IRTs, which were pointed North. These represented 3 % of the available data. Finally, the measured IRT temperatures were corrected to account for canopy emittance and the reflected sky radiation as described by Pinter et al. (2000). The sky radiation adjustment was computed based on air temperature and vapor pressure, determined from wet bulb measurements using aspirated psychrometers, applying the 8–14 µm equation derived by Idso (1981).

The total number of measured hourly data points that were compared to simulations totaled 5851.

In Maricopa FACE, Grignon and Braunschweig datasets, canopy height measurements were partly or completely missing. Canopy height was therefore estimated from empirical functions that relates the rate of canopy height increase to the cumulative degree day temperature (GDD, °Cd) as described in Supplementary Methods SM1. Soil water potential was estimated from volumetric soil water measurements
through water retention curves as described in Supplementary Methods SM2.

Finally, canopy temperature in the Maricopa FACE, Maricopa HSC and Braunschweig FACE experiments was measured with IRTs that pointed at the canopy surface from a specific distance and inclination angle. Therefore, the measured radiometric temperatures included both sunlit and shaded leaves, as well as soil surface before the canopy was closed. For these experiments, the simulated canopy temperature was calculated as the weighted temperature of sunlit and shaded leaves as seen by virtual IRTs which had the same height and inclination of those used in the experiments, as described in Supplementary Methods SM3.

### 2.3. Model evaluation metrics

We evaluated the agreement between simulated and measured variables with distance- and correlation-based metrics. We used the coefficient of determination ($R^2$), the root mean squared error (RMSE) and the normalized Nash-Sutcliffe modeling efficiency (nNSE) given by:

$$R^2 = \frac{\left[ \sum_{i=1}^{n} \left( (\mathbf{y}_{sim. i} - \bar{\mathbf{y}}_{sim}) (\mathbf{y}_{obs. i} - \bar{\mathbf{y}}_{obs}) \right) \right]^2}{\sum_{i=1}^{n} \left( \mathbf{y}_{sim. i} - \bar{\mathbf{y}}_{sim} \right)^2 \cdot \sum_{i=1}^{n} \left( \mathbf{y}_{obs. i} - \bar{\mathbf{y}}_{obs} \right)^2} \quad (2)$$

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n} (\mathbf{y}_{sim. i} - \mathbf{y}_{obs. i})^2}{n}} \quad (3)$$

$$nNSE = \left( 1 + \frac{\sum_{i=1}^{n} (\mathbf{y}_{sim. i} - \mathbf{y}_{obs. i})^2}{\sum_{i=1}^{n} (\mathbf{y}_{obs. i} - \bar{\mathbf{y}}_{obs. i})^2} \right)^{-1} \quad (4)$$

where $\mathbf{y}_{sim. i}$ and $\mathbf{y}_{obs. i}$ are the $i^{th}$ simulated and measured values, respectively, $\bar{\mathbf{y}}_{sim}$ and $\bar{\mathbf{y}}_{obs}$ are mean simulated and measured values, respectively, and $n$ is the total number of variable pairs.

We considered the model performance for simulating the energy balance terms as satisfactory when RMSE values were lower than 50 W m⁻² and $R^2$ higher than 0.55. These values were obtained from the data reported by Norman and Kustas (1995), Leuning et al. (1998), Kustas and Norman (1999), Arora (2003), Blonquist et al. (2010), Colaizzi et al. (2012) and Song et al. (2016). We selected these studies for establishing the performance criteria since they reported RMSE and $R^2$ values for at least two of the energy balance terms.

Model performance for simulating surface temperature was considered satisfactory when RMSE and $R^2$ values were respectively lower than 2.9 °C and 0.55 (values average from Kimball et al. 2015 and Webber et al., 2016, Webber et al., 2017, Webber et al., 2018).

To further investigate the sources of the model errors, the mean squared error was decomposed into its three components, respectively the squared bias (SB), the none-unity slope (NU) and the lack of correlation (LC) following Gauch et al. (2003):

$$SB = (\mathbf{y}_{sim} - \bar{\mathbf{y}}_{obs})^2 \quad (5)$$

$$NU = \left( 1 - \frac{\sum_{i=1}^{n} \left( (\mathbf{y}_{sim. i} - \bar{\mathbf{y}}_{sim}) (\mathbf{y}_{obs. i} - \bar{\mathbf{y}}_{obs}) \right)}{\sum_{i=1}^{n} (\mathbf{y}_{sim. i} - \bar{\mathbf{y}}_{sim})^2} \right)^2 \frac{\sum_{i=1}^{n} (\mathbf{y}_{sim. i} - \bar{\mathbf{y}}_{sim})^2}{n} \quad (6)$$

$$LC = (1 - R^2) \frac{\sum_{i=1}^{n} (\mathbf{y}_{obs. i} - \bar{\mathbf{y}}_{obs})^2}{n} \quad (7)$$

## 3. Results and discussion

### 3.1. Model behavior

Canopy representation did not affect the simulated energy balance terms at the canopy scale (Fig. 3). Identical total net radiation (Rn), soil heat flux (G), sensible heat flux (H⁰) and latent energy flux (λE⁰) were obtained with both Layered and BigLeaf canopy representations, which implies that the up scaling of resistances was appropriate.

Yet, distinguishing fluxes from sunlit and shaded leaf fractions resulted in slightly lower λE⁰ (0.6 %) and consequently higher H⁰ (4 %) fluxes compared to the case where a Lumped canopy was considered (Fig. 3, panels b vs. a, d vs. c, f vs. e, and h vs. g). This is due to the nonlinear response of the stomatal conductance to the incident solar irradiance (Eq. A52), combined with the non-linear response of the sunlit leaf fraction with the cumulative leaf area index with depth (Eqs. A57 and A58). The result of combining these two non-linear responses is that the equivalent surface conductance of sunlit and shaded leaf fractions (weighed by their respective leaf fractions) is systematically lower than that for their equivalent Lumped leaf (Fig. 4a). This behavior is consistent across a wide range of leaf area indexes and diffuse-to-total irradiance ratios (Fig. 4b), implying that the simulated Bowen ratio (β = H⁰/λE⁰) is lower when a Lumped canopy representation is considered compared to a Sunlit-Shaded one.

This behavior is consistent with results reported by Dai et al. (2004) from a study comparing energy balance simulations to measurements for two contrasting forest datasets, using BigLeaf Lumped and BigLeaf Sunlit-Shaded canopy representations, although stomatal conductance formalism was different in our study and in Dai et al. (2004) - the latter study considered a feedback control from photosynthesis). Conversely, Arora (2003) reported higher simulated canopy conductance rates for winter wheat using BigLeaf Sunlit-Shaded canopy representation compared to a BigLeaf Lumped one. The difference between their results and those of this study are likely due to different upscaling procedure. Arora (2003) obtained canopy surface conductance by directly applying the sub-leaf-based stomatal conductance model of Leuning (1995) to the canopy scale, assuming all leaves of sunlit or shaded canopy fractions to absorb identical solar irradiance rates. In MONTPEL and in Dai et al. (2004), leaf stomatal conductance is first calculated using leaf-level irradiance then up-scaled to the canopy using the integral process described by Raupach (1995).

An example of simulated irradiance absorption, leaf fraction and the corresponding leaf surface temperature profiles for well-watered and water-deficit conditions is shown in Fig. 5 for all possible canopy representations. The absorbed irradiance of Lumped leaf layers closely matched that of the sunlit fraction at all depths across the day (Fig. 5a, e, i, m and q). However, the Lumped leaf temperature was closer to sunlit leaf temperature only at the top of the canopy and progressively shifted towards shaded leaf temperature with canopy depth (Fig. 5c, g, k, o and s for WW and Fig. 5d, h, l, p and t for WD). This behavior has been reported by Leuning et al. (1995) and Wang and Leuning (1998) who indicated the shaded leaves to largely determine the energy fluxes at lower depths in the canopy.

The simulated sunlit leaf temperature increased with canopy depth (Fig. 5c, g, k, o and s for WW and Fig. 5d, h, l, p and t for WD), despite a decreasing radiative load per unit ground area (Fig. 5a, e, i, m and q). This behavior is expected and has been reported by Wang and Leuning (1998). It is due to the steep decrease of wind speed within the canopy (Eq. A42) which counterbalanced the decrease of available energy with depth (Supplementary Fig. S1).

The above example demonstrates the difficulty of determining a “mean” canopy temperature for use in temperature-dependent physiological processes. Inclán and Forkel (1995) highlighted that surface temperature prediction may lack accuracy with BigLeaf canopy representations when compared to a Layered representation. However, we postulate that this may become problematic only when the canopy is
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Daily time course of simulated wheat crop net radiation (Rn), crop sensible heat flux (Hn), latent energy flux (λE^n), and soil heat flux (G) simulated under well-watered (a to d) and water-deficit (e to h) conditions. The canopy was represented as one leaf (BigLeaf; a, b, e and f) or with discrete leaf layers (Layered; c, d, g and h) and considering Lumped leaves (Lumped; a, c, e and g) or Sunlit-Shaded leaf fractions (Sunlit-Shaded; b, d, f and h). For simulations with a Layered canopy (c, d, g, h) the canopy had four leaf layers of identical leaf area indexes of 1 m²leaf m⁻²ground. Simulations were performed for wheat (with nominal parameter values in Table S1) with weather conditions at Grignon, France on 29 June 2019 (maximum incident solar global radiation of 940 W m⁻²ground, minimum and maximum vapor pressure deficit of 0.26 and 3.97 kPa, respectively).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** (a) Leaf surface (stomatal) conductance as a function of the downward cumulative leaf area index for Lumped leaves, sunlit, and shaded leaf fractions. 'sunlit+shaded' denotes the aggregated conductance for both sunlit and shaded leaf fractions with depth; (b) ratio of the surface resistance simulated with Lumped and Sunlit-Shaded BigLeaf canopies versus the diffuse-to-total irradiance ratio for different values of leaf area index. Simulations were performed for wheat (with nominal parameter values in Table S1) with weather conditions at Grignon, France on 29 June 2019 (maximum incident solar global radiation of 940 W m⁻²ground, minimum and maximum vapor pressure deficit of 0.26 and 3.97 kPa, respectively).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Daily time course of simulated vertical profiles of absorbed irradiance, shaded leaf area fraction, and temperature for a canopy having a leaf area index of 4 m²/leaf m⁻², represented with four canopy representations. (a, e, i, m and q) Simulated absorbed irradiance, (b, f, j, n and r) shaded leaf surface area fraction, (c, g, k, o and s) temperatures under well-watered conditions and (d, h, l, p and t) temperature under water deficit conditions. Circles are simulated absorbed irradiance, fraction of shaded leaf, and canopy temperature profiles for *Layered* canopies, vertical dashed lines are simulated values for *BigLeaf* canopies, and squares are simulated values for the soil surface. For simulations with a *Layered* canopy, the canopy had four leaf layers of identical leaf area indexes of 1 m²/leaf m⁻². Simulations were for wheat (with nominal parameter values in Table S1) for weather conditions at Grignon, France on 29 June 2019 (maximum incident solar global radiation of 940 W m⁻², minimum and maximum vapor pressure deficit of 0.26 and 3.97 kPa, respectively). Simulations were for a silty soil with water potentials of 0 and ~3 MPa under well-watered and water deficit conditions, respectively.</figcaption></figure>
sparse and the fraction of sunlit leaves is high.

Finally, it can be shown from Fig. 5 that soil surface temperature, leaf surface temperature, and sunlit-to-shaded leaf surface temperature difference was higher under WD (Fig. 5d, h, l, p, and t) compared with WW (Fig. 5c, g, k, o, and s) conditions. It is noteworthy to highlight that the increase of soil temperature does not preferentially increase the temperature of the lowermost leaf layer, since the multi-component energy balance model does not consider the spatial arrangement of crop components. The resulting behavior is that the increase of the temperature of one component leads to an identical shift in the temperature of all other components (Fig. 6).

## 3.2. Sensitivity analysis

The simulated latent energy flux and leaf surface temperature were most sensitive to the maximum stomatal conductance ($g_s$, max) under almost all combinations of weather and soil water scenarios (Fig. 7). This sensitivity increased with soil water deficit.

The model showed a moderate sensitivity to the aerodynamic resistance parameters which were not affected by canopy representation. The empirical parameters $\eta$ (weight attributed to the free convection resistance to the overall aerodynamic resistance, Eq. A30), $C_d$ (drag coefficient), and $\xi$ (ratio of heat to momentum transfer roughness lengths) affected moderately the simulated outputs, and mostly under well-watered conditions.

The only parameters related to the canopy representation to which the simulated temperature and latent energy flux were sensitive were the leaf boundary-layer resistance parameter $\alpha$ (Eq. A41) and the soil surface resistance parameter $a_s$ (Eq. A59). Both parameters were nevertheless highly sensitive only under the LL weather scenario (LL: low radiation and water vapor deficit).

The parameter $\alpha$ had a large influence only on sunlit leaf surface temperature under low weather demand conditions and well-watered conditions. Its sensitivity was higher for the upper than for lower layers. This was probably due to the decreasing wind speed with depth into the canopy, which reduces the forced convection boundary layer

conductance and its mixing effect on leaf temperature. Moreover, irradiance flux density (irradiance per unit leaf surface area) is highest for sunlit leaf fractions, for which the convective mixing with wind plays a greater role in determining the temperature of the leaf surface. Model sensitivity to $\alpha$ decreased with water deficit. The shape parameter controlling the response of soil surface resistance to soil humidity ($a_s$) influenced the sensible heat flux mainly under well-watered conditions.

Overall, the sensitivity analysis shows that for weather conditions that are likely to prevail over the life cycle of a crop (i.e. all but the LL weather scenario), the maximum stomatal conductance is, by far, the most sensitive parameter. This result agrees with previous studies indicating the canopy surface conductance to have the strongest impact of simulated energy balance terms (e.g. Leuning et al., 1998; Arora, 2003; Lagos et al., 2013).

## 3.3. Model accuracy

### 3.3.1. Simulation of energy balance terms

The simulated crop energy balance terms (net radiation, latent, sensible and soil heat fluxes) agreed well with those measured for the Maricopa FACE experiment, but only satisfactorily for the sensible and ground fluxes (Fig. 8a to d). The RMSE were approximately 10% of the maximum measured values and all $R^2$ and $n$NSE were higher than 0.65 and 0.72, respectively, for all energy balance terms. The simulated net radiation was consistently overestimated (43 W m$^{-2}$ ground in average) compared to measurements (Fig. 8a). This may result from an underestimation of the surface albedo (canopy plus soil) as reported by Leuning et al. (1998) and Blonquist et al. (2010), due to the assumption of a constant leaf scattering property over the simulated period. Leuning et al. (1998) highlighted that such assumption may lead to underestimating surface albedo, which increases with the green area index and with the ratio of direct to global irradiance, resulting in overestimated irradiance absorption.

Blonquist et al. (2010) further reported that most of net radiation overestimation in Penman-Monteith equation may originate from underestimated net longwave losses when the latter are calculated following the empirical formula of Brunt (1932) as in the FAO-56 paper (Allen et al., 1998) which was considered in MONTPEL (Eq. A20). Blonquist et al. (2010) reported this empirical equation to typically underestimate the net longwave loss by approximately 20% on average, which may explain the overestimation trend of the simulated net radiation values presented in this paper.

A third reason for net radiation overestimation comes from the oversimplified equation used for estimating the soil heat flux (G) on hourly basis (Eq. A26) following Allen et al. (1998), i.e. as a fraction of the crop net radiation (Rn), without accounting for soil cover status. Under low leaf area index ($L_t$), MONTPEL strongly underestimated the amplitude of the diurnal variation of the simulated soil heat flux (Fig. 8k, 31 Jan. to 05 Feb.) and this trend was inversed for higher $L_t$ (Fig. 8k, 28 Feb. to 04 Mar.) This is consistent with the results obtained by Colaizzi et al. (2012) who reported the simulated G errors to increase under sparse canopy conditions. This highlights the limits of the model of Allen et al. (1998), which provides a robust estimation of the G/Rn fraction for closed canopies when a general description of the latter is needed (and Maricopa FACE data provide a good example when all measurements are pooled together, Supplementary Fig. S2). However, the high data discrepancy shown in Supplementary Fig. S2 highlights how the G/Rn fraction method is not suitable for hour-to-hour simulation accuracy and this is due to the omission of factors affecting this fraction that are not considered in this method.

MONTPEL tended to overestimate the latent energy flux but to underestimate sensible heat flux (Fig. 8b and c). Similar sensible to latent energy flux trade-offs were reported in previous studies where a two-source energy balance model based on or similar to Shuttleworth and Wallace (1985) was used (Kustas et al., 1996; Cammalleri et al., 2012;

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Effect of the increase of soil temperature with soil water deficit on canopy layers temperature. Simulations are for a Lumped leaf representation with four leaf layers of identical leaf area indexes of $1 \text{ m}^2_{\text{leaf}} \text{ m}^{-2}_{\text{ground}}$. Simulations were for wheat (with nominal parameter values in Table S1) for weather conditions at Grignon, France on 29 June 2019 at midday (incident solar global radiation of 933 W m$^{-2}$ ground, vapor pressure deficit of 3.09 kPa, air temperature of 32.3 °C). According to the multi-component approach, all components share the same vapor pressure deficit.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Extended Fourier Amplitude Sensitivity Test (eFAST) of model parameters sensitivity for latent energy flux and leaf temperature. The analysis was performed for four canopy representations (*BigLeaf Lumped*, *BigLeaf Sunlit-Shaded*, *Layered Lumped* and *Layered Sunlit-Shaded*) and four combinations of incident solar radiation (*R*<sub>inc, sw</sub>) and air vapor pressure deficit (*D*<sub>a</sub>): HH, high *R*<sub>inc, sw</sub> (775 W m<sup>-2</sup> ground) and high *D*<sub>a</sub> (3.54 kPa); LH, low *R*<sub>inc, sw</sub> (416 W m<sup>-2</sup> ground) and high *D*<sub>a</sub> (1.78 kPa); HL, high *R*<sub>inc, sw</sub> (914 W m<sup>-2</sup> ground) and low *D*<sub>a</sub> (0.57 kPa); LL, low *R*<sub>inc, sw</sub> (416 W m<sup>-2</sup> ground) and low *D*<sub>a</sub> (0.18 kPa) and three soil water conditions: (a) well-watered ψ<sub>soil</sub> = 0 MPa); (b) mild water deficit (ψ<sub>soil</sub> = -0.3 MPa) and; (c) severe water deficit (ψ<sub>soil</sub> = -3 MPa). Simulations were for wheat (with nominal parameter values in Table S1). Weather data were from Grignon, France (HH: 29 June 2019 at 14:00, LH: 14 October 2018 at 14:00, HL: 23 April 2018 at 11:00, and LL: 27 February 2018). Model parameters are grouped following the processes they regulate.</figcaption></figure>

Song et al., 2016; Xue et al., 2023). We were not able to correct this error without compromising the accuracy of canopy temperature predictions which suggested that this behavior was not related to a lack of appropriate parameterization. Rather, such behavior may result from the

conceptual structure of the multi-component energy balance framework itself.

Indeed, Penman-Monteith's multi-component approach relies on the assumption of negligible aerodynamic resistance inside the canopy
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Measured and simulated hourly canopy temperature and energy balance terms. (a) to (d) Simulated versus observed crop net radiation, sensibly heat flux, latent energy flux, and soil heat flux during daytime (light blue) and nighttime (dark blue) over the whole growing season for all experiments. Dashed lines are 1:1 relationships. (e) to (k) Measured (dots) and simulated (lines) leaf area index, canopy temperature, canopy temperature depression (canopy to air temperature difference), and energy balance terms for two periods of five consecutive days with low and high leaf area index values. Data are for the spring wheat cultivar Yecora Rojo (with parameter values in Table S1) grown in the field in the Maricopa FACE experiment during the 1995–1996 growing season. Simulations were performed using *BigLeaf* and *Sunlit-Shaded* canopy representations. *R*², coefficient of determination; *RMSE*, root mean squared error; *nNSE*, normalized Nash–Sutcliffe modeling efficiency.</figcaption></figure>
(Lhommé et al., 2013). While such assumption is reasonable over large areas with sparse canopies (conditions for which the multi-component approach was developed by Shuttleworth and Wallace, 1985), its use is less justified in dense canopies such as those of wheat were the leaf area index may reach as high as 12 $m_{leaf}^2 m_{ground}^{-2}$ (e.g. Moreau et al., 2012). Latent energy overestimation is frequently reported when multi-component energy balance models are used (e.g. Kustas et al., 1996; Colaizzi et al., 2012; Chen et al., 2021; Cuadra et al., 2021; Bao et al., 2022; Chen et al., 2022; Xue et al., 2023). Chen et al. (2022) reported that the two-source energy balance model of Shuttleworth and Wallace (1985) systematically overestimated the simulated latent energy flux of maize (Zea mays L.) canopies, regardless of the formalism used to describe canopy surface conductance. The authors furthermore reported the two-source energy balance to perform best for sparse canopies, while model accuracy decreased as the canopy densities increased. Chen et al. (2021) and Bao et al. (2022) reported similar results, which supports our hypothesis of overestimated latent energy fluxes with multi-component energy balance models. Further exploration for model improvement solutions is required but exceeds the scope of the present paper.

### 3.3.2. Simulation of canopy temperature

The simulated canopy temperature ($T_{canopy}$) and canopy-to-air temperature depression ($\Delta T_{canopy}$) agreed satisfactorily with measurements for all four datasets (Fig. 9). Note that “canopy temperature” refers here to the radiometric surface temperature, which is a measurable quantity, not the aerodynamic canopy temperature $T_m$ that is not directly measurable (see details in **Supplementary Methods SM3**). In all

experiments, $R^2$ was higher than 0.91 and 0.55 for $T_{canopy}$ and $\Delta T_{canopy}$, respectively. $nNSE$ was higher for $T_{canopy}$, ranging between 0.89 and 0.92, compared to $\Delta T_{canopy}$, which ranged between 0.24 and 0.66. Finally, RMSE was lower than 2.23 °C for all datasets, representing an error of less than 9% of the maximum value of $T_{canopy}$.

MONTPEL showed a slight tendency to overestimate canopy temperature for temperatures roughly greater than 15 °C (Fig. 9a, c and d). This tendency was further evident when the simulated sunlit and shaded leaves were compared to measurements (Fig. 10) showing an RMSE greater than 3.5 °C, although $R^2$ was greater than 0.7 and 0.37 for $T_{canopy}$ and $\Delta T_{canopy}$, respectively. This overestimation may be explained by the overestimated net canopy radiation as discussed in the previous section.

However, overall, canopy temperature simulations reported in this study outperformed previous results reported on Maricopa FACE and HSC experiments (for well-fertilized fields data). On Maricopa FACE, Webber et al. (2018) reported RMSE and $R^2$ for simulated $\Delta T_{canopy}$ within the range 1.36–2.94 °C and 0.03–0.52, respectively, obtained using seven crop models with different energy balance options (APSIM Nwheat, Ecosys, FASSET, SIMPLACE, SiriusQuality, Sirius, and STICS, refer to Webber et al., 2018 for model references). On Maricopa HSC, Webber et al. (2016) reported RMSE and $R^2$ for simulated $\Delta T_{canopy}$ within the range 0.3–4.1 °C and 0.016–0.3, respectively, using a modelling solution in the SIMPLACE framework. Later, Webber et al. (2017) simulated crop energy balance on Maricopa HSC using eight crop models with different energy balance options (HUMEWheat, Nwheat, FASSET, Hermes, Sirius2014, Sirius, SSM-Wheat, and SIMPLACE, refer to Webber et al., 2017 for model references). The authors reported RMSE and $R^2$ for simulated $\Delta T_{canopy}$ within the ranges of 2.9–6.7 °C and

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Simulated versus measured hourly canopy temperature (a to d) and canopy temperature depression (e to h) during daytime (light blue) and nighttime (dark blue) for the spring wheat cultivar Yecora Rojo grown in the field at Maricopa, Arizona, USA, in the FACE (a and e) and HSC (b and f) experiments, the winter wheat cultivar Trémie grown in the field at Grignon, France (c and g) and the winter wheat cultivar Batis grown in the field at Braunschweig, Germany. Simulations were performed using BigLeaf Sunlit-Shaded canopy representation for Maricopa FACE, Maricopa HSC and Braunschweig FACE experiments and Layered Sunlit-Shaded representation for Grignon experiment (only dataset with leaf temperature measurements over several canopy layers). $R^2$, coefficient of determination; RMSE, root mean squared error; $nNSE$, normalized Nash-Sutcliffe modeling efficiency. All simulations used the nominal parameter values in Table S1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Simulated versus observed daily maximum canopy temperature (a, b) and canopy temperature depression (c, d) for sunlit (a, c) and shaded (b, d) leaf fractions under well-watered (blue circles) and water deficit (orange circles) conditions for the spring wheat cultivar Yecora Rojo (with parameter values in Table S1) grown in the field at Maricopa, Arizona, USA, in the FACE experiment during the 1992–1993 growing season (Zadok stage from 25 to 93). Simulations were performed using BigLeaf and Sunlit-Shaded canopy representations. R², coefficient of determination; RMSE, root mean squared error; nNSE, normalized Nash–Sutcliffe modeling efficiency.</figcaption></figure>

0.016–0.3, respectively, on Maricopa HSC. Finally, Kimball et al. (2015) reported for Maricopa HSC RMSE for simulated *T*<sub>canopy</sub> between 2.4 and 7.1 °C, using a standalone energy balance model (i.e. crop state variables such as height or leaf area were inputs).

It is noteworthy to recall that the simulated canopy temperatures in this study were obtained with a standalone energy balance model, while those reported by Webber et al. (2016, 2017, 2018), Webber et al. 2017, Webber et al. 2018) were obtained with broader crop models which simulate a myriad of physiological processes, including leaf area, water use and growth dynamics, in addition to the crop energy balance. The performance metrics reported by the former studies may hence differ if energy balance models were to be used separately but such evaluation was not reported.

Finally, MONTPEL reproduced the leaf surface temperature with depth measured at the Grignon experiment as shown in Fig. 11 for two green area index (assumed equal to leaf area index) values of 0.19 and 1.34, respectively. Although an overall good agreement was obtained, greater differences between simulated and measured leaf temperature occurred during nighttime hours (Fig. 11a and b). Blonquist et al. (2010) argued that the underestimation of net radiation in crop energy balance models during nighttime was due to inadequate parameterization of the atmospheric emissivity in net longwave loss equation (Eq. A21), which may explain the underestimated nighttime canopy temperature values simulated. By fine-tuning the atmospheric emissivity model in this study, nighttime temperature errors were reduced (Supplementary Fig. S3). However, since the model presented this behavior only for the Grignon dataset, the same parameters were kept across all sites in order to allow sound comparison of model performance across all sites. It is

worthy to stress that the underestimation of canopy surface temperature during nighttime may also originate from other sources than over-estimated net longwave loss. For instance, Wei et al. (2023) reported that higher errors in the prediction of nighttime temperatures may be explained by inaccurate measurements due to dew formation during night, by conceptual limitations due to disregarding condensation effect on the energy balance (as in MONTPEL), or by extremely weak aerodynamic conductance simulations being concomitant with high long-wave loss rates.

### 3.3.3. Sources of error for canopy temperature

A significant linear relationship was obtained between almost all explored drivers (absorbed solar irradiance by leaves and soil surface, air temperature, canopy height, soil water potential, aerodynamic resistance, vapor pressure deficit, leaf area index, and net longwave radiation) and obtained canopy temperature error (Fig. 12a to j). This relationship was similar between daytime and nighttime simulations, except for wind speed, which had a stronger effect on temperature error for nighttime simulations (Fig. 12b).

Regression tree analysis (Fig. 12k) showed that the absorbed irradiance by leaves had by far the strongest influence on temperature error, accounting for almost 25 % of the total obtained error, followed by the vapor pressure deficit (15 %), and net longwave radiation (11 %). The aerodynamic resistance and its direct determinants (wind speed, canopy height, and air temperature) accounted for 30 %. The variables characterizing canopy structure (height and leaf area index) contributed to 11 % of the total error. Finally, the soil absorbed irradiance and water potential accounted for 14 % of the total error.
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Simulated and measured daily time course of vertical canopy temperature profiles for wheat canopies with a green area index of 0.19 $m_{leaf}^2$ $m_{ground}^{-2}$ (a, c, e, g, i and k) on 17 February 2012 and 1.34 $m_{leaf}^2$ $m_{ground}^{-2}$ (b, d, f, h, j and l) on 30 April 2012. The hourly simulated and measured leaf layer surface temperatures are shown for the whole day (a, b). Vertical temperature profiles are shown at 06:00 (c, d), 09:00 (e, f), 12:00 (g, h), 15:00 (i, j) and 18:00 (k, l). Data are for the winter wheat cultivar Trémie (with parameter values in Table S1) grown in the field at Grignon, France with a sowing density of 250 seeds m⁻². The measured temperature values (red circles) are shown for three independent replicates. Simulations were performed using Layered and Sunlit-Shaded canopy representation. Soil temperature was not measured in this experiment. Note that higher layer indexes are for upper leaves.</figcaption></figure>

The high temperature error attributed to irradiance inputs and net longwave radiation (45 %) suggests that a substantial improvement to model performance may be obtained with a more accurate estimation of the net crop radiation, as previously discussed in Section 3.3.1.

### 3.3.4. Impact of canopy representation and atmospheric stability correction on model performance

Canopy representation (Lumped vs. sunlit-Shaded) had an extremely weak effect on model performance for simulating the radiometric canopy temperature (Fig. 13a). Conversely, the correction for non-neutral atmospheric stability conditions improved the accuracy of simulated canopy temperature (Fig. 13b).

Our results agree with those obtained by Webber et al. (2017, 2018)

who showed that crop models that included the correction for atmospheric stability conditions simulated most accurately the canopy surface temperature. Further, we show herein that this improvement was higher under the semi-arid conditions of Maricopa than under the temperate conditions of Grignon and Braunschweig. The mean squared error (MSE) of simulated $T_{canopy}$ was reduced by 51 % and 41 % in the FACE and HSC experiments at Maricopa, respectively, compared with 23 % in Braunschweig FACE, while it increased by 7 % at Grignon. The reduction in MSE was mainly attributed to the reduction of the squared bias.

The most notable effect of the stability correction was obtained under unstable conditions for all experiments except Maricopa HSC (Fig. 13f and g). For the latter, stability correction was strongest under
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Analysis of the drivers of simulated canopy temperature error. Daytime (light blue dots, solid red lines) and nocturnal (dark blue dots, dashed red lines) differences between observed and simulated hourly canopy temperature versus (a) solar irradiance absorbed by canopies, (b) wind speed, (c) air temperature, (d) canopy height, (e) soil water potential, (f) absorbed solar irradiance by soils, (g) aerodynamic resistance, (h) vapor pressure deficit at source height, (i) green area index, and (j) net longwave radiation. Significant linear regressions are denoted by (*). (k) Contribution of each driver to the total error. This contribution was calculated from a regression tree of simulated temperature residues on the same data as for (a) to (j). Data are for the four field experiments used in this study: the spring wheat cultivar Yecora Rojo grown in the field at Maricopa, Arizona, USA, in the Maricopa FACE and Maricopa HSC experiments, the winter wheat cultivar Trémie grown in the field at Grignon, France, and the winter wheat cultivar Batis grown in the field at Braunschweig, Germany. Diurnal data are defined as those where the incident irradiance is higher than zero. All simulations used the nominal parameter values in Table S1.</figcaption></figure>

near neutral conditions, followed by, unstable then stable conditions. For Maricopa FACE and HSC, although more than 50 % of the simulated time steps fell under stable conditions (Fig. 13f) canopy temperature error reduction was highest under unstable conditions (Fig. 13g) which suggests a stronger modification of the forced aerodynamic resistance ($r_{a,0, \text{forced}}$, Eq. A29) under the latter conditions. This result contrasts with that reported by Tolk et al. (1995) who indicated that stability correction functions (four models used in their study) mostly impacted $r_{a,0, \text{forced}}$ under stable conditions. Tolk et al. (1995) indicated that wind speed higher than 3 m s⁻¹, which occurred in more than 50 % of measurements, significantly reduced the Richardson number and consequently strongly limited the stability correction effects on $r_{a,0, \text{forced}}$. In our study, wind speed in simulated Maricopa FACE and HSC data exceeded 3 m s⁻¹ in only 13 and 19 % of times, which may explain the

different results from those of Tolk et al. (1995). Our analysis further points to the importance of stability correction under low wind speed conditions.

## 4. Summary and conclusions

We presented MONTPEL, a multi-component crop energy balance model, which simulates hourly crop surface temperatures and latent energy fluxes using different canopy representation schemes. We designed MONTPEL so that the canopy can be represented by one leaf layer (BigLeaf) or by multiple leaf layers (Layered), each layer being represented by average (Lumped) leaves or by sunlit and shaded leaf fractions (Sunlit-Shaded) separately.

The analysis of MONTPEL behavior, using all the possible canopy
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Effect of atmospheric stability correction on model error for canopy temperature. (a) Mean squared error for canopy temperature in the Maricopa FACE, Maricopa HSC, Grignon, and Braunschweig FACE field experiments, for canopies represented with Lumped or Sunlit-Shaded leaves assuming "Neutral" or "Corrected" stability conditions for the aerodynamic resistance of the canopy boundary layer. (b) Effect of the canopy boundary layer stability correction on model lack of correlation (blue), non-unity slope (white) and squared bias (red) for simulated canopy with Sunlit-Shaded leaves, (c to d) measured air temperature, canopy temperature depression, and wind speed for each experiment. (e) Percentage of simulated time steps with unstable, neutral, and stable conditions and (f) reduction in absolute canopy temperature error with stability correction under each stability condition for each of the experiments for simulated canopy with Sunlit-Shaded leaves. All simulations used the nominal parameter values in Table S1.</figcaption></figure>

representations, under different weather and soil water conditions showed that it reproduced expected canopy temperature both at the whole canopy scale and at different depths throughout the canopy. The simulated energy balance terms were also coherent with expected results, showing the ratio of latent to sensible heat fluxes to reduce as soil dried down and inversely. The analysis also revealed an intrinsically unrealistic behavior of the multi-component model when multiple leaf layers were simulated, whereby increasing the soil surface temperature led to equally increasing all canopy component temperatures, regardless of their position in the canopy. Such a behavior is due to the absence of the impact of canopy structure on its interaction with the soil component.

Model behavior analysis also showed that the number of simulated leaf layers did not affect the simulated total canopy fluxes. Conversely, distinguishing sunlit from shaded canopy fractions systematically resulted in lower latent energy fluxes compared to *Lumped* canopy representation.

The sensitivity analysis, using the *BigLeaf* canopy representation,

revealed that MONTPEL was most sensitive to the maximal stomatal conductance ($g_s$, max) parameter. However, using a single $g_s$ max value across the simulated experiments yielded satisfactory results. This suggests that model's sensitivity to the temporal and site-to-site variability of $g_s$ max is weak.

The comparison between simulated and measured energy balance terms, using the *BigLeaf* canopy representation, highlighted the trend of latent energy fluxes overestimation using MONTPEL. Similar results were reported with multi-component models (Kustas et al., 1996; Chen et al., 2021; Cuadra et al., 2021; Bao et al., 2022; Chen et al., 2022; Xue et al., 2023). This suggests that such a behavior is related to the conceptual framework of multi-component models and more precisely to the assumptions of the absence of aerodynamic resistance inside the canopy. Notwithstanding this limit, MONTPEL demonstrated satisfactory performance across a wide range of weather and soil conditions. The simulated temperature and energy balance terms followed closely the measured values, with RMSE of energy balance terms lower than 10% of maximum measured values and $R^2$ higher than 0.65. The simulated
canopy surface temperature and temperature depression compared well to measurements, with R² higher than 0.91 and 0.55, respectively, and RMSE lower than 2.3 °C for both variables.

We demonstrated using BigLeaf canopy representation, that accounting for atmospheric stability greatly improved the model accuracy in simulating canopy temperature, in particular under well-watered and high evaporative demand conditions. This result supports earlier reports in literature calling for the inclusion of stability correction functions in crop energy balance models to improve the accuracy of predicted crop surface temperature (Webber et al., 2017, Webber et al., 2018). Better estimation of extreme temperature conditions will enable better simulation of the likely effects on crop physiology and consequent effects of future climate change.

## Data availability

MONTPEL was implemented as a stand-alone program written in Python™. The source code of the model is freely available at https://github.com/RamiALBASHA/crop_energy_balance.

The data used in this paper and the analysis source code are freely available at https://github.com/RamiALBASHA/energy_balance_model_paper.

All the original field data used in this study are available from the cited references.

## CRediT authorship contribution statement

**Rami Albasha:** Conceptualization, Methodology, Software, Validation, Investigation, Writing – original draft, Writing – review & editing.
**Loïc Manceau:** Software. **Heidi Webber:** Methodology, Writing – review & editing. **Michaël Chelle:** Investigation, Writing – review & editing. **Bruce Kimball:** Investigation, Writing – review & editing. **Pierre Martre:** Conceptualization, Investigation, Writing – review & editing, Supervision.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported by the EU Project Horizon 2020 Solutions for Improving Agroecosystem and Crop Efficiency (grant no. 727247 to Rami Albasha, Loïc Manceau, and Pierre Martre) and the French Research Agency (grant no. ANR-19-CE20–0024 to Rami Albasha).

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.agrformet.2024.110221.

## Appendix A. Details of the MONTPEL model

### A.1. Theoretical framework

The MONTPEL model is an n-component energy balance model developed following Lhomme et al. (2013) who generalized the bi-component energy-balance model proposed by Shuttleworth and Wallace (1985). A component designates a crop entity that has distinguished properties (e.g. surface temperature, surface resistance, boundary layer resistance). A component may therefore be a leaf layer, a fraction of a leaf layer (e.g. sunlit or shaded leaf fractions), or the underlying soil.

Sensible and latent energy fluxes are simulated by analogy to electrical circuits whereby all components are connected in parallel at the crop source height ($z_m$, m).

Fluxes emanating from individual components converge at the source height. The vapor pressure deficit ($D_m$, kPa) at the source height is linked to the vapor pressure deficit ($D_a$, kPa) at reference (i.e. measurement) height ($z_r$, m):

$$D_m = D_a + \frac{r_{a,0}}{\rho c_p} [sA - (s + \gamma) \lambda E^n] \quad (A1)$$

where $r_{a,0}$ (h m⁻¹) is the aerodynamic resistance between source height and reference height, $\rho$ (g m⁻³) is the dry air density, $c_p$ (W h g⁻¹ K⁻¹) is the specific heat capacity under a constant pressure, $s$ (kPa K⁻¹) is the slope of the curve of saturated vapor pressure deficit at air temperature, $\gamma$ (kPa K⁻¹) is the psychrometric constant, $A$ (W m⁻² ground) and $\lambda E^n$ (W m⁻² ground) are the available energy and latent energy flux of the whole crop (vegetation and soil), respectively.

Crop latent energy flux is the sum of latent energy fluxes emanating from all components ($\lambda E^n = \sum_{i=1}^n \lambda E_i$) which can be rearranged to give (Lhomme et al. 2013):

$$\lambda E^n = R_0 \lambda E_p \sum_{i=1}^{n} P_i + \left(\frac{s}{\gamma}\right) \sum_{i=1}^{n} P_i A_i r_{a, i} \quad (A2)$$

where $R_0$ (h m⁻¹) is the Lumped aerodynamic resistance, $\lambda E_p$ (W m⁻² ground) is the potential latent energy flux of the crop according to Penman (1948) (i.e. without considering the surface resistance of the crop), $P_i$ (m h⁻¹) is the reciprocal of the equivalent boundary and surface resistances ($R_i$, h m⁻¹) of the component i, $A_i$ (W m⁻² ground) is the available energy to the component i and $r_{a, i}$ (h m⁻¹) is the boundary layer resistance of the component i. $R_0$, $P_i$, $R_i$ and $\lambda E_p$ are given by:

$$R_0 = \left(1 + \frac{s}{\gamma}\right) r_{a, 0} \quad (A3)$$
$$ P_i = \frac{1}{R_i \left[ 1 + R_0 \sum_{j=1}^{n} \frac{1}{R_j} \right]} \qquad (A4) $$

$$ R_i = r_s \cdot i + \left( \nu_i + \frac{s}{\gamma} \right) r_a \cdot i \qquad (A5) $$

$$ \lambda E_p = \frac{sA + \frac{\rho c_p D_a}{r_a 0}}{s + \gamma} \qquad (A6) $$

where $r_s, i$ (h m⁻¹) is the surface resistance of the component i and $\nu_i$ (−) is a parameter taking the value of 1 for amphistomatal leaves (e.g. wheat, maize) and the value of 2 for hypostomatal leaves (e.g. grapevine).

Finally, the latent energy flux from each individual component i is related to the latent energy flux of the whole canopy (Lhomme et al., 2013) and is given by:

$$ \lambda E_i = \frac{R_0 (\lambda E_p - \lambda E^n) + \frac{s}{\gamma} r_{a,i} A_i}{R_i} \qquad (A7) $$

## A.2. Variables determination

### A.2.1. Available energy

The crop available energy (A, W m⁻² ground) is given by:

$$ A = Rn - G = \sum_{i=1}^{n-1} (R_{v, sw, i} + R_{v, lw, i}) + (R_{s, sw} + R_{s, lw} - G) \qquad (A8) $$

where $Rn$ (W m⁻² ground) is the crop net radiation, $R_{v, sw, i}$ and $R_{v, lw, i}$ (W m⁻² ground) are the absorbed shortwave irradiance and net longwave radiation of the vegetative component i, respectively, $R_{s, sw}$ and $R_{s, lw}$ (W m⁻² ground) are the absorbed shortwave irradiance and net longwave radiation of the soil component, respectively, and $G$ (W m⁻² ground) is the soil heat flux.

The equations used for calculating the net energy terms differ whether the component is a Lumped leaf layer, a “Sunlit” or “Shaded” leaf fraction layer or a soil component. Shortwave, longwave and soil heat flux equations are detailed in the following sections.

#### A.2.1.1. Absorbed shortwave irradiance by canopies and leaf layers.
The absorbed shortwave irradiance by a leaf layer component (represented by a Lumped leaf or by sunlit and shaded leaf fractions) is an input to MONTPEL. Yet, the model implicitly assumes that irradiance absorption by leaves is calculated following Goudriaan (1988) for sunlit and shaded fractions and its derived version by de Pury and Farquhar (1997) for Lumped leaves.

For Lumped canopy representation, the absorbed shortwave irradiance is given by (de Pury and Farquhar, 1997):

$$ R_{v, sw, lumped, l} = R_{inc, sw, b} (1 - \rho_b) (e^{(-k_b L_{u, l})} - e^{(-k_b L_{l, l})}) + R_{inc, sw, d} (1 - \rho_d) (e^{(-k_d L_{u, l})} - e^{(-k_d L_{l, l})}) \qquad (A9) $$

where $R_{inc, sw, b}$ (W m⁻² ground) is the incident direct beam shortwave irradiance, $R_{inc, sw, d}$ (W m⁻² ground) is the incident diffuse shortwave irradiance, $\rho_b$ (−) is the canopy direct beam irradiance reflectance, $\rho_d$ (−) is the canopy diffuse irradiance reflectance, $k_b$ (m² ground m⁻² leaf) is the extinction coefficient of direct beam irradiance, and $k_d$ (m² ground m⁻² leaf) is the extinction coefficient of diffuse irradiance. $L_{u, l}$ and $L_{l, l}$ (m² leaf m⁻² ground) are the cumulative leaf area indices at the top and bottom of the leaf layer l. For canopies BigLeaf canopies we have $L_{u, l} = 0$ and $L_{l, l} = L_t$ where $L_t$ (m² leaf m⁻² ground) is the total leaf area index of the canopy.

When sunlit and shaded leaf fractions are considered, the absorbed shortwave irradiance is given by (Goudriaan, 1988):

$$ R_{v, sw, sunlit, l} = R_{inc, sw, b} (1 - \sigma_s) \left( e^{(-k_b L_{u, l})} - e^{(-k_b L_{l, l})} \right) + R_{inc, sw, d} (1 - \rho_d) \frac{k_d}{k_d + k'_b} \left( e^{(-k_d + k'_b) L_{u, l})} - e^{(-k_d + k'_b) L_{l, l})} \right) + R_{inc, sw, d'} \left[ \frac{(1 - \rho_d) k_d}{k_b + k'_b} \left( e^{(-k_b + k'_b) L_{u, l})} - e^{(-k_b + k'_b) L_{l, l})} \right) - \frac{(1 - \sigma_s)}{2} \left( e^{(-2k'_b L_{u, l})} - e^{(-2k'_b L_{l, l})} \right) \right] \qquad (A10) $$
$$ \mathbf{R}_{\mathrm{v, sw, shaded, l}} = \mathbf{R}_{\mathrm{v, sw, lumped, l}} - \mathbf{R}_{\mathrm{v, sw, sunlit, l}} \quad (A11) $$

where $\sigma_s$ (–) is the leaf scattering coefficient (i.e. the sum of leaf reflectance and transmittance coefficients) and $k'_b$ ($m_{\mathrm{ground}}^2 m_{\mathrm{leaf}}^{-2}$) is the extinction coefficient of black leaves for direct beam irradiance, and $k_b$ ($m_{\mathrm{ground}}^2 m_{\mathrm{leaf}}^{-2}$) is the extinction coefficient for direct beam irradiance with scattering, and $k_d$ ($m_{\mathrm{ground}}^2 m_{\mathrm{leaf}}^{-2}$) is the extinction coefficient for diffuse irradiance.

$k'_b$ is calculated assuming an ellipsoidal leaf angle distribution (Campbell, 1986) and is given by:

$$ k'_{b} = C \frac{\sqrt{\chi^2 + \cot^2 \beta}}{\chi + 1.774(\chi + 1.182)^{-0.733}} \quad (A12) $$

where $C$ (–) is the leaf clumping index, $\beta$ (rad) is the direct beam inclination angle, $\chi$ is a shape parameter related to the average inclination angle of the green elements of the canopy ($\bar{\alpha}$, rad) through the following empirical equation (Campbell, 1990):

$$ \chi = \left( \frac{\bar{\alpha}}{9.65} \right)^{-0.6061} - 3 \quad (A13) $$

Following Liu et al. (2021), $\bar{\alpha}$ is assumed to be constant between crop emergence and the beginning of stem extension and increases linearly during the stem extension period until flag leaf ligule appearance:

$$ \bar{\alpha} = \begin{cases} \bar{\alpha}_{juv}, & T_t \le T_{t,GS30} \\ \bar{\alpha}_{juv} + (T_{t,GS30} - T_t) \frac{(\bar{\alpha}_{mat} - \bar{\alpha}_{juv})}{P \times (FLN - HS_{30})}, & T_{t,GS30} < T_t < T_{t,GS39} \\ \bar{\alpha}_{mat}, & T_t \ge T_{t,GS39} \end{cases} \quad (A14) $$

where, $\bar{\alpha}_{juv}$ and $\bar{\alpha}_{mat}$ (rad) are the average inclination angle of green elements of canopies during the vegetative phase (i.e. before beginning of stem extension) and after flag leaf ligule appearance, respectively, $T_t$ (°Cd) is thermal time since crop emergence, $T_{t,GS30}$ and $T_{t,GS39}$ (°Cd) are $T_t$ at the beginning of stem extension and flag leaf ligule appearance, respectively, $P$ (°Cd leaf⁻¹) is the phyllochron, FLN is the final leaf number, and $HS_{30}$ is the decimal number of visible mainstem leaves (Haun stage) at GS30.

$k_b$ is approximated by Goudriaan (1988) as:

$$ k_b = k'_b \sqrt{1 - \sigma_s} \quad (A15) $$

$k_d$ is derived assuming that diffuse irradiance decreases exponentially with canopy depth by integrating $k_b$ over the hemisphere:

$$ k_d = -\frac{1}{L_t} \ln \left( 1 - 2 \int_{0}^{\pi/2} (1 - \exp(-k_b L_t)) \times \cos\beta \times \sin\beta \times d\beta \right) \quad (A16) $$

Similarly to extinction coefficients, reflectance coefficients depend on leaves angle distribution and the inclination of solar beam and sky sectors, for direct-beam and diffuse irradiances, respectively. $\rho_b$ is calculated following Goudriaan (1988) and is given by:

$$ \rho_b = 1 - \exp\left(-\frac{2\rho_h k'_b}{1 + k'_b}\right) \quad (A17) $$

where $\rho_h$ (–) is the reflectance coefficient of a canopy having horizontal leaves, given by:

$$ \rho_h = \frac{1 - \sqrt{1 - \sigma_s}}{1 + \sqrt{1 - \sigma_s}} \quad (A18) $$

As for $k_d$, $\rho_d$ is calculated as the integral over the hemisphere of Eq. A16:

$$ \rho_d = 2 \int_0^{\pi/2} \left( 1 - \exp\left(-\frac{2\rho_h k'_b}{1 + k'_b}\right) \right) \times \cos\beta \times \sin\beta \times d\beta \quad (A19) $$

Equations Eq. A16 and Eq. A19 were integrated numerically with the Gauss Legendre Rule method.

### A.2.1.2. Longwave radiation fluxes.
The net longwave radiation flux at the top of the canopy ($R_{0, lw}$, W m⁻² ground) is calculated following Brunt (1932) as:

$$ \mathbf{R}_{0, lw} = - \sigma_B \mathbf{T}_a^4 (1 - \varepsilon_{\mathrm{atm}}) (0.34 - 0.14\sqrt{\varepsilon_a}) \quad (A20) $$

where $\sigma_B$ (W m⁻² K⁻⁴) is the Stefan-Boltzmann constant, $\varepsilon_a$ (kPa) is the air vapor pressure, and $\varepsilon_{\mathrm{atm}}$ (–) is the atmospheric emissivity given by (Brutsaert, 1975):

$$ \varepsilon_{\mathrm{atm}} = 1.24 \left( \frac{10 \varepsilon_a}{T_a} \right)^{1/7} \quad (A21) $$
$R_{v, lw}$ is partitioned among leaf layers by considering that energy loss is attenuated through the canopy as a function of the diffuse irradiance extinction coefficient of black leaves (Leuning et al., 1995).

For canopies represented by Lumped leaf layers, the net longwave loss ($R_{v, lw, lumped, 1}$, W m⁻² ground) writes:

$$R_{v, lw, lumped, 1} = R_{0, lw} (e^{(-k'_b L_{u, 1})} - e^{(-k'_b L_{t, 1})}) \quad (A22)$$

For canopies represented by sunlit and shaded leaf fractions, $R_{v, lw, lumped, 1}$ is given by:

$$R_{v, lw, sunlit, 1} = R_{0, lw} \frac{k'_d}{k'_d + k'_b} (e^{(-k'_d + k'_b)L_{u, 1}} - e^{(-k'_d + k'_b)L_{t, 1}}) \quad (A23)$$

$$R_{v, lw, lumped, 1} - R_{v, lw, sunlit, 1} \quad (A24)$$

Soil net longwave loss is the part of total loss that has not been attributed to leaves and is given by:

$$R_{s, lw} = R_{0, lw} e^{(-k'_d L_{t, 1})} \quad (A25)$$

### A.2.1.3. Soil heat flux. The net heat flux density into the soil (G, W m⁻² ground) is calculated following Allen et al. (1998) and is given by:

$$G = \begin{cases} 0.1 Rn, & (R_{inc, sw, b} + R_{inc, sw, d}) > 0 \\ 0.5 Rn, & (R_{inc, sw, b} + R_{inc, sw, d}) \le 0 \end{cases} \quad (A26)$$

where $Rn$ is given by:

$$Rn = \sum_{i=1}^{n-1} (R_{v, sw, i} + R_{v, lw, i}) + (R_{s, sw} + R_{s, lw}) \quad (A27)$$

### A.2.2. Aerodynamic resistance

Both free and forced (turbulent) convection contribute to momentum and heat transfer, and hence to the aerodynamic resistance ($r_{a,0}$, h m⁻¹) through the boundary layer above the canopy. Free convection generally dominates under low wind speed conditions, when the value of the Richardson number ($Ri, -$) differs from unity (Monteith and Unsworth, 2013). Webber et al. (2016) suggested to consider that free convection dominates when $Ri$ is lower than a threshold value ($Ri_{free}, -$), otherwise turbulent forced convection prevails. This approach can cause numerical instability since the model iterates over the aerodynamic resistance $r_{a,0}$ (see Section 2.1.3). Therefore, to ensure numerical convergence, the transition between free and forced convection is performed in MONTPEL using a sigmoidal function:

$$r_{a,0} = \frac{1}{\delta \frac{1}{r_{a,0, free}} + (1 - \delta) \frac{1}{r_{a,0, forced}}} \quad (A28)$$

with

$$r_{a,0, forced} = \frac{1}{k \nu} \left[ \ln \left( \frac{z_r - d}{z_{0, v}} \right) - \varphi_v \right] \quad (A29)$$

$$r_{a,0, free} = \frac{\rho c_p}{\eta |T_a - T_m|^{1/3}} \quad (A30)$$

$$\delta = \frac{1}{1 + e^{(Ri - Ri_{free})}} \quad (A31)$$

where $k$ (–) is the von Kármán constant, $ν$ (m h⁻¹) is the friction velocity, $z_r$ (m) is the reference height, $d$ (m) is the zero-plane displacement height, $z_{0,v}$ (m) is the roughness length for heat transfer, $\varphi_v$ (–) is the atmospheric stability correction function for heat transfer, $η$ (W K⁻⁴ m⁻²) is an empirical surface-characteristic shape parameter.

The friction velocity is given as:

$$\nu = \frac{ku}{\left[ \ln \left( \frac{z_r - d}{z_{0, u}} \right) - \varphi_u \right]} \quad (A32)$$

where $u$ (m h⁻¹) is the wind speed, $z_{0,u}$ (m) is the roughness length for momentum, and $\varphi_u$ (–) is the atmospheric stability correction function for momentum.

$d$, $z_{0,u}$ (Choudhury and Monteith, 1988), and $z_{0,v}$ (Monteith and Unsworth, 2013) are given by:

$$d = 1.1 z_h \ln(1 + (C_d L_t)^{0.25}) \quad (A33)$$
$$z_{0, u} = \min \left\{ \begin{array}{l} z_{0, \text{soil}} + 0.3 z_h \sqrt{C_d L_t} \\ 0.3 z_h \left(1 - \frac{d}{z_h}\right) \end{array} \right. \qquad (A34)$$

$$z_{0, v} = \xi z_{0, u} \qquad (A35)$$

where $C_d$ ($m^2$) is the drag coefficient, $z_h$ (m) is the average canopy height, $z_{0,soil}$ (m) is soil roughness length for momentum, and $\xi$ (—) is the ratio of heat to momentum transfer roughness lengths.

The formula to calculate Ri depends on the atmospheric stability condition which is assumed to be stable if the canopy is cooler than the surrounding air, otherwise unstable:

$$Ri = \begin{cases} \frac{z_r - d}{\xi_{MO} + 5(z_r - d)}, & H^n < 0 \ (\text{stable}) \\ \frac{z_r - d}{\xi_{MO}}, & H^n \ge 0 \ (\text{unstable}) \end{cases} \qquad (A36)$$

where $\xi_{MO}$ (m) is the Monin-Obukhov length given by:

$$\xi_{MO} = \frac{-\nu^3}{k \frac{g}{T_a} \frac{H^n}{\rho c_p}} \qquad (A37)$$

where $g$ (m h⁻²) is the gravitational acceleration.

Finally, the stability correction functions $\varphi_u$ and $\varphi_v$ depend on $Ri$ (Liu et al., 2007; Webber et al., 2016):

If $Ri < -0.8$ (strongly unstable) (a):

$$\varphi_u = 0$$

$$\varphi_v = 0$$

else if $Ri < -0.01$ (unstable) (b):

$$\varphi_u = 2\ln\left(\frac{1+x}{2}\right) + \ln\left(\frac{1+x^2}{2}\right) - 2\tan^{-1}(x) + \frac{\pi}{2}$$

$$\varphi_v = 2\ln\left(\frac{1+x^2}{2}\right)$$

with,

$$x = \left(1 - 16 \frac{z_r - d}{\xi_{MO}}\right)^{0.25} \qquad (A38)$$

else if $Ri < 0.2$ (stable) (c):

$$\varphi_u = -5 \left( \frac{z_r - d}{\xi_{MO}} \right)$$

$$\varphi_v = \varphi_u$$

else (d):

$$\varphi_u = 0$$

$$\varphi_v = 0$$

### A.2.3. Boundary layer resistances

The boundary layer resistance of a component i ($r_{a,i}$, h m⁻¹) is the reciprocal of the boundary layer conductance ($g_{b,i}$, m h⁻¹):

$$r_{a, i} = \begin{cases} \frac{\nu_i}{g_{b, i}}, & \text{for latent energy} \\ \frac{1}{g_{b, i}}, & \text{for sensible heat} \end{cases} \qquad (A39)$$

$g_{b,i}$ is calculated as the integral of leaf-based boundary layer conductance ($g_{b,L}$, m h⁻¹) over the leaf layer thickness $L_{i,L} - L_{u,L}$:

$$g_{b, i} = \int_{L_{u,i}}^{L_{i,i}} g_{b, L} \, dL = \int_{L_{u,i}}^{L_{i,i}} \left( g_{b, \text{free, L}} + g_{b, \text{forced, L}} \right) dL \qquad (A40)$$
where g<sub>b, free, L</sub> and g<sub>b, forced, L</sub> (m h<sup>-1</sup>) are the free and forced convection conductance, respectively.

Both free and forced convection conductance are considered in this model since wind speed decreases sharply with depth within the canopy and commonly reaches values between 0.1 and 0.5 m s<sup>-1</sup> at the bottom of a closed canopy (Monteith and Unsworth, 2013). Under such conditions, both processes contribute to heat transfer (Monteith and Unsworth, 2013).

g<sub>b, forced, L</sub> is calculated as a function of the leaf characteristic length (w, m) and wind speed in the direct vicinity of the canopy (u<sub>L</sub>, m h<sup>-1</sup>) and is given by (Jones, 1992):

$$
g_{b, \text{forced, L}} = \alpha \sqrt{\frac{u_L}{w}} \qquad (\text{A41})
$$

where α (m h<sup>-1/2</sup>) is an empirical shape parameter.

Wind speed u<sub>L</sub> is assumed to decrease exponentially with depth and is given by:

$$
u_L = u_h e^{(-k_u L)} \qquad (\text{A42})
$$

where k<sub>u</sub> (m<sup>2</sup> ground m<sup>-2</sup> leaf) is the wind speed extinction coefficient through the canopy, and u<sub>h</sub> (m h<sup>-1</sup>) is the wind speed at the top of the canopy given by:

$$
u_h = u \frac{\ln\left(\frac{z_h-d}{z_0, u}\right)}{\ln\left(\frac{z_r-d}{z_0, u}\right)} \qquad (\text{A43})
$$

g<sub>b, free, L</sub> is calculated following Leuning et al. (1995) as:

$$
g_{b, \text{free, L}} = 2 \frac{0.5 D_H Gr^{0.25}}{w} \qquad (\text{A44})
$$

where D<sub>H</sub> (m<sup>2</sup> h<sup>-1</sup>) is the heat molecular diffusivity of air, T<sub>s, L</sub> (K) is leaf temperature at depth L and Gr (—) is the Grashof number given by:

$$
Gr = 1.58 \cdot 10^8 |T_{s, L} - T_a| w^3 \qquad (A45)
$$

For a Lumped leaf layer spanning between L<sub>u, 1</sub> and L<sub>l, 1</sub>, the boundary layer resistance is given by:

$$
\begin{align}
g_{b, \text{l, lumped}} &= \int_{L_{u, 1}}^{L_{l, 1}} \left( \alpha \sqrt{\frac{u_h e^{(-k_u L)}}{w}} + \frac{D_H Gr^{0.25}}{w} \right) dL \tag{A46} \\
&= \frac{2\alpha}{k_u} \sqrt{\frac{u_h}{w}} \left( e^{-\frac{k_u}{2} L_{u, 1}} - e^{-\frac{k_u}{2} L_{l, 1}} \right) + \frac{D_H Gr^{0.25}}{w} (L_{l, 1} - L_{u, 1}) \nonumber
\end{align}
$$

The boundary layer conductance of Sunlit (g<sub>b, l, shaded</sub>) and Shaded (g<sub>b, l, shaded</sub>) leaf fractions are given by:

$$
\begin{align}
g_{b, \text{l, sunlit}} &= \int_{L_{u, 1}}^{L_{l, 1}} \left( \alpha \sqrt{\frac{u_h e^{(-k_u L)}}{w}} + \frac{D_H Gr^{0.25}}{w} \right) e^{(-k'_b L)} dL \nonumber \\
&= \frac{2\alpha}{k_u + 2k'_b} \sqrt{\frac{u_h}{w}} \left( e^{\left(-\left(\frac{k_u}{2} + k'_b\right)L_{u, 1}\right)} - e^{\left(-\left(\frac{k_u}{2} + k'_b\right)L_{l, 1}\right)} \right) + \frac{D_H Gr^{0.25}}{w} \frac{\left(e^{(-k'_b L_{u, 1})} - e^{(-k'_b L_{l, 1})}\right)}{k'_b} \tag{A47}
\end{align}
$$

$$
g_{b, \text{l, shaded}} = g_{b, \text{l, lumped}} - g_{b, \text{l, sunlit}} \qquad (\text{A48})
$$

For BigLeaf canopy representation, equations Eq. A46 and Eq. A48 are used with L<sub>u</sub> = 0 and L<sub>l</sub> = L<sub>t</sub>.

For the soil component, the aerodynamic resistance between the ground surface and the canopy source height (r<sub>a, soil</sub>, h m<sup>-1</sup>) is given by
(Choudhury and Monteith, 1988):

$$
r_{a, \text{soil}} = \frac{z_h e^{a_w}}{\alpha_w K(z_h)} \left[ e^{-\left(\frac{a_w z_{0, \text{soil}}}{z_h}\right)} - e^{-\left(\frac{a_w}{z_h}\left(d+z_{0, \text{u}}\right)\right)} \right] \quad (A49)
$$

where α<sub>w</sub> (—) is an empirical shape parameter, z<sub>0,soil</sub> (m) is the soil roughness length for momentum, and K(z<sub>h</sub>) (m<sup>2</sup> h<sup>-1</sup>) is the eddy diffusivity at canopy height z<sub>h</sub> given by (Choudhury and Monteith, 1988):

$$
K(z_h) = \frac{k^2 u (z_h - d)}{\ln\left(\frac{z_r - d}{z_0, u}\right)} \qquad (A50)
$$
### A.2.4. Surface resistances

Similarly, to the boundary layer resistance, the surface resistance of a leaf layer ($r_{s,l}$, h m⁻¹) is the reciprocal of its surface conductance ($g_{s,l}$, m h⁻¹):

$$r_{s,l} = \frac{\nu_{l}}{g_{s,l}} \qquad (A51)$$

$g_{s,l}$ is upscaled from the leaf to the layer or canopy scale assuming that leaf temperature inside each leaf component is constant (Leuning et al., 1995; Kelliher et al., 1995; Wang and Leuning, 1998, Irmak et al., 2008), that is to say, solar shortwave irradiance is the main driver of the variation of surface conductance inside the canopy and is given by:

$$g_{s,L} = g_{s,\text{res}} + g_{s,\text{max}} \frac{I_{\text{PAR}, L}}{(I_{\text{PAR}, L} + I_{\text{PAR}, 50}) \left(1 + \frac{D_L}{D_0}\right) \left(1 + \frac{\psi_{\text{soil}, 50}}{\psi_{\text{soil}, 50}}\right)^\beta} \qquad (A52)$$

where $g_{s,L}$ (m h⁻¹) is the leaf surface conductance at depth L, $g_{s,\text{res}}$ and $g_{s,\text{max}}$ (m h⁻¹) are the residual and maximum stomatal conductance, respectively, $I_{\text{PAR}, L}$ (W m⁻²) is the absorbed photosynthetically active radiation (PAR) per unit leaf surface area (the absorbed PAR is assumed equal to 0.48 of the absorbed shortwave irradiance), $I_{\text{PAR}, 50}$ (W m⁻²) is $I_{\text{PAR}, L}$ at which $g_{s,\text{max}}$ is reduced by 50%, $D_L$ (kPa) is the vapor pressure deficit of the leaf surface, $D_0$ (kPa) is an sensitivity parameter to $D_a$, $\psi_{\text{soil}, 50}$ (MPa) is the soil water potential, $\psi_{\text{soil}, 50}$ (MPa) is the soil water potential at which $g_{s,\text{max}}$ is reduced by 50%, and $\beta$ (−) is an empirical shape parameter.

In the MONTPEL model the absorbed PAR is modeled following Goudriaan (1988) model for sunlit and shaded leaf fractions and the derived version of the latter model by de Pury and Farquhar (1997) for Lumped leaf layers. For clarity and reproducibility, the equations of both models are recalled below.

For a Lumped leaf layer at depth L, the absorbed PAR ($I_{\text{PAR}, \text{lumped}, L}$, W m⁻²) is given by:

$$I_{\text{PAR, lumped}, L} = I_{\text{PAR, sunlit}, L} e^{-k_b L} + R_{\text{PAR, shaded}, L} (1 - e^{-k_b L}) \qquad (A53)$$

where $\varphi = e^{-k_b L}$ denotes the sunlit leaf fraction at the depth L.

The amount of PAR absorbed by sunlit ($I_{\text{PAR, sunlit}, L}$) and shaded ($I_{\text{PAR, shaded}, L}$) leaf fractions at depth L is given by:

$$I_{\text{PAR, sunlit}, L} = 0.48 \left[ \begin{array}{l} I_{\text{inc,sw, direct}} k_{\text{direct, black}} \times (1 - \sigma_s) \\ + I_{\text{inc,sw, diffuse}} (1 - \rho_d) k_d e^{-k_d L} \\ + I_{\text{inc,sw, direct}} \left( \begin{array}{l} (1 - \rho_b) k_b e^{-k_b L} \\ -(1 - \sigma_s) - k'_b e^{-k'_b L} \end{array} \right) \end{array} \right] \qquad (A54)$$

$$I_{\text{PAR, shaded}, L} = 0.48 \left[ \begin{array}{l} I_{\text{inc,sw, d}} (1 - \rho_d) k_d e^{-k_d L} \\ + I_{\text{inc,sw, b}} \left( \begin{array}{l} (1 - \rho_b) k_b e^{-k_b L} \\ -(1 - \sigma_s) - k'_b e^{-k'_b L} \end{array} \right) \end{array} \right] \qquad (A55)$$

On the scale of a leaf layer l spanning between depths $L_{u, 1}$ and $L_{l, 1}$ and having a Lumped leaf representation, the surface resistance ($g_{s, \text{lumped}, 1}$, m h⁻¹) is given by:

$$g_{s, l, \text{lumped}} = \int_{L_{u, 1}}^{L_{l, 1}} g_{s,L} dL = \int_{L_{u, 1}}^{L_{l, 1}} g_{s,\text{res}} dL + \int_{L_{u, 1}}^{L_{l, 1}} g_{s,\text{max}} \frac{I_{\text{PAR}, L}}{(I_{\text{PAR}, L} + I_{\text{PAR}, 50}) \left(1 + \frac{D_a}{D_0}\right) \left(1 + \frac{\psi_{\text{soil}, 50}}{\psi_{\text{soil}, 50}}\right)^\beta} dL \qquad (A56)$$

Similarly, the surface conductance of sunlit ($g_{s, l, \text{sunlit}}$, m h⁻¹) and shaded ($g_{s, l, \text{shaded}}$, m h⁻¹) leaf layer fractions of a layer l are given by:

$$g_{s, \text{sunlit}, l} = \int_{L_{u, 1}}^{L_{l, 1}} g_{s,L} e^{-k_b L} dL \qquad (A57)$$

$$g_{s, \text{shaded}, l} = \int_{L_{u, 1}}^{L_{l, 1}} g_{s,L} (1 - e^{-k_b L}) dL \qquad (A58)$$

The terms $e^{-k_b L}$ and $(1 - e^{-k_b L})$ refer to the surface fractions of sunlit and shaded leaves at depth L, respectively.

The integrals for calculating $g_{s,l, \text{lumped}}$, $g_{s,l, \text{sunlit}}$, and $g_{s,l, \text{shaded}}$ are highly nonlinear on L. They were solved numerically by fixing the value of dL as $\Delta L = 0.01 \text{ m}^2 \text{leaf} \text{m}^{-2} \text{ground}$ and summing up the resulting values over the layer thickness $L_{l, 1} - L_{u, 1}$.

Finally, the surface resistance of the soil component ($r_{s, n}$, h m⁻¹) is calculated following Sellers et al. (1992) and is given by:
$$r_{s,soil} = (e^{(a_s - b_s \theta/\theta_{sat})} / 3600) \qquad (A59)$$

where $a_s$ and $b_s$ ($\text{m}^{-1}$) are empirical parameters, and $\theta$ and $\theta_{sat}$ ($\text{m}^3\text{m}^{-3}$) are the actual and saturated soil volumetric water content, respectively. A.3. Resolution scheme

Equations Eq. A1 to Eq. A59 determine the energy balance of all components jointly, with complex interactions between the temperature of individual components, $T_s$, i, and aerodynamic temperature, $T_m$.

$T_s$, i determines the Grashof number (Gr, Eq. A45) and consequently the component boundary layer conductance under free convection (Eq. A44) and its bulk boundary layer resistance (Eqs. Eq. A39 and Eq. A40). The latter determines the latent heat flux of the component (Eq. A7) and the total heat flux $H^n$ (Eq. A2). $H^n$ is also directly determined by $T_m$ (Eq. A63) which also determines $r_{a,0}$, free (Eq. A30) and hence indirectly $r_{a,0}$ (Eq. A28 and Eq. A61) and $H^n$ (Eq. A63). $T_m$ also determines, via $r_{a,0}$, both $\xi_{MO}$ (Eq. A37) and Ri (Eq. A36). Finally, Ri determines the stability correction equations (Eq. A38) which determine in turn $r_{a,0}$, forced and hence $r_{a,0}$.

The complexity of the aforementioned interactions require an iterative solution with $T_s$, i as state variables. As described in the main text (Section 2.1.3), the iterative solution includes three levels of iterations. The first (inner) level (Fig. 2c) solves the energy balance for a fixed aerodynamic resistance value, the second (outer) level (Fig. 2a) solves energy balance with varying aerodynamic resistance values that are adjusted for atmospheric stability, and the third (intermediate) level (Fig. 2d) determines intermediate stability correction variables for a fixed source temperature value.

The inner iteration uses the temperatures of individual components as state variables. For each component i and each iteration j, using the temperature from the previous iteration ($T_s$, i, j-1) the latent energy fluxes $\lambda E_j^n$ and $\lambda E_{i,j}$ are calculated. Those fluxes are used in turn to estimate the new temperature value of the current iteration ($T_s$, i, j), given by:

$$T_{s, i, j} = T_{m, j} + \frac{r_{a, i}}{\rho c_p} (A_{i, j} - \lambda E_{i, j}) \qquad (A60)$$

$$T_{m, j} = T_a + \frac{r_{a, 0}}{\rho c_p} (A_j - \lambda E_j^n) \qquad (A61)$$

where $T_a$ and $T_m$ (K) are the air temperatures at reference (measurement) height and source temperature, respectively. If the difference $T_s$, i, j-1 - $T_s$, i, j between two consecutive estimates is lower than a predefined threshold (set to 0.02 K; Maes and Steppe, 2012), then iteration stops. Otherwise, a new value of surface temperature is calculated based on both $T_s$, i, j-1 and $T_s$, i, j, which is used for the next iteration:

$$T_{s, i, j} = T_{s, i, j-1} + \frac{1}{\omega_j} (T_{s, i, j} - T_{s, i, j-1}) \qquad (A62)$$

where $\omega_j$ (-) is a parameter that accounts for convergence “difficulty” ($\omega_j = 1$ for an explicit solution, $\omega_j = 2$ for a half implicit solution, can be set by the user). Higher values of $\omega_j$ increase the stability of the numerical solution at the expense of calculation speed.

The outer iteration uses the sensible heat flux ($H^n$, W m$^{-2}$ground) of the whole canopy as the state variable (Webber et al., 2016):

$$H^n = \frac{\rho c_p}{r_{a, 0}} (T_m - T_a) \qquad (A63)$$

For each iteration j, $r_{a, 0}$, j is calculated using the value of the source temperature of the previous iteration ($T_m$, j-1) and the corresponding sensible heat flux ($H_j^n$) is updated. The iteration continues until the difference between two consecutive estimations of sensible heat ($H_j^n - H_{j-1}^n$) drops below a predefined threshold (set to 0.01W m$^{-2}$ground).

Finally, in the intermediate iteration, both temperature and aerodynamic resistance values are fixed and the stability correction functions for momentum and heat (respectively $\varphi_u$ and $\varphi_v$) are used as state variables. For each iteration j, $\nu$, $\xi_{MO}$, and Ri are updated. The iteration stops when both values of $\varphi_u$ and $\varphi_v$ stabilize over two consecutive iterations (the acceptable error is set to 0.01). To improve numerical stability, a minimum value of 0.1 $\nu$ ($\text{m} \text{h}^{-1}$) is imposed following Webber et al. (2016).

## References

Allen, R., Pereira, L., Raes, D., Smith, M., 1998. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements. FAO Irrigation and Drainage Paper no. 56, Rome, Italy, p. 300.

Alves, I., Perrier, A., Pereira, L., 1998. Aerodynamic and surface resistances of complete cover crops: how good is the “Big Leaf”? Trans. Am. Soc. Agric. Eng. 41, 345–351.

Arora, V., 2003. Simulating energy and carbon fluxes over winter wheat using coupled land surface and terrestrial ecosystem models. Agric. For. Meteorol. 118, 21–47.

Bao, Y., Liu, T., Duan, L., Tong, X., Ji, H., Zhang, L., Singh, V., 2022. A comparative study of three stomatal conductance models for estimating evapotranspiration in a dune ecosystem in a semi-arid region. Sci. Total Environ. 802, 149937.

Beaudoin, N., Lecharpentier, P., Ripoche-Wachter, D., Strullu, L., Mary, B., Léonard, J., Launay, M., Justes, E. (Eds.), 2023. STICS Soil-Crop Model. Conceptual Framework, Equations and Uses. Editions Quæ, Versailles, France, p. 516.

Bernard, F., Chelle, M., Fortineau, A., Riahi El Kamel, O., Pincebourde, S., Sache, I., Suffert, F., 2022. Daily fluctuations in leaf temperature modulate the development of a foliar pathogen. Agric. For. Meteorol. 322, 109031.

Blonquist, J., Allen, R., Bugbee, B., 2010. An evaluation of the net radiation sub-model in the ASCE standardized reference evapotranspiration equation: implications for evapotranspiration prediction. Agric. Water. Manage. 97, 1026–1038.

Boote, K., Jones, J., White, J., Asseng, S., Lizaro, J., 2013. Putting mechanisms into crop production models. Plant, Cell Environ. 36, 1658–1672.

Brisson, N., Launay, M., Mary, B., Beaudoin, N., 2008. Conceptual Basis, Formalisations and Parameterisation of the Stics Crop Model. Editions QUAЕ, Versailles, France, p. 297.

Brunt, D., 1932. Notes on radiation in the atmosphere: I. Q. J. R. Meteorol. Soc. 58, 389–420.

Brutsaert, W., 1975. On a derivable formula for long-wave radiation from clear skies. Water. Resour. Res. 11, 742–744.

Cammalleri, C., Anderson, M., Ciraolo, G., D’Urso, G., Kustas, W., la Loggia, G., Minacapilli, M., 2012. Applications of a remote sensing-based two-source energy balance algorithm for mapping surface fluxes without in situ air temperature observations. Remote Sens. Environ. 124, 502–515.

Campbell, G., 1986. Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal inclination angle distribution. Agric. For. Meteorol. 36, 317–321.
Campbell, G., 1990. Derivation of an angle density function for canopies with ellipsoidal leaf angle distributions. Agric. For. Meteorol. 49, 173-176.

Chelle, M., 2005. Phylloxyle or the climate perceived by individual plant organs: what is it? How to model it? What for? New Phytol. 166, 781-790.

Chelle, M., Cellier, P., 2009. Horizontal variability in air temperature over time within a maize inter-row. Agric. For. Meteorol. 149, 1294-1300.

Chen, N., Li, X., Shi, H., Yan, J., Hu, Q., Zhang, Y., 2021. Modeling maize evapotranspiration and associated processes under biodegradable film mulching in an arid dripped field. Agric. For. Meteorol. 297, 108247.

Chen, X., Yu, L., Cui, N., Cai, H., Jiang, X., Liu, C., Shu, Z., Wu, Z., 2022. Modeling maize evapotranspiration using three types of canopy resistance models coupled with single-source and dual-source hypotheses—a comparative study in a semi-humid and drought-prone region. J. Hydrol. (Amst) 614, 128638.

Chenu, K., Porter, J., Martre, P., Basso, B., Chapman, S., Ewert, F., Bindi, M., Asseng, S., 2017. Contribution of crop models to adaptation in wheat. Trends. Plant Sci. 22, 472-490.

Choudhury, B., Monteith, J., 1988. A four-layer model for the heat budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.

Colaizzi, P., Kustas, W., Anderson, M., Agam, N., Tolk, J., Evett, S., Howell, T., Gowda, P., O'Shaughnessy, S., 2012. Two-source energy balance model estimates of evapotranspiration using component and composite surface temperatures. Adv. Water. Resour. 50, 134-151.

Cuadra, S., Kimball, B., Boote, K., Suyker, A., Pickering, N., 2021. Energy balance in the DSSAT-CSM-CROPGRO model. Agric. For. Meteorol. 297, 108241.

Dai, Y., Dickinson, R., Wang, Y., 2004. A two-big-leaf model for canopy temperature, photosynthesis, and stomatal conductance. J. Clim. 17, 2281-2299.

de Pury, D., Farquhar, G., 1997. Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models. Plant Cell Environ. 20, 537-557.

Ding, R., Kang, S., Du, T., Hao, X., Zhang, Y., 2014. Scaling up stomatal conductance from leaf to canopy using a dual-leaf model for estimating crop evapotranspiration. PLoS. One 9 (4), e95584.

Furon, A., Warland, J., Wagner-Riddle, C., 2007. Analysis of scaling-up resistances from leaf to canopy using numerical simulations. Agron. J. 99, 1483-1491.

Gauch, H., Gene Hwang, J., Fick, G., 2003. Model evaluation by comparison of model-based predictions and measured values. Agron. J. 95, 1442-1446.

Goudriaan, J., 1988. The bare bones of leaf-angle distribution in radiation models for canopy photosynthesis and energy exchange. Agric. For. Meteorol. 43, 155-169.

Grant, R., Kimball, B., Conley, M., White, J., Wall, G., Ottman, M., 2011. Controlled warming effects on wheat growth and yield: field measurements and modeling. Agron. J. 103, 1742-1754.

Guarin, J., Asseng, A., 2022. Improving wheat production and breeding strategies using crop models. In: Reynolds, M., Braun, H. (Eds.), Wheat Improvement: Food Security in a Changing Climate. Springer Nature Switzerland AG, pp. 573-591.

Holzworth, D., Huth, N., deVoil, P., Zurcher, E., Herrmann, N., McLean, G., Chenu, K., van Oosterom, E., Snow, V., Murphy, C., Moore, A., Brown, H., Whish, J., Verrall, S., Fainges, J., Bell, L., Peake, A., Poulton, P., Hochman, Z., Keating, B., 2014. APSIM – evolution towards a new generation of agricultural systems simulation. Environ. Model. Softw. 62, 327-350.

Idso, S., 1981. A set of equations for full spectrum and 8- to 14-µm and 10.5- to 12.5-µm thermal radiation from cloudless skies. Water. Resour. Res. 17, 295-304.

Inclán, M., Forkel, R., 1995. Comparison of energy fluxes calculated with the Penman-Monteith equation and the vegetation models SiB and Cupid. J. Hydrol. (Amst) 166, 193-211.

Irmaç, S., Mutibwa, D., Irmaç, A., Arkebauer, T., Weiss, A., Martin, D., Eisenhauer, D., 2008. On the scaling up leaf stomatal resistance to canopy resistance using photosynthetic photon flux density. Agric. For. Meteorol. 148, 1034-1044.

Jones, H., 1992. Plants and Microclimate: A Quantitative Approach to Environmental Plant Physiology, Second edition. Cambridge University Press, Cambridge, UK, p. 428.

Keating, B., Carberry, P., Hammer, G., Probert, M., Robertson, M., Holzworth, D., Huth, N., Hargreaves, J., Meinke, H., Hochman, Z., 2003. An overview of APSIM, a model designed for farming systems simulation. Eur. J. Agronomy 18, 267-288.

Kelliher, F., Leuning, R., Raupach, M., Schulze, E., 1995. Maximum conductances for evaporation from global vegetation types. Agric. For. Meteorol. 73, 1-16.

Kimball, B., White, J., Ottman, M., Wall, G., Bernacchi, C., Morgan, J., Smith, D., 2015. Predicting canopy temperatures and infrared heater energy requirements for warming field plots. Agron. J. 107, 121-141.

Kimball, B., Pinter, P., Lamorte, R., Leavitt, S., Hunsaker, D., Wall, G., Wechsung, F., Wechsung, G., Bloom, A., White, J., 2017. Data from the Arizona FACE (Free-Air CO₂ Enrichment) experiments on wheat at ample and limiting levels of water and nitrogen. Open Data J. Agric. Res. 3, 29-38.

Kjelgaard, J., Stockle, C., Evans, R., 1996. Accuracy of canopy temperature energy balance for determining daily evapotranspiration. Irrig. Sci. 16, 149-157.

Kustas, W., Humes, K., Norman, J., Moran, M., 1996. Single- and dual-source modeling of surface energy fluxes with radiometric surface temperature. J. Appl. Meteorol. Climatol. 35, 110-121.

Kustas, W., Norman, J., 1999. Evaluation of soil and vegetation heat flux predictions using a simple two-source model with radiometric temperatures for partial canopy cover. Agric. For. Meteorol. 94, 13-29.

Kuusk, A., 1991. The hot spot effect in plant canopy reflectance. In: Myneni, R., Ross, J. (Eds.), Photon-Vegetation Interactions. Springer, Berlin, Heidelberg, pp. 139-159.

Lagos, L., Martin, D., Verma, S., Irmaç, S., Irmaç, A., Eisenhauer, D., Suyker, A., 2013. Surface energy balance model of transpiration from variable canopy cover and evaporation from residue-covered or bare soil systems: model evaluation. Irrig. Sci. 31, 135-150.

Lawrence, D., Fisher, R., Koven, C., Olesen, K., Swenson, S., Bonan, G., Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., Kluzek, E., Lawrence, P., Li, F., Li, H., Lombardozzi, D., Riley, W., Sacks, W., Shi, M., Vertenstein, M., Zeng, X., 2019. The community land model version 5: description of new features, benchmarking, and impact of forcing uncertainty. J. Adv. Model. Earth. Syst. 11, 4245-4287.

Leuning, R., 1995. A critical appraisal of a combined stomatal photosynthesis model for C3 plants. Plant, Cell Environ. 18, 339-355.

Leuning, R., Kelliher, F., de Pury, D., Schulze, E., 1995. Leaf nitrogen, photosynthesis, conductance and transpiration: scaling from leaves to canopies. Plant, Cell Environ. 18, 1183-1200.

Leuning, R., Dunin, F., Wang, Y., 1998. A two-leaf model for canopy conductance, photosynthesis and partitioning of available energy. II. Comparison with measurements. Agric. For. Meteorol. 91, 113-125.

Lhomme, J., Montes, C., Jacob, F., Prévot, L., 2013. Evaporation from multi-component canopies: generalized formulations. J. Hydrol. (Amst) 486, 315-320.

Liu, S., Lu, L., Mao, D., Jia, L., 2007. Evaluating parameterizations of aerodynamic resistance to heat transfer using field measurements. Hydrol. Earth. Syst. Sci. 11, 769-783.

Liu, S., Baret, F., Abichou, M., Manceau, L., Andrieu, B., Weiss, M., Martre, P., 2021. Importance of the description of light interception in crop growth models. Plant Physiol. 186, 977-997.

Luo, X., Chen, J., Liu, J., Black, T., Croft, H., Staebler, R., He, L., Arain, M., Chen, B., Mo, G., Gonsamo, A., McCaughey, H., 2018. Comparison of Big-Leaf, Two-Big-Leaf, and Two-Leaf upscaling schemes for evapotranspiration estimation using coupled carbon-water modeling. J. Geophys. Res. 123, 207-225.

Maes, W., Steppe, K., 2012. Estimating evapotranspiration and drought stress with ground-based thermal remote sensing in agriculture: a review. J. Exp. Bot. 63, 4671-4712.

Maiorano, A., Martre, P., Asseng, S., Ewert, F., Müller, C., Rötter, R., Ruane, A., Semenov, M., Wallach, D., Wang, E., Alderman, P., Kassie, B., Biernath, C., Basso, B., Cammarano, D., Challinor, A., Doltra, J., Dumont, B., Rezaei, E., Zhu, Y., 2017. Crop model improvement reduces the uncertainty of the response to temperature of multi-model ensembles. Field. Crops. Res. 202, 5-20.

Manderscheid, R., Dier, M., Erbs, M., Hüther, L., Köhler, P., Luig, A., Oldenburg, E., Sickora, J., Weigel, H., Zörb, C., 2020. Experimental data from the Braunschweig FACE studies on wheat growth responses to elevated CO₂ in combination with nitrogen supply or infrared warming during grain filling. Open Data J. Agric. Res. 6, 28-33.

Martre, P., Jamieson, P., Semenov, M., Zyskowski, R., Porter, J., Triboï, E., 2006. Modelling protein content and composition in relation to crop nitrogen dynamics for wheat. Eur. J. Agron. 25, 138-154.

Martre, P., Kimball, B., Ottman, M., Wall, G., White, J., Asseng, S., Ewert, F., Cammarano, D., Maiorano, A., Aggarwal, P., Anothai, J., Basso, B., Biernath, C., Challinor, A., de Sanctis, G., Doltra, J., Dumont, B., Fereres, E., García-Vila, M., Zhu, Y., 2018. The Hot Serial Cereal Experiment for modeling wheat response to temperature: field experiments and AgMIP-Wheat multi-model simulations. Open Data J. Agric. Res. 4, 28-34.

Monteith, J., Unsworth, M., 2013. Principles of Environmental Physics. Plants, Animals, and the Atmosphere, Fourth edition. Elsevier Academic Press, p. 401.

Moreau, D., Allard, V., Gaju, O., Le Gouis, J., Foulkes, M.J., Martre, P., 2012. Acclimation of leaf nitrogen to vertical light gradient at anthesis in wheat is a whole-plant process that scales with the size of the canopy. Plant Physiol. 160, 1479-1490.

Norman, J., Kustas, W., Humes, K., 1995. Two source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature. Agric. For. Meteorol. 77, 263-293.

Ottman, M., Kimball, B., White, J., Wall, G., 2012. Wheat growth response to increased temperature from varied planting dates and supplemental infrared heating. Agron. J. 104, 7-16.

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825-2830.

Penman, H., 1948. Natural evaporation from open water, bare soil and grass. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 193, 120-145.

Pinter, P., Kimball, B., Wall, G., LaMorte, R., Hunsaker, D., Adamsen, F., Frumau, K., Vugs, H., Hendrey, G., Lewin, K., Nagy, J., Johnson, H., Wechsung, F., Leavitt, S., Thompson, T., Matthias, A., Brooks, T., 2000. Free-air CO₂ enrichment (FACE): blower effects on wheat canopy microclimate and plant development. Agric. For. Meteorol. 103, 319-333.

Raupach, M., 1995. Vegetation-atmosphere interaction and surface conductance at leaf, canopy and regional scales. Agric. For. Meteorol. 73, 151-179.

Rezaei, E., Webber, H., Gaiser, T., Naab, J., Ewert, F., 2015. Heat stress in cereals: mechanisms and modelling. Eur. J. Agron. 64, 98-113.

Rezaei, E., Webber, H., Asseng, S., Boote, K., Durand, J., Ewert, F., Martre, P., MacCarthy, D., 2023. Climate change impacts on crop yields. Nat. Rev. Earth Environ. 4, 831-846.

Rochette, P., Pattey, E., Desjardins, R., Dwyer, L., Stewart, D., Dubé, P., 1991. Estimation of maize (Zea mays L.) canopy conductance by scaling up leaf stomatal conductance. Agric. For. Meteorol. 54, 241-261.

Saltelli, A., Tarantola, S., Chan, K., 1999. A quantitative, model independent method for global sensitivity analysis of model output. Technometrics. 41, 39-56.

Sellers, P., Heiser, M., Hall, F., 1992. Relations between surface conductance and spectral vegetation indices at intermediate (100 m² to 15 km²) length scales. J. Geophys. Res. 97, 19033-19059.

Shuttleworth, W., 2007. Putting the "vap" into evaporation. Hydrol. Earth. Syst. Sci. 11, 210-244.
Shuttleworth, W., Wallace, J., 1985. Evaporation from sparse crops – an energy combination theory. Q. J. R. Meteorol. Soc. 111, 839–855.

Song, L., Liu, S., Kustas, W., Zhou, J., Xu, Z., Xia, T., Li, M., 2016. Application of remote sensing-based two-source energy balance model for mapping field surface fluxes with composite and component surface temperatures. Agric. For. Meteorol. 230–231, 8–19.

Stöckle, C., Donatelli, M., Nelson, R., 2003. CropSyst, a cropping systems simulation model. Eur. J. Agron. 18, 289–307.

Tolk, J., Howell, T., Steiner, J., Krieg, D., 1995. Aerodynamic characteristics of corn as determined by energy balance techniques. Agron. J. 87, 464–473.

van Zyl, W., de Jager, J., 1987. Accuracy of the Penman-Monteith equation adjusted for atmospheric stability. Agric. For. Meteorol. 41, 57–64.

Wall, G., Kimball, B., White, J., Ottman, M., 2011. Gas exchange and water relations of spring wheat under full-season infrared warming. Glob. Chang. Biol. 17, 2113–2133.

Wang, Y., Leuning, R., 1998. A two-leaf model for canopy conductance, photosynthesis and partitioning of available energy I: model description and comparison with a multi-layered model. Agric. For. Meteorol. 91, 89–111.

Wang, E., Brown, H., Rebetzke, G., Zhao, Z., Zheng, B., Chapman, S., 2019. Improving process-based crop models to better capture genotype × environment × management interactions. J. Exp. Bot. 70, 2389–2401.

Webber, H., Ewert, F., Kimball, B., Siebert, S., White, J., Wall, G., Gaiser, T., 2016. Simulating canopy temperature for modelling heat stress in cereals. Environ. Model. Softw. 77, 143–155.

Webber, H., Martre, P., Asseng, S., Kimball, B., White, J., Ottman, M., Wall, G., de Sanctis, G., Doltra, J., Grant, R., Kassie, B., Maiorano, A., Olesen, J., Ripoche, D.,

Rezaei, E., Semenov, M., Stratonovitch, P., Ewert, F., 2017. Canopy temperature for simulation of heat stress in irrigated wheat in a semi-arid environment: a multi-model comparison. Field. Crops. Res. 202, 21–35.

Webber, H., White, J., Kimball, B., Ewert, F., Asseng, S., Eyshi Rezaei, E., Pinter, P., Hatfield, J., Reynolds, M., Ababaei, B., Bindi, M., Doltra, J., Ferrise, R., Kage, H., Kassie, B., Kersebaum, K., Luig, A., Olesen, J., Semenov, M., Stratonovitch, P., Martre, P., 2018. Physical robustness of canopy temperature models for crop heat stress simulation across environments and production conditions. Field. Crops. Res. 216, 75–88.

Wei, G., Zhao, Z., Tian, W., Qiu, Z., Yang, Z., Wang, F., 2023. Two alternatives to the two-source energy balance evapotranspiration model. Water. Resour. Res. 59, e2022WR032958.

Williams, J., Jones, C., Kiniry, J., Spanel, D., 1989. The EPIC crop growth model. Trans. ASAE. Am. Soc. Agric. Eng. 32, 497–511.

Xue, K., Song, L., Xu, Y., Liu, S., Zhao, G., Tao, S., Magliulo, E., Manco, A., Liddell, M., Wohlfahrt, G., Varlagin, A., Montagnani, L., Woodgate, W., Loubet, B., Zhao, L., 2023. Estimating ecosystem evaporation and transpiration using a soil moisture coupled two-source energy balance model across FLUXNET sites. Agric. For. Meteorol. 337, 109513.

Yin, X., van Laar, H., 2005. Crop Systems Dynamics: An Ecophysiological Simulation Model For Genotype-By-Environment Interactions. Wageningen Academic Publishers, Wageningen, The Netherlands, p. 155.

Zhang, B., Liu, Y., Xu, D., Cai, J., Li, F., 2011. Evapotranspiration estimation based on scaling up from leaf stomatal conductance to canopy conductance. Agric. For. Meteorol. 151, 1086–1095.