How Does the Choice of the Lower Boundary Conditions in Large-Eddy Simulations Affect the Development of Dispersive Fluxes Near the Surface?

Large-eddy simulations (LES) are an important tool for investigating the longstanding energy-balance-closure problem, as they provide continuous, spatially-distributed information about turbulent flow at a high temporal resolution. Former LES studies reproduced an energy-balance gap similar to the observations in the field typically amounting to 10–30% for heights on the order of 100 m in convective boundary layers even above homogeneous surfaces. The underestimation is caused by dispersive fluxes associated with large-scale turbulent organized structures that are not captured by single-tower measurements. However, the gap typically vanishes near the surface, i.e. at typical eddy-covariance measurement heights below 20 m, contrary to the findings from field measurements. In this study, we aim to find a LES set-up that can represent the correct magnitude of the energy-balance gap close to the surface. Therefore, we use a nested two-way coupled LES, with a fine grid that allows us to resolve fluxes and atmospheric structures at typical eddy-covariance measurement heights of 20 m. Under different stability regimes we compare three different options for lower boundary conditions featuring grassland and forest surfaces, i.e. (1) prescribed surface fluxes, (2) a land-surface model, and (3) a land-surface model in combination with a resolved canopy. We show that the use of prescribed surface fluxes and a land-surface model yields similar dispersive heat fluxes that are very small near the vegetation top for both grassland and forest surfaces. However, with the resolved forest canopy, dispersive heat fluxes are clearly larger, which we explain by a clear impact of the resolved canopy on the relationship between variance and flux–variance similarity functions.


Introduction
The eddy-covariance (EC) method is well established and used worldwide in various networks for long-term measurements of energy and gas fluxes between ecosystems and the atmosphere (e.g., Baldocchi et al. 2001;Novick et al. 2018). However, the EC method often systematically underestimates these fluxes, which leads to a gap in the energy balance of about 10 to 30% (Hendricks-Franssen et al. 2010;Stoy et al. 2013;Soltani et al. 2018). This phenomenon has been widely discussed in the literature in the last few decades, and some possible causes, such as instrument and set-up errors in EC measurements (Laubach et al. 1994;Goulden et al. 1996;Kochendorfer et al. 2012;Nakai and Shimoyama 2012;Frank et al. 2013;, or measurement errors of other components of the energy balance, such as soil heat flux or net radiation (Liebethal et al. 2005;Kohsiek et al. 2007;Foken 2008), have already been excluded as a general problem across different sites (Foken 2008;Mauder et al. 2020).
An important reason for the underestimation of fluxes using single-tower measurements is the missed dispersive flux, i.e. the transport carried out by secondary circulations. These secondary circulations can be divided into two types, thermally-induced mesoscale circulations (TMC), which are generated by surface heterogeneity and are therefore spatially bound to the surface conditions (Foken 2008;Kenny et al. 2017;Bou-Zeid et al. 2020;Mauder et al. 2020), and slow-moving turbulent organized structures (TOS) that can develop even over homogeneous surfaces (Kanda et al. 2004;Inagaki et al. 2006). Both types of secondary circulations contribute to the vertical transport by a non-zero mean vertical velocity component, which cannot be captured by the EC method, since only the small-scale turbulent part of the flux (i.e. the temporal covariance-heat flux) is resolved. Field measurements have shown that secondary circulations reach well into the surface layer (Eder et al. 2015a) and contribute to the energy-balance gap (Eder et al. 2015b). In a large-eddy simulation (LES) study over a homogeneous forest where the canopy was explicitly resolved by the grid, Patton et al. (2016) showed that atmospheric structures that scale with the atmospheric boundary-layer height reach down to the canopy top and also occur in the understorey air space.
The use of long averaging intervals up to 24 h can reduce the energy-balance gap as they include the energy transport by TOS (Finnigan et al. 2003;Foken et al. 2006). With this approach, however, the temporal resolution is reduced. Thus, flux measurements on a half-hourly basis that are required for comparisons against numerical models are no longer available and diurnal variations are no longer detectable. Moreover, it is questionable whether stationarity can still be assumed, which is a prerequisite for applying Gaussian statistics, such as the calculation of a covariance .
Biomass heat storage is another major contributor to the energy-balance gap, especially in high vegetation such as forests (Lindroth et al. 2010;Leuning et al. 2012;Swenson et al. 2019). In this study, however, we do not address this factor.
In recent years, various experimental designs have been developed to systematically investigate the energy-balance-closure (EBC) problem in field measurements Foken et al. 2010;Butterworth et al. 2021). This requires a large number of spatiallydistributed EC measurements (Kanda et al. 2004;Xu et al. 2020), making the measurement campaigns accordingly complex and expensive.
Complementarily, LES are particularly well suited to investigating the influence of secondary circulations, as they provide continuous information with a high temporal and three-dimensional spatial resolution and can capture atmospheric motions on a wide range of scales (Inagaki et al. 2006;Schalkwijk et al. 2016). Furthermore, LES provide a controlled environment in which the boundary conditions are known, as opposed to realworld conditions, and virtual measurements without instrument errors (Inagaki et al. 2006;Schalkwijk et al. 2016;Sühring et al. 2018).
In recent years, computational resources have increased significantly, enabling LES with high spatial resolution. This, in turn, facilitates the investigation of turbulence structures and flow components near the surface, including the canopy layer (e.g., Kanani-Sühring and Raasch 2017;Kröniger et al. 2018). It is particularly important to investigate the EBC problem near the surface since EC measurements are typically taken in the surface layer at a height of 2-40 m, depending on the vegetation height (Hendricks-Franssen et al. 2010;Butterworth et al. 2021).
Several studies have already investigated the EBC problem using LES. Kanda et al. (2004), Steinfeld et al. (2007), and Huang et al. (2008) analyzed the height dependency of the imbalance and found that virtual tower measurements at higher altitudes underestimated sensible heat fluxes even if the surface itself was homogeneous and no TMCs were generated by heterogeneous surface heating (Inagaki et al. 2006). However, using a rather large vertical grid spacing of 25 m (Kanda et al. 2004) and 20 m (Huang et al. 2008), it was found that the energy-balance gap vanished near the surface. Moreover, using simulations of similar grid resolution, Steinfeld et al. (2007) showed that the imbalance close to the ground reduced to less than 5%, which does not match field observations. Hence, they concluded that TOS could not explain the magnitude of the energy-balance gap. Zhou et al. (2019) examined the magnitude of the energy-balance gap in relation to landscape heterogeneity and found it to be largest for heterogeneities on the scale of the boundary-layer height. However, Inagaki et al. (2006) studied the influence of surface heterogeneity on the imbalance and found the imbalance above a homogeneous surface to be on the same order of magnitude as above a heterogeneous surface. This was also supported by Margairaz et al. (2020), who showed that persistent structures in half-hourly averaged vertical wind velocity were as pronounced above homogeneous surfaces as over heterogeneous surfaces.
We hypothesize that the common use of prescribed surface fluxes (PSF) is one important reason why former LES studies underestimate the energy balance. With prescribed fluxes, the surface-atmosphere exchange is decoupled from the atmosphere (i.e. the atmosphere responds to the surface) while the surface does not respond to the atmosphere. With this mutual feedback missing, possibly important aspects such as self-reinforcement or weakening of secondary circulations are not accounted for in the simulations.
In the present study, we aim to discover the reasons why LES studies have so far been unable to reproduce near-surface EBC gaps that were similar to observations and if different lower boundary conditions have an influence on dispersive heat fluxes.
In order to answer our research question, we set up two branches of numerical experiments, one for short grass and one for tall vegetation. We performed LES with varying surface-boundary conditions (as illustrated in Fig. 1) and evaluated the effect of the surface-boundary condition with regard to the energy-balance closure for varying atmospheric stability.
For the short-grass experiments, we set up LES with an interactive land-surface model (LSM, Gehrke et al. 2020) where the surface-atmosphere exchange concerning water and heat exchange is modelled explicitly, and compared these against simulations with prescribed surface fluxes where the surface does not interact with the atmosphere at all. Maronga et al. (2020) describe in detail how near-surface air temperature feeds into the LSM.
For the tall-vegetation experiments, we set up simulations with an interactive LSM combined with a plant-canopy model (PCM) where trees are explicitly resolved by the numerical grid. Again, these were compared against simulations with prescribed surface fluxes with respect to the energy-balance closure. We aim to assess the importance of landsurface representation in LES for different kinds of land use.
Since we use an idealized, homogeneous set-up, we investigate the contribution of vertical transport by TOS only, hereafter referred to as the dispersive heat flux. We wish to test if the use of LSM and PCM lead to larger dispersive fluxes by adapting to changes in the atmosphere.
The next section describes our procedure divided into set-up of the different simulations (Sect. 2.1) and evaluation of the simulations, i.e. the calculation of the heat flows (Sect. 2.2). Section 3 presents the results, first considering the comparability of the respective PSF and LSM(+ PCM) studies (Sect. 3.1), before we examine the formation of TOS at different lower boundary conditions (Sect. 3.2), the resulting dispersive fluxes (Sect. 3.3) and finally the influences of the LSM in different set-ups (Sect. 3.4). The results are then analyzed and discussed against the background of other findings from the field in Sect. 4 and our findings are summarized in Sect. 5.

Methods
In this study, we compared different lower boundary conditions using PALM v6 ). We used a highly idealized set-up with homogeneous surfaces and cyclic lateral boundary conditions, while the stability classes were similar to previous LES studies, such as Patton et al. (2016) andDe Roo et al. (2018). To obtain a representative result, the simulations were carried out for various combinations of two different land-cover types, grassland (G) and forest (F), as well as for moderately unstable (MU), strongly unstable (SU), and free convective (FC) conditions. An overview of the combinations of lower boundary conditions and atmospheric stabilities is shown in Table 1.

Large-Eddy Simulation Set-up
We used the LES Model PALM, version 6, revision 4529 , for the numerical simulations. PALM solves the non-hydrostatic incompressible Boussinesq equations. For the subgrid model, the kinetic energy scheme of Deardorff (1980) Moeng and Wyngaard (1988) and Saiki et al. (2000) was used. The advection terms were discretized using a fifth-order scheme (Wicker and Skamarock 2002), and a third-order Runge-Kutta scheme by Williamson (1980) was used for the time integration.

modified by
To achieve a high grid resolution near the surface, we employed the vertical grid nesting technique available with PALM (Hellsten et al. 2021) where a child domain with smaller grid spacing but the same horizontal extent is placed within the parent domain with larger grid spacing, while these two domains interact with each other. The horizontal domain measured 7200 × 7200 m 2 , with the parent domain reaching up to 2400 m and the child domain up to 240 m. In the coarse grid, the horizontal and vertical grid spacings were 30 m and 20 m, respectively, resulting in (x, y, z) = 240 × 240 × 120 grid points, whereas, in the fine grid, the horizontal and vertical grid spacings were 6 m and 4 m, respectively, which yields (x, y, z) = 1200 × 1200 × 60 grid points. The parent and child domains both reached down to the lower boundary of the domain and simulations in the both domains were run parallelly with two-way nesting. The timestep was set to a constant value of 0.5 s. Each simulation consisted of 2 h of model spin-up time followed by a 4-h period during which data were captured. The latitude was set to 46 degrees north.
The initialization of the atmospheric conditions follows De . The initial horizontal velocity profile was vertically constant and homogeneous over the entire horizontal extent of the domain, with velocity in the x-direction and geostrophic wind speed varying among the simulations (see Table 2). The geostrophic velocity here featured a horizontal pressure gradient that is oriented in x-direction. Different geostrophic wind speeds are used for the three atmospheric stabilities as shown in Table 2.
The initial potential temperature at the surface was set to 295 K. Between 40 and 800 m, a vertical gradient of 3 × 10 -3 K m −1 was added and above, the gradient was 8 × 10 -3 K m −1 . The mixing ratio at the surface was set to 8 × 10 -3 kg kg −1 . Between 1000 and 1100 m, a vertical gradient of − 1 × 10 -5 m −1 was imposed, while below and above this area, no vertical gradient was applied. All profiles used for initialization (shown in Appendix, Fig. 11) were homogeneous over the entire horizontal extent of the domain. A stable inversion layer at the top of the domain ensured that the processes within the boundary layer were not affected by the vertical extent of the domain. The horizontal extent of the domain was at least seven times the boundary-layer depth for all atmospheric conditions. Randomly distributed perturbations were imposed on the horizontal velocity fields at the beginning of each simulation to initiate turbulence. At the lateral boundaries, cyclic conditions were applied. At the surface, we set an impermeable boundary with zero vertical velocity and imposed surface stress by applying Monin-Obukhov similarity theory (MOST) locally between the surface and the first vertical grid level. The resulting surface fluxes of horizontal momentum are then entered as lower boundary conditions via the subgrid-scale term. More precisely, surface horizontal momentum fluxes w ′ u ′ i were computed from (Maronga et al. 2015), where u∕ z is the vertical gradient of horizontal wind speed between the first prognostic grid level and the surface; u * is the friction velocity; = 0.4 is the von Kármán constant, and (z∕L) is the similarity function for momentum in the formulation of Businger-Dyer (see e.g. Panofsky and Dutton (1984)), with L being the Obukhov length. For the potential-temperature and mixing-ratio equations we also employed a fluxboundary condition at the surface, with fluxes either prescribed or computed by the land-surface model . For the boundary values of potential temperature and mixing ratio itself, we employed a zero-gradient Neumann boundary condition at the surface. This is to ensure that no flux contribution from the resolved-scale advection arises at the surface, leading to any double counting of the vertical transport.
At the domain top, we set zero-gradient Neumann conditions for the horizontal velocity components, reflecting the geostrophic wind in the upper part of the model domain. The vertical wind velocity at the domain top was set to zero to maintain continuity. Table 2 Settings for different combinations of land-cover type and atmospheric stabilities in LSM (simulation numbers 1-3), LSM + PCM (simulation numbers 4-6) and PSF (simulation numbers 7-12) simulations. The roughness lengths for momentum ( z 0 ) and heat ( z 0h ) are based on the vegetation types that are predefined in the PALM LSM , based on ECMWF-IFS classification Simulation number Zero-gradient Neumann conditions at the top boundary were also applied for potential temperature and mixing ratio. Moreover, as in De Roo et al. (2018), a vertical subsidencevelocity gradient of 4 × 10 -5 m s −1 m −1 up to 800 m and 2 × 10 -5 m s −1 m −1 between 800 and 1000 m was prescribed. By this, the boundary-layer depth during the analysis period was kept constant. At the surface, the subsidence velocity was zero.
In addition to varying geostrophic wind speeds, the three different atmospheric stabilities were defined by different prescribed surface fluxes, or incoming radiation where the LSM was used, as shown in Table 2. To ensure highest comparability between PSF and LSM(+ PCM) simulations, the resulting surface-heat fluxes for each atmospheric stability and vegetation type combination have to be similar in PSF and LSM(+ PCM) simulations, which is why we first ran the simulations with LSM(+ PCM) as a boundary condition and then used the resulting surface-heat fluxes in the PSF simulations as described below.

Set-up of the Land-Surface Model Simulations
In the LSM simulations (simulation numbers 1-3) over grassland, the vegetation type (VT) short grass (VT = 3) as specified in PALM Maronga et al. 2020) was used. The VT parameter provides predefined values of various vegetation parameters (e.g., leaf-area density, heat-capacity, canopy-specific resistance, aerodynamic roughness length z 0 , etc.) that determine the influence of vegetation on atmospheric processes that are shown in Table 2. For details of the default bulk parameters, we refer to Gehrke et al. (2020).
To ensure the highest comparability of simulations with the LSM to simulations with prescribed surface fluxes, a time-constant net radiation at the surface is necessary. Therefore, we used the PALM built-in clear-sky radiation model with a constant zenith angle and thus obtained a net radiation that remained almost constant over the 4 h of data acquisition. The zenith angles and resulting net radiation at the surface used for each stability and landcover combination can be found in Table 2. The zenith angles were chosen to give different net radiation at the surface for each atmospheric condition. For net radiation, we used 250 W m −2 (MU), 350 W m −2 (SU), and 450 W m −2 (FC) following De .
The LSM is coupled with a soil model for which we chose a medium-fine soil type . Since Liu and Shao (2013) have noted that a very thin top soil layer can lead to feedback effects as soil temperature and moisture can change due to short-term changes in the atmosphere directly above the soil, we have also adjusted the layer thickness in the soil model. The thickness of each soil layer, as well as initial soil temperature and moisture values, are shown in the Appendix, Table 6. The wilting point of the soil is defined at a soil moisture of 0.133 m 3 m −3 ) and the soil moisture in our simulations never fell below this value, therefore water availability for the plants was sufficient in all simulations.

Set-up of Land-Surface Model and Plant-Canopy Model Simulations
For the land-cover type forest, the land-surface model was additionally combined with PALM's embedded plant-canopy model (PCM) (Maronga et al. 2015) which explicitly considers the impact of grid-resolved vegetation on the momentum, potential temperature and heat equation (see simulations 4-6). The PCM was used with the horizontally homogeneous prescribed leaf-area density ( LAD ) profile shown in Fig. 2. The PCM follows Shaw and Schumann (1992) and Watanabe (2004), and adds a momentum sink. It was validated against wind-tunnel and lidar observations in Kanani et al. (2014). The PCM furthermore interacts with radiation (absorption, transmission, reflections) and provides volume sources of sensible and latent heat that enter the prognostic equations of potential temperature and mixing ratio, respectively (Krč et al. 2021). The PCM is currently not coupled with the LSM but assumes that water availability is always sufficient for transpiration. The shape of the LAD profile used in the fine grid (Fig. 2) is based on the plant-area density ( PAD ) profile used in Patton et al. (2016) but has a higher leaf-area index ( LAI ) of 5.95 m 2 m −2 .
Below the canopy, at the bottom boundary, we also used the LSM that provides a small part of the sensible and latent heat fluxes based on the radiation that penetrates the resolved canopy. Here, we used the vegetation type deciduous broadleaf forest (VT = 7), the characteristic parameters of which are also shown in Table 2. However, the roughness length is determined by the vegetation resolved in the PCM.

Set-up of Prescribed Surface-Flux Simulations
To set up PSF simulations comparable to the LSM simulations (simulation numbers 7-9), the sensible and latent surface-heat fluxes ( H s and E s ) resulting from the grassland LSM simulations were horizontally averaged over the entire domain and temporally averaged canopy top were averaged and usedover the four hours of data acquisition and then used as prescribed surface fluxes in the PSF simulations. For PSF simulations comparable to LSM + PCM simulations (simulation numbers 10-12), the fluxes at the as prescribed surface-heat fluxes (see Eqs. 2, 3). The exact approach for calculating the surface fluxes is described in Sect. 2.2. The prescribed H s and E s in each PSF simulation are shown in Table 2. The roughness lengths z 0 and z 0h were set according to the vegetation types used in the LSM(+ PCM) simulations.

Data Processing
30-min surface fluxes were calculated for each simulation and profiles of horizontal wind speed in the x -and y-directions ( u , v ), the component of vertical wind ( w ), and potential temperature ( ) were compared to estimate if the LSM(+ PCM) simulations are comparable to the respective PSF simulations, respectively. In this study, "surface" always refers to the canopy top unless stated otherwise, which means that the surface is at z b = 0 in PSF and LSM simulations, and at z b = 20 m in LSM + PCM simulations, where z b is the height above the bottom of the domain. Therefore, for the LSM + PCM simulations, the latent and sensible heat fluxes at z b = 20 m were used as surface fluxes.
The contributions of the 30-min individual flux components to the total energy fluxes originating from the vegetation surface H s and E s were calculated for all simulations. The surface-heat fluxes used in the PSF and LSM-only simulations were calculated directly by PALM as horizontal domain averages. For the simulations with PCM, we used the sum of the resolved heat fluxes, i.e. the temporal covariances w ′ ′ and w ′ q ′ , and the SGS fluxes w ′ ′ SGS and w ′ q ′ SGS provided by the SGS model at the canopy top, i.e. at z b = 20 m The overbar denotes temporal averaging over 30 min at each grid point. Furthermore, the covariances were spatially averaged in the x-and y-directions which is denoted by the angled brackets. To convert the heat fluxes from kinematic units to dynamic units (W m −2 ), they were multiplied with the latent heat of vaporization of air ( c p ) and density (ρ) and the latent heat of vaporization ( ), respectively.
In the surface layer, the total heat flux, i.e. the heat flux originating from the surface, is divided into individual flux components where H t and E t are the temporal covariance-heat fluxes, H d and E d are the dispersive heat fluxes, and S is the energy stored in the underlying air mass.
We calculated H t , E t , H d , and E d for each vertical grid level. H t and E t were calculated similarly to the surface fluxes in the PCM simulations (see Eqs. 2, 3) The resolved 30-min temporal covariances w ′ ′ and w ′ q ′ in Eqs. 2, 3 and 5, 6 were calculated from the temporal deviations from half-hourly averages of w , , and q , respectively: where nt is the number of timesteps in each 30-min averaging interval. To calculate w ′ q ′ , was replaced by q in Eq. 7. Those 30-min covariances where then spatially averaged over each x-and y-grid point with where nx and ny are the number of grid points in x and y directions. Again, was replaced by q to calculate ⟨w ′ q ′ ⟩ in Eq. 8. Similarly, the SGS fluxes calculated by PALM were spatially averaged over the entire domain, which is denoted by the angular brackets.
The dispersive fluxes were calculated at each grid level from the spatial covariance, where the local deviations of 30-min averages of vertical wind velocity w , potential temperature , and specific humidity q , marked by the star, from the spatial average were considered with Again, ⟨w * * ⟩ was calculated by replacing by q in Eq. 11. The resulting surface-heat fluxes, temporal covariance-heat fluxes and dispersive heat fluxes were temporally averaged over 30 min and horizontally averaged over the entire domain. Those values were then averaged over the entire data-capture period, i.e. 4 h.
We investigated the dispersive heat fluxes ( F d ), as these are typically considered as the energy-balance residual (e.g., Steinfeld et al. 2007). Because the absolute surface-heat fluxes of the compared simulations differ slightly, and because some energybalance terms like ground heat flux and energy stored within the canopy cannot be considered in the PSF simulations, the relative contributions of H d and E d to the total surface heat flux F s = H s + E s at the canopy top were investigated. The resulting F d values of the LSM(+ PCM) and corresponding PSF simulations were then compared at each grid level up to 100 m above the vegetation top.
We also analyzed xy cross-sections of 30-min averaged surface temperature ( s ) and vertical wind velocity ( w ) at different heights above the surface to see if the differences in F d are associated with the formation of mesoscale structures in the atmosphere and if these spatial structures are linked to surface properties. Furthermore, the Bowen ratio was evaluated to see whether the proportion of sensible and latent heat fluxes in the dispersive heat fluxes not captured by single-tower EC measurements in the field is the same as in the temporal covariance-heat fluxes.

Results
The focus in the data evaluation is on F d resulting from different lower boundary conditions since the main objective is how different surface-flux boundary conditions affect the surface energy balance and the occurrence of dispersive fluxes. However, before we analyze the impact of lower boundary conditions on F d , we evaluate the comparability of the LSM(+ PCM) simulations with the PSF simulations by comparing general boundary layer characteristics. Furthermore, xy cross-sections are presented to demonstrate that secondary circulations develop depending on the height and the lower boundary conditions.

Comparability of Simulations with Different Lower Boundary Conditions
In the simulations with prescribed surface-heat fluxes, the surface-heat fluxes are constant over the entire simulation period. To enhance comparability, the surface-heat fluxes in LSM(+ PCM) simulations should also be fairly constant in time. We therefore used constant sun angles, which lead to an almost constant net radiation. Figure 3 shows that the surface-heat fluxes are quasi-steady-state in the LSM simulations over grassland. Only when additionally using the PCM over forest does the temporal covariance-heat fluxes at the vegetation top fluctuate slightly, being almost constant during the first 2.5 h of data assimilation and then decreasing slightly. Overall, the surface-heat fluxes of LSM(+ PCM) simulations are constant within ± 12.2% and only deviate from the surface fluxes of the corresponding PSF simulations by 6.2% at maximum. Figure 4 shows profiles of horizontally averaged turbulence kinetic energy ( ⟨TKE⟩ ), ⟨ ⟩ and ⟨q⟩ . Here, the results of an LSM(+ PCM) and the respective PSF simulation are compared in each sub-panel. The profiles of ⟨TKE⟩ , ⟨ ⟩ and ⟨q⟩ are almost identical in the LSM(+ PCM) and respective PSF simulations.

Development of Turbulent Organized Structures for Different Lower Boundary Conditions
Figures 5, 6 and 7 show xy cross-sections averaged over the fifth 30-min interval of the data assimilation period for the three different atmospheric stabilities over forest with PSF (a-c) and LSM + PCM (d-f) as lower boundary conditions for three different heights above the vegetation top denoted by z v . Panels a, d in Figs. 5, 6 and 7 show the deviation from the spatially averaged temperature at the vegetation top ⟨ s ⟩ . In the MU case (Fig. 5), a striped pattern occurs at z v = 0 both with PSF ( Fig. 5a) and with LSM + PCM (Fig. 5d), but it is much more pronounced with LSM + PCM, which is confirmed by the higher spatial standard deviation in Table 3. In the SU case (Fig. 6), cold and warm areas also appear in the surface temperature, but here, the stripes are more fractured. Again, the cold and  Table 1) temporally averaged over the entire data assimilation period. The results of the PSF simulations are displayed as a black line, while the results of the corresponding LSM(+ PCM) simulations are displayed as red lines. Here, we show the output of the parent domain to cover a larger vertical extent warm patches are more pronounced with LSM + PCM. In the FC case (Fig. 7), the warm and cold areas appear in a cell-like structure and are, unlike in the first two cases, considerably more pronounced with PSF at the surface (Fig. 7a) than with LSM + PCM at z v = 0 (Fig. 7d). This is also reflected in the standard deviations of 0.119 K (PSF) and 0.096 K (LSM + PCM) shown in Table 3.  The structures in the vertical velocity (Figs. 5, 6, 7b-c, e-f) under different atmospheric conditions match the patterns in the surface temperature, respectively. Under all atmospheric conditions, updrafts and downdrafts in the vertical wind at z v = 4 m are significantly weaker in PSF boundary conditions (Figs. 5,6,7b) than in LSM + PCM boundary conditions (Figs. 5,6,7e) which is also reflected in the spatial standard deviations shown in Table 3. With greater distance from the surface, the patterns are still less distinct in the PSF simulations (Figs. 5, 6, 7c), but the difference to the LSM + PCM simulations (Figs. 5,6,7f) is reduced as shown, for example, at 20 m above the surface.
Panels g-j of Fig. 5, 6 and 7 additionally display profiles of ⟨w ′ w ′ ⟩ and ⟨ ′ ′ ⟩ for the same 30-min interval. The profiles show that ⟨w ′ w ′ ⟩ at the vegetation top is under all atmospheric conditions significantly larger in the LSM + PCM simulations (Figs. 5,6,7g) than in the PSF simulations (Figs. 5,6,7i) due to the resolved vegetation. However, the difference between PSF and LSM + PCM simulations becomes smaller with increasing height.
Due to the resolved vegetation, small differences between PSF (Figs. 5, 6, 7h) and LSM + PCM (Figs. 5, 6, 7j) simulations also arise for ′ ′ , especially directly above the vegetation top. These differences are, however, pronounced to varying degrees under different atmospheric conditions. While in the MU case, ′ ′ is larger for LSM + PCM (Fig. 5h) than for PSF (Fig. 5j), the profiles look almost the same above the canopy top with both lower boundary conditions in the FC case (Fig. 7). Table 4 shows that the Bowen ratios of temporal covariance ( Bo t ) and dispersive ( Bo d ) heat fluxes at 20 m above are almost equal and that they match the Bowen ratio of surfaceheat fluxes ( Bo s ). It furthermore shows that Bo values for LSM(+ PCM) and comparable PSF simulations agree quite well. Figure 8 shows comparisons of latent and sensible, as well as total dispersive fluxes ( F d ) that are normalized with F s between the simulations with land surface model ( F d,LSM ) ( Fig. 8a-c) or land-surface and plant-canopy model ( F d,LSM+PCM ) ( Fig. 8d-f) and the respective simulations with prescribed surface fluxes ( F d,PSF ) for each stability and VT combination.

Dispersive Fluxes for Different Lower Boundary Conditions
For both VTs, F d increases with atmospheric instability and with distance from the vegetation top. The Bowen ratios of F d ( Bo d ) nearly equal the Bowen ratio of F t ( Bo t ) in all simulations (see Table 4); F d is dominated by E d in all simulations (Fig. 8). In the  grassland case, F d,PSF is slightly smaller than F d,LSM under all atmospheric conditions as indicated by Fig. 8a, b. Figure 8d-f show the comparison between F d,PSF and F d,LSM+PCM over forest and demonstrate that the dispersive fluxes above the resolved vegetation behave differently than with PSF. In contrast to F d,PSF , F d,LSM+PCM does not approach zero close to the canopy top but remains significantly larger. F d,PSF accounts for 0.52 ± 0.01% (1.44 ± 0.0 W m −2 , MU case), 0.79 ± 0.02% (3.05 ± 0.0 W m −2 , SU case), and 1.46 ± 0.01% (7.27 ± 0.0 W m −2 , FC case) of F s,PSF 4 m above the canopy top, whereas F d,LSM+PCM contributes 1.14 ± 0.09% (3.23 ± 0.01 W m −2 , MU case), 2.76 ± 0.21% (11.23 ± 0.04 W m −2 , SU case), and 5.78 ± 0.45% (30.7 ± 0.17 W m −2 , FC case) to F s,LSM+PCM , respectively. In the FC case, this leads to an absolute difference of 23.42 W m −2 , though it must also be noted here that F d slightly differs between PSF and LSM + PCM simulations. Therefore F d,LSM+PCM is significantly larger near the surface than F d,PSF under all atmospheric conditions. However, with increasing distance from the surface, F d,LSM+PCM behaves very differently under changing atmospheric conditions. While F d,LSM+PCM under MU conditions is smaller than F d,PSF already 8 m above the vegetation top, it is never smaller than F d,PSF under SU conditions. Under FC conditions, F d,LSM+PCM becomes smaller than F d,PSF 20 m above the vegetation top. Figure 9 shows the differences between F d,PSF and F d,LSM or F d,LSM+PCM in more detail. It demonstrates that the differences between F d,PSF and F d,LSM increase linearly with height over grassland under all atmospheric conditions. The differences between F d,PSF and F d,LSM+PCM over forest are not linear due to the different behaviour of F d,LSM+PCM .
In the MU case, this leads to F d,LSM+PCM being smaller than F d,PSF already at 8 m above the canopy. However, since F d,LSM+PCM increases more strongly again from about 20 m above the canopy, F d,LSM+PCM and F d,PSF are equal again, at about 80 m above the vegetation top (Fig. 9a). In the SU case, the difference between F d,LSM+PCM and F d,PSF initially becomes smaller with increasing height, but always remains slightly larger and the difference increases, again, from about 30 m above the surface as shown in Fig. 9b. In the FC case, the difference between F d,LSM+PCM and F d,PSF close to the vegetation top is largest ( Fig. 9c) with 4.32 ± 0.47% (23.43 ± 0.17 W m −2 ). Here, again, F d,LSM+PCM approaches F d,PSF with increasing height but here, it becomes smaller than F d,PSF from 24 m above the canopy. Table 5 shows the contributions of different water vapour sources to the latent surface-heat flux directly at the bottom of the domain ( E s,b ) where the LSM is embedded. In the grassland simulations, the entire E s,b originates from the transpiration of plants ( E s,veg ). Here, no water evaporates directly from the soil ( E s,soil ) as the plant coverage is set to 100%. For VT 7 (broadleaf forest) that is used for the LSM + PCM simulations, the share of E s,soil in E s,b amounts to less than 10%. Furthermore, the contributions from liquid water on plants ( E s,liq ) are negative here, which is defined as the condensation of water vapour. Overall, E s,b in the simulations with PCM are very small compared to the total latent heat flux at the canopy top, meaning that the main portion of the latent heat flux originates from the crown space of the resolved plants.

Comparison of Vegetation Types in Land-Surface Model and Plant-Canopy Model Simulations and Surface-Heat Fluxes from the Land-Surface Model in General
In comparison to E s,b at the bottom of the domain shown in Table 5, the latent heat fluxes at the canopy top ( E s,v ) are significantly larger with 243.13 ± 4.46 W m −2 (MU), 318.50 ± 12.81 W m −2 (SU), and 398.74 ± 26.56 W m −2 (FC). In the LSM + PCM simulations, the sensible surface-heat flux at the bottom of the domain ( H s,b ) shown in Table 5 are negative, indicating that the canopy heats up more than the shaded ground, causing warm air masses to sink towards the ground. At the canopy top, however, the sensible heat fluxes ( H s,v ) are similar to the prescribed sensible heat fluxes shown in Table 2 with 39.63 ± 1.89 W m −2 (MU), 88.20 ± 5.39 W m −2 (SU), and 132.37 ± 11.56 W m −2 (FC).

Discussion
The investigation of F d resulting from different lower boundary conditions showed that the differences between LSM and PSF close to the vegetation top are very small for simulations over grassland. The PSF lower boundary condition results in slightly larger F d with increasing distance from the canopy top than when using the LSM.
The comparison between the LSM + PCM and the PSF over the forest shows varying results depending on which height is considered. Close to the surface, the LSM + PCM boundary condition yields significantly larger F d values, and this difference increases with increasing instability. This agress well with the developing TOS, i.e. roll-like structures in MU and SU cases, and cell-like structures in the FC case, being much more pronounced near the vegetation top in LSM + PCM simulations than in PSF and LSM-only simulations. This indicates that the magnitude of dispersive fluxes is correlated with the strength of TOS when averaged over 30 min.
One reason for the undeveloped structures close to the surface despite high heterogeneity in surface temperature in PSF simulations might be the fact that structures at the lowest grid levels are only partially resolved in LES and significant parts of the vertical transport are carried out by the SGS model (Bou-Zeid et al. 2005). When using the PCM, the canopy top is at the fifth grid level, so structures directly above the trees can be resolved well. Other studies also found that resolved canopies enhance the turbulence, e.g. by causing ejection sweep events (Finnigan et al. 2009), which explains why the dispersive fluxes behave so differently with the PCM than with PSF or the LSM only. However, in our study, resolved vegetation affects dispersive fluxes only directly above the vegetation top and no consistent positive effect on dispersive fluxes is observed as distance increases. Patton et al. (2016) observed that the structures form already within the trunk layer, vanish at the treetop, and are more pronounced again above the forest. We do not observe this effect in our simulations. Instead, the structures in our simulations become constantly more pronounced with increasing height. This discrepancy might be explained by the different canopy structure of Patton et al. (2016), where the crown layer has a lower plant area density.
To investigate how the height-dependent differences in dispersive heat fluxes over the forest in Fig. 9 can be explained and whether the different behaviour of the variance profiles in Figs. 5, 6 and 7, especially for w , provided an explanation, we investigated flux-variance similarity functions for u , w , and variances. The similarity function for u ( u ) was calculated based on Panofsky et al. (1977) The variance similarity functions for w ( w ) and ( ) were calculated similarly by replacing u ′ u ′ with w ′ w ′ and ′ ′ , respectively, and u * with the temperature scale T * . The differences in between PSF and LSM + PCM simulations were very small and are hence not shown. However, we found significant differences between PSF and LSM + PCM simulations in the u and w profiles that are shown in Fig. 10. We furthermore compared the profiles obtained from the simulations to theoretical variance similarity functions based on (Foken et al. 2004) using 1 3 with and with Figure 10 shows a clear dependence of the u and w profiles on stability, as does the proportion of dispersive heat fluxes in the total surface-heat fluxes in Fig. 8. For all atmospheric stabilities, u and w are larger in the LSM + PCM simulations than in the PSF simulations directly above the canopy top. This is because the profiles in the PSF simulations are forced to zero at the canopy top. Within the canopy, u and w are highly variable and depend on the vegetation and canopy structure at each site (Rannik et al. 2003) and we do not expect the model to provide realistic values in this range. However, a clear effect of the resolved vegetation is that u and w are not forced to zero at the vegetation top, which explains why the resulting dispersive heat fluxes near the canopy top are larger in the LSM + PCM simulations than in the PSF simulations (see Fig. 9). The fact that w becomes larger in the PSF simulations than in the PCM + LSM simulations from a height of about 20 m above the canopy top also explains the negative difference further away from the canopy top in Fig. 9. The bulbous shape of u , and w increasing with altitude suggest that secondary circulations are forming that provide greater variance in the vertical wind at higher altitudes and greater variance in the horizontal wind at lower altitudes. Both effects are very pronounced in the FC case, which causes the profiles derived from the simulations to deviate strongly from the profiles calculated using the theoretical functions, which only include the purely turbulent fluxes. This indicates more intensive secondary circulations for the FC case, which is why the dispersive heat fluxes are also larger here than for the less unstable cases.
Whether the use of PSF or LSM(+ PCM) as lower boundary condition yields larger dispersive heat fluxes, depends on the land cover type, the measurement height and atmospheric stability. For more realistic simulations over heterogeneous surfaces of energy-balance closure field experiments, it should therefore be considered which type of vegetation is predominant. Moreover, the PCM can be used only for those areas where the vegetation is tall enough to be resolved by the grid. In most cases, depending on the grid spacing, only a forest can be resolved using the PCM and the grassland would be simulated by using the LSM or the PSF method. Assuming a typical EC set-up is roughly mounted 4 m above the canopy for grasslands and 12 m above the canopy for forests, the use of the LSM instead of PSF as the lower boundary condition for grasslands would decrease the contribution of F d to F s by only − 0.08 ± 0.06% (− 0.16 ± 0.0 W m −2 ) for the MU case, − 0.09 ± 0.03 (− 0.24 ± 0.0 W m −2 ) for the SU case, and − 0.41 ± 0.08% (− 1.05 ± 0.0 W m −2 ) for the FC case. Over the forest, the use of LSM + PCM instead of PSF would decrease the dispersive heat fluxes by − 0.71 ± 0.19% (− 1.86 ± 0.01 W m −2 ) in the MU case, but increase them by 0.32 ± 0.14% (1.89 ± 0.04 W m −2 ) and 1.0 ± 0.22% (7.67 ± 0.2 W m −2 ) for SU and FC cases, respectively. Additional factors should also be considered when choosing the lower boundary condition, such as the fact that the flow modification by forest edges cannot be represented by PSF or LSM. These effects are not investigated in this study, although they can have a strong influence on temporal covariance and dispersive heat fluxes (Kanani-Sühring and Raasch 2017;Kenny et al. 2017).
Furthermore, forests typically have a roughness length of > 1 m (Foken 2017). However, z 0 is limited to a quarter of the vertical grid spacing in PALM because MOST is not valid within the roughness sublayer (Basu and Lacser 2017). Therefore, forests cannot be represented correctly by the LSM with vertical grid spacings < 4 m. To estimate the roughness length resulting from the resolved canopy in the LSM + PCM simulations over forest, we have performed an additional simulation with near neutral atmospheric stratification. We fitted where is the von Kármán constant (0.4) and d is the zero-plane displacement to the horizontal wind profile and found z 0 to be 2.2 m. Thus, z 0 is almost a factor of 10 larger than in the PSF simulations and much closer to z 0 determined from various field measurements over forests (Reithmaier et al. 2006).
Finally, when using the PCM, we found fluxes from the ground to be very small. Most of the incoming shortwave radiation is converted into longwave radiation and heat fluxes at the top of the tree canopy. Due to the shading by the trees, almost no energy reaches the ground where it is fed into the LSM and coupled soil model. This also means that the choice of vegetation type in LSM in combination with the PCM plays hardly any role.
As mentioned earlier, biomass heat storage was found to play a major role in energybalance closure, especially in forests (Lindroth et al. 2010;Swenson et al. 2019). However, the simulations with PSF or LSM only do not provide information on biomass storage or processes within the canopy and in the understory air space. To answer our research question about the influence of lower boundary conditions on dispersive heat fluxes, we have therefore limited ourselves to the area above the vegetation surface, i.e. dispersive flux components relative to the total heat flux at the canopy top. Therefore, we cannot draw any conclusions on the influence of biomass heat storage.

Conclusion
To answer whether the choice of lower boundary conditions has an influence on the resulting dispersive heat fluxes in LES, we compared two different boundary conditions over grass and forest, respectively, using different atmospheric conditions. We found that especially the use of the PCM, i.e. resolving the vegetation in the model, indeed has a significant influence on the dispersive fluxes. The comparison between PSF and LSM for grassland shows no significant differences near the vegetation top for all atmospheric conditions with respect to the proportion of dispersive heat fluxes to the total surface-heat flux. The use of LSM + PCM for forests causes significantly larger dispersive fluxes than PSF near the surface regardless of the stability regime. The more unstable the atmospheric conditions are, the greater this difference. At larger measurement heights, however, the effect is reversed depending on the atmospheric stability. We showed that resolving the vegetation in the PCM also has a clear impact on the relationship between variance and flux-variance similarity functions, which is related to the relative importance of dispersive heat fluxes. The use of the PCM is therefore recommended for low measurement heights within 20 m above the vegetation top, which is the case for many field measurement set-ups. Nevertheless, these simulated dispersive fluxes are much smaller than the EBC gap of realworld EC measurements. Despite this discrepancy, the PCM also has other advantages, for example, it allows for surface types with a high roughness, such as forests, to be correctly represented even at high grid resolutions on the order of metres. In the case of heterogeneous surfaces, resolved forest edges can also have an additional effect on dispersive fluxes which was, however, not investigated in this study, where only homogeneous conditions were analyzed. Figure 11 shows horizontally-averaged profiles for , q , and w that were used to initiate all simulations. Within the 2-h spin-up time, the boundary layer develops so that the profiles evolve into those shown in Fig. 4.  Table 6 shows the initial set-up of the soil model that is part of the LSM. The soil temperature and moisture were applied in each layer throughout the entire horizontal extent of the domain. This set-up was used in all simulations where LSM(+ PCM) was employed as lower boundary condition.