Incorporating the effect of heterogeneous surface heating into a semi-empirical model of the surface energy balance closure

It was discovered several decades ago that eddy covariance measurements systematically underestimate sensible and latent heat fluxes, creating an imbalance in the surface energy budget. Since then, many studies have addressed this problem and proposed a variety of solutions to the problem, including improvements to instruments and correction methods applied during data postprocessing. However, none of these measures have led to the complete closure of the energy balance gap. The leading hypothesis is that not only surface-attached turbulent eddies but also sub-mesoscale atmospheric circulations contribute to the transport of energy in the atmospheric boundary layer, and the contribution from organized motions has been grossly neglected. The problem arises because the transport of energy through these secondary circulations cannot be captured by the standard eddy covariance method given the relatively short averaging periods of time (~30 minutes) used to compute statistics. There are various approaches to adjust the measured heat fluxes by attributing the missing energy to the sensible and latent heat flux in different proportions. However, few correction methods are based on the processes causing the energy balance gap. Several studies have shown that the magnitude of the energy balance gap depends on the atmospheric stability and the heterogeneity scale of the landscape around the measurement site. Based on this, the energy balance gap within the surface layer has already been modelled as a function of a nonlocal atmospheric stability parameter by performing a large-eddy simulation study with idealized homogeneous surfaces. We have further developed this approach by including thermal surface heterogeneity in addition to atmospheric stability in the parameterization. Specifically, we incorporated a thermal heterogeneity parameter that was shown to relate to the magnitude of the energy balance gap. For this purpose, we use a Large-Eddy Simulation dataset of 28 simulations with seven different atmospheric conditions and three heterogeneous surfaces with different heterogeneity scales as well as one homogeneous surface. The newly developed model captures very well the variability in the magnitude of the energy balance gap under different conditions. The model covers a wide range of both atmospheric stabilities and landscape heterogeneity scales and is well suited for application to eddy covariance measurements since all necessary information can be modelled or obtained from a few additional measurements.

It was discovered several decades ago that eddy covariance measurements systematically underestimate sensible and latent heat fluxes, creating an imbalance in the surface energy budget. Since then, many studies have addressed this problem and proposed a variety of solutions to the problem, including improvements to instruments and correction methods applied during data postprocessing. However, none of these measures have led to the complete closure of the energy balance gap. The leading hypothesis is that not only surfaceattached turbulent eddies but also sub-mesoscale atmospheric circulations contribute to the transport of energy in the atmospheric boundary layer, and the contribution from organized motions has been grossly neglected. The problem arises because the transport of energy through these secondary circulations cannot be captured by the standard eddy covariance method given the relatively short averaging periods of time (~30 minutes) used to compute statistics. There are various approaches to adjust the measured heat fluxes by attributing the missing energy to the sensible and latent heat flux in different proportions. However, few correction methods are based on the processes causing the energy balance gap. Several studies have shown that the magnitude of the energy balance gap depends on the atmospheric stability and the heterogeneity scale of the landscape around the measurement site. Based on this, the energy balance gap within the surface layer has already been modelled as a function of a nonlocal atmospheric stability parameter by performing a large-eddy simulation study with idealized homogeneous surfaces. We have further developed this approach by including thermal surface heterogeneity in addition to atmospheric stability in the parameterization. Specifically, we incorporated a thermal heterogeneity parameter that was shown to relate to the magnitude of the energy balance gap. For this purpose, we use a Large-Eddy Simulation dataset of 28 simulations with seven different atmospheric conditions and three heterogeneous surfaces with different heterogeneity scales as well as one homogeneous surface. The newly developed model captures very well the variability in the magnitude of the energy balance gap under different conditions. The model covers a wide range of both atmospheric stabilities and landscape heterogeneity scales and is well suited Introduction 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.

Theory
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.

The atmospheric stability dependent energy balance gap model of De Roo et al. [58]
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 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].

The combination of the atmospheric stability and thermal heterogeneity parameters into a new model
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. Parametrization of the imbalance with respect to atmospheric stability and surface heterogeneity (new model). 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.

Dataset and study cases
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.

Reference models
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). Results are then represented in Fig 2B as 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.

New model
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 Table 2. Overview over characteristic variables that are relevant for this study, including geostrophic wind speed U g , boundary layer height z i , the friction velocity u�, the Deardorff velocity w�, the atmospheric stability parameters u�/w� and -z i /L, the heterogeneity parameter H, and the energy imbalance I for each simulation. 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 Fig 5 where xy-cross-sections of the 30-min averaged vertical wind speed w for combinations of U g = 1, 2, 3, 4 m s -1 and L h = 200, 400 m are displayed. While in the case of U g = 1 m s -1 (Fig 5A and 5E) there are large cellular circulations taking place, they disappear with increasing wind speed to give place to the formation of roll-type turbulent structures for U g > 3 m s -1 (Fig 5C, 5D, 5G and 5H). The structures resulting from U g = 2 m s -1 (Fig 5B and 5F) are neither cellular nor roll-like and are therefore excluded from the analysis. Fitting F 3 to the two datasets, we receive the following scaling functions:

Name U g (m s -1 ) z i (m) u� w� u�/w� -z i /L H I (%)
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.

Discussion
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 Kármá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  Table 1) and the respective scaling functions F 3,c or F 3,r (Eqs 13-14). The blue line shows the fitted scaling function F 2,N (Eq 15). The scaling function derived by De Roo et al. [58] is shown in grey for comparison.
https://doi.org/10.1371/journal.pone.0268097.g006 PLOS ONE 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%.

Conclusion
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.