Modeling of land-surface interactions in the PALM model system 6.0: Land surface model description, first evaluation, and sensitivity to model parameters

In this paper the land-surface model embedded in the PALM model system is described and evaluated against in-situ measurement data in Cabauw. For this, two consecutive clear-sky days are simulated and the components of surface energy balance, as well as near-surface potential temperature, humidity and horizontal wind speed are compared against observation data. For the simulated period, components of the energy balance agree well during dayand nighttime, and also the daytime Bowen ratio agrees fairly well compared to the observations. Although the model simulates a significantly more stably-stratified 5 nocturnal boundary layer compared to the observation, near-surface potential temperature and humidity agree fairly well during day. Moreover, we performed a sensitivity study in order to investigate how much the model results depend on land-surface and soil specifications, as well as atmospheric initial conditions. By this, we find that a false estimation of the leaf area index, the albedo, or the initial humidity causes a serious misrepresentation of the daytime turbulent sensible and latent heat fluxes. During night, the boundary-layer characteristics are mostly affected by grid size, surface roughness, and the applied radiation 10 schemes.


Introduction
The land surface influences atmospheric dynamics significantly through the exchange of energy, mass, and momentum. Therefore, an accurate representation of surface-atmosphere interactions is essential for any numerical modeling of the atmospheric boundary layer (Garratt, 1993;Betts et al., 1996). More specifically, surface roughness as well as sensible and latent heat fluxes 15 at the surface act as the lower boundary condition for the momentum, temperature and humidity equations in the atmosphere, respectively. If this information is not available, it is necessary to parameterize the land-surface processes with a Land-Surface Model (LSM). In simple terms, the input to an LSM is the type of surface, vegetation and soil, as well as the radiative forcing.
Based on that, LSMs solve the surface energy budget equation by means of a set of prognostic equations for the surface temperature and compute soil moisture and temperature in a multi-layer soil model. Nevertheless, the results strongly depend on 20 the input data (e.g. Avissar and Pielke, 1989).
LSMs are required in various situations: Most often, observational data of the sensible and latent heat fluxes are unavailable or do not adequately reflect the complexity of the landscape. E.g. in case of forecasts, it is trivial that a prognostic approach has details, see also Viterbo and Beljaars (1995); van den Hurk et al. (2000); Balsamo et al. (2009Balsamo et al. ( , 2011 and literature referenced therein.

Energy balance solver
Within the LSM, the energy balance is calculated as 70 where C 0 and T 0 are the heat capacity and radiative temperature of the surface, respectively. Note that C 0 is zero by default in case of surfaces covered by vegetation or water surfaces, where it is assumed that T 0 is the temperature of a skin layer covering the surface that does not have a heat capacity. In all other cases (i.e., pavements and bare soils), no skin layer is assumed (see below). R n , H, LE, and G are the net radiation, sensible heat flux, latent heat flux, and ground (soil) heat flux at the surface, respectively. R n is defined positive downwards whereas H, LE, and G are defined positive away from the surface. R n is 75 defined through the sum of the radiative fluxes: where SW ↓ , SW ↑ , LW ↓ , and LW ↑ are the shortwave incoming (downward), shortwave outgoing (upward), longwave incoming (downward), and longwave outgoing (upward) flux, respectively. The radiation components are defined positive according to their direction (SW ↓ and LW ↓ positive downwards; SW ↑ and LW ↑ positive upwards). The radiative fluxes are provided by 80 one of the available radiation schemes in PALM (for details, see Maronga et al., 2020).

Parameterization of fluxes
The turbulent heat fluxes both are parameterized using a resistance parameterization. H is calculated as where ρ is the density of the air, c p = 1005 J kg −1 K −1 is the specific heat at constant pressure, and r a is the aerodynamic 85 resistance. θ 0 and θ mo are the potential temperature at the surface and at a fixed height within the atmospheric surface layer (at height z mo , usually at height of the first atmospheric grid level, i.e., z mo = 0.5∆ z where ∆ z is the vertical grid spacing), respectively. The potential temperature is linked to the actual temperature via the Exner function, viz.
with pressure p and the gas constant for dry air R d . r a is calculated via Monin-Obukhov similarity theory (MOST) as where u * and θ * are the friction velocity and characteristic temperature scale, respectively, and which are calculated locally based on MOST. Note that the values for u * and θ * from the previous time step are used in Eq. (5). The roughness lengths are individually set for momentum, heat, and moisture (see Table 1). Note that for water surfaces, a Charnock parameterization can be switched on for taking into account the effect of subgrid-scale wave motions through the roughness lengths. For details 95 on the particular implementation of MOST and the Charnock parameterization in PALM, see Maronga et al. (2020). Note that r a is calculated based on u * and θ * values from the current time step to calculate H at the prognostic time step.
The ground heat flux, G, is parameterized after (Duynkerke, 1999) as with Λ being the total thermal conductivity between skin layer and the uppermost soil layer. T soil,1 is the temperature of the 100 uppermost soil layer (calculated at the center of the layer). Λ is calculated via a resistance approach as a combination of the conductivity between the canopy and the soil-top (Λ skin , constant value) and the conductivity of the top half of the uppermost soil layer (Λ soil ): When no skin layer is used (i.e. in case of pavements and bare soils), Λ reduces to the heat conductivity of the uppermost soil 105 layer, viz.
with λ T,pave being the thermal conductivity of the pavement and ∆z soil,1 being the depth of the uppermost soil layer. In this case, it is assumed that the soil temperature is constant within the uppermost 25 % of the top soil layer and equals the radiative temperature at the surface. C 0 is then set to a non-zero value according to the material properties and the layer thickness.

110
The total latent heat flux, LE, is parameterized as r a + r s (q v,mo − q v,sat (T 0 )) .
Here, l v = 2.5 × 10 6 J kg −1 is the latent heat of vaporization, r s is the total surface resistance, q mo is the water vapor mixing ratio at height z mo , and q sat is the water vapor mixing ratio at saturation at the surface, which is a function of T 0 . In practice, parameterized as where r c is the canopy resistance. Analogous the bare soil fraction evaporation (LE soil ) is calculated via with r soil being the soil resistance. The liquid water reservoir evaporation (LE liq ) is given by i.e., only the aerodynamic resistance exists for liquid water. The total evapotranspiration is then given by a combination of the three individual components by (see Viterbo and Beljaars, 1995) Here, c veg and c liq are the fractions of the surface covered with vegetation and liquid water, respectively. Liquid water from 125 precipitation can be stored on the vegetation and bare soil. Note that for paved and water surfaces both LE veg and LE soil are set to zero and the only possible source of evaporation is the liquid water reservoir.
All equations above are solved locally for each surface element of the model grid.

Liquid water reservoir
In order to account for the evaporation of liquid water on plants on impervious surfaces, an additional equation is solved for 130 the liquid water reservoir: where m liq and ρ l are the water column on the surface and the density of water, respectively. The maximum amount of water that can be stored on plants is calculated via

Calculation of resistances
The resistances are calculated separately for bare soil and vegetation following Jarvis (1976). The canopy resistance, r c , is calculated as with r c,min being a minimum canopy resistance. f 1 − f 3 are correction functions depending on LAI, the incoming shortwave 145 radiation (SW ↓ ) and the water-vapor pressure deficit (e def = e sat − e, with e sat and e being the water-vapor pressure at saturation and the current water-vapor pressure, respectively). The layer-averaged volumetric soil moisture content ( m) is given bỹ where N is the number of soil layers, R fr,k is the root fraction in layer k, m soil,k is the volumetric soil moisture content in 150 layer k, and m wilt is the permanent wilting point.
The correction functions f 1 and f 2 read which accounts for the reaction of plants to sunlight (opening/closing stomatas); the reaction of plants to water availability in the soil is considered via with m fc being the soil moisture at field capacity. Furthermore, a correction for the water vapor pressure deficit is given by where g D is a correction factor that is used for high vegetation.
The soil resistance (r soil ) is calculated as where r soil,min is the minimum soil resistance. The correction function f 4 is given by with m min being a minimum soil moisture for the soil matrix based on the wilting point and the residual moisture m res , calculated as Note that the total surface resistance (r s , cf. Eq. (9)) is calculated as a diagnostic quantity from LE after the energy balance is solved.

Soil model
The layer. Viterbo and Beljaars (1995) and Balsamo et al. (2009) give more details

Soil heat transport
The Fourier law of diffusion reads with (ρC) soil and λ T being the volumetric heat capacity and the thermal conductivity of the soil layer in question, respectively. 180 λ T is calculated as with λ T,sat , λ T,dry , and Ke being the thermal conductivity of saturated soil, of dry soil and the Kersten number, respectively.
λ T,sat is given by

185
Here, λ T,sm is the thermal conductivity of the soil matrix and λ m is the heat conductivity of water. The Kersten number (Ke) is calculated as At the bottom boundary a fixed deep soil temperature T deep is prescribed (Dirichlet conditions), which is a plausible assumption for short term simulations covering only a few days.

Soil moisture transport
The vertical transport of water within the soil matrix is calculated using Richards' equation, viz.
where λ m , γ, and S m are the hydraulic diffusion coefficient, hydraulic conductivity, and a sink term due to root extraction, respectively. The hydraulic diffusion coefficient is calculated after Clapp and Hornberger (1978) as with b = 6.04 being a fixed parameter, γ sat being the hydraulic conductivity at saturation, and Ψ sat = −338 m being the soil matrix potential at saturation. The hydraulic conductivity (γ) is calculated after van Genuchten (1980) (as in HTESSEL): Here, α, n, and l are van Genuchten coefficients that depend on the soil type (see Table 2). h is the pressure head, which is The root extraction of water from the respective soil layer S m,k is calculated as follows: where m total is the total water content of the soil,

Treatment of pavements
Pavements are treated as a common soil (allowing varying number and depths of the pavement layers) but with physical 215 properties of the pavement material. The pavement layer is impermeable to water and prohibits the vertical transport of soil moisture. Soil layers are placed below the pavement layers.

Treatment of water bodies
For water surfaces, PALM currently only allows for prescribing a bulk water temperature. The energy balance is then solved as for land surfaces, but without evapotranspiration from vegetation and bare soil (see above). A skin layer is adopted so that 220 C 0 = 0 and Λ = 1 × 10 11 in order to calculate the heat flux into the water body.

Numerical methods
In order to solve for the energy balance for the surface temperature (T 0 ), Eq. (1) is first linearized around T 0 at the current time step and then discretized in time using PALM's default Runge-Kutta third-order time stepping scheme. In this way, an iterative procedure to solve the energy balance is avoided and the prognostic equation then reads where t is the time index, ∆t is the current time step. A and B are coefficients given by and
site at Cabauw. The period was chosen because the forcing from the surface was dominant and larger-scale advection played a minor role. We used direct measurements of temperature, humidity and wind and derived observations of sensible, latent and ground heat flux as well as net radiation. CESAR features a 213 m high measurement mast with instruments in 1. 5, 10, 20, 40, 80, 140 and 200 m (Bosveld, 2020b). For temperature, humidity and wind, 10-minute averages are derived with a measurement accuracy for temperature and humidity of 0.1 K and 3.5 %, respectively (Meijer, 2000). Soil temperature is observed at a depth The surface soil heat flux is derived from soil temperature measurements by means of a Fourier extrapolation. Net radiation is derived from the budget of the four radiation components (Eq. 2). Turbulent fluxes of sensible and latent heat are derived by 245 means of the eddy-covariance (EC) method. In the process of calculating EC-fluxes from raw turbulence data it is unavoidable that low frequency contributions to the flux are getting lost. For the current Cabauw data, 10 minute intervals are used with simple subtraction of the mean from the raw turbulent timeseries; afterwards a low frequency correction is applied based on the spectra by Kaimal et al. (1972), taking into account wind speed and stability (Bosveld, 2020a). The method of low-frequency loss correction assumes that all turbulence characteristics follow surface-layer scaling. However, this is not always true, as 250 for example horizontal advection by organized turbulent structures (Eder et al., 2015;De Roo and Mauder, 2018) may add further low-frequency contributions which is not accounted for in the surface-layer scaling. For further information on the instrumentation of the Cabauw site we refer to Bosveld (2020b).
Cabauw is located in the western part of the Netherlands and surrounded by mostly meadows with ditches passing through as well as villages, orchards and lines of trees. The CESAR tower itself is installed over an area of short grass that is kept at 255 a height of about 8 cm. The immediate surroundings of the measuring tower are free of significant heterogeneities for a few hundreds of meters. During the simulation period from 5 th to 6 th of May, 2008 the prevailing wind direction is from south-east and the 10 m average wind speed ranges from 2 to 6 ms −1 . The convective boundary layer reaches a height of around 2 km.
The groundwater level is in 1.3 m depth. The profiles of temperature θ, humidity q and horizontal wind v h retrieved from radiosounding are shown in Fig. 2 together with tower measurements from all available levels (see details in Sect. 5). The soil 260 temperature T soil and soil moisture m soil at 2008-05-05 05:00 UTC is depicted in Fig. 1 (see Sect. 4).

Simulation setup
To evaluate a land-surface parameterization scheme the relevant vegetation and soil information is required. With regards to this, the CESAR site is well described in literature (e.g. Beljaars and Bosveld, 1997). The reference simulation (hereafter referred to as case REF) was set up as a best guess with the Cabauw land surface parameters according to Beljaars 265 and Bosveld (1997), Ek and Holtslag (2004) and Maronga and Bosveld (2017). The vegetation type of the surface is 'short grass' with some modifications. The roughness length for momentum is set to 0.15 m, which is representative for a few kilometers of upstream terrain from the Cabauw tower and the roughness length for temperature is set to 2.35 × 10 −5 m (Ek and Holtslag, 2004). The leaf area index is 1.7 m 2 m −2 , the minimum canopy resistance 110 s m −1 and we chose a heat conductivity between skin layer and soil of 4 W m 2 K −1 . The heat capacity of the skin layer is set to 0 J m −2 K −1 ,  Table 2), the layer between 24-60 cm is the same as the uppermost layer except for the porosity which had to be increased due to observed large values of soil moisture, and between 60-225 cm the parameters are set to organic soil according to Table 2. In all simulations, the land surface and soil parameters are homogeneous over the model domain. This means that buildings of the 280 small town of Lopik, west of the CESAR tower, are neglected, as well as the small ditches which cross the observation site.
The root fraction and initial soil profiles of temperature and moisture are shown in Fig. 1. In 3 cm depth begins a layer of relatively high root density R dens down to 18 cm followed by a layer of relatively low root density down to 60 cm depth.
According to Jager et al. (1976), no roots are found near the surface (< 3 cm) and in the deep soil layers (> 60 cm). The initial soil temperature and moisture profiles are taken from measurement data at 2008-05-05 05:00 UTC (shown in Fig. 1). A 285 summary of the land surface scheme configuration and its used values are listed in Table 4.
Case REF is driven by external forcing of incoming short-and longwave radiation taken from the Cabauw measurements.
In this way, the effects of high altitude aerosols, moisture and clouds are included. Accordingly, the degree of freedom is reduced and we can focus on parameters of the LSM rather than additional uncertainties of a radiation model. Nonetheless, we performed sensitivity tests using the Rapid Radiation Transfer Model for Global Models (RRTMG, Clough et al., 2005) as 290 well as a clear-sky radiation parameterization, which are described in detail in Maronga et al. (2020). The longwave outgoing radiation of the surface is calculated from the skin-layer temperature using Stefan-Boltzmann law. The long-and shortwave albedos of diffusive radiation are set to 0.34 and 0.14, respectively, to fit the dominating grassland. Albedos of direct radiation are calculated according to Briegleb (1992) considering a weak solar zenith-angle dependence as such that their direct values equal the diffusive ones at 60 • .   All sensitivity simulations are based on the setup of the reference simulation and only differ by the respective parameter to be analyzed. An overview of the sensitivity simulations and their well-defined change compared to case REF is given in Table 5.

Boundary-layer profiles
At first, we will look at the vertical profiles of the radiosonde and the tower measurements. Figure 2 shows vertical profiles of θ, q and v h indicating the evolution of the boundary layer during the considered period of time in Cabauw at the times of the radiosonde ascents. Both nights show a stable nocturnal boundary layer before sunrise (at 05:00 UTC). A nighttime 310 low-level jet is observed in the horizontal wind. Profiles of potential temperature and water vapor mixing ratio show a wellmixed convective boundary layer at 10:00 UTC and 16:00 UTC of either day. Over the course of the two days, the depth of the shortwave albedo (at 60 • ) of 0.24

ALBE_44
shortwave albedo (at 60 • ) of 0.44 ADV_tq advection of T and q in all heights according to mean change in radiosounding data between 2.5 km and 4 km convective boundary layer increases. The profiles of v h show a mismatch between the radiosounding measurements and the tower data. This is explained by different analysis timescales. One the one hand, the radiosonde records instantaneous values with a sampling frequency of 0.1 Hz, which results in five recording heights below 300 m. On the other hand, data from the cup 315 anemometers at the mast is temporally averaged over 10 min, based on a 3 s running mean calculated with an update frequency of 4 Hz (Bosveld, 2020b). Thus, the mean tower profiles can be over-or underestimated by the instantaneous values of the radiosonde. Note that we initialize the simulations with data from the radiosonde ascent, because it reaches high enough, but the comparison later is made against tower data due to its higher temporal resolution. shows no good agreement between model and observations. At 16:00 UTC of the first day, the simulation slightly overestimates the boundary layer depth despite a fairly good agreement of model and observations regarding mixed-layer temperature and 325 humidity values. The wind profiles, however, agree well. At night (2008-05-06 05:00 UTC), the near-surface temperature is significantly lower (out of range at ca. 281 K) than measured. Another difference is that the nocturnal boundary layer is too shallow in the simulation. At the same time, the near-surface humidity shows small differences between model and reality. The simulated horizontal wind speed also depicts a low-level jet, but compared with observations it occurs closer to the surface, in accordance with the simulated nocturnal boundary-layer depth. Even though the LES cannot reproduce the nocturnal situation 330 very well, the mixed-layer quickly develops in the next morning (2008-05-06 10:00 UTC) which is in agreement with findings of van Stratum and Stevens (2015). At 10:00 UTC of the second day the simulation indicates a warmer and deeper boundary layer compared to the observation, which could be caused by advection processes in reality modifying the residual layer and thus the boundary-layer evolution during the morning transition. The wind profiles agree fairly well. The temperature profile at 2008-05-06 16:00 UTC shows that the simulation underestimates the mixed-layer temperature but instead has a higher 335 boundary-layer top, suggesting that the energy input into the boundary layer is similar, but distributed over a deeper layer.
The humidity profiles agree well in the lowest 1000 m, but deviate above in accordance with the differences in boundary-layer depth. The horizontal wind is slightly overestimated. In general, the simulated profiles are much more constant with height in the well-mixed layer, because they depict domain averages, as opposed to the local measurements, so that a direct comparison is inherently improper. Above the boundary layer, in the free atmosphere, synoptic-scale processes dominate in reality. Since 340 we did not consider these in our simulations, profiles may deviate. Given the fact that the surface forcing in this particular case was the dominant forcing for the development of the boundary layer, this deficiency should not compromise the present evaluation study. included and can also be seen in the other surface fluxes (cf. Fig. 4, 5, 6). Figure 3 also depicts the effect of different landsurface properties and simulation setups on the surface radiation. Please note that we will only highlight the cases which lead to the largest differences compared to case REF for the respective variable. In Fig. 3, these are mostly the changes to the albedo and the radiation models as well as cases LAI_05 and HUMID_sat. The most obvious differences are seen in case of a changed

Evaluation of energy-balance components
albedo. An increase of the albedo (case ALBE_44) leads to a decrease in R n and vice versa (case ALBE_24) with maximum 355 deviation to case REF of about ±100 W m −2 at noon, indicating that mismatches in the estimated albedo cause significant errors on the surface net radiation during daytime. However, also the spread of the non-highlighted simulations (gray) is about 50 W m −2 at noon, meaning that also variations in the surface parameters (e.g. emissivity, roughness length) affect the surface net radiation significantly. Using the RRTMG, radiative fluxes are calculated for each horizontal grid box of the LES based on time and geographical location as well as air pressure, and local profiles of temperature and humidity (Clough et al., 2005) 360 instead of being taken from the measurements. One key feature of the RRTMG, which is neither included in observations nor in the clear-sky model, is the direct cooling of air due to longwave radiative flux divergence, particularly during nighttime. This can result in cooling rates in the order of 0.1 to 0.3 K h −1 (see, e.g. Maronga and Bosveld, 2017, their Fig. 1). In case RRTMG, the simulated R n approaches that of the observations during day and night, though it slightly underestimates the surface net radiation at noon by about 30 W m −2 . However, we have to note that we use default input of the RRTMG with 365 standard profiles of water vapor, other trace gases, and aerosol concentration above model top, which do not necessarily reflect the real conditions at Cabauw during the simulated time period and may impact the simulated surface net radiation, too. Using the clear-sky model, the simulated R n becomes even better during day, but also marginally overestimates net radiation during night by 10 to 20 W m −2 . Case LAI_05 has a smaller R n compared to case REF, because a smaller LAI significantly decreases LE and increases H, which leads to higher surface temperatures during day resulting in higher LW ↑ . In case of a saturated 370 humidity mixing-ratio (case HUMID_sat) LE, H and G are significantly altered in a way that higher surface temperature is simulated during day and night. This again leads to higher LW ↑ and hence lower R n . Except for ALBE_24 the simulations tend to underestimate the surface net radiation during daytime while the simulated surface net radiation tend to overestimate the observed one during night.    (Fig. 4). Moreover, case HUMID_sat stands out because it shows a negative LE during night, which is explained by dew formation in the model.     Table 5.
Only case REF and some relevant cases are highlighted, all others are shown in gray.
lighted. From Fig. 4 and Fig. 5 we can calculate the daytime Bowen ratio β 0 = H LE (not shown). We find that the difference between the Bowen ratio of the simulations and that observed at Cabauw (β 0,sim − β 0,obs ) ranges from −0.1 to +0.2, except for cases HUMID_sat and LAI_05, which show significantly lower Bowen ratios compared to observations. As opposed to other models, which were shown to systematically overestimate the Bowen ratio measured in Cabauw during summer (Schulz et al., 1998;Sheppard and Wild, 2002), we cannot identify a systematic bias for Bowen ratio of PALM's LSM by means of the 400 analyzed period. Figure 6 shows the timeseries of the ground heat flux, which reaches 60 W m −2 during daytime in the observations. With respect to the amplitude of G, case REF shows good agreement with the observation. However, the shape of the timeseries is discernibly different. At daytime, OBS shows a more sinusoidal shape, whereas the simulations show a more humped-shaped diurnal variation. We attribute this to the method to derive the observed ground heat flux, which involves the average soil heat 405 flux in 10 cm and 5 cm depth and the soil temperature difference between 0 and 2 cm (Bosveld, 2020b). In the model, on the other hand, the ground heat flux is parameterized according to Eq. (6) and thus only the temperature gradient between the surface and 0.5 cm are taken into account. Hence, the simulated G resembles the diurnal cycle of the surface temperature, while the observed G correlates with the soil temperature. Observations and case REF agree well at noon with values G around 55 W m −2 , whereas in the afternoon at 15:00 UTC case REF is almost 10 W m −2 higher than the observations. It strikes that all 410 simulations have a fairly small spread at that point in time. During the day, changing the heat conductivity between skin layer and uppermost soil layer has the largest impact on the ground heat flux, with smaller and larger G observed in case COND_2 and COND_6, respectively. This can be attributed to the linear relationship between G and Λ in Eq. (6). Also, the simulation cases with a different LAI show a relatively large deviation from case REF at noon, with increased (decreased) G for LAI_05 (LAI_3). For example, with a smaller leaf area and thus less transpiration the available energy is preferentially partitioned 415 into G and H (see Fig. 4) rather than into LE. At noon, the spread among all simulations is about 40 W m −2 , whereas the non-highlighted cases show a maximum deviation from REF of no more than ±10 W m −2 . In the night, cases HUMID_sat and RRTMG stand out. In both cases the controlling variable is the (near-)surface temperature (see Fig. 9). In case RRTMG, the surface cools down while the soil temperature does not change in the same amount (see Fig. 8), which leads to a strong ground heat flux directed towards the surface. In case HUMID_sat, the surface stays relatively warm, therefore it results in a 420 small negative G. Except for cases HUMID_sat and COND_2, the model tends to simulate a stronger upward directed ground heat flux at night. Compared to the observations, case COND_2 overestimates G during night whereas it underestimates G at noon, which points to a stability dependence of the conductivity λ. Even though this is technically realized in the code, the standard values for short grass do not differ between stable and unstable conditions due to a lack of knowledge of the correct relationship between λ and stability. As cases COND_2 and COND_6 do not have a significant impact on any of the other 425 variables, we can infer that uncertainties in modeling G is relatively small. Nonetheless, please note that we only analyzed a short period of two days. For longer simulation times covering multiple days, e.g. for heat-wave scenarios where heat storage within the soil becomes important, this might become relevant again. Furthermore, it shall be mentioned that the ground heat flux is small compared to the surface net radiation as well as the surface sensible and latent heat fluxes.
Similar to many other sites, also eddy-covariance measurements at the Cabauw site suffer from the well-known problem 430 of energy-balance non-closure (de Roode et al., 2010). To date, the leading hypothesis is that low-frequency contributions are inherently not captured by the eddy-covariance method, leading to the situation that the sum of the surface heat fluxes underestimates the available energy by 10 to 30 % (Foken et al., 2011). Figure 7 shows the timeseries of the residual as well as the individual differences between the simulated (case REF) and observed energy balance components. Note that the shown residual is effectively the residual of the observations, because the energy balance in the model is zero at all times which indicates that the simulated H and LE are overestimated compared to the observations. In reverse, however, this may also suggest that the observed H and LE are underestimated, even though we explicitly note that we cannot know from the 440 simulations how the missing energy is partitioned onto the individual energy-balance components in reality as both, observation and simulation may contain a bias. Taking this into account, we summarize that the four components of the energy balance are represented reasonably well by the LES-LSM interface. Also, the daytime Bowen ratio agrees fairly wells with that observed at Cabauw. By contrast, climate models have often shown to overestimate the Bowen ratio in summertime (Schulz et al., 1998;Sheppard and Wild, 2002).  that there are doubts about the quality of the CESAR soil temperature observations, because of problems during calibration (Bosveld, 2020b) and the speculation that the sensors had sunken deeper into the soil over time (Bosveld, 2020a). Hence, observations cannot be used as a reliable reference to evaluate the model. Nevertheless, we will discuss the variation among the different simulations. In the late afternoon of the first day (18:00 UTC), case REF shows a continuous decrease with depth 450 from 288 K close to the surface to 283.5 K in −50 cm. During night, the soil begins to cool down starting at the surface, where it reaches 285 K at 04:00 UTC. The maximum soil temperature is now found in a depth of −15 cm. As soon as the net radiation becomes positive (at 08:00 UTC), the soil temperature close to the surface follows such that the maximum is again found close to the surface. The parameters with the largest spread -up to 4 K during day and 2 K during night -are, like for G, the conductivity and atmospheric humidity. The non-highlighted cases have a spread of about 1 K indication that e.g. changing 455 radiation or vegetation parameters have only minor effect on the soil temperature. A higher conductivity (COND_6) leads to higher soil temperatures, especially in the late afternoon (18:00 UTC), because the heat of the surface net radiation is more easily conducted into the soil (cf. Fig. 6); correspondingly COND_2 has smaller T soil . During night, the lower atmosphere in case HUMID_sat does not cool down as much as case REF (cf. Fig. 9), therefore the soil remains relatively warm, too.

Evaluation of atmospheric quantities 460
Next, we will evaluate the quantities characterizing the lower atmosphere, that is temperature, humidity, and wind speed. First, we will discuss the nighttime situation followed by the daytime situation. Finally, results from the sensitivity study will be discussed. Figure 9 shows the vertical profiles of potential temperature during night. At 18:00 UTC, the lowest levels start to

08
:00 UTC, observations shows an already vertically well-mixed lower boundary layer, while in the simulations a stable layer at about 100 m is still present and gets eroded due to surface heating and mixing processes. This is consistent with Fig. 11 470 which indicates a faster temperature increase at z = 40 m compared to the observation.
One reason for the misrepresentation of the nocturnal stable layer could be too coarse grid spacing. However, compared to case REF, case Dz_2 with smaller grid spacing is even more stable during the night hours, having a lower boundary layer depth and cooler near-surface temperature in agreement with van Stratum and Stevens (2015). This suggests that the coarse resolution might not explain the misrepresentation of the stable boundary layer, nonetheless with only one grid sensitivity simulation, a 475 non-linear relationship cannot be fully ruled out. Another parameter affecting the diffusivity is the subgrid-scale scheme of the LES. Yet, the subgrid-scale scheme of Dai et al. (2020), which is more diffusive in the middle of the stable boundary layer but less diffusive towards the surface, does not change the stability of the simulation (not shown). Alternatively, too low diffusivity could be due to low wind speed, however the simulated wind speed in Fig. 13  Furthermore, compared to the observation, the simulations indicate a less negative H (see Fig. 4), though we note that the 485 simulated H shows a domain-averaged value while the observed value is taken from a point measurement with more temporal fluctuation. Nevertheless, according to the less negative H, the near-surface layer should be less stable in the simulations, which, however, is in contrast to the profiles shown in Fig. 11, which would suggest that the simulated H should show larger negative values.  occurs in a height above the tower. According to the radiosounding measurements at 2008-05-06 05:00 UTC (Fig. 2), the maximum wind speed occurs in approximately 500 m. In the LES (case REF), however, the low-level jet has its maximum between 60 m and 70 m height (cf. Fig. 2, 10). In the morning, when convective mixing sets in, the wind profiles become more 495 well-mixed. At 08:00 UTC, we already find a mixed-layer in reality. By contrast, the stable layer has not been fully eroded in the LES, therefore the wind speed in the residual layer (above ca. 100 m) is still higher than that of the lower 100 m. Even though case RRTMG shows a large influence on the whole temperature profile (Fig. 9), its wind profile only deviates from case REF above the low-level jet and during the morning transition where the stable layer is already further eroded. Likewise, case HUMID_sat does not deviate much in Fig. 10. Conversely, case ROUGH_001 deviates significantly in the wind profiles but 500 does not show significant changes in the temperature profiles (Fig. 9). Only case Dz_2 shows a strong interconnection between temperature and wind profiles.  The CESAR tower samples temperature, humidity and horizontal wind speed in 10, 40, and 80 m height (wind speed not available at 1.5 m). We chose to compare timeseries of the 40 m level, because on the one hand, 10 m lies between the first and second grid level and on the other hand, 80 m might already be above the nocturnal boundary layer (see discussion of Fig. 9).
505 Figure 11 shows the diurnal cycle of temperature at z = 40 m which, at first sight, agrees fairly well between observations and case REF. The temperature amplitude is around 10 K and between 16:00 and 17:00 UTC peak temperatures of 295 K and 296 K are reached on the first and second day, respectively. Significant differences are, however, found during night as well as during the second day: During night, case REF becomes more stable compared to the observations (Fig. 9), hence, the temperature in 40 m height is lower compared to the observations, as discussed earlier. The following morning, the boundary layer is heated 510 up quickly between roughly 06:00 UTC and 08:00 UTC in the simulation, whereas in reality temperature increases slower.
Again, the reason is that case REF develops a much shallower stable boundary layer (with a residual layer above, cf. Fig. 2).
Once the stable layer is eroded up to 40 m the temperature can increase rapidly at that level as it is heated from the surface as well as from above, where warmer air from the residual layer is mixed into the shallow layer. In reality, where the stable layer is much deeper, the air in 40 m height is first only heated from the surface and not by entrainment. The sensitivity study shows 515 that during night case HUMID_sat is in agreement with the observations at that height (cf. Fig. 9). This is because the humidity reaches saturation and condenses on the vegetation as dew, and heat from condensation is released to the air by different partitioning of the surface fluxes. However, dew formation was not observed at the CESAR site as shown by timeseries of LE ( Fig. 5), thus humidity is not suitable to explain the misrepresentation of near-surface temperature. In Fig. 11, also case Dz_2 seems to agree with observations. However, as shown in Fig. 9, this is only true for the depicted height of 40 m, e.g. below 20 m, 520 the nocturnal air temperature is lower compared to case REF. Even though case Dz_2 has a much higher resolution, we cannot be certain that turbulence in the nocturnal boundary layer is sufficiently resolved and hence the resolution is still a possible explanation for the misrepresentation of the stable boundary layer in the LES. Another parameter to point out in Fig. 11   cooling of the air volume at night but similar boundary-layer heating. This misrepresentation propagates into the next day. Figure 12 shows that the observed humidity in 40 m has only small diurnal variations. The reference simulation mostly agrees with the observations. The small decrease in the measurement data around 09:00 UTC of the first day can be reproduced by case ADV_tq, which includes advection tendencies of θ and q according to the mean change well above the boundary layer. However, during the morning transition of the second day (around 08:00 UTC) the humidity in the simulations first rises 530 and then drops while near-surface air is mixed with dry air from above. By contrast, this is not observed in the measurement data, most likely because no strong stable layer had developed during nighttime. to higher wind speed close to the surface. Reducing the grid size (case DZ_2) has the same effect mostly shown at night.

Summary
In this paper we gave a description of the land-surface model embedded into the PALM model system which is applied to model the surface energy balance at vegetated or paved land surfaces as well as at water surfaces. We evaluated the land-surface model 550 implementation against in-situ observations of the energy-balance components as well as near-surface wind-, temperature-, and humidity-profiles taken at Cabauw over a quasi homogeneous flat grass site (Monna and Bosveld, 2013) for two consecutive diurnal cycles. A sensitivity study showed the relative importance of the choice of land-surface input parameters and thereby gives valuable reference for the user.
The diurnal cycles of surface latent and sensible heat flux are well represented. Even though the model seems to overestimate 555 the fluxes, the differences can be explained by the non-closure of observed energy balance, i.e. missing energy in H and LE observed at the CESAR site. During daytime the simulated Bowen ratio agrees reasonably well with the observed one, whereas climate models overestimate the summer Bowen ratios observed in Cabauw (e.g. Schulz et al., 1998;Sheppard and Wild, 2002). The diurnal cycle of the modeled ground heat flux agrees with the observations, though the modeled flux overestimates the observed flux during the morning and evening hours. During nighttime, the modeled ground heat flux shows slightly larger 560 negative values compared to the observation. Due to its relatively small contribution compared to the surface net radiation or the surface latent and sensible heat fluxes, these mismatches do not affect boundary-layer development too much, if only one or two days are simulated. However, for longer simulation periods heat storage in the soil may become an important factor, e.g.
when heat waves built-up over days (Miralles et al., 2014). During daytime the model tends to underestimate the surface net radiation, while during nighttime the surface net radiation shows slightly less negative values compared to the observations, 565 indicating underestimated longwave radiative cooling. The near-surface temperature matches well with the observed one at the first simulated day. During nighttime, however, the near-surface temperature is underestimated in the model, the nocturnal boundary layer is too stable and too shallow compared to the observations. The diurnal cycle of the near-surface wind is well represented, though the low-level jet in the model occurs much closer to the surface.
In addition to the evaluation against observations, we carried out a comprehensive sensitivity study. Land-surface parameters 570 and the initial state of the atmosphere were varied within a typical range for the respective quantity. The net radiation is significantly influence by the albedo, the radiation model as well as many of the land surface parameters. The ground heat flux, though not as important as the other energy-balance components, is mostly influenced by the soil conductivity. The distribution of the available energy into surface sensible and latent heat fluxes depends mostly on the leaf area index as well as the initial atmospheric humidity. Within the investigated range of LAI values for short grass, differences of up to 50 % are possible.

575
Moderately less relative deviation is found for the range of completely dry to saturated initial atmospheric humidity. Potential temperature is influenced by a change in LAI of ±2 m 2 m −2 as much as it is influenced by an initial temperature deviation in the whole atmosphere of ±1 K. While most of the sensitivity studies show the most significant difference during day, the choice of initial humidity, grid size, roughness length and radiation model plays an important role at night. Overall, we could not identify a single parameter as being the most sensitive in all quantities at the same time. In fact, different parameters become 580 relevant if different quantities are analyzed.
In order to evaluate and possibly improve land-surface schemes also for different types of surfaces, e.g. pavements or different vegetation coverage, reliable measurements of the energy-balance terms at different surfaces types are required.
However, eddy-covariance measurements often suffer from the well-known problem of energy-balance non-closure, where the sum of surface sensible, latent, and ground heat flux underestimates the surface net radiation by about 10 to 30 % (Foken et al.,585 2011). Besides measurement uncertainties and footprint biases, one leading hypothesis is that low-frequency contributions from organized turbulent structures or surface heterogeneity-induced circulations are inherently not captured by the eddy-covariance method (e.g. Finnigan et al., 2003;Foken, 2008;Eder et al., 2015). As the modeled energy balance is closed by definition, it is thus difficult to draw final conclusions that may point directions for further improvements of the land-surface parameterizations as both, the model as well as the observation may contain a bias.

590
As the description of the land-surface model embedded in PALM only reflects its current state, a short outlook into future development is given below. Until now, the LSM implementation does not incorporate a tile approach (as the embedded building-surface model (Resler et al., 2017) does), such that a land-surface grid cell in PALM is either classified as water or as pavement-/vegetation-covered. Particularly for coarser grids, however, patchy landscapes such as e.g. with small ditches like in Cabauw, might be filtered, meaning that the relative contributions of surface types and thus the area-averaged energy-balance 595 terms become a function of the horizontal grid size. In order to avoid this, one of our next steps will be to implement a tile approach into the land-surface model. Furthermore, the current LSM implementation does not include heat storage within water bodies. However, especially for multi-day simulations, e.g., for heat-wave scenarios, heat storage in water bodies might become important to accurately represent the cool-air production in urban environments. Hence, we plan to improve the representation of water surfaces by implementing a lake parameterization (e.g. Mironov et al., 2010). Further lines of future development will 600 be the implementation of a snow parameterization as well as a parameterization to consider phase transitions to also consider frozen soil. In this study we only considered a homogeneously flat surface. Nonetheless, heterogeneous surfaces as well as step-like orography is implemented in the LSM. In addition, we plan to implement an immersed boundary method (Mason and Sykes, 1978) where elevation changes can also be represented by slanted surfaces.
Code and data availability. The PALM model system is freely available from http://palm-model.org and distributed under the GNU Gen-