Validation of the Dynamic Core of the PALM Model System 6.0 in Urban Environments: LES and Wind-tunnel Experiments

We present the capability of PALM 6.0, the latest version of the PALM model system, to simulate neutrally stratified urban boundary layers. The studied situation includes a real-case building setup of the HafenCity area in Hamburg, Germany. Simulation results are validated against wind-tunnel measurements of the same building layout utilizing PALM’s virtual measurement module. The comparison reveals an overall very high agreement between simulation results and wind-tunnel measurements not only for mean wind speed and direction but also for turbulence statistics. However, differences between 5 measurements and simulation arise within close vicinity of surfaces where the resolution prevents good representation of the inertial sub-range of turbulence. In the end, we give guidance how these differences can be reduced using already implemented features of PALM.


Introduction
The PALM model system version 6.0 is the latest version of the computational fluid-dynamics (CFD) model PALM, a Fortran 10 based code to simulate atmospheric and oceanic boundary layers. Version 6.0 was developed within the scope of the Urban Climate Under Change ([UC] 2 ) framework funded by the German Federal Ministry of Education and Research (Scherer et al., 2019;Maronga et al., 2019). The aim of the [UC] 2 project is to develop a fully functional urban climate model capable of simulating the urban canopy layer from city scale down to building scale with grid sizes down to 1 m or less. A detailed description of the model system is given by Maronga et al. (2015Maronga et al. ( , 2020. PALM has already been applied to a variety of 15 studies within the area of urban boundary-layer research (e.g., Letzel et al., 2008;Park et al., 2012;Kanda et al., 2013;Kurppa et al., 2018;Wang and Ng, 2018;Paas et al., 2020). Built upon PALM version 4.0, the latest version 6.0 contains many new features and improvements of already existing components of the model system. One of the most impacting changes is the new treatment of surfaces within PALM. While PALM did not distinguish between different surface types within former versions, it is now possible to directly specify a surface type to each individual solid surface within a model domain via the 20 land-surface model (Maronga et al., 2020) or the building surface model (Resler et al., 2017;Maronga et al., 2020). Also, a fully three-dimensional obstacle representation is possible while former versions allowed only a 2.5-dimensional representation of obstacles (no overhanging structures like bridges or gates). These additions, however, required extensive re-coding of the former version PALM 4.0 which affected also the dynamic core of the model.
Whereas former versions of PALM were already validated against wind-tunnel measurements, real-world measurements, 25 and other CFD codes (Letzel et al., 2008;Razak et al., 2013;Park et al., 2015b;Gronemeier and Sühring, 2019;Paas et al., 2020), the drastic changes of PALM's code base require a new validation from scratch. A sufficient validation is inevitable to ensure confidence in the results of the PALM model system as it is also the case for every other CFD code (Blocken, 2015;Oberkampf et al., 2004).
Because of the high complexity of PALM, validating the model is a very longsome and costly exercise and a complete 30 validation of all model components would easily go beyond the scope of any article. Within this study, we therefore focus on the validation of the model's flow dynamics which make up the core of the model system and build the foundation for all other features within PALM. In order to isolate the pure dynamics from all other code parts, PALM is operated in a pure dynamic mode, i.e. all thermal effects (temperature and humidity distribution, radiation, surface albedo, heat capacity, etc) are switched off. The simulation results can then be validated using wind-tunnel measurements (Leitl and Schatzmann, 2010) which are 35 recorded in a similar setup as the simulation data. While it is virtually impossible to neglect temperature or humidity effects in real-world measurements, wind-tunnel experiments can provide exactly the same idealized conditions as used in our idealized simulation. Also, other difficulties such as additional non-resolved obstacles like trees or sub-grid features on building walls existing in the real world can make a comparison with real-world measurements troublesome (Paas et al., 2020). Paas et al. (2020) compared PALM simulations to measurements of a mobile measurement platform. Although overall good agreement 40 was found between PALM and the measurements, some non-resolved obstacles like trees complicated the comparison at several points and led to differences in results. Hence, we decided to validate PALM against an idealized wind-tunnel experiment for this study.
A realistic building setup, in this case the HafenCity area of Hamburg, Germany, is chosen for validation. A real-case building setup has the advantage over more idealized, i.e. a single-cube, cases that a variety of different building configurations, also 45 including more or less solitary buildings, can be covered in a single simulation. Likewise, it may show the capability of PALM to correctly reproduce a complex realistic wind distribution.
To further improve the quality of the validation, a blind test was conducted, i.e. both experiments (CFD and wind-tunnel) were made independently from one another. Only information about the building layout and the inlet wind profile were shared between both experiments. Such a blind test gives high value to the validation as no tweaking of either experiment is possible 50 to tune the results towards each other. Although, there are methods to adjust CFD results to better match to measurements (e.g., Blocken et al., 2007), these adjustments depend on the individual case and need to be re-calculated for each new studied situation. Within this study, we chose not to utilize any adjustments to the results in order to show PALM's capability to simulate a certain situation without any customized correction of the results. Also, often, measurements might even be unavailable to calculate such adjustment factors.  Figure 1 shows a photograph 60 from within the wind tunnel for reference. For each wind tunnel campaign, a neutrally stratified model boundary layer flow is generated by a carefully optimized combination of turbulence generators at the inlet of the test section, and a compatible floor roughness. For the present study a boundary layer flow was modelled to match full scale conditions for a typical urban boundary layer measured at a 280 m tall tower in Billwerder, Hamburg. The mean wind profile can be described by a logarithmic wind profile with a roughness length z 0 = (0.66 ± 0.22) m and by a power law with a profile exponent α = 0.21 ± 0.02. The 65 approaching flow profile is depicted in Fig. 2 and was modelled at a wind direction of 110 • .
A 2D Laser-Doppler-Anemometry (LDA) System was used to measure component-resolved flow data at sampling rates of 200 Hz-800 Hz (model scale), resolving even small-scale turbulence in time. At each measurement location a 3 minutes time series was recorded which corresponds to a length of about 25 hours at full scale. The reference wind speed was permanently monitored close to the tunnel inlet through Prandtl tube measurements. An area of 2.6 km 2 covering the HafenCity of Hamburg, 70 Germany, was modelled at a scale of 1/500 within the wind tunnel (see Fig. 1). Measurements were taken at 25 different locations within the building setup as shown in Fig. 3.

PALM simulation
The PALM Model System 6.0, revision 3921, was used to conduct the simulation for this study. PALM was operated using a fifth-order advection scheme after Wicker and Skamarock (2002) in combination with a third-order Runge-Kutta timestepping 75 scheme after (Williamson, 1980). At this point, we skip a detailed description of the PALM model. A detailed description is provided by Maronga et al. (2015Maronga et al. ( , 2020. For the PALM simulation, real-world scales were used. The simulation domain covered an area of 6000 m by 2880 m horizontally and 601 m vertically at a spacial resolution of ∆x = ∆y = ∆z = 1 m in each direction. This resulted in about 10.4 × 10 9 grid points for the used staggered Arakawa C grid (Harlow and Welch, 1965;Arakawa and Lamb, 1977). The area of interest, i.e. the HafenCity area, was situated downstream of the simulation domain.

80
The model domain was oriented so that the mean flow direction was aligned with x direction. With a mean wind direction of 110 • , the model was hence rotated counter-clockwise by 200 • .
The building layout as used in PALM is depicted in Fig. 4. In PALM, topography is considered using the mask method (Briscolini and Santangelo, 1989) where a grid volume is either 100 % fluid or 100 % obstacle. In combination with PALM's rectilinear grid, this can cause buildings not aligned with the grid to appear differently, more brick-like, than they were within 85 the wind tunnel.
The basic setup for this study is based on the settings used in the former study of Letzel et al. (2012).
A heterogeneous building setup usually requires a non-cyclic boundary condition along the mean flow direction to ensure that building-generated turbulence is not cycled over and over the analysis area which otherwise might influence the results.
However, tests with non-cyclic boundary conditions along the mean flow showed that extremely long simulation times would be 90 required to generate a stationary state. Hence, we used cyclic boundary conditions instead, which reduced the required CPUtime significantly. The domain was extended in mean flow direction to allow the building-generated turbulence to dissipate before the flow hits the target area again due to the cyclic conditions. As the simulation was aimed at a pure neutral case and without releasing any trace gases or alike within the city area, there was no disadvantage in using cyclic boundary conditions instead. After a simulation time of 7 h steady-state conditions were reached.

95
At the top boundary, a fixed wind speed of 1.54 m s −1 was considered according to measurements within the wind tunnel at the corresponding height. A constant flux layer was assumed between the surface and the first computational grid level to calculate the surface shear stress, using a roughness length of z 0 = 1 m. Due to the staggered grid, the first computational level is positioned 0.5∆z above the surface.
Within the wind-tunnel experiment, the modelled roughness length was z 0 = (0.66 ± 0.22) m within the upstream area (cf. Sect. 2.1). With the chosen resolution of 1 m this roughness length cannot be used in the PALM simulation, because it would be larger than the assumed height of the constant flux layer (which is 0.5 m) and hence would cause numerical problems. We used a roughness length significantly smaller than half the grid spacing instead (i.e. z 0 0.5 m), with a value of z 0 = 0.1 m at each surface. Other comparable city simulations prove to give reasonable results using such value (e.g., Park et al., 2015a;Gronemeier et al., 2017;Paas et al., 2020).

105
To match the conditions within the wind tunnel, a strictly neutral atmosphere was considered with temperature being constant over time. Also, the Coriolis force was neglected.
In the past, streak-like structures were reported for LES of neutral flows using cyclic boundary conditions (Munters et al., 2016). To avoid such numerical artifacts, a shifting method was used according to Munters et al. (2016). The flow was shifted by 300 m in y direction, i.e. perpendicular to the mean wind direction, before entering the domain at the left boundary.

110
The wind field was initialized using a turbulent wind field from a precursor simulation via the cyclic-fill method (Maronga et al., 2015). The setup of the precursor simulation was similar to the main simulation but with a reduced domain size of 600 m by 600 m in horizontal direction. To initialize the precursor simulation, the approaching wind profile measured in the wind tunnel was used as shown in Fig. 2.
The total simulation time of the main simulation was 12 h from which the first 7 h were required to reach a steady state of 115 the simulation. The final 5 h were used for the analysis. Figure 2 shows the mean wind profile of the flow approaching the building area during the analysis time. Differences between the approaching flow within the wind tunnel and PALM are due to differences in the roughness of the upstream area. While roughness elements were present within the wind-tunnel experiment, there were no such elements present within the PALM simulation and surface roughness was considered only via the roughness length used within the constant-flux layer assumption 120 at the lower-most grid level. Due to the mismatch in z 0 as discussed above, the smoother surface within PALM causes higher wind speed close to the ground and lower speed at greater heights. Antoniou et al. (2017) also reported a similar behaviour in their comparison of LES and wind-tunnel results. Future simulations should hence try to include roughness elements similar to those of the wind-tunnel to overcome this discrepancy.

Measurement stations 125
Within the wind-tunnel experiments, wind speed was measured at certain measurement stations within the building array. The locations of which are shown in Fig. 3. To be able to mirror the measurements as best as possible, the virtual measurement module of PALM was used (Maronga et al., 2020). This module allows to define several virtual measurement stations within the model domain via geographical coordinates. The model domain itself then needs to be geo-referenced in order to identify the grid points closest to the measurement location. Referencing is done by assigning geographical coordinates and orientation 130 to the lower left corner of the domain.
When mapping the measurement stations onto the PALM grid, there were two difficulties: First, there was not always a grid point at the exact same location as the measurement points and therefore no data was available at the exact same position; second, the topography in close vicinity of a measurement point might have been slightly different due to topography representation used in PALM (cf. Sect. 2.2). To overcome these two issues, virtual measurements not only from the closest grid 135 point to a measurement position were saved but also values from the neighbouring grid points. In post-processing, the area of each measurement station was analyzed and a grid point selected which best fitted the wind-tunnel measurements.
At each measurement station presented in this study, nine profiles were recorded with a sampling rate between 2 Hz-2.5 Hz (measurements recorded during each time step).

PALM simulation
The PALM simulation required a spin-up time of 7 h as can be seen by the time series of the domain-averaged kinetic energy E = 0.5 √ u 2 + v 2 + w 2 and the friction velocity u * (see Fig. 5). Both quantities stabilized after 7 h at around 0.8856 m 2 s −2 (E) and 0.0594 ms −1 (u * ). Therefore, only data from the last 5 h of the simulation were used for the following validation.
The horizontally and time-averaged vertical profile of the stream-wise component of the vertical momentum flux wu is 145 shown in Fig. 6. The vertical momentum flux wu can be split into a resolved and a subgrid-scale (SGS) part which is parameterized via an SGS model. The higher the resolved part, the less the SGS model contributes to the flux therefore indicating that the flux and hence the turbulence causing the flux is well resolved. The ratio of the total, i.e. resolved plus SGS part, momentum flux and its resolved part is close to 1 revealing that turbulence is properly resolved within the simulation domain (see Fig. 6). At the surface, turbulence is less resolved due to the fact that turbulent structures tend to become smaller the closer 150 they get to the surface and cannot be resolved by the grid spacing any more. However, the ratio between resolved and total wu stays above 0.9 except for the lowest two grid levels.
To get an impression of the turbulent structures, Fig. 7 shows a snapshot of the magnitude of the rotation of the velocity field as a measure of turbulence. One can clearly identify strong turbulent features (yellow and red structures) within the vicinity of buildings while only weak turbulence is present above smooth surfaces.  The general lower wind speed recorded within PALM can be partly explained by the difference in measurement height.
Within PALM, measurements were located 0.5 m lower than in the wind-tunnel experiment due to the staggered Arakawa C At measurement heights 50 m and 75 m (not shown), PALM reports only 2 % and less than 1 % lower wind speed, respectively, than the wind tunnel which supports the idea that a mismatch in z 0 might be responsible for the differences.
The large discrepancy between PALM and the wind-tunnel results measured at station 13 can be explained by differences within the building setup. Although the building data was checked thoroughly before the simulation, the windward building at station 13 has different heights within the wind-tunnel experiment and the PALM simulation. While in the wind tunnel the 190 building height is about 6 m, in PALM, it is set at 24 m (cf. Figs. 1 and 4). This leads to significant differences in the flow around the building and hence causes different results at station 13. Wind direction and wind speed match again between PALM and wind-tunnel measurements at heights above 30 m, i.e. above the building height (not shown).
The building layout was again thoroughly checked after discovering this flaw but no further deviations were detected within the vicinity of any of the other measurement stations. A new simulation with an updated building height would yield different 195 results around the changed building including station 13. However, given the fact that the local building configuration was found to dominate the wind flow at the other measurement locations, it is very unlikely that the results would change in a revised simulation. Therefore, repeating the simulation with an updated building layout around station 13 would not change the overall outcome of the study.
In the following, we limit the discussion to three stations: 4, 11, and 7 which are chosen to represent a good, an average  Measurements present lower wind speed within the building canopy and higher wind speed above the canopy in PALM compared to wind-tunnel measurements. Within the canopy, the flow is more turbulent in PALM as I u and I v are higher than within the wind tunnel. This in turn leads to a reduced mean wind speed. The wind direction also differs slightly within the The rooftop of the surrounding buildings is at 26 m and 36 m height (cf. Fig. 4). Differences might therefore be related to rooftop vortices being of slightly other shape within PALM compared to the wind tunnel. Results from stations 5, 6, 10, 15, and 17 are comparable to those of station 11.
Measurements at station 4 (Fig. 12) show a very close agreement between wind-tunnel and PALM results. The v profile matches very well to that of the wind-tunnel measurement while the u values tend to be slightly lower than the wind-tunnel measurement. This results in a good agreement of wind direction and horizontal wind speed as well. PALM does not only produce similar mean values as the wind tunnel experiment at station 4 but also the variance of wind speed and turbulence intensity 225 matches well between wind-tunnel and PALM results. A similar agreement between PALM and wind-tunnel measurements can be observed at the majority of stations, i.e. at stations 1-3, 8, 9, 12, 14, 16, 18, 19, 22-25. A further insight of the simulated turbulence distribution is derived by analyzing the spectral energy density S. Figure 13 shows the dimensionless energy spectra for station 4 (left panel) and 11 (right panel) at different heights z. Spectra measured at station 7 are comparable to those of station 11 and are therefore not shown. Note that the covered range of frequencies 230 differ between PALM and wind-tunnel measurements as the sampling rate of the measurements and the measured time interval vary between the PALM and wind-tunnel experiment. However, results of both experiments still overlap over a large range of frequencies.
In general, spectra for both wind-vector components coincide to a high degree between PALM and wind-tunnel measurements at all heights. The inertial range of the turbulence spectrum is clearly visible within both experiments at 75 m height 235 (above the canopy layer) at both stations ( Figs. 13a and b). The normalized energy spectrum decays with roughly f − 2 3 following Kolmogorov's theory. At high frequencies, spectra of the PALM measurements strongly decay which is related to numerical dissipation and is a typical behaviour of LES models using high-order differencing schemes (e.g., Glendening and Haack, 2001;Kitamura and Nishizawa, 2019).
At rooftop height (Figs. 13c and d), PALM's spectra are shifted towards higher frequencies compared to those of the wind 240 tunnel at the same height. This indicates that PALM simulates turbulence of smaller scales at these points. The effect is more pronounced for station 4 than for station 11 and station 7 (not shown). The smaller turbulence within PALM might be due to a higher roughness length z 0 of the rooftops which in turn might either be due to the specified z 0 being of too high value or due to the increased roughness of the buildings introduced by the brick-like representation of the buildings in PALM (cf. Sect. 2.2). Around station 7 and 11, building roofs are not flat but possess small additional structures which might give a higher 245 roughness as a result of the grid-following topography structure in PALM compared to the building model within the windtunnel experiment. This discrepancy then also leads to the difference already mentioned in the profiles of station 11 (Fig. 11).
At station 4, this effect is not as pronounced because the closest building to station 4 has a flat roof.
Spectra close to the surface again agree better between PALM and the wind tunnel measurements. Though, due to the limited measurement frequency and very small turbulent structures at the surface, the inertial range is not covered by the wind-250 tunnel measurements at 3 m height. This height corresponds to the third grid level above the surface of the PALM simulation and therefore it can be truly expected that the inertial range is only poorly resolved. Comparing the measured spectra to the theoretical decay of f S ∝ f − 2 3 , the inertial range is indeed hardly represented within the data.
Results within the vicinity of the measurement stations could be improved by utilizing PALM's self-nesting feature. This allows to use a higher grid resolution within specific areas of the model domain. We recommend that future simulations should 255 try using this feature for areas requiring high resolution.

Conclusions
In this study, we analyzed PALM's capability to simulate a complex flow field within a realistic urban building array. Simulation results were compared with measurements done at the EWTL facility at University of Hamburg, Germany. The aim was to validate the dynamic core of the newest version 6.0 of PALM which underwent significant code-changes in recent model 260 development.
The comparison of PALM results with the wind tunnel data proves that the model can correctly simulate a neutrally stratified urban boundary layer produced by a realistic complex building array. Measurements were compared at several different positions throughout the building array at non-obstructed locations, at the windward and leeward site of buildings as well as within street canyons and at intersections. Overall, the PALM results highly agreed with the corresponding wind-tunnel measurements 265 in regards to wind speed and direction as well as turbulence intensity. Differences could be observed close to surfaces where PALM reported slightly lower wind speed compared to wind-tunnel measurements. Such discrepancies were also reported recently by Paas et al. (2020) when comparing PALM simulations to real-world measurements. These differences can to large extent be ascribed to deficiencies in the model setup, namely the choice of a relatively large and homogeneous roughness length in the whole model domain, which creates excessive turbulence and a reduced mean wind in the vicinity of buildings.

270
Our choice of z 0 = 0.1 m is also not consistent with the recommendations given by Basu and Lacser (2017) who recommend that the first grid level should at least be at 50z 0 . Given that the first computational grid level is set at 0.5 m, z 0 should have been no larger than 0.01 m which could have also better resembled the smooth surface of the buildings within the wind-tunnel experiment. An improved simulation setup would have to set a significantly smaller z 0 within the building array to improve the wind field close to the surface. To also resemble a rough surface within the upstream area, roughness elements similar to those 275 from the wind-tunnel experiment should then be simulated as well. Moreover, changes in the measurement heights between wind tunnel and LES led to systematic differences, particularly near the surface where the vertical wind gradients are large.
When planning a simulation, it must be taken into consideration, that surfaces are always aligned with the rectilinear grid within PALM. If the alignment of walls and the simulation grid do not match, the resulting brick-like representation can lead to a higher surface roughness and hence to higher turbulence presumably altering the results and leading to differences to other 280 measurements. To overcome this error or at least reduce its effect, a higher grid resolution should be used which will then minimize the size of artificially created steps for topography elements non-aligned to the numerical grid. In order to achieve this without increasing the computational costs significantly, PALM's ability of grid nesting should be used.
In the present study, we also experienced that input data must always be checked with very high caution. Especially large building data sets might contain errors and false building heights or missing/displaced buildings, which are more difficult 285 to spot than in setups with a limited number of buildings. This is, of course, of utmost importance for the area-of-interest, but also the upwind region requires proper verification as it directly affects the analysis area. Also, when comparing to other experiments like real-world or wind-tunnel measurements, positioning of the measurements must be thoroughly checked as mentioned by Paas et al. (2020). This is also true for positioning virtual measurements within the PALM domain. At positions with highly complex wind fields, it can make a large difference for the results if measurement positions are off by only a single 290 grid point (as can be seen, e.g., by the range of profiles for station 11, Fig. 11). This of course depends strongly on the grid spacing used and will be most relevant when using relatively coarse grids.
This study focused only on a single but the most essential part of PALM, the dynamic core. However, validation of the entire model requires additional studies focusing on the other model parts like the radiation module, the chemistry module or the land-surface module to mention only a few. Some of these are already validated (Resler et al., 2017;Kurppa et al., 2019;295 Fröhlich and Matzarakis, 2019), others will follow in future publications.
Code and data availability. The    Figure 11. Profiles of mean horizontal wind speed Uhor, wind direction, wind components u and v, and turbulence intensity Iu and Iv along x and y, respectively, at measurement station 11. Error bars denote the standard deviation of the respective quantity. Note that z = 0 denotes street-level height.