Assessing the Departures from the Energy- and Flux-Budget (EFB) Model in Heterogeneous and Urbanized Environment for Stable Atmospheric Stratification

Turbulence closure schemes, besides their intrinsic theoretical importance, represent a fundamental component in the atmospheric numerical models. Among his numerous and diverse scientific contributions, Prof. Sergej S. Zilitinkevich, with his coauthors, elaborated a turbulence closure model for stably-stratified geophysical flows, the Energy and Flux Budget (EFB) model. This closure has been verified and applied on many different experimental datasets and case studies, for steady state and homogeneous conditions. Having available observational datasets for urban and suburban sites in different cities in Italy, we investigate the deviation of the observations of turbulent kinetic energy and momentum flux from the EFB turbulence closure model in heterogeneous conditions. This allows addressing and interpreting the features that induce such deviation between the model and the observations. The EFB model is then revisited including residual terms that can account for the non-stationarity and heterogeneity of the considered cases. The correction with the residual terms leads to improve the agreement between the theoretical formulations and the observed behaviour for the turbulent kinetic energy shares and for the vertical momentum flux.


Introduction
Since the seminal works by Obukhov (1946) and Monin and Obukhov (1954) and the many field experiments, carried out mainly in homogeneous terrain (see for instance Foken 2006), most of the interpretation of the data and the following atmospheric and dispersion modelling has been based on the Monin-Obukhov similarity theory (MOST). This has been leading to semi-empirical formulations for the non-dimensional gradients of mean velocity and temperature, and for the second-order moments of the fluctuations of velocity and scalars, as functions of the stability parameter ζ = z/L, where z is the height and L is the Obukhov length. The widely accepted log-linear formulation for the mean velocity U and potential temperature Θ in the stably stratified surface layer leads to a critical value of the flux Richardson number Ri f , beyond which turbulence is supposed to be damped by buoyancy. As far the observations lead to investigate in detail stable conditions, it turns out to be clear that turbulence is present beyond the critical Richardson number, and the mean profiles depart from the simple log-linear relation (see among others Cheng and Brutsaert 2005;Yagüe et al. 2006;Grachev et al. 2007). Moreover, it is evident that the spectra of the velocity components and of the scalars in almost all the actual realizations of the planetary boundary layer (PBL) have a low frequency part containing much more energy with respect to the 'ideal' Kansas spectra (Kaimal et al. 1972;Kaimal and Finnigan 1994;Mortarini and Anfossi 2015;Mortarini et al. 2016). This aspect is related to different phenomena, like waves or geometrical effects or unsteadiness of the mesoscale forcing ). Gravity waves have been early recognized, being the spectra characterized by a gap separating the waves from the small scale turbulence as discussed, for instance, by Bolgiano (1962) and Finnigan and Einaudi (1981) for the atmosphere, and Baumert and Peters (2009) for the ocean. Yet, in general submeso motions (Mahrt 2014;Cava et al. 2019;Mortarini et al. 2019) do not exhibit such a gap, leading to an intrinsic difficulty in the definition of low and high frequency ranges.
As for the revision of the MOST, since long time the need of modifying the universal similarity functions when using them in heterogeneous conditions has been discussed and remarked (Foken 2006;Wilson 2008). The related research in the field of urban boundary layer is extensive and advanced since years, important and recognized achievements have been obtained, also in the frame of international projects (see Fisher et al. 2006;Kanda 2007;Baklanov et al. 2008;De Ridder 2010;Barlow 2014;Theeuwes et al. 2019;Schmutz and Vogt 2019). Clearly, in complex geometries and urban environments, and within any kind of canopy, the applicability of the MOST finds its limitation in the roughness sublayer, adapted formulations are to be designed and the need to take into account modification of the spectra cannot be overlooked.
Facing this complex frame, different approaches to the description of the PBL turbulence have been developed. From one hand, extensions of similarity have been proposed in order to take into account the presence of submeso motions, going beyond MOST by using different relations for specific cases (Sun et al. 2012) or adding new parameters (Stiperski et al. 2021). In particular, assessing the applicability of the MOST theory in complex terrain and urban environments, identifying the departures of its formulations and investigating possible corrections, is since long an established field of research (Rotach 1999;Roth 2000;Moraes et al. 2005;Rotach et al. 2005;Fisher et al. 2006;Christen et al. 2007Christen et al. , 2009Mahrt 2010;Mortarini et al. 2013;Trini Castelli and Falabino 2013;Falabino and Trini Castelli 2017;Srivastava et al. 2020). On the other hand, the equations for the ensemble averaged second order moments (SOM) have been used as a frame to understand and model the behaviour of stably stratified PBL. Different solutions of the SOM steady state and homogeneous equations result for different parameterizations of the third order moments (in particular, the terms containing the covariances between velocity or temperature and pressure) and of the dissipation. Mellor and Yamada (1982) and followers (Nakanishi and Niino 2009;Cheng et al. 2020) focused on the budget equations, and the last authors were able to overcome the limitation due to the critical Richardson number (which is present in the original formulation of Mellor and Yamada 1982) modifying some parameterizations, related to the choice of the mixing length and/or the parameterization of the pressure-temperature covariance. The parameterization of the inter-component energy exchange terms represents a critical issue in turbulence closure models. Generalized forms of the traditional "return-to-isotropy" hypothesis by Rotta (1951) have been proposed for both stably stratified and convective turbulence, as in Zilitinkevich et al. (2007Zilitinkevich et al. ( , 2013; Kleeorin et al. (2021); Rogachevskii et al. (2022). In the turbulence closure models, the particular form of the inter-component energy exchange term does not affect the equations for the other turbulent parameters than the turbulent kinetic energy shares. Therefore, it is of utmost importance to evaluate and verify their parameterization, and possible corrections, for different conditions.
The original and fundamental contribution by Zilitinkevich was to recognize the importance of the turbulent potential energy (TPE) E P , which increases as the turbulent kinetic energy (TKE) E K decreases due to the buoyancy sink. In this perspective, no critical gradient Richardson number exists. Moreover, the experimental evidence of the increase of the turbulent kinetic energy associated with the horizontal components of the velocity as stability increases has been considered as a necessity for the correct representation of turbulence for large stability. The Energy and Flux Budget (EFB) model (Zilitinkevich et al. 2007(Zilitinkevich et al. , 2013Kleeorin et al. 2021) uses the equations for TKE and TPE, together with those for the fluxes, with some novel parameterizations to cope with the experimental evidences. In particular, its steady state and homogeneous formulation (i.e., the SOM equations without time derivative and divergence of third order moments) gives rise to a description of the stable boundary layer (SBL) which is local (in the meaning of Nieuwstadt 1984Nieuwstadt , 1985, without critical Ri, and is free from any similarity hypothesis. Moreover, the EFB model contains some empirical constants which have to be determined according to the observations: as far as the observations are representative of an ensemble average and satisfy the steady state and homogeneous conditions, the model can be adapted on such observations. The EFB model has been evaluated and verified over a number of datasets collected or elaborated in homogeneous and steady conditions, from field experiments, to data from numerical simulations with Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) models (see Zilitinkevich et al. 2013). In this work we use the model as a reference paradigm to investigate non-ideal SBLs in complex environments, dealing with observed data gathered in suburban and urban conditions. The rationale of this study is twofold: to assess the deviations of the EFB model formulations in heterogeneous conditions, then to propose an approach to adapt them to such cases. First, the EFB functions for the TKE shares are estimated on an empirical basis using the observations at the different measuring sites, then comparing them to the original formulations in Zilitinkevich et al. (2013). This comparison provides hints to interpret the differences related to the interaction of the flow with urban geometries of various degrees of complexity. Then, to account for the possible transport and redistribution of turbulence due to the heterogeneity, residual terms are added to the budget equations for the horizontal and vertical TKE "components", E x , E y , E z , as defined in the EFB model from the diagonal Reynolds stresses, recalculating the TKE shares and vertical momentum flux accordingly. The measurement sites and the main characteristics of the observed datasets are described in Sect. 2, whereas how the raw data have been selected and elaborated for the present work is presented in Sect. 3 together with their overall analysis. In Sect. 4 the first approach used to apply the EFB model in the non-homogeneous cases considered in this study (Sect. 4.1), then the modification proposed to adapt the EFB formulations to such heterogeneous conditions (Sect. 4.2) are detailed. In Sect. 5 the results of the comparison among the observations, the two approaches and the original EFB model is presented and discussed. Conclusions are drawn in Sect. 6.

The Measurement Sites
Data measured by sonic anemometers at four suburban and urban sites in three cities, in the north-west (Torino), Centre (Roma) and south-east (Lecce) of Italy, were collected and elaborated for the present study.
Torino city (220-m average altitude above sea level) lies at the north-west edge of the Po Valley in an area characterized by complex topography, surrounded by the Alps at the northern, western and southern sides (crest line about 100-km distance) and by a chain of hills at the eastern side (maximum altitude is about 700 m). The Torino measurement site (45 • 1 4 N, 7 • 38 34 E, 243 m a.s.l.), hosted at the CNR research area, is located in the southern outskirts of the city in a heterogeneous and mixed geometry and represents a typical meteorological suburban site. The height of buildings in the area ranges from 30 m at about 150 m in the north-north-east direction, down to 18 m and 4 m at about 70 m to 90 m distance in the other directions. The general conditions of the Po basin, characterised by prevailing low wind, distinguish the climatology of the site. Given the heterogeneity of the site, the roughness z 0 and displacement height z d largely vary depending on the wind direction. Based on different morphometric methods and considering all wind directions together, z 0 was found to range between 0.5 and 1 m, z d between 0.7 and 4.2 m. During a field campaign (Urban Turbulence Project, UTP, Ferrero et al. 2009;Trini Castelli et al. 2012) lasting from January 18 2007 to March 19 2008, a mast was placed on a lawn in flat terrain, surrounded by buildings and some patches of open and grassy fields. The mast was equipped with two Gill Solent 1012R2 anemometers, at 5 m and 9 m height, and one Gill Solent 1012R2A, at 25 m, recording data at 21 Hz frequency. The raw data were synchronized, archived and interpolated to obtain a dataset with a 20-Hz frequency. The percentages of low-wind occurrences estimated on hourly averages, considering wind speed lower than 1.5 m s −1 , are 92% at 5 m, 86% at 9 m and 60% at 25 m. From previous analyses (Trini Castelli et al. 2014), on the basis of two different morphometric methods and a micrometeorological approach it was estimated that the observations at 5 m and 9 m are characteristic of the roughness sublayer, while the data measured at 25 m are representative of the inertial sublayer, where the MOST holds.
The Metropolitan City of Roma is located in the central-western part of Italy. It is the largest city and most populated municipality in Italy, with more than four millions inhabitants and high-density urbanised areas around the city center. The whole area is characterised by complex orography. The Mountains of Tolfa lie to the north, while the Apennines mountain chain extends on the eastern side across the Tiber Valley. From the southeast, a valley separates the Apennines from the volcanic area of Alban Hills. The western side is characterised by low land extending for 30 km to the Tyrrhenian coastal area. Rome has a typically Mediterranean climate characterised by mild winter and relatively hot summer seasons, with monthly minimum and maximum temperatures ranging from 3.1−4.4 • C (December-February) to 26−29.6 • C (June-August). The low-level, mesoscale circulation is mainly influenced by geographic effects generating two different and alternative patterns, namely sea/land breezes and drainage flows along the Tiber valley (Petenko et al. 2011). During daytime, the typical wind direction close to the coast is from the west−south-west sector, shifting towards the north-east sector as it approaches the mountains. Two measurement sites are considered in this study, belonging to the network of the Regional Environmental Protection Agency of Lazio Region (ARPA), located in Roma city centre (Roma BNC, 41 • 54 33.541 N, 12 • 29 47.544 E, 72 m a.s.l.) and in its SW suburbs (in Tor Vergata area, Roma TVG, 41 • 50 30.170 N, 12 • 38 51.320 E, 104 m a.s.l.). Roma BNC is an urban station, located on ARPA Lazio headquarter roof, a six-story building in the historical centre of Roma. In Roma BNC the roughness z 0 and displacement height z d strictly depend on wind directions, reflecting terrain and buildings heterogeneity around the measurements site, and vary in the range 1.04−1.92 m and 11.99−17.77 m, respectively. Roma TVG station is hosted by CNR-ISAC (Gobbi et al. 2019) and is approximately 15 km far from Roma BNC, in a semi-rural area close to suburban agglomerations, where the roughness z 0 varies between 0.18 and 0.27 m . In Roma TVG station a meteorological mast is equipped with a number of meteorological and micrometeorological sensors: among them, a three-axial ultrasonic anemometer/thermometer (USA 1, Metek, GmbH) operating at 10 Hz is deployed at its 10m top. Roma BNC station is equipped as Roma TVG one, apart from the mast height, that is 6 m above the 23-m rooftop where it is installed. At that height sonic anemometer measurements are routinely acquired at a frequency of 10 Hz. For the present study, the data are available from February 1 2013 to December 31 2014 in Roma TVG and from June 20 2013 to December 31 2014 in Roma BNC.
Lecce is a main city of the Salento peninsula in the SE of Italy, positioned between the Adriatic and Ionio Seas. No relevant orography is present in the area, yet the Apennines in the north-west act with a blocking effect so that precipitations are generally scarce. Clear skies and strong insulation usually characterize its climatology, however, during the night moisture contribution derives from the surrounding seas and the pronounced diurnal thermal/wind cycle above the ground. The most frequent wind directions on this site result to be around north and south. The first is associated to the enhanced effect of the Otranto Channel over the wind speed in anticyclonic conditions (north-west) and to cold outbreaks from the Balkans (north−north-east), while the south direction is typically associated to incoming cyclonic activity (south) and sea breezes (east-south-east). Temperature may rise up to more than 40 • C during the warm and dry season, from April to September, and keeps being mild also in autumn and winter, hardly reaching values below zero. The CNR-ISAC micrometeorological station (40 • 19 59.69 N, 18 • 7 0.86 E, 21 m a.s.l.), hosted at the Salento University campus, is located in a dismissed stone cave in a suburban area 3.5 km SW of the city, characterized by typical local vegetation -shrubs, pines and olive trees -and several buildings. Most buildings and trees have heights between 10 and 15 m. The roughness length z 0 and displacement height z d were estimated to be about 0.5 m and 7.5 m respectively (Martano 2000). The station is equipped with a 6-elements telescopic mast, on which a Gill-Solent 1012R2 ultrasonic anemometer, with a 21-Hz measuring frequency, is installed at a 14-m height above the main surface (street level). The station is mainly devoted to long-term measurements and data are automatically processed as 30-minute averages and stored and available in a web database (Martano et al. 2013(Martano et al. , 2015. For this study, 21-Hz frequency data from available raw data files from the sonic anemometer were used and processed in the same way as for the other sites. These data files include the almost complete months of January, June, July, August, September and October 2019 and part of May 2019. To depict the main topographical characteristics of the sites, in Fig. 1 maps of the imperviousness for the three urbanized areas are reported, together with the plots of the Crosswind Integrated Flux-Footprint (CIFF), estimated on the basis of the work of Kljun et al. (2015). In the Online Resource document, pictures of the sites and their areas and the distributions of the upwind distances where the maximum of the CIFF sets (CIFF-xmax) and within which the 90% of the sources influence the measurement of the turbulent flow (CIFF-x90), are provided.
The presence of artificially sealed surface can be estimated by satellite measurements of Impervious Surface Area (IMP), made available on the Copernicus Land Monitoring Service website. In more detail, the 2015 imperviousness densities satellite data (with a spatial resolution of 20 m) were used in this study, calculating the IMP percentage around each of the three measurement sites as the average of the values within a circumference of radius r = 1400 m, i.e. the value estimated in Cecilia et al. (2021) to quantify the influence of IMP on local temperature in Rome. For reference, the IMP of residential areas ranges from about 20% in sparsely built areas rising up to 60% in compact high-rise zones, based on the Local Climate Zones (LCZ) classification by Stewart and Oke (2012). The differences among the four sites can be appreciated in Fig. 1, clearly showing that Roma BNC is located fully inside a compact urban area, while Roma TVG and Torino sites are in heterogeneous suburban areas and the surrounding of Lecce site is mostly semi-rural.
The CIFF values (Table 1) reveal and confirm the specific characteristics of the three sites. Lecce and Roma TVG sites, the least heterogeneous, show very similar behaviours. They have analogous values of their CIFF-xmax and CIFF-x90 and can be considered the measurement sites that may best match the paradigm expressed by the EFB model, since their measured fluxes are expected to be representative of a rather homogeneous condition. Roma BNC location feels the effect of a dense yet homogeneous urban structure, its CIFF-xmax and CIFF-x90 indicate at the same time that the closest inhomogeneities do not affect it and that the source area is relatively large. The fluxes measured at Roma BNC are thus expected to represent a blending resulting from the interaction of the flow with the urban structures in a large portion of the city centre. Rather diverse is the situation in Torino, where the three levels of measurements are capturing different characteristics of the incoming flow, proved also by the CIFF-xmax and CIFF-x90 values. These last imply that source areas of different extension are affecting the measurements at the three heights and reveal the complex and composite structure of the local flow and turbulence. This reflects the fact that the two lower levels lie inside the roughness sublayer, while the highest one is characterising the flow in the inertial sublayer. The Torino datasets may be expected to be the farthest ones from the EFBmodel paradigm. The six sets of measurements gathered at the four sites are representing different conditions, all of them typical and frequent in complex and heterogeneous locations. This variety favours a thorough assessment of the EFB model, its related deviations and its possible adaptation to account for the heterogeneity.

Elaboration and Analysis of the Observations
For the present study, no detrending was applied to the available raw data and a 1-hour averaging time was considered. To calculate the averaged data, we retained the time intervals where the percentage of valid data is greater than 75%. The observed data were rotated into the local streamline reference system by applying a double rotation (McMillen 1988). To assess the EFB model over the four suburban and urban sites, the data selection was performed on the basis of the following criteria. In order to discard non-stationary transition periods, data relative to the night-time periods in a time range from one hour after sunset and one hour  More, night-time data relative to the stability parameter ζ < 0 have been discarded. The selection yielded to a number of 2567, 2736 and 2021 hly-averaged data for Torino site at 5, 9 and 25 m respectively, 4188 and 1025 data for Roma TVG and Roma BNC, 703 observed data for Lecce site. This approach allows investigating the possible departures of the experimental data from the theoretical EFB curves in relation to the occurrence of heterogeneous conditions, both on the horizontal direction (surface heterogeneity of the sites) and on the vertical direction (third-order moments divergence).
To provide a general description of the dynamical features of the selected data at the three sites, in Fig. 2 the wind-charts of the hourly-averaged sets are plotted, in Fig. 3 the distribution of the hourly-averaged TKE values is reported.
The wind-charts for Torino site show two prevailing wind directions, from south and southsouth-east and from north and north-north-east, and they reflect the climatological means of the area even for the selected neutral and stable sub-set of data. The dominance of low-wind conditions is clearly visible at all levels, the strongest speeds recorded at 25-m height are related to winds originated in the main valley of the area, the Susa Valley, and blowing from north-north-west to west-north-west. The TKE values distribute in a limited range and tend to increase with the height, not showing particular features that could be induced by the surrounding buildings. Under stable conditions the prevalent wind regime at Roma BNC is north-east to south-west, as a result of the combined effect of the complex orography of the city and the proximity of the stations to the narrow valley of Tiber, that splits the city in two.
The suburban site Roma TVG shows a wind direction distribution varying from east to south, which reflects the presence of the Albani Hills about 15 km south-east of the station (Ciardini et al. 2019). Compared to Roma TVG, Roma BNC TKE distribution has a heavier right tail that can be attributed to the urban heat island effect, which results in higher minimum values Fig. 3 Distribution of the TKE hourly-averaged measured values for the selected stable data: a-c Torino at 5, 9 and 25 m; d Roma BNC, e Roma TVG, f Lecce of both the buoyant and the mechanical production/loss terms. The TKE distribution for the Lecce dataset is similar to that of Torino (Fig. 3), in some sense in between the distributions at 9 m and 25 m of the Torino site. The suburban characteristics of the two sites are similar in some aspects, with some arboreal vegetation and low buildings around and comparable estimated values of z 0 . Moreover, the wind speed distribution for the Lecce data is somehow close to that of Torino (Fig. 2), in spite of the fact that the Lecce site is climatically windy if compared with the Torino site. Indeed, the choice of the stable condition for this study tends to increase the number of weaker wind conditions in the selected data for the Lecce site, as the anticyclonic channeling wind generally shows a diurnal cycle, increasing in daytime and decreasing during the night.

Revisiting the Energy-and Flux-Budget Model for Disturbed Atmospheric Surface Layer
The two approaches adopted to assess the departure from the EFB model when applying it in the urbanized area considered in this study are detailed hereafter.

Investigating the Deviation from the EFB Model of Observations in Heterogeneous Conditions
As first step, we investigated the applicability of the EFB model in the heterogeneous conditions characterising the four measurement sites, to examine the possible departure between the observations and the EFB formulations. Following (Zilitinkevich et al. (2013), Z2013hereafter) discussion on the inconsistency of the Rotta (1951) hypothesis on the return-to-isotropy, we assessed the inter-component exchange of turbulent kinetic energy E K estimating the energy shares eqs. (51) in Z2013], recalled hereafter: In order to compare the TKE shares at the different sites with the plot in Fig. 3 of Z2013, hereafter they are written through their dependency on the ζ = z/L parameter in place of the Richardson flux Ri f , using Z2013 Eq. (71), which reads: where R ∞ = 0.25 is defined in Z2013 as the upper limit for the Ri f attainable in steady-state regime of turbulence. For the parameter ζ the Obukhov length is defined as: being τ i j the Reynold stresses (i,j=x,y,z) and where the von Karman constant k is not included. This leads to the following formulations for the A i as functions of the stability parameter ζ : A y = C r 3(1 + C r ) Since the Richardson number can be estimated only for the Torino site where a profile of data is available, this approach allowed us estimating the energy shares at all sites by using the ζ stability parameter. However, in principle the model gives a relationship between Ri f and ζ tied to the existence of a log-linear velocity profile, condition that is generally not met in urban context. Thus, first we need to postulate the validity of such approximation when examining the results, then we consider a ζ -independent representation of the data to overcome it, looking at the normalized momentum flux ( τ xz E K ) 2 , being τ xz = uw , as function of the vertical TKE share. The energy shares asymptotic values for A z are given by Eqs. (52) and (53) in Z2013 [see also Appendix 1, Eqs. from (50) to (55)], so that the following formulations for the inter-component energy exchange constants can be derived: The expressions for C 1 and C 2 can be derived either by Eq. (1) or Eq.
(2) is chosen, they read: Empirical values for the constants (see Table 3 in Sect. 5) have been calculated for each dataset from the corresponding asymptotic values of the A i functions, estimated as the median of the series of the observed data at the two extremes of the ζ range and reported in Table 2, together with the original asymptotic values determined in Z2013. Both the original values for the constants as estimated in Z2013 -C r =1.5, C 0 =0.125, C 1 =0.5, C 2 =0.72 -and the values empirically determined at each site, have been used to draw the theoretical curves for the A i functions (see Fig. 5 and related discussion in Sect. 5), as done also in Kleeorin et al. (2021). As remarked above, since the EFB model determines the relationship between Ri f and ζ for log-linear velocity profiles, in general giving the results as function of ζ is affected by an approximation which effects are difficult to disentangle. Instead, when Ri cannot be estimated, a formulation of the normalized momentum flux as a function of the vertical share of the TKE can be used to evaluate the departure of the observations from the reference paradigm. In fact, following Z2013 the vertical momentum flux normalized with the TKE, E K , is among the parameters characterizing the turbulence state of the surface layer.

Revising the Budget Equations in the Energy-and Flux-Budget Model to Account for the Heterogeneity
As second step, in order to account for the unsteadiness of TKE and/or the horizontal transport due to the inhomogeneity of the TKE field and/or vertical divergence of third order terms, the EFB model has been revisited including residual terms R χ (χ = x, y, z, K , τ ) in the budget equations. Referring to Z2013, the budget equations for the diagonal terms of the Reynold stress, the TKE budget equation together with the budget equation for the τ xz Reynold stress can be written as: ( where R K = R x + R y + R z , and ∂ V /∂z = 0 was assumed, so that the shear production acts only on the x component of the velocity. The R χ terms at the r.h.s. are the sum of the total derivative of the momentum and the divergence of the third order terms which represent the velocity component covariances and pressure-velocity covariances. The terms E χ /t T parameterise the dissipation of the velocity components and of the TKE.
We use the following definitions: From Eqs. (16) and (17) we have, respectively: In a sense, our approach is analogous to the work by Chamecki et al. (2018), who focused on a reduced form of the TKE budget equation, where all terms potentially producing a local imbalance between production and dissipation are lumped together into a residual term. Chamecki et al. (2018) referred to Eq. (16) as the reduced budget equation, in which the causes of local imbalance cannot be distinguished. A positive value of R K is associated to regions in which the production term is larger than the dissipation term (see also Chamecki et al. 2020). Considering the diagonal Reynold stresses, when R x = R y = R z = 0 the boundary layer is thus in a state of local balance between production and dissipation of TKE. In this sense, in their turn the R χ (χ = x, y, z) can be considered as local imbalance terms for the boundary layer. Notice that Chamecki et al. (2018) normalize the TKE reduced budget equation using the dissipation term . Being the dissipation parameterized, in the following we shall make the residuals non-dimensional using the shear production term τ S. The formulations (49) and (50) in Z2013 are used respectively for the correlations between the fluctuations of pressure and velocity shear, Q i j , and for the part of the TKE participating in the inter-component energy exchange, E ⇔ . For reference, in Appendix 1 the Z2013 equations (20) and (21) are reported in their components in Eqs. (39) and (40), (49) and (50) are reported as (45), (46), (47) and (48).
Combining Eqs. (13), (14), (15) and (16), the modified version of Eqs. (1), (2) and (3) for the shares A x , A y and A z , can be derived as follows: where A Rx , A Ry and A Rz indicate the modified TKE shares and P χ = R χ τ S (χ = x, y, z, K ) are the normalized residual terms. Notice that the residual terms, P x , P y and P z , are not constant but depend on Ri f . When writing the system for the asymptotic values of the A i , including the P χ (χ = x, y, z, K ) residual terms, since the three equations of both systems sum up to 1, only two of the three equations are independent. Thus, two equations are available for the three unknown residual terms for the two asymptotic conditions at neutrality, Ri f = 0, and at large stability, Ri f = R ∞ . This implies that assumptions are needed to determine the three residuals. To simplify and reduce the possible arbitrariness, we combine the two horizontal shares of the TKE and the two horizontal residual terms as: where A H is the total horizontal share of TKE and R H is the residual term of the horizontal components of the TKE.
The asymptotic values for the TKE shares at neutrality, Ri f = 0, and at large stability, Ri f = R ∞ , can then be evaluated from Eq. (24) as: for Ri f = 0, and for Ri f = R ∞ .
To circumscribe and shed light on the effect of adding the residual terms to the EFB model budget equations, we chose to set the values C r = 1.5, C 0 = 0.125, C 1 = 0.5 and C 2 = 0.72 as in Z2013. The values determined by Z2013 can be considered universal for stationary and horizontally homogeneous turbulence. The Eqs. (27), (28) can be used to find an expression of the residual terms as functions of the TKE shares' asymptotic values: It is straightforward noticing that even considering the four C i constants as known, the system of equations for the residual terms (29) remains undetermined, since P z depends on P H both at Ri f = 0 and at Ri f = R ∞ . Also, notice that choosing the value P H = 1 for all stability range would lead to P z = 0, and to fixed values of A 0 z and A ∞ z , which contradicts the behaviour of the data, thus the condition P H = 1 has to be observed.
With respect to the approach used in the first step of the analysis, in this way the adaptation of the EFB model formulations to heterogeneous conditions moves from the estimation of the local empirical constants to the inclusion of residual terms that are based on mathematical functions, as explained in the following, thus offering a more rigorous and generalized methodology.
The asymptotic values for P H have to be assigned based on assumptions for known cases. In the following analysis two cases are considered. In the first case P 0 H and P ∞ H are assumed to be null, P 0 H = P ∞ H = 0, in the second case P 0 H = 0 and P ∞ H = 0 are considered and their values are estimated empirically. For each site, P 0 H and P ∞ H are chosen considering a combination of P x and P y asymptotic values that, when used in Eqs. (22) and (23), return the most reasonable agreement with the observations. Once both values of P H are given, the system (29) is closed and provides the vertical residual terms P z for neutral and very stable conditions as functions of the measured values of A 0 z and A ∞ z . All values are reported in Table 4.
As cited before, the residual values introduced in the budget equations for the diagonal Reynolds stresses and the TKE are clearly not constant and depend on the stability through the flux Richardson number Ri f , as can be seen in Eqs. (22), (23) and (24). Therefore, it is necessary to elaborate functions that provide their dependence on stability through Ri f . Based on our analysis, the following analytical formulation is proposed: in which the parameters α 1 and β 1 are chosen to satisfy the asymptotic constraints. The behaviour of Eq. (30) for a generic residual term P χ is depicted in Fig. 4a for three different values of the exponent n (n = 1, 2, 3). To draw the curves, the values P 0 χ = 0.1 and P ∞ χ = −0.02 were chosen consistent with the ones determined for P z in the Torino dataset in the case that both P z and P H were considered. The effect in taking into account the residual terms on the parameterization of the vertical share of the TKE are shown in Fig. 4c i.e.
Since τ E K > 0 and c < 0, the solution with the minus in front of the square root must be discarded. Equation (34) where γ determines the values of the momentum flux residual for Ri f = 0. Figure 4b shows the behaviour of Eq. (35) for γ = −0.15 (yellow curve), γ = −0.1 (black curve) and γ = −0.05 (light blue curve). As for the asymptotic values of the normalized horizontal residual terms of TKE shares, for each site, γ values are chosen to return the most reasonable agreement with the observations, when used in Eq. (34). These values are γ = −0.08, −0.01, −0.02, −0.1 for Torino, Roma BNC, Roma TVG and Lecce, respectively. In panel d of Fig. 4 the influence of the residual terms of the TKE shares and momentum flux residuals is shown. The dashed line refers to the Z2013 behaviour [Eq. (57)], the continuous colored lines (black, yellow, light blue) curves represent Eq. (34) where the corresponding R χ (shown in panel a with respectively the same colour) and b (shown in panel b) were substituted, the grey dotted line refers to Eq. (34) where A z takes into account the residual term (panel a) but the momentum residual term, b, is neglected. Notably, the inclusion of the momentum residual term changes the concavity of the curve predicted by Z2013 relationship. It is crucial to notice that while the reliability of plotting on ζ relies on Eq. (4), depicting the normalized momentum flux versus the vertical TKE share is independent from the choice of the stability parameter. The whole behaviour of the normalized residual terms for the TKE shares (P K , P x , P y , P z ) and for the momentum flux (b) as a function of ζ , at the different sites and for each choice of P H , is depicted in Fig. 7 of the Online Resource document.
In Appendix 2 the Eqs. (22), (23) and (24) are rewritten making explicit the terms in the TKE shares' modified equations, including the residuals, which differ with respect to the original Z2013 formulations.

Results and Discussion
As for the first step in our approach, we compare the Z2013 asymptotic values for the TKE shares with the values estimated from the hourly-averaged datasets (Table 2). For the horizontal components we notice a general differentiation between the less heterogeneous sites, Lecce, Roma TVG and Torino 25 m, and the values estimated for Torino 5 m, fully in the roughness sublayer, and for the urban site of Roma BNC, as can be seen for A 0 x and A 0 y . The vertical component is always lower for A 0 z while keeps very similar to the Z2013 A ∞ z value. The asymptotic values of Lecce dataset are the closest to the Z2013 values, followed by those of Roma TVG. This suggests that the departure from the asymptotic values as estimated in Z2013, and consequently from the related empirical constants, may be ascribed to the degree of heterogeneity characterizing the different sites.
In particular, the inter-component energy exchange constant C 0 , determining the TKE vertical share A z , takes negative values for all datasets, as seen in Table 3. We verified that this leads to an increase in the vertical component at the expenses of the horizontal transverse component. The Lecce set of constants are the least different with respect to the Z2013 set. It might be useful to recall the role of the inter-component energy exchange in the EFB model framework. C r [Eq. (9)] determines the vertical share of TKE in neutral conditions (ζ = 0), while C 0 [Eq. (10)] determines the vertical share of TKE in the very stable limit and, consequently, the decrease of A z with ζ , the larger C 0 the smaller A z . Interestingly, due to the presence of the first term on the r.h.s of Eq. (6), the longitudinal share of TKE, A x , is less sensitive than the transverse component, A y , to changes in C 0 and the energy redistribution towards the A z is mainly done at A y expenses, as found in our case, this being one of the main differences between the EFB and the Rotta model (Bou-Zeid et al. 2018). C r and C 0 regulate how energy is distributed between the TKE horizontal and vertical components, thus they also regulate their ratio. Mortarini et al. (2019) identified that for A z /(A x + A y ) < 0.1 the SBL dynamics is characterized by intense submeso activity and wind meandering takes place. Assuming R ∞ = 0.25 and the C r and C 0 estimated in Z2013 the meandering threshold is Ri f ∼ 0.2. Further, assuming R ∞ = 0.25 and imposing Eq. (8) > 0, the constraint C 0 > 1 2 1 − C −1 r is found. The C 1 and C 2 constants [Eqs. (11), (12)] distribute the horizontal share of the TKE between A x and A y , with the first accounting for the neutral limit and the second for the very stable one. If C 2 satisfies the condition (formula): horizontal isotropy, A x = A y = 1 2 (1 − A z ) is achieved (Kleeorin et al. 2021). The four inter-component energy exchange constants make the EFB model versatile and able to adapt to different asymptotic values of the TKE components. Figure 5 compares the A i functions as estimated from the local asymptotic values with the original ones from Z2013, and it highlights the different behaviour of the A i in the different sites. Looking at the functions estimated from the datasets, we find that (i) as in Z2013 in neutral stratification A 0 z is essentially smaller than A 0 y ; (ii) with strengthening stability not always A y increases and A x decreases, so there is not a general tendency to horizontal isotropy: this occurs in Roma TVG and Lecce, but in the other sites there is a mild trend to increase for A x in the strongly stable range together with an increase also for A y , or a rather constant A y behaviour as in Torino 5 m and 9 m; (iii) the vertical energy share A z decreases with increasing ζ and at ζ > 1 levels off at a similar small value as in Z2013: in general the local A z tends to be lower than Z2013 one in the neutral and weakly stable conditions, while they meet in strongly stable range, ζ > 1.
Overall, the findings in Kleeorin et al. (2021), that isotropy is achieved on the horizontal plane at large stratification, is not met for all datasets. The empirically estimated asymptotic values, and in general the trend of the local A i functions, remark the differences due to the distinctive features of the sites. The EFB model represents a steady homogeneous case, , dashed for A y and dotted for A z : black lines indicate the original Z2013 EFB curves; blue, orange and magenta indicate respectively the A x , A y and A z from observations at the site, with the corresponding point in the ζ class. The coloured ribbons cover the areas between the 5 th and 95 th percentiles from previous analysis it is shown that it may work well in less 'disturbed' sites, while the departures from it at the different sites reveal the fingerprint of the local heterogeneity.
As for the second step in the analysis, Figs. 6, 7 and 8 show the shares A z , A x and A y respectively, allowing a detailed comparison of the data with the original EFB model and the extended model including the residuals. Note that the residuals for the extended model are computed using the asymptotic values of the shares (for negligible and large stability) so that the trend at intermediate values of ζ is defined by the choice of the exponent n in Eq. (30). Concerning the choice of the residual, the extended model was applied with P H = 0 and with P H = 0, using the values reported in Table 4.
As far as the vertical share A z is concerned, overall no major differences result using the original model with modified constants (Fig. 5) and the extended one ( Fig. 6; see also  (37) and (38)]. In the Roma TVG site the transition from near neutral to very stable values of A z occurs at values of ζ smaller (approximately one order of magnitude) than those typical of the other sites and predicted by the EFB model, and is coupled with a large variability. Correspondingly, an increase of the cross-stream share is observed. Thus, the transfer of TKE from the vertical direction to the horizontal plane occurs at smaller values of ζ than in the  Fig. 7, the choices about the residuals are more critical than for the vertical share. Note that the choice P H = 0 (continuous coloured lines) turns out to give a wrong behaviour in the Torino site, predicting a decrease with increasing stability, while the data show an opposite trend. The data are quite satisfactorily described allowing the horizontal residual to vary. Large differences occur at large stability for the Torino site, suggesting that the effect of the residual terms is more important under stable conditions.
As far as the cross-stream share is concerned, Fig. 8, it may be noted that the choice P H = 0 approximates well the original Z2013 model. The data are better represented allowing P H to be different from zero and thus possibly varying with stability. This in general holds both for small and large stability conditions, stressing again the importance to account for the residual term on the horizontal component share of TKE.
Looking at Table 4, when the horizontal residual P H is null, the vertical one P z is always positive, both in neutral and stable conditions. The non-zero P 0 H values are always negative and the P 0 z always positive in the neutral limit. In the strongly stable limit P ∞ H changes its sign to positive in all the sites but two, Roma TVG and Roma BNC, and the sign of P ∞ z keeps being always the opposite to P ∞ H one. When considering the total P K = P H + P z residual, for the analysed datasets both in neutral and strongly stable conditions mostly a positive value is found. Therefore, the residual acts as a further dissipation term, subtracting energy from the turbulent system. Instead, P K negative value means that there is an input of turbulent kinetic energy which sums to production by shear and destruction by buoyancy. Such case occurs for Torino 5 m in neutral conditions (P 0 K = −0.02), for Roma BNC (P ∞ K = −0.07) and Roma TVG (P ∞ K = −0.02) in the strongly stable limit. For Torino 5 m the negative P K value is due to a higher absolute value for the horizontal residual than for the vertical, indicating a different redistribution of the energy in the roughness sublayer with respect to the other levels above it. As for the two sites which exhibit a negative P K value in stable conditions, we notice that Roma BNC is located in a fully urbanized environment, whereas Roma TVG is a suburban site but affected by winds blowing across groups of buildings for . For these three cases, the different contribution of the P K residual, bringing an input to the turbulent kinetic energy, can be due to the effect of the built environment. In other words, the interaction of the incoming flow with the obstacles upstream the measuring site lead to increasing turbulence.
Further insight in the features of the turbulence is derived by the analysis of the vertical momentum flux normalized with the TKE. In Fig. 9 the normalized momentum flux is plotted vs ζ . In near neutral conditions the data exhibit values not too different from the reference one (the black line, i.e. the value suggested by Z2013), but the ratio approaches zero at large stability, much smaller than the reference value, and decreases quickly at intermediate ζ values. The extended model with the residuals on the TKE shares but without a residual on the vertical momentum flux (i.e. b = 0) is plotted as a gray line, and is qualitatively and quantitatively far from the data. In other words, it is necessary to introduce a residual also on the momentum flux equation. The colored lines show the results according to different choices for P H but always with b = 0. It results that the data are well reproduced; in particular, the fast decrease in the 0.1 < ζ < 1 and the very low values at large stability are well gathered.
The normalized momentum flux is reported as function of the vertical share in Fig. 10. The introduction of a residual on the τ xz budget, equivalent to b = 0, (the coloured lines) with P H = 0 and P x = 0.1P H ; points: median of the observed data in the corresponding ζ bin. The coloured ribbons cover the areas between the 5 th and 95 th percentiles of the observed data is essential to modify the curvature of the extended model computed line and to cope with the observations. Remember that the left side of the x-axis corresponds to stable conditions, the right one to near-neutral ones. The stable cases exhibit very low normalized momentum flux values (as noted above) for all the sites; approaching neutrality the different sites show different behaviours, ranging from values higher than the Z2013 ones for Lecce, Torino 9 m and 25 m, the less affected by urban effects, to smaller values, evident at Roma BNC for example. Thus, the presence of obstacles leads to normalized momentum fluxes smaller than the ones observed in almost homogeneous sites, suggesting that the mechanically induced turbulence due to the obstacles increases the TKE but has a minor effect on the vertical momentum flux.
It is worth noticing that the values of the constants shown in Table 3, which are different from the EFB original ones, can be interpreted in terms of the presence of residuals. For instance, the value of C r [see Eq. (9)] can be derived from Eq. (27) and reads: The first term on the r.h.s., coincident to the original EFB expression as in Eq. (9), decreases as A 0 z decreases, and the second one is positive with the actual values of the residuals. Thus a value of A 0 z smaller than the EFB one (as occurs in our data sets, see Table 2) leads to a C r value smaller than the EFB one, when the residuals are not accounted for.
An analogous calculation can be made for the constant C 0 deriving its expression from Eq. (28) and comparing it with Eq. (10). It reads: An example of the direct comparison between the curves of the vertical TKE shares A z calculated by the original formulations and using the local constants [as in Fig. 5, Eq.(3)] or keeping the original Z2013 constants but introducing the residual terms in the budget equations [see Eq. (24)], together with the original EFB function, is reported in Fig. 8 of the Online Resource document. By summarizing, forcing the steady homogeneous version of the model to the actual data leads to changes in the constants, that compensate neglecting of the residuals. This hold essentially for the velocity shares, not for the vertical momentum The points refer to the median of the observed data in the corresponding ζ bin. The coloured ribbons cover the areas between the 5 th and 95 th percentiles of the observed data flux. In other words, the steady homogeneous version of the model is suitable to cope with different data sets, at the price of its universality.

Conclusions
The Energy-and Flux-Budget turbulence closure models, formulated by Zilitinkevich and coauthors in the seminal paper dated 2013, constitute a hierarchy of models applicable to stably-stratified geophysical flows. An extension of the EFB model to convective condition is under development (online documentation can be found in https://arxiv.org/abs/ 2112.14121v1). In steady state and homogeneous conditions this closure gives, among others, relationships describing the turbulent kinetic energy shares and the normalized vertical momentum flux as function of stability, expressed either by Richardson numbers or Obukhov length. These relationships derive from the budget equations for the variance of the velocity components and the vertical momentum flux. In the present work, by adding to each budget equation a residual term which accounts for unsteadiness and heterogeneity, new relationships have been derived, which cope with the effects of transfer and redistribution of the turbulent kinetic energy and momentum flux. The analysis has been based on observational datasets collected in four suburban and urban sites in Italy. First, the original functions describing the TKE shares as formulated in the EFB model have been used, estimating their empirical constants from the observed data. This approach allowed assessing the deviation between the original EFB model formulations and the shape they assume when following the observations. The agreement turns out to be fair, especially as far the vertical share is concerned, and highlights the flexibility of the EFB model, if we give up the universality. However large qualitative and quantitative discrepancies are observed for instance for the normalized momentum flux.
After this initial analysis, the relationships derived from the budget equations including residual terms have been applied. In this case, the original constants of the EFB model have been kept. The discrepancies deriving from the actual asymptotic values of the shares (different from site to site, and different from the original ones used to estimate the EFB constants) have been interpreted as due to the presence of non-zero residuals, which are thus computed (with a certain degree of arbitrariness) from the data. The agreement between the theoretical formulation and the observed behaviour for the turbulent kinetic energy shares and for the vertical momentum flux improves after the correction by the residuals. In particular, the normalized momentum flux dependence from the vertical share catches the characteristic concavity of the curve as seen from the observed data, which is missed using the original equations.
The analysis highlights the importance of accounting for the transfer among the different turbulent kinetic energy components when dealing with heterogeneity. This aspect was already considered by Wyngaard and Coté (1971) in the homogeneous conditions of the well-known Kansas experimental campaign. They observed that in the homogeneous case, the transport (in the present terms, the residual) was of minor importance in the turbulent kinetic energy budget.
The main lesson we got from the present analysis is that understanding the behaviour of the turbulent stable boundary layer requires to go beyond simple and well established schemes, a fundamental legacy of Prof. Sergej S. Zilitinkevich.