Understanding how energy in the form of sensible and latent heat is exchanged between the biosphere and the atmosphere is of great importance in different fields. For example, it is critical for weather forecasting and climate modelling [1][2][3][4], understanding of CO 2 sequestration by plants [5][6][7][8][9], and developing management recommendations for pastures, croplands, and forestry that enable efficient use of water resources [10][11][12][13][14][15].
Traditionally, the eddy covariance (EC) method is the approach used to measure the momentum, energy, and mass fluxes between the earth's surface and the atmosphere [16,17]. This is the most direct and non-destructive method to quantify momentum, energy, and mass fluxes between a given ecosystem and the atmospheric boundary layer (ABL) [18]. Nonetheless, it has been repeatedly found that when experimentally measuring the surface energy balance (SEB, balance between the energy reaching/leaving the ground surface through net radiation and the corresponding ground and turbulent fluxes, storage, and metabolic terms) there is a 10-30% imbalance, resulting in an energy gap in the SEB [19][20][21].
Multiple possible causes for this gap have been investigated over the years. Some of them are instrument errors, including the systematic error of sonic anemometers or humidity measurements [22][23][24][25][26][27], systematic error in measurements of other SEB components like soil heat flux or radiation [28][29][30], footprint mismatch [31], and heat storage in tall vegetation canopies [32][33][34][35]. These error sources have been progressively addressed by improving the measurement techniques and the development of correction methods that can be applied during data post-processing [36][37][38]. Nonetheless, besides all these significant efforts there remains an important SEB gap [37,39,40].
The EC method relies on very high temporal resolution measurements of the three-dimensional wind speed and any other additional scalar of interest. For example, if one is interested in measuring the sensible heat flux, then the temperature would be the additional scalar of interest. The sensible heat flux is then calculated as the covariance of the vertical wind speed and the temperature around the average for a defined time period of typically 30 minutes [17,41]. As a result, the EC method can only capture the turbulent contribution of the energy fluxes, where the turbulent fluctuations are defined around the adopted averaging period [42]. Initially, this was thought to be sufficient given that almost all atmospheric transport in the boundary layer is considered of turbulent nature [43,44].
However, more recently, it was found that, under certain atmospheric conditions, a significant part of the energy is not transported by turbulent motions but rather large-scale persistent atmospheric circulations that contribute to the vertical mean wind and reach far into the atmospheric surface layer [45,46]. This transport of sensible and latent heat by secondary circulations can significantly contribute to the SEB non-closure [30,40]. However, because it is expressed through a mean advective transport in the differential equations, it can only be captured through a spatial array of sensors and not by single-tower EC measurements [47][48][49]. At present, it is possible to differentiate between two types of secondary circulations. The first type are the so-called thermally-induced mesoscale circulations (TMCs) which result from heterogeneous surface forcing and are therefore spatially bound [30,40,50,51]. The second type are slow-moving turbulent organized structures (TOSs) that develop randomly even over homogeneous surfaces [52,53] and can translate with time. Recent data analysis has shown that extending the averaging period to several days instead of 30 minutes could almost close the energy balance gap for some sites [37,54] but not for all [41]. This could be explained because the TMCs are bound to the surface and thus do not move over time [40,53,55]. Moreover, such long averaging periods typically violate the stationarity requirement that has to be fulfilled to calculate a covariance [37,40].
Multiple approaches to correct for the SEB non-closure have been developed already, e.g. by extending the averaging period [37,41,54] applying the Bowen ratio of the measured turbulent fluxes to the missing dispersive fluxes [38], attributing the entire residual to the sensible [56] or latent [57] heat flux, or modelling the energy balance gap [58][59][60]. While some of these correction methods have proven to improve SEB closure [61][62][63], these models do not consider the factors and processes that cause the SEB gap. Some approaches consider the influence of atmospheric stability or heterogeneity in surface roughness [59,60], but they do not take into account the influence of thermal surface heterogeneity.
We hypothesize that it is possible to overcome the SEB non-closure problem by considering both, the influence of thermal landscape heterogeneity, and the effect of atmospheric stability. Our study expands beyond the earlier LES works of De Roo et al. [58], and Margairaz et al. [64]. Specifically, we use the correction method developed by De Roo et al. [58] that models the SEB non-closure as a function of the atmospheric stability factor u�/w� (here w� is the Deardorff velocity) and take it one step further by including the effect of landscape heterogeneity. For this second step, we use the thermal heterogeneity parameter defined in Margairaz et al. [64]. The use of the LES technology is ideally suited to investigate the influence of atmospheric stability and surface heterogeneity on the SEB gap because it allows the control of both, the atmospheric conditions, and the surface characteristics. This facilitates the development of idealized analysis that can later shed light on the datasets of more complex field experiments [53,65,66]. Furthermore, LESs provide information on the structure of the atmospheric flow as a function of time, and the contribution of turbulent and advective transport of latent and sensible heat fluxes at each point in space [53,65].
This paper is organized as follows. In the next section, we provide a brief overview of former LES-based energy balance closure approaches, the two studies by De Roo et al. [58] and Margairaz et al. [64], and the theory underlying our new model. This is followed by a description of the dataset and study cases. Afterwards, we present the resulting reference models and our new model, which are then further discussed. The last section provides a short summary of our findings.
Several field studies have investigated EC measurements at multiple sites to understand the systematic behavior of the SEB closure, and have found relations with surface inhomogeneity [60,[67][68][69][70], friction velocity u� [19,20,71,72], and atmospheric stability [19,20,72,73]. Also, large-eddy simulation (LES) studies confirm the dependence of the SEB gap with surface heterogeneity [74], u� [65] and atmospheric stability [59,75]. The relation between the SEB gap and surface heterogeneity can be explained as follows: the patches in heterogeneous surfaces heat up differently, which favors the formation of TMCs in addition to TOSs, with the amplitude and size of the individual surfaces conditioning how strongly these TMCs will be [53,66,76,77]. There is also a causal relation between the SEB gap and atmospheric stability: a large horizontal geostrophic wind speed, i.e., neutral to stable atmospheric stratification, results in enhanced horizontal mixing, which is why the influence of TOSs and TMCs on the measured flux is less pronounced than under free convective conditions [65,78].
At present, there exists only a reduced set of approaches to model the SEB closure based on the underlying processes by considering the factors that determine the magnitude of the energy balance gap such as atmospheric stability or surface heterogeneity. One of them is the model of Huang et al. [59] that depends on u� and w�, the measurement height z, and the atmospheric boundary layer height z i . This model is applicable to 30-min flux measurements, but it was only developed for homogeneous surfaces, and only heights between 0.3 and 0.5 z/z i were considered, so it is not applicable to typical EC measurement heights within the surface layer [58].
Another model is the one of Panin and Bernhofer [60]. They developed a heterogeneitydependent energy balance gap parametrization that depends on changes in surface roughness and a corresponding heterogeneity length scale. However, this model does not include the effect of thermal heterogeneity [45,53]. Furthermore, it does not account for the effect of changing atmospheric conditions [30,50] and only provides the average energy balance closure for a site. As a result, it is rarely applicable to 30-min flux measurements.
De Roo et al. [58] developed a parametrization for the SEB gap within the surface layer that results from the energy transport by TOSs. They use the so-called imbalance (I) as a suitable measure of the missing part of the energy fluxes, i.e., the advective and dispersive fluxes that do not contribute to the Reynolds flux [58,59]. It is based on the flux balance ratio that is computed as the Reynolds sensible heat flux H divided by the total available sensible heat flux, which equals the surface flux H s at the bottom of the domain, and defined as
Following the findings of Huang et al. [59], De Roo et al. [58] assumed that the underestimation of the heat fluxes, i.e. the imbalance I can be described by a function of the non-dimensional scaling parameter u�/ w�, as well as a function of the measurement height z relative to the boundary layer height z i :
To determine the shape of functions F 1 and F 2 , they developed a LES dataset of ABL flow over idealized homogeneous surfaces using PALM [79]. They considered nine combinations of atmospheric stability and Bowen ratios (Bo) with a vertical grid spacing of only 2 m to investigate the energy imbalance at a height of 0.04 z i , i.e. within the atmospheric surface layer where most EC stations around the world are employed [19,80].
They found that combining two sets of scaling functions described well the imbalances in the sensible and latent heat fluxes. Specifically, they found that the sensible heat flux imbalance within the surface layer can be described with
and
The thermal heterogeneity parameter of Margairaz et al. [64] As part of the Idealized Planar Array study for Quantifying Spatial heterogeneity (IPAQS) [70], Margairaz et al. [81] developed a set of idealized LES of convective boundary layers over homogenously rough surfaces with embedded thermal heterogeneities of different scales. In their work, a wide range of mean geostrophic wind was implemented to vary the flow characteristics from inertia driven to buoyancy dominated. The goal of the study was to determine under what flow conditions TMCs are formed and to unravel the relation between the surface heterogeneity length scales and the dynamic length scales characterizing the TMCs. In their work, the authors show how TMCs express through mean advective transport of heat, which when unresolved either due to coarse numerical grid resolution or coarse experimental distribution of sensors can then be equivalently expressed through dispersive fluxes [81]. Furthermore, in their work, a scaling analysis between the vertical mean momentum equation and the continuity equation lead the authors to a non-dimensional parameter, referred to therein as heterogeneity parameter, that was shown to scale well with the contribution of dispersive fluxes when normalized by the turbulent fluxes [64].
The thermal heterogeneity parameter developed therein not only depends on the horizontal heterogeneity length scale L h but also on the length scale characteristic of the TMCs, L d which also depends on buoyancy and the mean horizontal wind speed. Specifically, the thermal heterogeneity parameter was defined as
where T s is the surface temperature and ΔT is the amplitude of the surface temperature heterogeneities, calculated from the absolute deviations (indicated by the vertical bars) of T s from the averaged T s following
The angular brackets denote horizontal averaging over the entire domain and the overbars denote temporal averaging over 30 minutes. Interestingly, the heterogeneity parameter can also be interpreted as a modified Richardson number, representing a balance between the mean buoyancy forces developed by the thermal surface heterogeneities, and the inertia forces represented by the geostrophic wind that tend to blend the surface effects.
In this work, we will revisit the scaling relation from De Roo et al. [58] developed for homogeneous surfaces, and the one from Margairaz et al. [64] for heterogeneous surfaces, and we will illustrate how they complement each other and can be generalized to a single relation valid for both, TMCs and TOSs. Results from the work presented herein will therefore lead to a generalization of the correction scaling relation for the closure of the SEB presented initially in De Roo et al. [58].
To our knowledge, no existing approach considers both the influence of atmospheric stability and thermal surface heterogeneity on the magnitude of the SEB gap. The SEB model based on the atmospheric stability of De Roo et al. [58] and the thermal heterogeneity parameter developed by Margairaz et al. [64] proved to capture the changes in the magnitude of the SEB very well. We hypothesize, that combining their findings in one model will lead to a very powerful tool to parameterize the SEB gap in EC measurements. This new model could then be applied to various combinations of atmospheric stability and surface heterogeneity found in numerous eddy covariance measurements worldwide.
In this work, the focus is placed on the atmospheric surface layer (ASL) because eddycovariance measurements are typically carried out close to the ground, within the surface layer [19,80]. Correspondingly, the analysis is carried out at the height of z = 0.04 z i , which corresponds to 52-59 m above the surface in our simulations. We calculate the imbalance ratio as defined in De Roo et al. [58] following Eq 1. Specifically, the turbulent flux, H, is calculated using the 30-min averaged values of vertical wind speed w and temperature θ, as well as the 30-min averaged temporal covariance of w and θ and the subgrid-scale contribution H sgs ,
The overbars indicate temporal averaging and the angled brackets denote horizontal averaging over the entire extent of the domain. In contrast to De Roo et al. [58], we therefore use the horizontally averaged imbalance instead of the local one. The sensible surface heat flux at the ground H s corresponds to H at the lowest grid point (dz/2).
To parametrize the imbalance, we first produce a set of reference models by adapting the existing model of De Roo et al. [58] to each heterogeneity scale in our dataset as described in the following subsection. Then, we proceed with developing the new model by including another scaling function that accounts for the influence of heterogeneity.
Parametrization of the imbalance with respect to atmospheric stability (reference models). First, we adapt the existing model of De Roo et al. [58] to each of the datasets to obtain a benchmark for our new model. This results in four F 1 scaling functions that are similar to the scaling function presented in De Roo et al. [58], but represent one heterogeneity case, respectively. Following De Roo et al. [58], we factorize the imbalance following Eq 1, assuming that the imbalance can be described by two scaling functions that are functions of the stability parameter u�/w� and the normalized measurement height z/z i . Based on the findings by De Roo et al. [58], we first assume that F 1 is an exponential function of the form F 1 = a exp(b u�/ w�) + c, and F 2 is a linear function of the form F 2 = i z/z i + j, where a, b, c, i, j, are fitting constants. Thus, we first fit F 1 to each of the simulation sets, individually, and later, we fit all of them onto a single F 2 function to observe their collapse on a unique curve.
For this analysis, we calculate the friction velocity u� and the Deardorff velocity w� directly using the 30-min averaged covariances as it would be done with experimental data obtained from EC systems. Thus, we calculate u� following
where u and v are the horizontal wind speeds in x-and y-direction, and w� following
where g is the gravitational acceleration (9.81 m s -2 ). Here, z i is determined as the height at which the total sensible heat flux crosses the zero value prior to reaching the capping inversion.
The resulting set of four F 1 scaling functions for each of the datasets and one F 2 scaling function for all of the datasets is then used as a benchmark for our new model and referred hereafter as reference models.
To consider the effect of surface heterogeneity, we assume that instead of describing the imbalance with a different scaling function F 1 for each set of simulations, it is possible to use the scaling function that describes the imbalance in the simulations with a homogeneous surface, F 1,HM , and add another scaling function, F 3 , that accounts for the heterogeneity:
where H is the thermal heterogeneity parameter introduced in Margairaz et al. [64] (Eq 4).
After analyzing the relationship between I normalized with F 1,HM , we assume that the relationship between I/F 1,HM and H is of linear nature and fit I/F 1,HM to F 3 = m H + n. Once F 3 is found, we proceed to identify the new F 2 similarly to the previous section.
The data used in this study was originally developed in the computational work of Margairaz et al. [64,81]. They used the pseudo-spectral LES approach that was first introduced by Moeng [82] and Albertson and Parlange [83] and further developed by Bou-Zeid et al. [84], Calaf et al. [85], and Margairaz et al. [86]. The data consists of a set of numerical simulations of a characteristic ABL developed over a homogeneously rough and flat surface. The simulations represent an idealized dry ABL, forced through a geostrophic wind at the top with Coriolis force, and an imposed surface temperature at the bottom of the domain.
Study cases include a set of simulations with homogenous surface temperature (referred hereafter as HM) and a second set with heterogeneous surface temperature distributions (referred hereafter as HT). In both sets, the geostrophic wind is varied between 1 m s -1 to 15 m s -1 . For the set of heterogeneous surface temperature conditions, the corresponding length scale of the characteristic surface heterogeneities is also varied, considering cases with 800 m, 400 m, and 200 m patches (referred hereafter as HT200, HT400, HT800, see Fig 1). In this case, the surface temperature variations are randomly distributed following a gaussian distribution with a standard deviation of ± 5 K and mean temperature equal to that of the homogenous cases, namely 290 K.
In all cases the surface temperature is initialized at a temperature of 5 K higher than the air temperature to promote the development of a convective boundary layer. All simulations have a domain size (l x , l y , l z ) = (2π, 2π, 2) km with a horizontal grid-spacing of Δx = Δy = 24.5 m and a vertical grid-spacing of Δz = 7.8 m, resulting in (N x , N y , N z ) = (256, 256, 256) grid points. At the bottom boundary, the surface heat flux is computed from the imposed surface temperature θ s using Monin-Obukhov similarity theory. In all cases, the initial boundary layer height z i was set to 1000 m by applying a capping inversion of 0.012 K m -1 . While θ s remains stable over the entire simulation time, the air temperature increases over time, leading to slightly less unstable atmospheric conditions over time. However, this effect was found to be negligible over the short duration of the simulation [64].
In total, 28 simulations were performed with different atmospheric conditions, controlled by seven different geostrophic wind speeds (i.e. U g = 1, 2, 3, 4, 5, 6, 9, 15 m s -1 ) for each set of homogeneous and heterogeneous surface conditions. In the simulations, the Coriolis parameter was set to 10 -4 Hz, representative of a latitude of 43.3˚N. Also, the roughness length was set to 0.1 m for all simulations, and the used thermal roughness was set to 1/10 z 0 following [87]. More details on the numerical simulations can be found in the original work of Margairaz et al. [64].
For the analysis presented in this work, we use statistics over a 30-minute interval recorded after 4 hours of spin-up time.
As described previously, we first fitted the exponential function F 1 to each one of the simulation cases resulting in four different sets of parameters, shown in Table 1 for the scaling function F 1,h :
Note that for each simulation case, there exist seven data points corresponding to the changes in geostrophic forcing and hence to different thermal stratification. These four different fits describe the imbalance ratio for each surface heterogeneity condition as a function of the non-dimensional term u�/w�.
The values calculated for u�/w� are shown in Table 2 where all relevant parameters characterizing the simulations are summarized. Fig 2A shows that these fitted functions for I collapse into the same value of roughly 6% of H s under less unstable conditions (u�/w� > 0.4). Only for HT800, the imbalance (I) settles at around 8% of H s the total available flux for the weaker unstable conditions. Alternatively, the imbalance increases with increasing instability, with the weakest increase found in the homogeneous surface cases and stronger increases with heterogeneous surfaces. The increase also depends on the patch size, being strongest with the largest patch size.
We then normalized the imbalance ratios (i.e. Eq 1, also vertical axis in Fig 2A ) with the four different scaling functions for the respective simulations (Eq 11, Table 1
where i R is 20.05 and j R is 0.157.
Fig 3 also shows the normalized imbalances, but in this case, the scaling function that was derived for the homogeneous simulations (F 1,HM , Eq 11, Table 2) was used for all simulations.
Here, the profiles don't collapse into a single curve, but instead present a data spread, with the largest deviation found once again in the L h = 800 m configuration. Next, we investigate whether these deviations can be reduced if the imbalance (I) is normalized by F 1,HM and represented as a function of the thermal heterogeneity parameter (H). In this case, Fig 4 shows that two different linear relationships can be differentiated for those cases with weak geostrophic forcing (U g = 1 m s -1 ) and those with a more moderate or stronger wind (U g � 3 m s -1 ). We find those two groups to correspond to the formation of cellular and roll-like secondary circulations. This is shown in Fitting F 3 to the two datasets, we receive the following scaling functions:
with m c = 0.018 and n c = 0.973for u�/w� < 0.1, which is valid for all simulations where cellular structures develop (U g = 1 m s -1 ), and
with m r = 0.116 and n r = 1.07 for u�/w� > 0.14, which is valid for all simulations where roll-like structures develop (U g � 3 m s -1 ). The fit for the very unstable simulations (u�/w� < 0.1) describes the normalized imbalance with a very high R 2 of 0.996. For the corresponding fit to the less unstable conditions (u�/w� > 0.14), the R 2 value is slightly lower with 0.841. When normalizing the imbalance additionally with F 3,c or F 3,r , respectively, the vertical profiles of imbalance collapse similar to when they are normalized with different F 1 scaling functions for each stability, as shown in Fig 6 . In this case, the remaining imbalance can be described following Eq (15):
with i N = 20.2 and j N = 0.153. All characteristic variables that are relevant for our simulations are summarized in Table 2. 2).
The reference models derived by fitting one curve for each heterogeneity scale are a more direct way to parametrize the energy imbalance than the new model, as they rely on fewer assumptions. Specifically, they are tailored to each heterogeneity scale and do not rely on the additional assumption that the magnitude of the SEB gap relates to the heterogeneity scale. However, it is not practical to use as a correction method to real measurements because it is only applicable to the discrete heterogeneity scales covered in this study. This means, that for each study case, there is a need to re-derive the corresponding scaling relation. Alternatively, using the new method proposed in this work that adds a third scaling function to parametrize the imbalance as a function of the heterogeneity parameter to account for the surface characteristics which facilitates the generalization of the correction method to EC towers surrounded by landscapes featuring any characteristic heterogeneity scale. However, because our dataset only covered heterogeneity length scales up to L h = 800 m, which corresponds to 0.57 z i on average, it is questionable whether the resulting scaling function F 3 would hold for larger heterogeneity length scales. Zhou et al. [76] investigated the relation between the scale of surface heterogeneity and the SEB gap and found that the SEB gap increases with heterogeneity length scale, reaching its maximum when the heterogeneity length scale is of the order of the boundary layer height, and decreases again, with even larger heterogeneity length scales. Our results confirm that the imbalance increases with the heterogeneity scale, especially under very unstable conditions (Fig 2A, Table 2), at least up to L h = 0.57 z i which is the heterogeneity scale our study is limited to.
While the new model is very flexible regarding the landscape heterogeneity scale, it is not applicable to all atmospheric conditions. This is because we were unable to define the scaling function F 3 for the atmospheric conditions that cause sub-mesoscale circulations to form neither uniquely cellular nor roll-shaped. While Margairaz et al. [81] found that different geostrophic forcing leads to clearly cellular or roll-like structures using the roll factor defined by Salesky et al. [88], there is a transition zone in which the structures could not be clearly assigned to a cell or roll regime. In our analysis, we therefore excluded the simulations with U g = 2 m s -1 . Several studies found the transition from cellular to roll-like structures to be rather sharp, occurring somewhere between -z i /L = 4.5 and -z i /L = 45 [89] or -z i /L = 8 and -z i /L = 65 [90], or at around -z i /L = 25. Other studies have found the transition to occur more gradually with transitional structures or co-existing rolls and cells for -z i /L = 14.1 [91], or for -z i /L < 21 [92]. For better comparison, we converted u�/w� to -z i /L for our simulations using
where κ is the von Ka ´rma ´n constant (0.4) [93]. The resulting -zi/L values are shown in Table 2. For the simulations with U g = 2 m s -1 , -z i /L varies between 156.17 and 338.84, indicating that the transition to clearly roll-like structures occurs at larger -z i /L values than reported by other studies. The model presented in this study can be applied to correct field measurements under unstable and free convective atmospheric conditions with u�/w� < = 0.1 (or -z i /L > = 400) using F 3,c or u�/w� > = 0.14 (or -z i /L < = 145) using F 3,r .
To apply the correction method, a certain amount of information on the atmospheric conditions and the surrounding landscape is required. The atmospheric conditions are considered in F 1 , using u�/w� which can be calculated from the EC measurements similar to Eqs 8-9. F 2 is a function of z/z i which means that z i needs to be known, which cannot be derived from EC measurements, only. Mauder et al. [61] already tested the correction method proposed by De Roo et al. [58] using ceilometer measurements of z i . For one site where no ceilometer measurements were available, they followed the method of Batchvarova and Gryning [94] to calculate z i using radiosonde measurements of the morning temperature gradient. They found the correction method leading to a good energy balance closure, even though the radiosonde measurements were taken at a distance of 170 km. Finally, the characteristic heterogeneity parameter can be derived using remote sensing methods or already available land cover maps [60]. At permanent measurement sites with continuous flux measurements, the temperature amplitude can be derived by performing ground-based measurements of surface temperature over the different landcover types surrounding the tower. In extensive measurement campaigns, additional airborne measurements can provide information on the temperature amplitude [80]. If additional measurements are too costly, however, it is also possible to model the surface temperature using the radiation measurements and landcover characteristics [95][96][97]. What is clear from these results, is that for accurate SEB studies, the use of single point measurements is not sufficient, but obtaining spatial information of the surroundings as well as from the flow is proven to be critical. This is a strong motivation for a paradigm change in the standard single point EC measurement approaches.
To compare the performance of our newly developed model with the reference models and with the parametrization developed by De Roo et al. 2018, we computed the dispersive flux H d using the scaling functions derived by the different approaches with
The share of H and H d in H s , i.e. the total available heat flux, is shown in Fig 7 . Without any correction, H is on average 90.24 ± 4.77% of the total heat flux at 0.04 z/z i . With H d calculated using the reference models based on De Roo et al. [58], we obtain H + H d,R = 99.49 ± 0.86%. Reaching nearly 100% means that the energy balance gap is almost closed. At the same time, the standard deviation becomes significantly smaller, indicating that the method captured the deviations in the energy balance gap well. The use of our newly developed model for imbalance calculation gives similar results with H + H d,N = 99.53 ± 0.87%. This shows that the newly developed model, which is much more flexible in its application to measurements, achieves just as good results as the reference models.
Using the scaling functions defined by De Roo et al. [58] (Eqs 3-4) results in H + H d,DR = 101.28 ± 4.2%. This shows that the method of De Roo et al. [58] generally works well with our data set, but it slightly overestimates the energy balance gap on average. This correction method has already been tested on EC measurements by Mauder et al. [61] who also found the method to yield good results. Furthermore, it does not capture the deviation of the imbalance due to heterogeneity as shown in Fig 7, which is also reflected in the almost unchanged standard deviation. This is to be expected since this method was developed for homogeneous surfaces only. However, we do not recommend combining the scaling functions defined in De Roo et al. [58] and F 3 derived in this study to address the effect of the heterogeneity as it leads to a clear overcorrection with H + H d,DR,N = 106.34 ± 2.41%.
We extended the energy balance gap correction method initially developed by De Roo et al. [58] taking into account the effects of spatial surface heterogeneity onto the atmospheric flow. We compared our new model to the reference models that are based on the already existing approach. The use of the reference models resulted in sets of two scaling functions for different heterogeneity scales, respectively. This approach is the more direct way to determine the imbalance and produces very good results. However, those sets of scaling functions are restricted to the distinct heterogeneity scales used in this study, which is why this approach is not transferable to all characteristic continuously distributed heterogeneity scales of the landscape surrounding an EC system, i.e. an area of about 20 × 20 km [20,60]. Our new model proved to yield similar results and its application to real-world EC tower sites is very flexible, since a third scaling function characterizing the influence of heterogeneity was introduced. Therefore, this correction method can be used for a wide range of characteristic heterogeneity scales of a landscape surrounding an EC tower. To apply the correction method, the atmospheric stability parameter u�/w�, the boundary layer height z i , the heterogeneity scale L h , and the amplitude of the surface temperature ΔT need to be known that can be either calculated from the EC measurements together with nearby operational radiosonde measurements or by using a ceilometer, and remotely-sensed land-surface-temperature data products.
Conceptualization: Luise Wanner, Matthias Mauder.
PLOS ONE | https://doi.org/10.1371/journal.pone.0268097 June 1, 2022