Accounting for Turbulence-Induced Canopy Heat Transfer in the Simulation of Sensible Heat Flux in SEBS Model

Abstract

Surface turbulent heat fluxes are crucial for monitoring drought, heat waves, urban heat islands, agricultural water management, and other hydrological applications. Energy Balance Models (EBMs) are widely used to simulate surface heat fluxes from a combination of remote sensing-derived variables and meteorological data. Single-source EBMs, in particular, are preferred in mapping surface turbulent heat fluxes due to their relative simplicity. However, most single-source EBMs suffer from uncertainties inherent to the parameter $k{B}^{-1}$, which is used to account for differences in the source of heat and the sink of momentum when representing aerodynamic resistance in single-source EBMs. For instance, the parameterization of $k{B}^{-1}$ in the commonly used single-source Surface Energy Balance System (SEBS) model uses a constant value of the foliage heat transfer coefficient ($C_t$), in the parameterization of the vegetation component of $k{B}^{-1}$ ($k{B}_v{}^{-1}$). Thus, SEBS ignores the effect of turbulence on canopy heat transfer. As a result, SEBS has been found to greatly underestimate sensible heat flux in tall forest canopies, where turbulence is a key contributor to canopy heat transfer. This study presents a revised parameterization of $k{B}_v{}^{-1}$ for the SEBS model. A physically based formulation of $C_t$, which considers the effect of turbulence on $C_t$, is used in deriving the revised parameterization. Simulation results across 15 eddy covariance (EC) flux tower sites show that the revised parameterization significantly reduces the underestimation of sensible heat flux compared to the original parameterization under tall forest canopies. The revised parameterization is relatively simple and does not require additional information on canopy structure compared to some more complex parameterizations proposed in the literature. As such, the revised parameterization is suitable for mapping surface turbulent heat fluxes, especially under tall forest canopies.

1. Introduction

Quantifying surface turbulent heat fluxes is crucial in the determination of evaporation and transpiration, drought monitoring, heat wave prediction, agricultural water management, and monitoring urban heat islands, among other applications. Direct observations of turbulent heat fluxes are highly scarce and are only limited to point scale hindering their utility in the applications above. Energy balance models (EBMs), such as SEBS [1], SEBAL [2], METRIC [3], and TSEBS [4,5], have risen to prominence in the mapping of surface heat fluxes. Their popularity is mainly due to their ability to capture spatial variability in fluxes related to variations in topography, land use, and land cover, since they rely on spatially varying remote sensing-derived inputs [6]. However, especially under heterogeneous landscapes, the accuracy of surface turbulent heat fluxes simulated by most EBMs mainly depends on the parameterization of the aerodynamic resistance of the underlying surface [7,8].

EBMs are broadly classified into either single or dual source models depending on how the aerodynamic resistance is parameterized [6]. Dual source models such as TSEBS separate aerodynamic resistance into two components: one resistance emanating from the bare surfaces below the canopy and the other emanating from the vegetation components of the canopy. On the other hand, in single-source models, e.g., SEBS, SEBAL, and METRIC, the bare-soil and vegetation resistances are lumped into a single bulk aerodynamic resistance. Using a single bulk aerodynamic resistance necessitates using the parameter $k{B}^{-1}$ proposed by Owen and Thomson [9], to account for the differences in the source of heat and sink of momentum in the canopy [10]. The $k{B}^{-1}$ parameter has been subject to extensive research and has been found to be highly variable both in space and time [11,12,13,14,15]. Besides, using surface radiometric temperature instead of surface aerodynamic temperature, as applied in single-source EBMs, introduces a radiometric component to $k{B}^{-1}$, further exacerbating the uncertainties associated with the parameter [15]. Although dual-source EBMs are more physically consistent with interactions between the land surface and the atmosphere compared to single-source EBMs, they are complex and require extra ancillary data [7] and calibration [16] compared to single-source EBMs. Thus, despite their relatively lower accuracy, single-source EBMs are often preferred over dual-source EBMs in operational applications.

SEBS is one of the most widely used single-source EBMs for simulating surface turbulent heat fluxes. Although SEBS has been shown to simulate heat fluxes at a reasonable accuracy [1], some studies have reported an underestimation of sensible heat flux and, subsequently, overestimation of latent heat flux with the model [10,16,17,18,19,20,21]. The underestimation of sensible heat flux with SEBS is mainly attributed to the parameterization of $k{B}^{-1}$ in the model. Most of the modifications to the parameterization of $k{B}^{-1}$ in SEBS proposed by various studies [18,20,21] are meant to alleviate the underestimation of sensible heat flux in arid and semi-arid areas. However, the study by Chen et al. [19] indicates that the model also significantly underestimates sensible heat flux in tall forest canopies. Chen et al. [19] attributed the behavior of SEBS in tall forest canopies to the use of a constant canopy heat transfer coefficient ($C_t$) in the parameterization of the vegetation-related component of $k{B}^{-1}$ ($k{B}_v{}^{-1}$). In their study, they propose a variable $C_t$ as opposed to the constant value of 0.01 adopted in the original parameterization of $k{B}_v{}^{-1}$ in SEBS. They showed that their parameterization significantly reduced the underestimation of sensible heat flux observed under tall vegetation canopies using the SEBS model. The parameterization proposed by Chen et al. [19], however, is rather complex and requires extra land cover-specific parameters and information on leaf area density profiles, which can only be obtained from lidar sensors, limiting its applicability for operational mapping of surface heat fluxes.

In this study, an improved parameterization of the parameter $k{B}_v{}^{-1}$ in SEBS suitable for use in tall forest canopies is derived. The proposed $k{B}_v{}^{-1}$ parameterization does not assume a constant value of $C_t$, but rather takes into account the influence of turbulence on canopy heat transfer through a more physically based formulation of $C_t$ proposed by Brutsaert [22]. A physically sound representation of canopy heat transfer is crucial to the ability of SEBS to simulate reliable surface turbulent heat fluxes, especially under tall forest canopies.

2. Materials and Methods

2.1. Data

This study used Eddy Covariance (EC) data to derive the required input data for SEBS and the validation dataset. The dataset was accessed from the Protocol for the Analysis of Land Surface Models (PALS) Land Surface Model Benchmarking Evaluation Project (PLUMBER) data repository (https://dapds00.nci.org.au/thredds/catalog/ks32/CLEX_Data/PLUMBER2/v1-0/catalog.html, accessed on 27 July 2022). PLUMBER is a land surface and ecosystem models intercomparison and evaluation project across 170 flux tower sites across the globe [23,24]. The PLUMBER dataset offers a high-quality filtered and gap-filled EC dataset from the original La Thuile, FLUXNET2015 [25], and OzFlux [26] EC datasets. Acquisitions of EC variables are usually done at high temporal resolution and averaged to half-hourly or hourly fluxes. The half-hourly averaged variables were used in this study. Besides fluxes of energy (sensible and latent heat fluxes), scalars (concentrations of CO2, CH4, water vapor), and various meteorological variables, the PLUMBER dataset includes leaf area index (LAI) products derived from Moderate Resolution Imaging Spectroradiometer (MODIS) and PROBA sensors for each of the sites. The MODIS-derived LAI product was used in this study. The MODIS LAI provided with the PLUMBER dataset is the 8-day temporal and 500 m spatial resolution MCD15A2H product. According to Ukkola [24], additional processing has been carried out on the LAI products provided with the PLUMBER EC flux dataset. The additional processing mainly involves gap filling and temporal smoothing to match the spatial-temporal variations in vegetation dynamics at each EC site. Besides, the gap-filled smoothened time series LAI data have been interpolated to match the 30-min temporal resolution of the EC flux dataset for seamless ingestion into climate models. A detailed description of the additional processing carried out on the PLUMBER dataset and the algorithms employed for each respective component of the PLUMBER dataset is presented in Ukkola [24].

SEBS requires inputs of remote sensing-derived surface variables and meteorological data. Meteorological data can either be from observations or simulations by weather and climate models. The various inputs necessary for simulating surface heat fluxes using the SEBS model are shown in

Table 1. Required inputs for the SEBS model.

Variable	Meteorological	Remote Sensing
Land surface temperature (LST)		yes
Leaf Area Index (LAI)		yes
Fraction vegetation cover		yes
Land Surface albedo		yes
Land surface emissivity		yes
Vegetation roughness/Canopy height		yes
Down-welling solar radiation	yes	yes
Down-welling longwave radiation	yes	yes
Air temperature	yes
Air pressure	yes
Humidity	yes
Wind speed	yes

.

The meteorological variables listed in Table 1 were obtained from the flux tower measurements. The temporal resolution of the meteorological inputs used in this study was 30 min. In addition to the down-welling solar and longwave radiation, the up-welling solar and longwave radiation measured at the flux tower sites were used to calculate the radiation budget. Since the up-welling solar radiation was available, land surface albedo was not used as input into the model. The measured canopy height (

Table 2. Summarized description of the eddy covariance (EC) station sites used in this study.

Country	Location	Station Code	Longitude	Latitude	Vegetation	Canopy Height [m]	Study Period
Germany	Hainich	DE-Hai	10.4522	51.0792	DBF	33.0	[2007, 2008]
Denmark	Soroe	DK-Sor	11.6446	55.4859	DBF	25.0	[2011, 2012]
USA	Morgan Monroe State Forest	US-MMS	−86.4131	39.3232	DBF	27.0	[2009, 2010]
Australia	Cumberland Plains	AU-Cum	150.7236	−33.6152	EBF	23.0	[2016, 2017]
Australia	Robson Creek	AU-Rob	145.6301	−17.1175	EBF	44.0	[2015, 2016]
France	Puechabon	FR-Pue	3.5957	43.7413	EBF	6.5	[2010, 2011]
Netherlands	Loobos	NL-Loo	5.7436	52.1666	ENF	15.5	[2009, 2010]
France	Le Bray	FR-LBr	−0.7693	44.7171	ENF	20.0	[2005, 2006]
Germany	Oberbärenburg	DE-Obe	13.7213	50.7867	ENF	19.0	[2012, 2013]
USA	Walnut Gulch Kendall Grasslands	US-Wkg	−109.9419	31.7365	Grassland	1.0	[2011, 2012]
Australia	Emerald	AU-Emr	148.4746	−23.8587	Grassland	2.0	[2012, 2013]
Australia	Riggs Creek	AU-Rig	145.5759	−36.6499	Grassland	0.4	[2015, 2016]
Australia	Alice Springs	AU-ASM	133.2490	−22.2830	Savannah	6.5	[2016, 2017]
Australia	Calperum	AU-Cpr	140.5891	−34.0021	Savannah	4.0	[2013, 2014]
Australia	Great Western Woodlands	AU-GWW	120.6541	−30.1913	Savannah	18.0	[2014, 2015]

) at each site was used as input into the model to describe the vegetation roughness.

A total of 15 stations, each with a continuous 2-year time series, were used in this study. The stations are spread across five major vegetation types, namely, evergreen broad-leaf forests (EBF), deciduous broad-leaf forests (DBF), evergreen needle-leaf forests (ENF), savannah, and grasslands. Three stations were selected for each vegetation type. The selected stations were the ones where the canopy height is not expected to vary significantly across the study period, since the PLUMBER dataset provides a constant value for the canopy height throughout the year for most of the sites.

The remaining remote sensing variables were derived directly or indirectly from LAI, as shown in Equations (1)–(3). Besides being used in the computation of the surface radiation budget, the measured longwave radiation components were used to derive high temporal resolution land surface temperature (LST), which is a key input in SEBS, by inverting the Stefan–Boltzmann formula (Equation (1)):

$${\mathrm{Lw}}_{\mathrm{u}}=\left(1-\mathsf{\epsilon }\right){\mathrm{Lw}}_{\mathrm{d}}+\left({\mathsf{\epsilon }\mathsf{\sigma }\mathrm{LST}}^4\right)$$

(1)

where $L{w}_u$ is the up-welling longwave radiation, $L{w}_d$ is the down-welling longwave radiation, $\epsilon$ is the land surface emissivity (Equation (2)) which is derived based on the method of Sobrino, Jiménez, and Paolini [27], and $\sigma$ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴).

$$\mathsf{\epsilon }={\mathsf{\epsilon }}_{\mathrm{v}}{f}_c+{ \mathsf{\epsilon }}_{\mathrm{s}}{f}_s+\left(1-{\mathsf{\epsilon }}_{\mathrm{s}}\right)f_s{\mathrm{F}\mathsf{\epsilon }}_{\mathrm{v}}$$

(2)

where ${\epsilon }_v$ is the emissivity of vegetation (taken as 0.99), ${\epsilon }_s$ is the bare-soil emissivity (taken as 0.96), $F$ is the shape factor, and $f_c$ is the fraction of vegetation cover, which is calculated from LAI using Equation (3), and $f_s$ is its complement.

$$f_c=1-\exp \left(-0.5LAI\right)$$

(3)

The energy balance closure corrected version of the sensible heat flux observed at EC sites was used to validate model estimates of sensible heat flux since most of the EC tower sites are often characterized by a lack of energy balance closure [25,26,28]. The evaluation metrics considered for the validation are the root mean square error (RMSE), coefficient of determination (R²), and bias.

2.2. SEBS Model Description

Similar to other EBMs, SEBS solves for the energy balance at the surface based on Equation (4). The solution to the energy balance equation in SEBS is bounded within dry and wet limiting boundary conditions, as described in Su [1]. The main outputs of SEBS are the sensible heat flux, latent heat flux, and evaporative fraction.

$$R_n=H+LE+G_0$$

(4)

where $R_n$ is the net radiation (W m⁻²), $H$ is the sensible heat flux (W m⁻²), $LE$ is the latent heat flux (W m⁻²), and $G_0$ is the ground heat flux (W m⁻²).

The most critical step in SEBS is the determination of $H$. In SEBS, $H$ is calculated based on the Monin–Obukhov similarity theory (MOST), as formulated in Equations (5) and (6). MOST applies surface-layer approximations to provide a robust scaling of fluxes between the surface and the atmosphere [29,30]:

$$H=\frac{\rho c_{p}u{k}^2\left({\theta }_s-{\theta }_a\right)}{\left[ln\left(\frac{z-d}{z_{0m}}\right)-{\Psi }_m\left(\frac{z-d}{L}\right)+{\mathsf{\Psi }}_m\left(\frac{z_{0m}}{L}\right)\right]\left[ln\left(\frac{z-d}{z_{0h}}\right)-{\Psi }_h\left(\frac{z-d}{L}\right)+{\mathsf{\Psi }}_h\left(\frac{z_{0h}}{L}\right)\right]}$$

(5)

where $\rho$ is the density of air (kg m⁻³), $c_p$ is the specific heat capacity of air (J kg⁻¹ K⁻¹), u is the wind speed (m s⁻¹), k is the von Kármán constant (0.4), ${\theta }_s$ is the potential temperature at the surface (K), ${\theta }_a$ is the potential air temperature (K), z is the reference height (m), d is the displacement height (m), $z_{0m}$ is the roughness length for momentum transfer (m), $z_{0h}$ is the roughness length for heat transfer (m), ${\Psi }_m$ and ${\Psi }_h$ are the stability correction functions for momentum and heat, respectively, and L is the Obukhov length (m), see Equation (6):

$$L=-\frac{\rho c_p{u}_{*}{}^3{\theta }_a}{kgH}$$

(6)

where $u_{*}$ is the frictional velocity (m s⁻¹) and g is the acceleration due to gravity (m s⁻²).

In SEBS, $z_{0m}$ and $d$, two critical parameters that describe the influence of canopy structure on $H$ when using the MOST approach, are calculated using Equations (7) and (8):

$$d=h_c\left(1-\frac{1}{2n_{ec}}\left[1-exp\left(-2n_{ec}\right)\right]\right)$$

(7)

where $h_c$ is the canopy height (m) and $n_{ec}$ is the within-canopy wind profile extinction coefficient, which is formulated as $n_{ec}=\frac{\left(C_{d}LAI\right)}{\left(2\left[\raisebox{1ex}{$u_{*}{}^2$}\!\left/ \!\raisebox{-1ex}{$u{\left(h\right)}^2$}\right.\right]\right)}$, with $C_d$ being the foliage drag coefficient (0.2), and $\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.$ is the representation of the surface drag coefficient, approximated as $\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.=C_1-C_{2}exp\left(-C_3{C}_{d}LAI\right)$, where $C_1=0.32$, $C_2=0.264$ and $C_3=15.1$.

$$z_{0m}=h_c\left(1-\left(\frac{d}{h_c}\right)\exp \left[\frac{-k}{\left[\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.\right]}\right]\right)$$

(8)

SEBS uses Equation (9) to relate $z_{0h}$ and $z_{0m}$ through the parameter kB⁻¹:

$$z_{0h}=\frac{z_{0m}}{exp\left(k{B}^{-1}\right)}$$

(9)

where $B^{-1}$ is the inverse of the Stanton number.

2.2.1. Parameterization of kB⁻¹ in SEBS and Its Shortcomings

In SEBS, $k{B}^{-1}$ is parameterized by combining the full canopy kB⁻¹ model of Choudhury and Monteith [31], with the bare-soil $k{B}^{-1}$ model of Brutsaert [32], using a quadratic weighting function based on the fractional vegetation cover, as shown in Equation (10):

$$k{B}^{-1}=k{B}_v{}^{-1}{f}_c{}^2+2f_c{f}_{s}k{B}_m{}^{-1}+k{B}_s{}^{-1}{f}_s{}^2$$

(10)

where $k{B}_v{}^{-1}$ is the full canopy $k{B}^{-1}$ (Equation (11)), $k{B}_m{}^{-1}$ is used to account for interactions in mixed surfaces (Equation (12)), and $k{B}_s{}^{-1}$ is the $k{B}^{-1}$ for bare soil (Equation (13)).

$$k{B}_v{}^{-1}=\frac{k{C}_d}{4C_t\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.\left(1-e^{-\frac{n_{ec}}{2}}\right)}$$

(11)

where $C_t$ is the heat transfer coefficient of the leaf.

$$k{B}_m{}^{-1}=\frac{k\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.\frac{z_{0m}}{h}}{C_t{}^{*}}$$

(12)

where $C_t{}^{*}$ is the heat transfer coefficient of the soil, which is given by $C_t{}^{*}=P{r}^{-\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}R{e}_s{}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$, where $Pr$ is the Prandtl number and $R{e}_s$ is the roughness Reynolds number of bare soil, calculated from $u_{*}$, with the roughness height of the soil ($h_s$) and the kinematic viscosity of air ($\nu$) calculated as $R{e}_s=\raisebox{1ex}{$h_s{u}_{*}$}\!\left/ \!\raisebox{-1ex}{$\nu $}\right.$.

$$k{B}_s{}^{-1}=2.46{\left(R{e}_s\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}-\ln \left[7.4\right]$$

(13)

In SEBS, $C_t$ in Equation (11) is held constant and assumes a value of 0.01. However, $C_t$ will be affected by, among others, the response of the leaves to water availability and the prevailing environmental conditions [19,22]. As such, a constant value of $C_t$ may not adequately describe heat transfer within the canopy under varying environmental conditions. Thus, a realistic representation of heat transfer dynamics in the canopy is needed to avoid uncertainties related to the parameterization of $k{B}_v{}^{-1}$ propagating into the simulated sensible heat flux.

2.2.2. Revisions to the Parameterization of $k{B}_v{}^{-1}$ in SEBS

In this section, a revised $k{B}_v{}^{-1}$ which does not explicitly consider $C_t$, is derived. According to Brutsaert [22], the foliage heat transfer coefficient can be parameterized using Equation (14):

$$C_t=C_L{P}_r{}^{-m}R{e}_f{}^{-n}$$

(14)

where $P_r$ is the Prandtl number, $R{e}_f$ is the local foliage Reynolds number calculated as $R{e}_f=\raisebox{1ex}{$l\times u\left(z\right)$}\!\left/ \!\raisebox{-1ex}{$\nu $}\right.$, where $l$ is the characteristic leaf length scale (m) and $u\left(z\right)$ is the wind speed at a given height inside the canopy (ms⁻¹), and $C_L$, m, and n are parameters that depend on the canopy structure and turbulence.

Substituting $C_t$ from Equation (14) into Equation (11) yields Equation (15):

$$k{B}_v{}^{-1}=\frac{k{C}_d}{4C_L{P}_r{}^{-m}R{e}_f{}^{-n}\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.\left(1-e^{-\frac{n_{ec}}{2}}\right)}$$

(15)

Brutsaert [22] suggested that $C_L={\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.}^{0.5}$ if the Reynolds number is parameterized as a function of wind speed. On the other hand, various studies suggest using a value of 0.67 for $m$ [19,22]. Wang et al. [33] showed that for water flow through a vegetated channel, $C_d$ and the Reynolds number ($Re$) can be related through $C_d=R{e}^{-n}$. Assuming that the same relationship is valid for the typical range of Reynolds numbers associated with wind flow over canopies, and subsequently substituting it into Equation (15), the new parameterization of $k{B}_v{}^{-1}$ that does not explicitly consider $C_t$ can be derived using Equation (16):

$$k{B}_v{}^{-1}=\frac{k}{4P_r{}^{-0.67}{\left(\raisebox{1ex}{$u_{*}$}\!\left/ \!\raisebox{-1ex}{$u\left(h\right)$}\right.\right)}^{\left(3/2\right)}\left(1-e^{-\frac{n_{ec}}{2}}\right)}$$

(16)

3. Results

3.1. Savannah Sites

The time-averaged diurnal variations of sensible heat flux at savannah sites, shown in Figure 1, indicate that both the original and revised parameterizations consistently underestimate the observed sensible heat flux for most of the months at both the AU-ASM and AU-Cpr sites. At the AU-GWW site, however, the original parameterization matches the observed sensible heat flux better than the revised parameterization, while both parameterizations overestimate the observed sensible heat flux.

In Figure 2, it is observed that the performance of both parameterizations at the AU-ASM and AU-Cpr is quite similar. However, the revised parameterization performs slightly better than the original parameterization at these two sites. The revised parameterization shows a slightly less underestimation of the observed sensible heat flux and marginally better RMSE values. On the other hand, at the AU-GWW site, the revised parameterization tends to overestimate the observed sensible heat flux, resulting in slightly higher RMSE values than the original parameterization.

Figure 3 shows the temporal variations in the observed versus the simulated half-hourly sensible heat flux at the AU-GWW site for five days. The trend observed in Figure 1 and Figure 2 is also evident in the half-hourly time series. Both parameterizations tend to overestimate the observed sensible heat flux, mainly during midday to afternoon hours. The revised parameterization shows a relatively higher overestimation of peak sensible heat flux than the original parameterization. Both parameterizations perform relatively well during periods characterized by low sensible heat flux, especially at night. A minimal difference is observed between the two parameterizations during periods of low sensible heat flux.

3.2. Deciduous Broad-Leaf Forest Sites

Figure 4 shows time-averaged diurnal variations in sensible heat flux obtained using the two parameterizations versus the observations at the three DBF sites. The revised parameterization matches the observations at DE-Hai and US-MMS better than the original parameterization. The original parameterization shows a significant underestimation of the observed sensible heat flux at both sites during months with relatively high sensible heat fluxes. The underestimation is more pronounced during daytime peak hours of sensible heat flux. However, at the DK-Sor site, the original parameterization matches the observations better than the revised parameterization. The revised parameterization overestimates the observed sensible heat flux at the DK-Sor site.

The statistical metrics shown in Figure 5 prove that the revised parameterization generally performs better than the original parameterization at the DBF sites. The revised parameterization outperforms the original parameterization in terms of the R² at all three DBF sites. The revised parameterization is less biased and shows lower RMSE values at DE-Hai and US-MMS sites. However, the revised parameterization is more biased at the DK-Sor site and shows higher RMSE values than the original parameterization. The revised parameterization shows a better spread of scatter points along the one-to-one line than the original parameterization at all three sites.

The half-hourly time series for five days for the US-MMS site shown in Figure 6 indicate that, similar to the savannah sites, more uncertainties exist in the simulation of peak sensible heat fluxes than low sensible heat fluxes. The original parameterization largely underestimates peak sensible heat fluxes. On the other hand, the revised parameterization matches the observed peak sensible heat flux at the US-MMS site during the five days fairly well.

3.3. Evergreen Broad-Leaf Forest Sites

The time-averaged diurnal variations in sensible heat flux shown in Figure 7 indicate better performance by the revised parameterization at EBF sites than by the original parameterization. The original parameterization at all three sites shows a significant underestimation of the observed sensible heat flux. On the other hand, the revised parameterization matches the observed sensible heat flux variations, especially at the AU-Cum and AU-Rob sites, fairly well. At the FR-Pue site, the revised parameterization still underestimates the observed sensible heat flux but to a lesser extent than the original parameterization.

The variations in sensible heat flux observed in Figure 7 are also supported by the statistical results shown in Figure 8. The revised parameterization outperforms the original parameterization based on all three statistical metrics considered for all three EBF sites. In particular, a significant improvement is observed with regard to the RMSE and bias at the AU-Cum and AU-Rob sites.

Figure 9 shows a half-hourly time series of the observed versus the simulated sensible heat flux by the two parameterizations at the FR-Pue site over five days. Although the revised parameterization simulates a relatively higher peak sensible heat flux than the original parameterization, both parameterizations greatly underestimate the observed peak sensible heat flux at this site. However, similar to the previously presented sites (Figure 3 and Figure 6), during hours of low heat fluxes, both parameterizations match the observed sensible heat flux fairly well.

3.4. Evergreen Needle Leaf Forest Sites

Figure 10 shows the time-averaged diurnal variations of sensible heat flux across the three ENF sites. A similar trend to the one observed in most of the DBF and EBF sites is also observed at the ENF sites. Significant differences in performance between the two parameterizations are observed in months with relatively high observed sensible heat flux, especially at peak hours during the day. There are minimal differences in the performance of both parameterizations during observed low sensible heat flux periods.

Statistically, it is observed in Figure 11 that the revised parameterization outperforms the original parameterization at all the ENF sites based on all three evaluation metrics. Significantly higher R² values and lower RMSE and bias values are observed across the three ENF sites.

The half-hourly time series over 5 days for the NL-Loo site shown in Figure 12 further affirms the superiority of the revised parameterization over the original parameterization in simulating peak sensible heat flux, more so in tall forest canopies.

3.5. Grassland Sites

Figure 13 shows the time-averaged diurnal variations in the modeled versus the observed sensible heat fluxes across the grassland sites. Unlike in the forest sites (DBF, EBF, and ENF), no major differences are observed between the sensible heat fluxes simulated by the two parameterizations at grassland sites. The sensible heat flux simulated by the two parameterizations matches the observed sensible heat flux at the three sites relatively well compared to the forest sites. The revised parameterization slightly overestimates sensible heat flux during some months (January, February, and August) at the AU-Emr site. On the other hand, the original parametrization shows a slight underestimation of the observed sensible heat flux at the AU-Emr, specifically in March, April, May, and November. Both parameterizations underestimate the observed sensible heat flux from November to March at the AU-Rig site, albeit with a marginally less underestimation with the revised parameterization.

The scatterplots and the statics shown in Figure 14 indicate that the difference in the performance between the two parameterizations at the grassland sites is almost insignificant. Based on the RMSE and bias, it appears the revised parameterization shows a marginally improved performance than the original parameterization, especially at the AU-Emr and AU-Rig sites.

The trend observed in the half-hourly time series for the AU-Emr grassland site in Figure 15 is more comparable to the trend observed in the time series for the AU-GWW savannah site (Figure 3) than the forest sites (Figure 6, Figure 9 and Figure 12). Peak sensible heat flux is less underestimated by the original parameterization at the AU-Emr site compared to forest sites. On the other hand, similar to the savannah site in Figure 3, the revised parameterization tends to slightly overestimate the observed sensible heat flux at the AU-Emr grassland site.

4. Discussion

One of the shortcomings of SEBS reported in various studies is the underestimation of sensible heat flux, especially in arid and semi-arid environments [10,16,17,20,21]. The results shown in Figure 1, Figure 2 and Figure 3 for the Savanna sites agree with the findings presented in these studies. Both the original and revised parameterization underestimate the observed sensible heat flux at the savannah sites apart from at the AU-GWW sites. The relatively similar performance observed with both parameterizations at savannah sites is because the two models differ in the parameterization of $k{B}_v{}^{-1}$ only. Most Savanna environments are characterized by low vegetation coverage most of the year. Thus, the contribution of $k{B}_s{}^{-1}$ to the overall $k{B}^{-1}$ as parameterized in Equation (10) will dominate, resulting in the relatively low differences in the performance between the two parameterizations observed at savannah sites. The parameterization of $k{B}_s{}^{-1}$ as adopted in SEBS model is arguably the main contributor to the underestimation of sensible heat flux observed with the SEBS model in arid and semi-arid environments [20,21].

However, as observed in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12, the SEBS model based on the original parameterization significantly underestimates sensible heat flux at DBF, EBF, and ENF sites. These sites are forest sites with characteristically tall canopies and are generally not water limited. Chen et al. [19] also reported a significant underestimation of sensible heat flux with SEBS under tall forest canopies. As pointed out in Brutsaert [22], heat transfer within the canopy is highly dependent on turbulence. Tall forest canopies will induce a significantly greater proportion of large eddies compared to short canopies. Large eddies are associated with greater turbulence and are, thus, very efficient in the transport of heat fluxes between the surface and the atmosphere. Since the original parameterization of $k{B}_v{}^{-1}$ in SEBS assumes a constant value of $C_t$, it implies that canopy structure dependent variations in heat transfer such as the ones related to large eddy transport are largely ignored with the original parameterization. The better match between the observed and the sensible heat flux simulated based on the revised parameterization compared to the original parameterization at the forest sites (Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12), especially during periods of peak sensible heat flux, is mainly related to a better ability of the new parameterization to account for canopy structure-dependent variations in turbulence compared to the original parameterization. During periods of low sensible heat flux, the performance of both parameterizations is quite similar irrespective of the vegetation type. Thus, the differences observed in the performance between the original and the revised parameterization can be mainly attributed to turbulence-induced canopy heat transfer.

The relatively less unbiased sensible heat flux simulated by the SEBS model based on the original parameterization at grassland sites (Figure 13 and Figure 14) in comparison to the other vegetation types indicates that the constant value of $C_t$ (0.01) adopted in SEBS might be more suited for short smoother vegetation canopies. In such canopies, turbulence-induced heat transfer will be significantly less compared to taller canopies, and thus, $C_t$ might not vary considerably. The relatively good performance observed with the original parameterization despite a constant value of $C_t$ might be because the development and validation of the original $k{B}^{-1}$ parameterization used in SEBS was mainly carried out on grass and cotton sites [1]. Figure 13 and Figure 14 show that the performance of the two parameterizations at the grassland sites is quite similar. However, although not significant during peak sensible heat flux periods, the original parameterization still underestimates the observed sensible heat flux more than the revised parameterization. The differences in the two parameterizations during periods of peak sensible heat flux at grassland sites imply that although a constant $C_t$ value of 0.01 seems adequate in describing heat transfer in grassland sites, it might fail for such sites during extreme conditions, e.g., during heatwaves when turbulence is more pronounced.

Although the SEBS model based on the revised parameterization simulates less unbiased sensible heat flux than the SEBS model based on the original parameterization, some biases persist. Some of the biases might be inherent to the parameterization of $k{B}_m{}^{-1}$ and $k{B}_s{}^{-1}$. In addition, as argued in various studies, e.g., in Garrat [30] and Arnqvist and Bergström [29], MOST-based models tend to under-predict momentum flux, leading to an underestimation of sensible heat flux in tall canopies. Under tall canopies, the roughness elements create a roughness sublayer in which the flow deviates from the surface-layer approximations employed in MOST. Another source of uncertainties in the simulated sensible heat flux may be related to uncertainties in the LST data. The field of view of the radiometers that were used to derive LST using Equation (1) is generally smaller than the footprint of the flux measurements. Hence, in situations where the surface is inhomogeneous, the LST from the radiometers might differ significantly from the average LST within the footprint of the flux measurements. Besides, the land surface emissivity used to derive the LST will be highly influenced by the ability of the LAI product to capture the temporal dynamics in vegetation cover. Furthermore, the bare-soil emissivity value of 0.96 adopted in Equation (2) will vary depending on the soil characteristics and soil moisture at each respective site. It should also be noted that a constant value for canopy height was assumed for each site across the respective study period. Assuming a constant value for canopy height will introduce uncertainties in the estimation of sensible heat flux in case the canopy height deviates from the assumed value since $d$, $z_{0m}$, and, consequently, $z_{0h}$ are highly dependent on canopy height.

Apart from uncertainties related to the inadequate description of the physical processes by the two parameterizations, it should be noted that the validation data have inherent uncertainties as well. A lack of energy balance closure generally characterizes most flux towers. Lack of energy balance closure at flux tower sites emanates from a mismatch between footprints of the various sensors, energy storage within the canopy, and advection, among others. Although the sensible heat flux used for validation has been corrected for energy balance closure, uncertainties might exist in the corrected sensible heat flux, since most correction approaches assume the available energy is correctly measured [24,26,28]. This assumption might not be valid, especially under inhomogeneous sites, since the footprints of the radiation and ground heat flux sensors are comparatively smaller than the fetch of the flux sensors.

5. Conclusions

This study presents a revised parameterization of $k{B}_v{}^{-1}$ that does not assume a constant value of $C_t$, as opposed to the original parameterization used in the SEBS model. The revision incorporates a physically based formulation of $C_t$ that includes the influence of turbulence on canopy heat transfer.

Evaluation of the revised parameterization across 15 EC flux tower sites indicates that the parameterization greatly reduces the underestimation of sensible heat flux at tall forest canopies compared to the original parameterization used in SEBS. It is concluded that a constant value of $C_t$, as employed in the original parameterization used in SEBS, is inadequate to capture turbulence-induced canopy heat transfer in tall forest canopies.

The generally better performance observed with the revised parameterization compared to the original parameterization indicates that it is more suited for simulating turbulent surface heat fluxes under tall forest canopies. Besides, the revised parameterization is simple and does not necessitate additional information on the canopy structure beyond the one needed in the original parameterization. As such, it is more suited for the operational mapping of turbulent surface fluxes compared to more complex parameterizations proposed in the literature.