Explaining the convector effect in canopy turbulence by means of large-eddy simulation

Semi-arid forests are found to sustain a massive sensible heat flux in spite of having a low surface to air temperature difference by lowering the aerodynamic resistance to heat transfer (rH) – a property called the “canopy convector effect” (CCE). In this work large-eddy simulations are used to demonstrate that the CCE appears more generally in canopy turbulence. It is indeed a generic feature of canopy turbulence: rH of a canopy is found to reduce with increasing unstable stratification, which effectively increases the aerodynamic roughness for the same physical roughness of the canopy. This relation offers a sufficient condition to construct a general description of the CCE. In addition, we review existing parameterizations for rH from the evapotranspiration literature and test to what extent they are able to capture the CCE, thereby exploring the possibility of an improved parameterization.


Introduction
Understanding the role of turbulence in interactions between vegetation canopies and the atmosphere is crucial for interpreting momentum and scalar fluxes above vegetation.This is relevant for a number of practical applications, such as regional and global weather and climate modeling, energy balance closure studies, and development of forest management strategies.Measurement campaign networks such as FLUXNET monitor carbon, water, and energy fluxes on a long-term basis for this same reason (Baldocchi et al., 2001) to study how different ecosystems in-teract with the atmosphere and influence local and global weather and climate.One such measurement campaign (Rotenberg and Yakir, 2011) focused on semi-arid ecosystems, specifically the Yatir forest situated in the Negev desert in Israel, to study the survival and productivity of the pine forest in spite of the high radiation load and suppressed latent heat flux.An important outcome of this campaign was the concept of the "canopy convector effect" (CCE) introduced by Rotenberg and Yakir (2010), hereafter called RY10.To quote RY10, "With suppressed latent heat flux (LE) because of lack of water, the forest is transformed into an effective "convector" that exploits the low tree density and open canopy and, consequently, high canopy-atmosphere aerodynamic coupling".RY10 ascribed the origin of the CCE to the roughness difference between desert and forest.However, in the present work, we demonstrate that the canopy convector effect appears more generally in canopy turbulence.In fact, we show that the CCE is also a generic artifact of homogeneous canopy turbulence by using large-eddy simulation (LES).In doing so, the canopy aerodynamic resistance to heat transfer (r H ) is revisited.The canopy aerodynamic resistance is a concept borrowed from the evapotranspiration literature where it represents the resistance between the idealized "big-leaf" (a reduced-order representation of the fully heterogeneous 3-D canopy) and the atmosphere for heat or vapor transfer (Monteith, 1973;Foken et al., 1995;Alves et al., 1998;Monteith and Unsworth, 2007).The Penman-Monteith equation to calculate evapotranspiration requires parameterization of the aerodynamic resistance which requires information on roughness lengths for heat and momentum and stability (Penman, 1948;Allen et al., 1998;Cleverly et al., 2013).r H parameterizations are also used in global climate models to describe the canopy-atmosphere interaction at the canopy surface layer (Walko et al., 2000).Thus better parameterizations of r H are of fundamental importance in modeling canopy level fluxes of heat and water vapor which can be used in assessing impacts of climate change, disturbance effects such as vegetation thinning and forest fires, as well as in developing forest management strategies.
We investigate whether the existing parameterizations of the canopy aerodynamic resistance exhibit the CCE, and we identify uncertainties in their application.As the CCE is the crucial mechanism that ensures the survival of the Yatir forest, an improved physical understanding of the CCE is of primordial importance when considering large-scale afforestation in semi-arid regions.
2 Background and theory

The canopy convector effect and aerodynamic resistance
As mentioned earlier, the canopy convector effect was introduced by Rotenberg andYakir (2010, 2011) while studying the interaction of vegetation cover with the surface radiation balance for the Yatir forest.The annual average incoming solar radiation in the Yatir forest is about 238 W m −2 , comparable to that in the Sahara, but the net radiation (R n ) is about 35 % higher than the Sahara (RY10) due to the lower albedo of the forest.However, both remote sensing and local measurements indicated that the surface temperature of the forest canopy in Yatir is lower than the surface temperature of the nearby non-forested area -on annual average by about 5 K.This is striking, as firstly, the lower albedo (by 0.1) of the forest than that of the surrounding shrubland translates into an approximately 24 W m −2 increase in radiation load on the forest canopy.Secondly, the cooler canopy surface suppresses the upwelling longwave radiation, resulting in an additional increase in radiation load by about 25 W m −2 .The combined annual increase in radiation load by about 50 W m −2 associated with the Yatir afforestation in the Negev is quite high and is comparable to the net annual radiation difference between the Sahara and Denmark, for example (RY10).Thirdly, the latent heat flux of evapotranspiration (LE), the obvious cooling and energy dissipation mechanism in temperate forests, is not an option since water is virtually unavailable for about 7 months a year.Thus sensible heat flux (H ) is the only major heat dissipation route, translating into a Bowen ratio (H /LE) as high as 20 or more -unlike temperate forests with a Bowen ratio ≈ 1.
In the Yatir forest, the entire net solar radiation flux (up to 800 W m −2 ) is equilibrated by a massive sensible heat flux (H ) of similar magnitude.Note that this high H cannot be explained by the difference between surface and air temperature ( T = T s − T a ) as the canopy surface is cooler than the surrounding desert surface in this case, but the air temperatures above desert and forest canopy are similar.To expound this apparent contradiction of larger sensible heat flux for smaller T , it is important to recall that, when adopting the simplified big-leaf representation of the forest as a single surface, where ρ and C p are the density and specific heat capacity of air, respectively, T a is air temperature, T s is canopy surface temperature, and r H is the apparent canopy aerodynamic resistance to heat transfer (the word "apparent" is used to indicate that this property is a construct of the formulation and not a direct physical property).Hence the large H is not explained by the temperature difference ( T ) but by a decreased r H . Thus the semi-arid forest with its low tree density and large surface area becomes an efficient low aerodynamic resistance "convector" that is well coupled to the atmosphere above (Rotenberg andYakir, 2010, 2011).This "canopy convector effect" (CCE) is adequate enough to support the massive sensible heat flux larger than the surrounding Negev desert, still maintaining a relatively cooler (than the desert) surface temperature (of the canopy top).It is worth noting here that Eq. ( 1) offers a very simplistic description of the complex mixing process in the surface layer; however, it should be interpreted as a zeroth-order representation of the corresponding processes.RY10 identified the difference of roughness between the desert and forest as the underlying mechanism of the CCE by arguing that r H ∝ 1/PAI where PAI denotes the plant area index.However, in this work, we attempt to identify a more fundamental mechanism behind the CCE which is more strongly connected to the feature of canopy turbulence.Therefore we hypothesize that even with the same physical roughness, variation of the aerodynamic roughness is a sufficient condition for displaying the CCE.This difference of aerodynamic roughness for the same physical roughness (of the same vegetation canopy) can be generated by changing the intensity of atmospheric stratification (Zilitinkevich et al., 2008).Thus observing the variation of the canopy aerodynamic resistance to heat transfer (r H ) with atmospheric instability is a sufficient condition to demonstrate the generality of the CCE.To be more precise, if r H is found to decrease with increasing unstable stratification, that would exhibit the fact that canopy turbulence effectively reduces the aerodynamic resistance to cope with heat stressed environments; i.e., the canopy convector effect would manifest itself.Therefore, to summarize the canopy convector effect in simpler terms, it can be mentioned that the darker and colder canopy surface reduces albedo, which leaves more of the incoming energy on the canopy surface.However, the organization of these dark leaf surfaces is such that they are spread Hydrol.Earth Syst.Sci., 21, 2987Sci., 21, -3000, 2017 www.hydrol-earth-syst-sci.net/21/2987/2017/ over a relatively thick canopy depth (relative to grassland or shrubland, where all leaves are condensed in a much thinner layer).Because canopy in dry forests is sparse, wind can easily penetrate it and can easily exchange heat with the leaf surfaces.Therefore, forests would have intrinsically lower aerodynamic resistance to heat transfer than shorter biomes because of the higher roughness.Moreover, the same forest (with the same physical roughness) could have higher aerodynamic roughness and consequently lower aerodynamic resistance to heat transfer for more heat stressed conditions.Given Eq. (1), that would mean higher heat flux.Thus while the CCE would always be present in a forest compared to a grass or shrubland because of the obvious roughness difference, we establish that the CCE can also be present within the same forest for different conditions of heat stress -which is a more subtle point and will be further discussed in the following sections by using large-eddy simulation (LES).LES provides a useful and meanwhile standard tool for studying canopy turbulence under different conditions of atmospheric stratification.A recent publication by Patton et al. (2015) studied the influence of different atmospheric instability classes on coupled boundary layer-canopy turbulence.In this work, those same instability classes are simulated to put our hypothesis to the test.

Parameterizations for canopy aerodynamic resistance to heat transfer
Apart from the LES outcomes, it is also important to study whether the existing parameterizations of r H can exhibit the CCE.Parameterizations of r H in the literature use Monin-Obukhov similarity theory (MOST) extensively.MOST can provide corrections for the vertical profile of the mean longitudinal velocity u and potential temperature (T a − T s ) under thermal stratification, which deviates from the traditional log-law under neutral conditions.Thus under MOST, with the assumption that the vegetation is low, dense, and horizontally homogeneous, and where u * is the friction velocity, κ is the von Kármán constant, z is the height from the ground, d is the zero-plane displacement height, often approximated as (2/3)h c as per the literature (Seginer, 1974;Shuttleworth and Gurney, 1990;Alves et al., 1998;Maurer et al., 2013Maurer et al., , 2015)), and ζ = (z − d)/L is called the stability parameter.L is the Obukhov length, computed as where g = 9.81 m s −2 , the gravitational acceleration.w T is the sensible heat flux -assumed to be constant in the surface layer (Foken, 2006).Negative ζ indicates unstable stratification and thus ζ decreases with increasing instability.z 0m and z 0h are the characteristic roughness lengths for momentum and heat transfer, respectively.ζ 0m = z 0m /L and ζ 0h = z 0h /L are the stability parameters associated with roughness lengths.Pr 0 = K m /K h is the turbulent Prandtl number where K m and K h are eddy diffusivities of momentum and heat, respectively.T * is a characteristic temperature scale, obtained from H and the characteristic velocity scale, i.e., Combining Eqs. ( 1), ( 2), (3), and ( 5), one can write where ψ m and ψ h are the integral stability correction functions for momentum and heat, respectively.Following Liu et al. (2007), they can be parameterized for unstable conditions as (Dyer and Hicks, 1970;Paulson, 1970;Dyer, 1974;Garratt, 1977;Webb, 1982) where . Different values for the parameters γ m and γ h are reported in the literature, and the ones suggested by Paulson (1970) are used, i.e., γ m = γ h = 16.This formulation for r H given by Eq. ( 6) with some approximations (ζ 0m = ζ 0h = 0) was first used by Thom (1975) and is called the "reference parameterization" (Liu et al., 2007).
The full form of Eq. ( 6) was used by Yang et al. (2001) with their only approximation being Pr 0 = 1.Several other studies also used semi-empirical and empirical parameterizations and included the bulk Richardson number Ri B (Monteith, 1973) given by with U the horizontal wind speed at the height that corresponds to the T a measurement.Liu et al. (2007) compiled different parameterizations of r H which we will test in the context of the canopy convector effect against our LES output.2007).These parameterizations based on MOST (Thom, 1975;Yang et al., 2001), empirical (E) (Verma et al., 1976;Hatfield et al., 1983;Mahrt and Ek, 1984;Xie, 1988) and semi-empirical (SE) (Choudhury et al., 1986;Viney, 1991) assumptions can be classified into two categories.Formulations by Thom (1975), Choudhury et al. (1986), Yang et al. (2001) and Viney (1991) have assumed z 0m = z 0h , which should be a more realistic assumption.On the other hand, formulations by Verma et al. (1976), Hatfield et al. (1983), Mahrt and Ek (1984) and Xie (1988) assumed z 0m = z 0h .Different parameters used in the empirical formulations are also listed in Table 1.One important point to note is that only the formulation by Yang et al. (2001) uses the stability parameters associated with the roughness lengths ζ 0m and ζ 0h .Also note that all pa-rameterizations assume a turbulent Prandtl number of unity; i.e., the diffusivities for momentum and heat are assumed to be the same.We shall later discuss the consequence of letting this parameter vary.Another important approximation necessary to evaluate all formulations in Table 1 is a prescription for the roughness lengths z 0m and z 0h .Effects of different roughness lengths will be investigated in the following section.However, a relation between the two roughness lengths (κB −1 = ln(z 0m /z 0h )) was proposed by Owen and Thomson (1963) and Chamberlain (1968), where κB −1 is called an "excess resistance parameter".Yang et al. (2001) suggested an average value of κB −1 = 2.0 (Liu et al., 2007), which will be used throughout this work.
Before moving on to the usage of LES, it warrants mentioning that the entire roughness length formulation (Eqs.2-8 and Table 1) is based on different variants of analytic approximation approaches to reduce the complexity of flow in and above the forest canopy to a 2-D surface equivalent.It is widely accepted that the MOST approach is not completely accurate close to the canopy (Foken, 2006).It was proposed that a mixing length driven approach can be applied (Harman and Finnigan, 2007).Nonetheless, large-scale models, which cannot vertically resolve the canopies, still use MOST, and it has been demonstrated to be relatively accurate.Thus from an operational perspective, the present formulation revisits the current leading approach for simplification of the physics in a parameterized way that can be used by coarse-resolution models.

Methodology
The PALM large-eddy simulation model (Raasch and Schröter, 2001;Maronga et al., 2015) is used to investigate this generic nature of the canopy convector effect.The representation of the canopy in the LES follows the standard distributed drag parameterization (Shaw and Schumann, 1992;Watanabe, 2004;Patton et al., 2015) by adding an additional term in the momentum budget equations as where a is a one-sided frontal plant area density (PAD), C d is a dimensionless drag coefficient assumed to be 0.3 (Katul et al., 2004;Banerjee et al., 2013), |u| is the wind speed, and u i is the corresponding velocity component (i = 1, 2, 3, i.e., u, v, and w).The effect of the canopy on the subgrid-scale (SGS) turbulence is accounted for by adding a sink term to the prognostic equation for the SGS turbulent kinetic energy (e) as F = −2C d a |u|e.For closure of the SGS covariance terms, PALM uses the 1.5 order closure developed by Deardorff (1980) as modified by Moeng and Wyngaard (1988) and Saiki et al. (2000), which assumes a gradient-diffusion parameterization.The diffusivities associated with this gradient diffusion are parameterized using the subgrid-scale turbulent kinetic energy (SGS-TKE) and include a prognostic equation for the SGS-TKE.This SGS-TKE scheme after Deardorff (1980) is deemed to be an improvement over the more traditional Smagorinsky (1963) parameterization since the SGS-TKE allows for a much better estimation for the velocity scale corresponding to the subgrid-scale fluctuations (Maronga et al., 2015).Further details of the LES model can be found in the literature and are not discussed here (Shaw and Schumann, 1992;Watanabe, 2004;Maronga et al., 2015;Patton et al., 2015).For our simulation, the number of grid points in the x, y, and z directions are 320, 320, and 640, respectively, with grid resolutions of 3.91, 3.91, and 1.95 m in the respective directions.Each simulation has a simulated time of 10 000 s with a time step of 0.1 s, while the outputs of the first 6400 s are discarded before achieving computational quasi-equilibrium.The canopy height (h c ) is taken as 35.0 m with a plant area index (PAI) of 5.0.It is important to note that Rotenberg and Yakir (2011) reported an effective PAI of about 5-6 for heat exchange for the Yatir forest.This makes our PAI similar to a recent simulation study of Dias-Junior et al. (2015).In fact, as we already simulate a homogeneous canopy to show that the CCE appears more generically above vegetation canopies, we have decided to tailor our simulations following the examples of Patton et al. (2015) and Dias-Junior et al. (2015) in order to allow a better comparison of the LES data.The vertical distribution of plant area density (a) follows the probability density function (pdf) of a Beta distribution as described in Markkanen et al. (2003) and the parameters α and β controlling the vertical distribution of foliage are set as 3.0 and 2.0, respectively, to simulate a PAD distribution similar to Dias-Junior et al. (2015).The parameters to drive the simulations for five different instability classes, namely, near neutral (NN), weakly unstable (WU), moderately unstable (MU), strongly unstable (SU), and free convection (FC), are similar to those of Patton et al. (2015) and are presented in Table 2.Note that the canopy convector effect as a general phenomenon should not depend on water content in the soil-plant-atmosphere continuum and, moreover, the PALM-LES does not take into account any physiological processes which normally happen with a larger timescale.Nevertheless, instead of simulating a specific dry water free environment, some moisture at the lower surface is provided and the boundary conditions for surface moisture content are taken as similar to the simulations of Dias-Junior et al. (2015) as well.The initial conditions of the potential temperature (and moisture) profile are also taken as similar to Dias-Junior et al. (2015).PALM's canopy module allows sensible heat flux input at the canopy top only, and the sensible heat flux is attenuated exponentially due to the decay of the incoming energy by absorption and reflection by the leaves.Thus the ground surface heat flux would be different from Patton et al. (2015).Another important point to note is that instead of lowering the wind speeds while maintaining similar sensible heat fluxes, the different stability classes can also be achieved by maintaining the same wind speed and ramping up the surface sensible heat fluxes.However, this should not affect the generic feature of the CCE as discussed at the end of Sect.2.1.Further details and boundary conditions about the LES are discussed in Appendix B.
It is worth highlighting again here that the large-eddy simulations have been conducted with an explicit 3-D canopy.This means that the surface assumptions are not needed to develop a revised approximation approach for the surface equivalence that accounts for the forest density effects.Only the outcomes of the LES are parameterized in a way that will allow resolving of the canopy convector effect even in largescale models.
Table 2. Parameters to drive the simulations for five different instability classes, namely, near neutral (NN), weakly unstable (WU), moderately unstable (MU), strongly unstable (SU), and free convection (FC), are similar to Patton et al. (2015).U g and V g denote geostrophic wind speeds, w T toc denotes canopy top surface sensible heat flux, T s denotes ground surface potential temperature, and q s denotes specific humidity at the ground surface.

Stability class
(U g , V g ) (m s −1 ) w T s (K m s −1 ) T s (K) q s (g g −1 ) Near neutral (NN) 20, 0 0.18 307.7 0.02 Weakly unstable (WU) 10, 0 0.18 307.7 0.02 Moderately unstable (MU) 5, 0 0.18 307.7 0.02 Strongly unstable (SU) 2, 0 0.18 307.7 0.02 Free convection (FC) 0, 0 0.18 307.7 0.02  (2015).It is interesting to observe that the magnitude of velocity, the velocity fluctuations, and the turbulent intensity decrease gradually from the near neutral to free convective conditions, i.e., with increasing instability.The potential temperature also reduces with increasing instability at all heights.On the other hand, the sensible heat flux appears to increase with increasing stability, especially more above the forest (z/h c = 1 indicates the canopy top).These results are physically consistent.The near neutral case is dominated by mechanical shear driven turbulence -given by the highest mean velocity.The free convection case is fully buoyancy driven and the motion is fully upwards -as is evident from the near zero mean hori- zontal velocity.For the same reasons, the turbulent intensity and friction velocity follow the same pattern.The strongly unstable cases have the highest heat fluxes, which is also physically consistent.Note that the canopy top sensible heat flux is similar to the imposed value of 0.18 K m s −1 that was used to drive the simulations.
Figure 2 shows temporally and spatially averaged vertical profiles for the different stability conditions with the same color coding as Fig. 1.We investigate the vertical profile of  1).As is evident from panel (c), the aerodynamic resistance reduces with increasing instability, confirming the hypothesis constructed earlier and thus clearly demonstrating the canopy convector effect (CCE).As noted by Zilitinkevich et al. ( 2008), with increasing instability, "convective updraughts developing at side walls of roughness elements extend upwards and provide extra resistances to the mean flow.Then the mean flow interacts with both solid obstacles and their virtual extensions (updraughts), which results in the increased roughness length".This increased roughness can be recognized as the aerodynamic roughness.For the same physical roughness of the canopy, an increase in instability increases this aerodynamic roughness and, in turn, reduces r H .The low aerodynamic resistance effectively allows larger eddies to form above the forest canopy which are more efficient to dissipate the sensible heat by promoting buoyancy.This description refers to a more general phenomenon as opposed to the description by Rotenberg and Yakir (2011) which identifies the higher physical roughness of the canopy compared to the desert and is thus a more site specific description.Nevertheless, it is acknowledged that the more generic description presented here can be reconciled with the explanation from Rotenberg and Yakir (2011) by noting that increased physical roughness can also result in increased aerodynamic roughness.Also incidentally, RY10 reported a value of r H ≈ 16 for the Yatir forest which is of similar order of magnitude as what is found in panel (c) of Fig. 2. One important point to note in Fig. 1 is the magnitude of the Prandtl number, which is almost fixed to about 0.335 above the canopy.This can be reconciled with the theoretical prediction of the variation of Pr 0 with stability by Li et al. (2015).For stability ranges 1 ≤ −ζ ≤ 10, Pr 0 is also estimated to be approximately 0.33, consistent with the stability ranges plotted in Fig. 2. The variation of Pr 0 with stability is discussed further in Appendix A.

Testing different parameterizations
It is interesting to study whether the different parameterizations capture the correct behavior of r H at different heights across stability.To compute r H variations, the LES generated profiles of mean velocity u, sensible heat flux, air temperature, and Prandtl number (thus the diffusivities) are used where all of them have z variations.The friction velocity u * and the roughness lengths are fixed.Figure 3 plots the variation of r H with height as obtained from the LES (black "+" markers), and the predicted r H from different parameterizations for the different stability cases -near neutral (column 1), weakly unstable (column 2), mildly unstable (column 3), and strongly unstable (column 4).The top row compares the parameterizations by Thom (1975)  Figure 3. Variations of r H with height across stability ranges and comparisons with different parameterization schemes as described in Table 1. and where d s and z 0ms are the stability dependent zero-plane displacement length and roughness lengths for momentum, respectively, d and z 0m being their neutral counterparts.d = (2/3)h c can be assumed as usual.The neutral z 0m can be assumed to be related to LAI as given by Shuttleworth and Gurney (1990).According to the relation used by Shuttleworth and Gurney (1990), for an LAI of 5, z 0m = 0.12 h c can be obtained (which is almost constant for a wide range of canopy drag coefficients and LAI).Moreover, if one uses the correct stability dependent Prandtl number Pr 0 (ζ ) instead of setting it to unity, an improved parameterization based on Yang et al. (2001) can be written as Note that z 0ms and z 0hs are still related by the same relation κB −1 = ln(z 0ms /z 0hs ) with κB −1 = 2.0 as discussed earlier and ζ = (z − d s )/L and ζ 0ms = z 0ms /L, ζ 0hs = z 0hs /L.If a Prandtl number of unity is still assumed but the roughness lengths are assumed to be varying with stability as given by Eq. ( 13) with a neutral value of z 0m = 0.2 h c , the formulation by Yang et al. (2001) is found to display the correct behavior of canopy convector effect with stability as observed in Fig. 6.Panel (a) shows the variation of r H with height across stability according to the improved formulation as given by Eq. ( 13).Panel (b) shows similar variations of r H computed from the LES repeated again for comparison.The profile for the near neutral case crosses over the more highly unstable cases at heights around 6 h c ; however, the general behavior of the CCE is captured well.On the other hand, if the full complexity of Eq. ( 13) is used including a stability dependent Prandtl number (discussed in Appendix A) but using the canopy top surface value of the sensible heat flux for all computations involved, the variation of modeled r H is shown in panel (a) of Fig. 7. r H computed from the LES results using the surface value of the heat flux is shown in panel (b).This assumption of a constant sensible heat flux in the canopy sub layer or the atmospheric surface layer is a more realistic one than a monotonically reducing sensible heat flux with height as shown in panel (e) of Fig. 1.In fact, the surface layer is defined as a constant flux layer (Stull, 2012).The reducing flux profiles in LES are common features of large-eddy simulations since the top boundary of the LES domain is assumed stress free (Shaw and Schumann, 1992).Figure 7 correctly captures the order of magnitude of the r H observed from the simulations, and also captures the CCE correctly.However, it is acknowledged that the exact profiles of the observed r H can not be captured.However, these different comparisons highlight the uncertainties involved in the parameterization of r H .   13), but using the top-of-canopy surface flux throughout all heights; (b) similar variations of r H computed from the LES using the top-of-canopy surface flux assumed to be constant in the surface layer.

Conclusion
The canopy aerodynamic resistance is a concept borrowed from the evapotranspiration literature where it represents the resistance between the idealized "big-leaf" (a reduced-order representation of the fully heterogeneous 3-D canopy) and the atmosphere for heat or vapor transfer (Alves et al., 1998).In semi-arid ecosystems, vegetation canopies maintain a relatively cool surface temperature in spite of the high sensible heat flux by reducing the canopy aerodynamic resistance to heat transfer (r H ) -a phenomenon named the "canopy convector effect" by Rotenberg and Yakir (2010).In the present work, a large-eddy simulation is used to examine this canopy convector effect and, in the process, several existing parameterizations for r H are examined.The objectives behind this exploration are 2-fold.The first one is to investigate whether the existing parameterizations exhibit the canopy convector effect and the second one is to identify the uncertainties associated with these different parameterizations since they are applied in different climate models often under conditions of thermal stratification.As illustrated by the LES results, r H above the canopy is found to reduce systematically as the strength of unstable stratification increases.This is deemed to be the core feature of the canopy convector effect, since with increasing instability, more convective updraughts enhance the roughness over the canopy elements that the mean flow encounters.The height variation of r H is also found to have a highly nonlinear profile; thus, any model prescribing a parameterization for r H needs to employ considerable caution regarding the height it is prescribed.Existing parameterizations of r H employ either Monin-Obukhov similarity theory (MOST) or Richardson number based empirical or semi-empirical formulations to account for thermal stratification.However, most of them are found to be unable to de-scribe the correct trend of the CCE.Among different formulations, the one by Yang et al. (2001) is found to be the most promising candidate.This parameterization employs MOST, and accounts for stability parameters associated with roughness lengths for momentum and heat transfer.It is found that a stability dependent zero-plane displacement height as well as stability dependent roughness lengths for momentum and heat transfer can improve its performance.Moreover, if the surface layer or the canopy sublayer is assumed to have a constant sensible heat flux equal to the flux at the canopy top, and a stability dependent Prandtl number is used, the performance improves further.These assumptions also lead to a less nonlinear height variation.These explorations highlight the uncertainties associated with the parameterizations of r H .One possible major source of uncertainty is the usage of Monin-Obukhov similarity theory in the canopy sublayer (CSL) (up to 3 h c to 6 h c ) since it is not expected to perform in the CSL (Kaimal and Finnigan, 1994).Nevertheless, MOST formulations are found to outperform other semiempirical formulations using Richardson numbers.Thus future research work will involve studying these uncertainties of r H parameterizations in regional and global climate models.The consequence of this CCE for local circulation, atmospheric moisture, and tree physiology will also be investigated, extending the preliminary study of Eder et al. (2015).However, the fact that the CCE is a more generic feature of canopy turbulence provides hope that the afforestation of an area larger than the Yatir forest would also be able to cope with a high-radiation load under water scarcity in semi-arid climates.

Figure 1 .
Figure 1.Summary statistics of five LES simulations showing the variations between different stability classes in increasing order of instability -from near neutral to free convection color coded as indicated in the legend.

Figure 2 .
Figure 2. Aerodynamic resistance to the heat transfer exhibiting the canopy convector effect.(a) Difference between surface and air temperature (T s − T a ); (b) stability parameter ζ ; (c) canopy aerodynamic resistance (r H ).
r H in order to assess the uncertainty that arises from varying the reference height for the air temperature under varying stability.The temperature of the canopy top is taken as the surface temperature (T s ) and thus results are shown from above the canopy top, i.e., z/h c = 1.Panel (a) shows the difference of surface and air temperature (T s − T a (z)).Panel (b) shows the stability parameter ζ at every level computed as ζ = (z−d)/L as explained in Sect. 2. Panel (c) plots canopy aerodynamic resistance to heat transfer (r H ) at every level computed from Eq. ( (blue line),Yang et al. (2001) (red line),Choudhury et al. (1986) (black line), and Viney (1991) (pink line), which assume z 0m = z 0h .The bottom panel compares the parameterizations by Verma www.hydrol-earth-syst-sci.net/21/2987/2017/ Hydrol.Earth Syst.Sci., 21, 2987-3000, 2017

Figure 4 .
Figure4.Variations of r H with height for different stability classes computed for each parameterization scheme as described in Table1.

Figure 5 .
Figure 5. Variations of r H as given by the parameterization of Yang et al. (2001) with height across stability ranges and a wide range of z 0m .Black "+" markers indicate the observed r H from LES at any particular stability state.

Figure 6 .
Figure 6.(a) The variation of r H with height across stability according to the improved formulation as given by Eq. (13); (b) similar variations of r H computed from the LES repeated again for comparison.

Figure 7 .
Figure 7. (a) The variation of r H with height across stability according to the improved formulation as modeled by Eq. (13), but using the top-of-canopy surface flux throughout all heights; (b) similar variations of r H computed from the LES using the top-of-canopy surface flux assumed to be constant in the surface layer.

Table 1 .
Table 1 lists the details of the different parameterizations as compiled by Liu et al.Different parameterizations of r H as compiled by Liu et al. (2007).