Turbulent transport of energy across a forest and a semiarid shrubland

The role of secondary circulations has recently been studied in the context of well-defined surface heterogeneity in a semiarid ecosystem where it was found that energy balance closure over a desert–forest system and the structure of the boundary layer was impacted by advection and flux divergence. As a part of the CliFF (“Climate feedbacks and benefits of semi-arid forests”, a collaboration between KIT, Germany, and the Weizmann Institute, Israel) campaign, we studied the boundary layer dynamics and turbulent transport of energy corresponding to this effect in Yatir Forest situated in the Negev Desert in Israel. The forest surrounded by small shrubs presents a distinct feature of surface heterogeneity, allowing us to study the differences between their interactions with the atmosphere above by conducting measurements with two eddy covariance (EC) stations and two Doppler lidars. As expected, the turbulence intensity and vertical fluxes of momentum and sensible heat are found to be higher above the forest compared to the shrubland. Turbulent statistics indicative of nonlocal motions are also found to differ over the forest and shrubland and also display a strong diurnal cycle. The production of turbulent kinetic energy (TKE) over the forest is strongly mechanical, while buoyancy effects generate most of the TKE over the shrubland. Overall TKE production is much higher above the forest compared to the shrubland. The forest is also found to be more efficient in dissipating TKE. The TKE budget appears to be balanced on average both for the forest and shrubland, although the imbalance of the TKE budget, which includes the role of TKE transport, is found to be quite different in terms of diurnal cycles for the forest and shrubland. The difference in turbulent quantities and the relationships between the components of TKE budget are used to infer the characteristics of the turbulent transport of energy between the desert and the forest.

1 Introduction

Understanding the interaction between vegetation canopies and atmosphere is a crucial component in the quantification of biosphere–atmosphere exchange of heat, carbon dioxide, water and trace gas fluxes. It is also important for the development of numerical weather and climate models where the fluxes in the canopy surface layer (CSL) and the atmospheric surface layer (ASL) are parameterized through bulk exchange coefficients of momentum and scalar. However, idealizations of the forest canopies as horizontally homogeneous momentum sinks and scalar sources introduces uncertainties in flux estimations and estimating diffusion coefficients. The presence of heterogeneities such as roughness transitions, complex topography and mesoscale circulations are common sources of such uncertainties that give rise to nonlocal motions and secondary circulations. These secondary circulations not only occur in forests but are also generic characteristics of boundary layer flows over natural and man-made landscapes with discongruity of land use types, surface moisture, temperature, etc. (Eder et al., 2015; Higgins et al., 2013). Different types of land cover such as agricultural lands or urban areas can affect local energy balance closure and the structure of the overlying boundary layer as well as cloud formation and regional weather (Eder et al., 2015; Fuentes et al., 2016). Strong differences in surface properties and large swaths of such surface patches are known to induce secondary circulations (Banerjee et al., 2013; Courault et al., 2007; Dalu and Pielke, 1993; Dixon et al., 2013; Garcia-Carreras et al., 2010; Kang and Lenschow, 2014; Mahfouf et al., 1987; Raupach and Finnigan, 1995; Sühring and Raasch, 2013; Van Heerwaarden et al., 2014; van Heerwaarden and Guerau de Arellano, 2008). Recent works by Mauder et al. (2007), Stoy et al. (2013) and Eder et al. (2014) have suggested that non-closure of the energy balance is also related to advection and flux divergence due to secondary circulations (Foken, 2008; Kanda et al., 2004). The non-closure of the energy balance refers to the fact that the available energy R_n−G is often higher than the turbulent energy H+LE at micrometeorological sites, where R_n is net radiation, G is soil heat flux, H is sensible heat flux and LE is latent heat flux. Thus, it is established that studies involving surface heterogeneities such as a difference in roughness characteristics and albedo are crucial for the advancements of our understanding of biosphere–atmosphere interaction since the quasi-universal scaling laws of turbulent moments and simple parametrizations of exchange coefficients are disturbed and rendered nonoperational.

Several studies have attempted to study the nature of turbulence across a roughness transition such as a grassland and a forest canopy by means of experimental and numerical methods (Banerjee et al., 2013; Belcher et al., 2003; Cassiani et al., 2008; Chatziefstratiou et al., 2014; Dalpe and Masson, 2009; Detto et al., 2008; Dupont and Brunet, 2009; Fesquet et al., 2009; Gavrilov et al., 2010, 2011; Huang et al., 2011; Irvine et al., 1997; Kanani-Sühring and Raasch, 2015; Kröniger et al., 2017; Li et al., 1990; Markfort et al., 2014; Peltola, 1996; Queck et al., 2016; Rominger and Nepf, 2011; Schlegel et al., 2012; Yang et al., 2006) and documented several length scales associated with the roughness transitions, recirculation zones and the nature of the turbulent momentum budget. However, all of these studies are concerned with the flow adjustment in the immediate vicinity of the roughness transition (edges or gaps). Eder et al. (2015) have studied the dynamics of the convective boundary layer over a well-defined surface heterogeneity – namely Yatir Forest and the shrubland surrounding it, which are located in the northern part of the Negev Desert in Israel. Eddy covariance (EC) and Doppler lidar measurements were conducted by Eder et al. (2015) at two sites approximately 6.5 km apart: one in the forest and one in the desert. The forest has a darker surface and consequently lower albedo (12.5 %) than the desert (33.7 %). Moreover, the higher surface roughness of the forest results in higher turbulence intensity, which leads to more efficient heat transfer above the forest, a phenomenon called canopy convector effect (Banerjee et al., 2017a; Rotenberg and Yakir, 2011). The region being very dry, there is very little latent heat flux (Bowen ratio >10 over the summer), resulting in a spatial difference in surface buoyancy flux of 220–290 W m⁻² between the desert and the forest. Furthermore, the length scale of surface heterogeneities (6–10 km) is larger than the minimal length scale needed for the development of secondary circulations: $L_{\mathrm{rau}}=C_{\mathrm{Rau}}U/w_{*}\approx \mathrm{2}$–5 km (Eder et al., 2015; Raupach and Finnigan, 1995), where U is mean wind speed, w_* is the convective velocity scale and C_Rau=0.8 is an empirical parameter, so that it is possible for secondary circulations to develop.

The present work is an attempt to examine this hypothesis of secondary circulations in more detail. We use eddy covariance and Doppler lidar measurements at two sites 4.3 km apart over the shrubland and Yatir Forest, where the shrubland is upwind of the forest in the path of the principal wind direction (during the summer, there exists a heat-induced low-pressure system to the east, resulting in the main wind direction from the northwest). We investigate the individual components of the turbulent kinetic energy budget, as well as the nature of advection and turbulent transport over the forest and the desert and determine if there is a relationship between them. Not many instances were found in the literature where the nature of turbulent transport was studied across large-scale surface roughness heterogeneities, except for Nadeau et al. (2011) and Yue et al. (2015). However, Yue et al. (2015) only studied turbulent production and the turbulent velocity fluctuations in the presence of a complex topography – so the nature of turbulent transport via secondary circulations was not highlighted. Nadeau et al. (2011) studied the decay of turbulence over different land surface types. Hence, the difference in turbulence production and simultaneous transport across different land use types was not studied, which determines the scope of the current work.

2 Method

2.1 Theory

The turbulent kinetic energy (TKE) budget is given by Stull (2012) without invoking any special assumption:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial e}{\partial t}}+U_j{\displaystyle \frac{\partial e}{\partial x_j}}={\mathit{\delta }}_{i\mathrm{3}}{\displaystyle \frac{g}{T}}\left(\stackrel{\mathrm{‾}}{u_i^{\prime }{T}^{\prime }}\right)-\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}{\displaystyle \frac{\partial U_i}{\partial x_j}}-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u_j^{\prime }e}\right)}{\partial x_j}}{\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u_i^{\prime }{p}^{\prime }}\right)}{\partial x_i}}-\mathit{ϵ},\end{array}$$

where i and j are the usual tensor indices, which can take the values of 1, 2 and 3 to indicate x, y and z directions, respectively, and δ_i3 is the Kronecker delta. $e=(\mathrm{1}/\mathrm{2})({\mathit{\sigma }}_u^{\mathrm{2}}+{\mathit{\sigma }}_v^{\mathrm{2}}+{\mathit{\sigma }}_w^{\mathrm{2}})=(\mathrm{1}/\mathrm{2})(\stackrel{\mathrm{‾}}{u^{\prime \mathrm{2}}}+\stackrel{\mathrm{‾}}{v^{\prime \mathrm{2}}}+\stackrel{\mathrm{‾}}{w^{\prime \mathrm{2}}})$ is the TKE, U denotes mean longitudinal velocity; u^′, v^′ and w^′ denote the fluctuations from mean for the longitudinal, transverse and vertical velocity components; g is acceleration due to gravity; T denotes mean potential temperature; T^′ is the potential temperature fluctuation; p^′ is the dynamic pressure perturbation; ρ is density of air. The first term on the left-hand side (LHS) denotes storage or TKE tendency. The second term on the LHS indicates advection of TKE by mean wind flow. The first term on the right-hand side (RHS) denotes buoyant production/destruction of TKE. The second term on the RHS denotes mechanical/shear production of TKE. The third term on RHS denotes turbulent transport of TKE and can also be called turbulent flux divergence. The fourth term on RHS denotes transport of TKE by pressure velocity correlation. ϵ is the dissipation of TKE.

Expanding the equations in terms of x, y and z coordinates, the full TKE budget can be written as Eq. (A1) as shown in Appendix A. Since it is difficult to keep track of the full equation due to the large number of terms, it would be easier to use a simple form of the TKE budget (Stull, 2012)

$$\begin{array}{}\text{(2)}& \mathrm{0}=-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}z}}+{\displaystyle \frac{g}{T}}\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}-\mathit{ϵ}-\text{Imbalance},\end{array}$$

where the “imbalance” is defined in Eq. (A2). Note that $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$ denote vertical momentum flux and sensible heat flux, respectively. Also note that if the term imbalance is set to zero, one recovers the TKE budget for an idealized surface layer where the coordinate system is aligned with the mean wind, and a planar, homogeneous flow with zero subsidence is assumed. Since our objective in the current problem is to study the effect of heterogeneity, we cannot make these assumptions. Moreover, we are also constrained by being able to measure only at two single points in space quite far apart. Single point eddy covariance measurements cannot compute spatial gradients, and the pressure perturbations are not measured either. Thus, explicit computations of the imbalance terms are not possible. Due to the three-dimensional nature of the problem, it is also difficult to anticipate what degrees of assumptions are sufficient, so that some of the terms can be ignored safely.

Under these constraints, a strategy is needed to evaluate the TKE budget. The dominant mechanical production term, the buoyant production/destruction term and the dissipation term will be evaluated directly from the data. The residual of the TKE budget will be described as the imbalance as per Eq. (3) which would contain the effects of advection and transport terms. The advantage of using this strategy is that since the original TKE budget equation has to be closed, the errors in computing the production and dissipation terms can also be assumed to be inside the imbalance term.

$$\begin{array}{}\text{(3)}& \text{Imbalance}=-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}z}}+{\displaystyle \frac{g}{T}}\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}-\mathit{ϵ}.\end{array}$$

To compute the mechanical production term, we momentarily assume that the TKE budget is well balanced and Monin–Obukhov similarity theory (MOST) (Monin and Obukhov, 1954) is valid (Banerjee et al., 2016; Kaimal and Finnigan, 1994; Li et al., 2016). This allows us to write

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}z}}={\mathit{\varphi }}_{\mathrm{m}}\left(\mathit{\zeta }\right){\displaystyle \frac{u_{*}}{\mathit{\kappa }z}},\end{array}$$

where ϕ_m is the stability correction function for momentum which varies with the stability parameter $\mathit{\zeta }=(z-d)/L$ and κ=0.4, the von Kármán constant. $u_{*}=\sqrt{{\left|\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\right|}^{\mathrm{2}}+{\left|\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}\right|}^{\mathrm{2}}}$ is the friction velocity, z is the measurement height, and $L=-u_{*}^{\mathrm{3}}/\left(\mathit{\kappa }\right(g/T\left)\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}\right)$ is the Obukhov length; d is zero plane displacement height, taken as 2∕3 of canopy height. The standard MOST scaling relations for ϕ_m are used, i.e., ${\mathit{\varphi }}_{\mathrm{m}}=\mathrm{0.74}+\mathrm{4.7}\mathit{\zeta }$ for stable (ζ>0) and ${\mathit{\varphi }}_{\mathrm{m}}=(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{-\mathrm{1}/\mathrm{4}}$ for unstable (ζ<0) stratification (Businger et al., 1971; Dyer, 1974).

Equation (4) allows us to compute the mechanical production term in Eq. (2) as

$$\begin{array}{}\text{(5)}& -\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}z}}={\mathit{\varphi }}_{\mathrm{m}}{\displaystyle \frac{u_{*}^{\mathrm{3}}}{\left(\mathit{\kappa }z\right)}}.\end{array}$$

The buoyancy term can directly be computed from the EC measurements as well. To compute the dissipation term ϵ, we use the scaling relation of second-order structure function $D_{uu}=\stackrel{\mathrm{‾}}{{\left[u(x+r)-u\left(x\right)\right]}^{\mathrm{2}}}$ in the inertial subrange (Banerjee et al., 2015, 2016; Li et al., 2016; Salesky et al., 2013)

$$\begin{array}{}\text{(6)}& D_{uu}\left(r\right)=C_u{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{r}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where C_u≈2 (Stull, 2012) and r is the spatial lag in the longitudinal direction, which can be computed by multiplying the sampling time interval with the mean longitudinal velocity, assuming that Taylor's frozen turbulence hypothesis is valid (r=|u|Δt). The range of r where this relation is valid is found to be between 0.2 to 2 m, and ϵ is found by regression of Eq. (6). Note that the computation of ϵ is independent of any assumptions used to compute the production terms.

2.2 Research site

Figure 1 Map of Yatir Forest in Israel and locations of the measurement stations. Insets: snapshots of measurement setups. Bottom panel: topography map of Yatir Forest from Google maps. The blue arrow indicates north.

The measurements were conducted in Yatir Forest and the surrounding shrubland in Israel between 18 and 30 August 2015 as part of the “Climate feedbacks and benefits of semi-arid forests” (CliFF) campaign, a joint collaboration between Karlsruhe Institute of Technology (KIT), Germany, and the Weizmann Institute, Israel. Figure 1 gives an idea of the locations of the EC towers. Tower 1 (lat 31.375728, long 35.024262) was located in the semiarid shrubland 620 m above sea level and tower 2 (lat 31.345315, long 35.052224) was located inside the forest 660 m above sea level. The linear distance between the two locations was measured to be 4.3 km, and as can be observed from Fig. 1, there is a distinct surface heterogeneity between the two sites. The climate of the area is in between Mediterranean and semiarid, with a mean annual precipitation of about 285 mm (Eder et al., 2015). Note that the measurement sites reported in this work are different from those in Eder et al. (2015). The trees in the forest were mostly Aleppo pine (Pinus halepensis), with an average height of 10 m with negligible height variation. The surrounding land was sparsely populated by small shrubs, and in the dry season, when the measurements were conducted, was mostly free of vegetation. Thus, it is referred to as “desert” for easy distinction (Eder et al., 2015). The measurement height for the forest was 19 m above ground (9 m above the canopy height). Note that with this height selection, the measurements were conducted above the roughness sub-layer, which ends at approximately 2 times the canopy height (Harman and Finnigan, 2007). A mast was used over the desert and the measurement height was 9 m until 23 August, after which it was changed to 15 m for the remaining period. In this zone of the atmospheric surface layer, the longitudinal and crosswise velocity variances decrease logarithmically with height and the vertical velocity variance shows independence from height (Banerjee and Katul, 2013a; Marusic et al., 2013; Perry and Chong, 1982; Townsend, 1976). High-frequency turbulent data were collected at 20 Hz and 30 min averaging periods were used for both sites. After conducting quality control of the data following Eder et al. (2015), a planar fit coordinate rotation is applied to the velocity components since the data are collected on a sloped ground. The coordinate rotation following Wilczak et al. (2001) ensures that the cross stream velocity component v is zero and corrects the tilting of the anemometer with respect to the local streamlines. Moreover, a different set of coordinate rotation is applied for the desert data after 23 August.

In addition, two Doppler lidars were used at the two locations which measured vertical velocities (Brugger et al., 2018; Kröniger et al., 2018). The Doppler lidars used were StreamLine systems from HaloPhotonics. They were operated in a vertical stare mode most of the time (interrupted every half hour for less than 90 s). Technical specifications and instrument settings of the Doppler lidars are given in Table 1. The Doppler lidar at tower 1 was not working from 19 August 2015, 15:00 UTC, until 21 August 2015, 10:30 UTC and very briefly on the 23 August 2015 around 10:00 UTC due to power cuts.

Table 1 Instrument specification and settings of the Doppler lidars. From top to bottom: serial number of the forest and desert lidar, pulse length of the laser pulse at full width at half maximum, range gate length, pulse repetition frequency, number of averaged pulses for a backscatter coefficient profile, and the wavelength of the emitted laser pulse (short wavelength infrared).

Figure 2 Time series of half-hourly averages of mean speed (ms⁻¹), mean vertical velocity (ms⁻¹), friction velocity (ms⁻¹) and mean potential temperature (K) for the measurement period. The black line indicates desert, and the red line indicates forest.

Figure 3 Time series of longitudinal velocity variance (m² s⁻²), vertical velocity variance (m² s⁻²), momentum flux (m² s⁻²) and sensible heat flux (K ms⁻¹) for the measurement period. The black line indicates desert, and the red line indicates forest.

3 Results and discussion

3.1 Time series of turbulence statistics

Time series of mean speed (ms⁻¹), mean vertical velocity (W, ms⁻¹) (after applying coordinate rotation), friction velocity (u_*, ms⁻¹) and mean near-surface air (potential) temperature (T, K) for the measurement period are shown in Fig. 2. Figure 3 shows time series of longitudinal velocity variance ($\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$, m² s⁻²), vertical velocity variance ($\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$, m² s⁻²), momentum flux ($\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$, m² s⁻²) and sensible heat flux ($\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$, K ms⁻¹). The black line indicates desert, and the red line indicates forest. As noted, the desert is associated with a higher wind speed because of a lower amount of friction on the desert surface. The higher vertical velocity over the desert indicates the presence of stronger updrafts, which would be explained by higher buoyancy-driven turbulence. The friction velocity (u_*) over the forest is much higher compared to the desert, especially in the daytime, which is expected because of higher surface roughness over the forest. u_*, above both the forest and the desert, shows a strong diurnal cycle. However, there seems to be a prominent increase in u_* over the desert after 23 August. This can be attributed to the raising of the tower height. Moreover, the gentle topography around the desert could result in the strong vertical updrafts above the desert. Interestingly, the near-surface air temperatures over both the forest and the desert show a strong diurnal cycle and their differences are about 5 K on average during daytime and almost zero at night.

The longitudinal velocity variance $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$ over the forest and the desert show similar variations over time. The vertical velocity variance $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$ over the forest is higher than its desert counterpart; however, after 23 August, the levels of $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$ over the desert increase as well and become similar to the forest. This is due to changing the tower height. As the vertical profiles of $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$ are different between the desert and the forest (due to roughness length differences), the observed differences between $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$ are a function of observation height. At 15 m above the desert and 19 m above the forest floor, high enough to be in the “constant flux layer”, the vertical profiles of TKE ($\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}+\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$) converge. However, when observed at a lower elevation and below the constant flux layer, the data show clear differences in $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$.

The vertical momentum flux $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ over the forest is much higher compared to the desert, which is also expected because of the higher surface roughness of the forest, making it a much more efficient momentum sink compared to the desert. Note that the shear transport of momentum flux is still much more effective over the forest compared to the desert because of roughness effects even though the mean quantities can be higher over the desert. The sensible heat flux $\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$ over the forest is also higher, as discussed before, due to the canopy convector effect.

Figure 4 Time series of mechanical production of TKE (m² s⁻³), buoyant production of TKE (m² s⁻³), full TKE production (m² s⁻³), dissipation of TKE (m² s⁻³) and imbalance of TKE (m² s⁻³). The black line indicates desert, and the red line indicates forest.

3.2 Nature of TKE budget

Figure 4 shows the time series of the components of the TKE budget as discussed in Sect. 2.1. The first row shows mechanical production of TKE (P_Mech, m² s⁻³); the second row shows buoyant production of TKE (P_Buoy, m² s⁻³); the third row shows full TKE production (P_TKE, m² s⁻³), which is the sum of mechanical and buoyant TKE production. The fourth row shows dissipation of TKE (ϵ, m² s⁻³) and the fifth row shows an imbalance of TKE (Imb, m² s⁻³). The black line indicates desert, and the red line indicates forest. As noted in Fig. 4, the production of turbulence is mostly by mechanical or shear forcing because of the roughness of the forest, whereas mechanical production of TKE over the desert is very small and does not have a strong diurnal cycle like the forest, although it increases slightly after 23 August. On the other hand, TKE production over the desert is mostly carried by buoyancy. Buoyant TKE production is slightly larger over the forest. The buoyant TKE production over the desert is also higher after 23 August. Given the moderate temperature difference between the desert and the forest, the difference in their corresponding buoyant TKE production is interesting. It also indicates that mechanical forcing and not buoyancy makes a difference (mechanical production is higher by approximately 1 order of magnitude than buoyant production) in the turbulence generation over the desert and the forest. The diurnal cycle of the TKE dissipation ϵ is interesting as well. The dissipation of TKE seems to be higher above the forest as well compared to the desert.

A smaller TKE dissipation is recorded when the measurement location is further from the ground and above the roughness sub-layer. One strong argument for observed changes after 23 August being tower-height effects rather than a change in any large-scale forcing is that changes in the desert are observed only after the 23 August, while the forest observations maintain rather consistent dynamics.

The diurnal cycles of the TKE imbalance computed by Eq. (3) are also very interesting. The imbalance over the forest is often positive over the daytime, while over the desert it is often negative, highlighting the difference in turbulent transport and advection over the two different regimes. Also note that the positive imbalance for the forest and negative imbalance for the desert almost have a phase (anti-)synchronization, indicting that the turbulence above the forest and the desert are responsive to one another and that they are part of a coupled system, indicating again the role of the secondary circulations.

3.3 Transport of TKE over the desert and the forest

Figure 5 is used to better understand the nature of turbulent transport between the desert and the forest. Panel (a) depicts the TKE imbalance over the desert vs. the net production of TKE over the forest. As observed, there is a significant correlation (0.5) between them, indicating that the advection and transport of TKE by flux divergence and pressure fluctuations reach downstream by means of the secondary circulations and produce TKE over the forest. On the other hand, the reverse is not true, as observed in panel (b) of Fig. 5. There is little correlation between the imbalance of TKE over the forest and the production of TKE over the desert (0.14). As observed in panel (c), the production over the desert is also well correlated with the production over the forest (0.3) as both the desert and the forest are subject to the same forcing. However, the TKE production over the desert is not that well correlated with the TKE imbalance over the desert as seen in panel (d). Thus, while there should be some cross correlation in panel (a) because of desert production, that is not the only effect. The nonlocal large-scale motions contribute to the transport over the desert (without significantly altering TKE production over the desert) which in turn cause TKE production above the forest because of the higher mechanical forcing.

Thus, it can be stated that at least in the canopy sub-layer and in the atmospheric surface layer, the effects of secondary circulations are transported from over the desert towards the forest following the background wind direction, and it is not the other way around. It is worth noting here that the term “secondary circulation” has been used somewhat loosely here and contain the effects of horizontal transport as well, since partitioning the imbalance term is not possible within the scope of this campaign. In the case of transport from the forest towards the desert, it is more likely that horizontal advection is the main mechanism. The nature of the full extent of the secondary circulation are a part of a much larger flow pattern and are not fully captured by the eddy covariance towers, which only capture the fine-scale turbulence. To reveal the full nature of the secondary circulations, one can look at lidar observations as shown in Appendix B (Brugger et al., 2018) as well as large eddy simulations (Kröniger et al., 2018).

Figure 5 (a) TKE imbalance for the desert vs. TKE production for the forest. (b) TKE imbalance for the forest vs. TKE production for the desert. (c) TKE production for the desert vs. TKE production for the forest. (d) TKE production for the desert vs. TKE imbalance for the desert. Significance: 0.05 level.

3.4 Effect of nonlocal motions

Figure 6 shows the time series of the triple moments $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{u}^{\prime }}$, $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{w}^{\prime }}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{T}^{\prime }}$ in the first three rows. The vertical velocity skewness term $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{w}^{\prime }}$ (second row) is of importance as it appears in the transport term of the TKE budget (Eq. 2) and is a measure of non-Gaussian turbulence, which indicates the presence of nonlocal coherent motions such as sweeps and ejections. Note that the vertical velocity skewness is often negative above the canopy, which is consistent with the generic feature of canopy turbulence (Chamecki, 2013; Dias-Junior et al., 2015; Kaimal and Finnigan, 1994). The daytime vertical velocity skewness over the desert is often positive, indicating again the presence of nonlocal coherent structures active over the desert. The measure of skewness increases over the desert after 23 August, which is also due to the height change. The other two terms $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{u}^{\prime }}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{T}^{\prime }}$ are also associated with turbulent transport of momentum and heat as evident from their respective budget equations (Banerjee et al., 2017b; Cava et al., 2006; Katul et al., 2013; Raupach et al., 1986; Zhuang and Amiro, 1994).

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{w^{\prime }{u}^{\prime }}}{\partial t}}=\mathrm{0}=-\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{U}}{\partial z}}-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{u}^{\prime }}\right)}{\partial z}}+{\displaystyle \frac{g}{\stackrel{\mathrm{‾}}{T}}}\stackrel{\mathrm{‾}}{u^{\prime }{T}^{\prime }}{\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left(\stackrel{\mathrm{‾}}{u^{\prime }{\displaystyle \frac{\partial p^{\prime }}{\partial z}}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}}{\partial t}}=\mathrm{0}=-\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{T}}{\partial z}}-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{T}^{\prime }}\right)}{\partial z}}+{\displaystyle \frac{g}{\stackrel{\mathrm{‾}}{T}}}\stackrel{\mathrm{‾}}{T^{\prime }{T}^{\prime }}{\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left(\stackrel{\mathrm{‾}}{T^{\prime }{\displaystyle \frac{\partial p^{\prime }}{\partial z}}}\right).\end{array}$$

Figure 6 Top three panels: time series of triple moments $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{u}^{\prime }}$, $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{w}^{\prime }}$ (m³ s⁻³) and $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{T}^{\prime }}$ (K m² s⁻²). Bottom two panels show the integral timescales of horizontal (In_u) and vertical velocities (In_w) in seconds.

Moreover, the triple moments have been shown to be directly correlated with the relative contributions of nonlocal events such as sweeps and ejections (Banerjee et al., 2017b; Cava et al., 2006; Katul et al., 2013; Nakagawa and Nezu, 1977; Raupach et al., 1986). Note that momentum transport term $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{u}^{\prime }}$ is also opposite in sign for the desert and the forest, and it shows a strong diurnal cycle. After 23 August, an increase in momentum transport is noted for the desert. However, the diurnal cycle of the heat transport term $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }{T}^{\prime }}$ is not as strong as its momentum counterpart, but it is often found to be larger over the desert compared to the forest, consistent with the findings from the TKE budget that show heat is transported from over the desert towards the forest. It is, however, important to note that the structures are representative of the fine-scale turbulence and not directly representative of the large-scale circulation structures spanning the whole boundary layer. The fourth and fifth rows of Fig. 6 show the time series of the integral timescale of horizontal (In_u) and vertical (In_w) velocity components in seconds. In_u and In_w for every half hour time period are computed by integrating the normalized autocorrelation function of u and w until the first zero crossing (Kaimal and Finnigan, 1994). They can be interpreted as the characteristic timescale of the most energetic eddies in each direction. As noted in Fig. 6, timescales in the horizontal directions are larger compared to the vertical direction. More interesting is the observation that the integral timescales for the eddies above the desert are larger than those above the forest, which also increase after 23 August. This is another indicator of buoyant production of turbulence, which generates larger eddies than shear production.

4 Conclusions

We studied the nature of turbulent transport over a well-defined surface heterogeneity, comprising a desert and forest in the Yatir semiarid area in Israel. Eddy covariance and Doppler lidar measurements were conducted for 12 days between 18 and 31 August 2015 over two locations in the forest and the shrubland (referred to as “desert” because of the almost complete lack of vegetation during the observation period). Earlier campaigns in this area focused on energy balance closure and hypothesized that there are secondary circulations because of surface heterogeneity. The present work was aimed to study the nature of turbulent transport over the forest and the desert in more detail to address the following questions:

1. How does Yatir Forest affect the boundary layer dynamics such as eddy size distribution, boundary layer height and diurnal variations in turbulent statistics and fluxes compared to the surrounding desert?

2. Can the existence of secondary circulation be confirmed?

3. Is there any horizontal energy transport between the forest and the desert and how does it vary with time?

To answer the abovementioned questions, we computed half hour average turbulent statistics for both the desert and the forest and looked at their diurnal variations. We also computed individual components of the turbulent kinetic energy (TKE) budget and argued that the turbulent transport of energy should be contained in the imbalance of the TKE budget, which consists of the effects of advection, transport by turbulent flux divergence and pressure velocity interactions, since we could not compute those terms explicitly. Moreover, we also computed triple moments, which are associated with nonlocal motions and coherent structures, and integral timescales, which are associated with the most energetic eddies. The findings to the questions are listed below.

1. The forest is found to be associated with a higher level of turbulent intensity because of higher roughness although the desert had higher mean speeds and vertical updrafts, possibly due to the presence of secondary circulations. Gentle topography around the desert might contribute to the updrafts over the desert as well. The higher roughness of the desert is also responsible for higher wind speeds above the desert. There is little air temperature difference between the desert and the forest, although the mean velocities and temperature have strong diurnal cycles. Momentum and heat flux are also found to be stronger above the forest. The presence of the secondary circulation enhances the turbulent fluxes as well as the turbulent intensity above the desert.

2. The role of secondary circulations can be better understood once the components of the TKE budget are studied. Over the forest, the production of turbulence is mechanical, while over the desert, TKE production is mostly carried by buoyancy. The forest is more efficient in dissipating TKE as well. The imbalance of TKE is taken as the indicator of TKE transport and is found to vary diurnally almost anti-synchronously over the desert and the forest, confirming the role of a secondary circulation. The TKE budget is closed better over the forest compared to the desert. Turbulent triple moments, which are indicators of nonlocal motions and coherent structures, also show strong variability over the desert and are opposite in signs also confirming the role of secondary circulations. The integral timescales are found to be greater over the desert compared to the forest. This suggests that the secondary circulations that transport energy are more active over the desert – however, they cannot produce much turbulence over the desert since they only rely on buoyancy-driven turbulence as mechanical forcing is missing over the desert. This is also highlighted by the fact that mean velocities are higher above the desert while turbulent fluctuations are higher above the forest.

3. To elucidate the role of horizontal transport between the desert and the forest, we studied the correlation between the TKE imbalance over the desert and the TKE production over the forest. The moderately high correlation suggests that the secondary circulation is transported from over the desert towards the forest, enhancing TKE production over the forest, at least in the canopy sub-layer and the atmospheric surface layer. The low correlation between the TKE imbalance over the forest and TKE production over the desert confirms the directionality of this horizontal exchange, which is from the desert towards the forest and not the other way around.

To summarize, we have examined the existence and role of secondary circulations that exists because of large-scale surface heterogeneities and possible due to some topography effects between the desert and the forest by looking at proxy quantities computed from turbulence measurements. Although the campaign was conducted at a particular site, the conclusions drawn are fairly general and can be extended to other scenarios involving surface heterogeneities, such as urban landscapes and agricultural fields. Future work will attempt to highlight a more spatially detailed picture of the turbulent structure under the interesting scenario of secondary circulations and horizontal energy transport.

Appendix A: Full form of the TKE budget

$$\begin{array}{ll}\text{(A1)}& {\displaystyle }& {\displaystyle \frac{\partial e}{\partial t}}+U{\displaystyle \frac{\partial e}{\partial x}}+V{\displaystyle \frac{\partial e}{\partial y}}+W{\displaystyle \frac{\partial e}{\partial z}}={\displaystyle \frac{g}{T}}\left(\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}\right){\displaystyle }& {\displaystyle }-\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial U}{\partial x}}-\stackrel{\mathrm{‾}}{v^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial V}{\partial x}}-\stackrel{\mathrm{‾}}{w^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial W}{\partial x}}\\ {\displaystyle }& {\displaystyle }-\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial U}{\partial y}}-\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial V}{\partial y}}-\stackrel{\mathrm{‾}}{w^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial W}{\partial y}}\\ {\displaystyle }& {\displaystyle }-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial U}{\partial z}}-\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial V}{\partial z}}-\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial W}{\partial z}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u^{\prime }e}\right)}{\partial x}}-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{v^{\prime }e}\right)}{\partial y}}-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }e}\right)}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u^{\prime }{p}^{\prime }}\right)}{\partial x}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{v^{\prime }{p}^{\prime }}\right)}{\partial y}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }{p}^{\prime }}\right)}{\partial z}}-\mathit{ϵ}.\end{array}$$

Thus, to be consistent, with Eq. (2), all the terms in Eq. (A1) that cannot be evaluated using one-point measurements can be clubbed in the imbalance term, which can be described by

$$\begin{array}{ll}\text{(A2)}& {\displaystyle }& {\displaystyle }\mathrm{Imbalance}={\displaystyle \frac{\partial e}{\partial t}}+U{\displaystyle \frac{\partial e}{\partial x}}+V{\displaystyle \frac{\partial e}{\partial y}}+W{\displaystyle \frac{\partial e}{\partial z}}{\displaystyle }& {\displaystyle }+\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial U}{\partial x}}+\stackrel{\mathrm{‾}}{v^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial V}{\partial x}}+\stackrel{\mathrm{‾}}{w^{\prime }{u}^{\prime }}{\displaystyle \frac{\partial W}{\partial x}}\\ {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial U}{\partial y}}+\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial V}{\partial y}}+\stackrel{\mathrm{‾}}{w^{\prime }{v}^{\prime }}{\displaystyle \frac{\partial W}{\partial y}}\\ {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial V}{\partial z}}+\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial W}{\partial z}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u^{\prime }e}\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{v^{\prime }e}\right)}{\partial y}}+{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }e}\right)}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u^{\prime }{p}^{\prime }}\right)}{\partial x}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{v^{\prime }{p}^{\prime }}\right)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{w^{\prime }{p}^{\prime }}\right)}{\partial z}}.\end{array}$$

Figure A1 Vertical mean velocity profile averaged from 18 to 29 August (only times with both instruments simultaneously online and the nearest three range gates are discarded). Left to right: 4 h window centered on noon, daytime (sunrise to sundown) and nighttime (sundown to sunrise). The forest is shown as a solid red line, the desert as solid black line and a vertical line at $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$ as a black dashed line. Note that near the surface, the desert always has larger w, but only during the noontime with the updrafts of the forest is there a change in sign.

Thus, if no assumptions or idealizations are invoked, the imbalance of the commonly used operational TKE budget (Eq. 2) consists of TKE tendency, advection, shear production, TKE flux divergence and pressure velocity interactions. Using an array of sonics in each direction will enable determination of all these terms. However, as evident from the myriad of terms contributing to the imbalance, it is difficult to determine what degree of assumptions of homogeneity in which direction are sufficient so that certain terms can be ignored. Thus, unless all terms in Eq. (A2) can be determined, it is easier to stick to the most idealized form of Eq. (2) and treat all other terms as imbalances. Future work will try to determine the partitioning of advection, flux divergence and the other shear production terms contributing to TKE budget imbalances in the presence of heterogeneities.

Appendix B: Further evidence of secondary circulation

Figure A1 shows mean vertical velocity W above the forest and the desert averaged over all observations using the Doppler lidars. Secondary circulation cannot be thought of as a single large rotational system spanning the desert and the forest; rather, it is a much more complex and three-dimensional structure. Close to the surface layer and the canopy sub-layer, the transport of energy is indeed from the desert to the forest (Fig. 5). Further, we observe that the desert has more updrafts and the forest has more downdrafts close to the surface. However, as we go up above roughly 100 m, this behavior flips. Lastly, Kröniger et al. (2018) found in his simulations that large rotational systems developed at specific locations connected to surface features. Therefore, we conclude that the bulk transport in the convective mixed layer by a secondary circulation is from the forest to the desert, but it is advected with the mean wind and heavily influenced by surface features on a smaller scale than the forest itself.