the Creative Commons Attribution 4.0 License.
Special issue: The PALM model system 6.0 for atmospheric and oceanic boundarylayer...
Model description paper 03 Sep 2021
Model description paper  03 Sep 2021
Mesoscale nesting interface of the PALM model system 6.0
 ^{1}Deutscher Wetterdienst, Offenbach, Germany
 ^{2}Institute of Meteorology and Climatology, Leibniz University Hannover, Hanover, Germany
 ^{1}Deutscher Wetterdienst, Offenbach, Germany
 ^{2}Institute of Meteorology and Climatology, Leibniz University Hannover, Hanover, Germany
Correspondence: Eckhard Kadasch (eckhard.kadasch@gmail.com)
Hide author detailsCorrespondence: Eckhard Kadasch (eckhard.kadasch@gmail.com)
In this paper, we present a newly developed mesoscale nesting interface for the PALM model system 6.0, which enables PALM to simulate the atmospheric boundary layer under spatially heterogeneous and nonstationary synoptic conditions. The implemented nesting interface, which is currently tailored to the mesoscale model COSMO, consists of two major parts: (i) the preprocessor INIFOR (initialization and forcing), which provides initial and timedependent boundary conditions from mesoscale model output, and (ii) PALM's internal routines for reading the provided forcing data and superimposing synthetic turbulence to accelerate the transition to a fully developed turbulent atmospheric boundary layer.
We describe in detail the conversion between the sets of prognostic variables, transformations between model coordinate systems, as well as data interpolation onto PALM's grid, which are carried out by INIFOR. Furthermore, we describe PALM's internal usage of the provided forcing data, which, besides the temporal interpolation of boundary conditions and removal of any residual divergence, includes the generation of stabilitydependent synthetic turbulence at the inflow boundaries in order to accelerate the transition from the turbulencefree mesoscale solution to a resolved turbulent flow. We demonstrate and evaluate the nesting interface by means of a semiidealized benchmark case. We carried out a largeeddy simulation (LES) of an evolving convective boundary layer on a clearsky spring day. Besides verifying that changes in the inflow conditions enter into and successively propagate through the PALM domain, we focus our analysis on the effectiveness of the synthetic turbulence generation. By analysing various turbulence statistics, we show that the inflow in the present case is fully adjusted after having propagated for about two to three eddyturnover times downstream, which corresponds well to other stateoftheart methods for turbulence generation. Furthermore, we observe that numerical artefacts in the form of gridscale convective structures in the mesoscale model enter the PALM domain, biasing the location of the turbulent up and downdrafts in the LES.
With these findings presented, we aim to verify the mesoscale nesting approach implemented in PALM, point out specific shortcomings, and build a baseline for future improvements and developments.
The simulation of urban flows under realistic conditions poses a multiscale problem where evolving synoptic scales interact with building and streetsize scales. While the continuing growth of available computational resources enables largeeddy simulation (LES) to be applied to more and more realistic scenarios at regional scales (Schalkwijk et al., 2015), it still remains unfeasible to simulate mesoscale processes at resolutions fine enough to represent smallscale turbulence generated by urban structures. To overcome this hurdle and consider synoptically evolving conditions in highresolution LES models, various concepts with different degrees of idealization have been developed to couple LES models to largerscale models. To date, modellers either employ cyclic boundary conditions and add largescale forcing and nudging terms to the prognostic equations (e.g. Heinze et al., 2017), or they may employ gridnesting approaches (e.g. Mirocha et al., 2014) with timedependent in and outflow boundary conditions.
Both approaches face particular challenges, mainly linked to the representation of the turbulent flow at the domain boundaries, requiring large buffer zones to move boundary effects away from the region of interest. In the first approach, periodic domain boundaries allow unrealistic flow feedbacks due to reentering flow structures caused by complex terrain, urban surfaces, or other surface heterogeneity. Furthermore, feedbacks are not limited to the velocity field. When anthropogenic heat or chemical compounds are emitted, unrealistic thermodynamic conditions or concentrations would reenter the model domain on the opposite boundary modifying the upstream conditions for the urban environment, which in turn may bias the distribution of heat and mass concentrations. Here, buffer zones help to move the affected flow region outwards (Letzel et al., 2012; Maronga and Raasch, 2013). Schalkwijk et al. (2015) used a hybrid nesting approach to minimize scalar and mean flow feedbacks from reentering wakes. They used cyclic boundary conditions in order to retain turbulent fluctuation across the domain but relax horizontal velocities, temperature, and specific humidity towards the parent mesoscale model in a region close to the LES domain boundaries. However, since the relaxation only shifts the mean state towards the parent model, turbulent wakes generated by local surface heterogeneities like orography, buildings, etc. may reenter on opposite boundaries nevertheless.
Alternatively, gridnesting approaches can be employed, which realize a oneway coupling via timedependent inflow and outflow boundary conditions derived from a largerscale parent model. In mesoscale models with horizontal grid spacings on the order of 𝒪(1 km), however, the turbulent exchange of momentum, heat, and water vapour is parameterized so that the their prognostic fields and derived LES boundary conditions lack turbulent fluctuations. For proper representation of the turbulent flow in the atmospheric boundary layer within the domain of interest, the incoming flow field should be fully spatially developed; i.e. it should not depend on the distance to the inflow boundary layer any more (Lee et al., 2019). This requires buffer zones at the inflow boundaries where turbulence can spatially develop. Mirocha et al. (2014) showed that without adding any perturbations it takes a fetch length of several tens of kilometres to obtain fully spatially developed turbulence, meaning that significant parts of the computational resources are only spend on the buffer zones.
To reduce the required fetch length, various approaches to generate turbulent inflow conditions exist; for a comprehensive overview about existing methods, we refer to Wu (2017). For simulations of atmospheric boundarylayer flows, turbulence recycling approaches are often used (e.g. Park et al., 2015; Munters et al., 2016; Gronemeier et al., 2017). For simulations with one defined inflow boundary, PALM offers a turbulence recycling method according to Kataoka and Mizuno (2002), where a turbulent signal is read from a defined plane of the model grid and imposed onto the stationary mean profiles at the inflow boundary. To apply this approach, the flow conditions within the recycling plane should be statistically homogeneous, in order to avoid feedbacks between the turbulent signal at the recycling plane and the inflow boundary (Munters et al., 2016). In simulations with realistic land surface distributions, complex terrain or buildings present, however, statistically homogeneous turbulence at the recycling plane cannot be guaranteed without adding large buffer zones. Moreover, due to the evolving boundary conditions accompanied with changing inflow boundaries, the location of the recycling plane may change and it is not clear what happens, e.g. for diagonal inflows, making the turbulence recycling difficult to apply for evolving inflow conditions.
In contrast to recycling methods, volumeforcing approaches do not necessarily require homogeneous inflow conditions. To accelerate the spatial development of a turbulent flow, MuñozEsparza et al. (2014) developed and implemented the socalled cellperturbation method into the Weather Research and Forecasting Model (WRF, Skamarock et al., 2008), where boxlike perturbations are added onto the potential temperature within a defined region close to the inflow boundary. This was further developed by MuñozEsparza et al. (2015) and MuñozEsparza and Kosović (2018) by scaling the thermal perturbation amplitude depending on the stability regime. With this approach, MuñozEsparza and Kosović (2018) could significantly reduce the required fetch length to obtain fully adjusted turbulence, even under neutral and stable stratifications. The cellperturbation method has shown promising results when applied in nested WRFLES simulations of a full diurnal cycle for a realworld setup (MuñozEsparza et al., 2017), as well as in simulations for ocean–island interactions (Jähn et al., 2016) and coastal sea breeze events (Jiang et al., 2017). Furthermore, prescribing WRF output data as boundary conditions in a PALM simulation, Lee et al. (2019) demonstrated the ability of the cellperturbation method in a densely builtup urban environment, where the required buffer zones could be significantly reduced compared to a nonperturbed simulation. To test a more direct method of turbulence generation, Mazzaro et al. (2019) extended the original cellperturbation method by adding scaled perturbations onto the velocity components at the near inflow region. This approach showed a comparable performance compared to the original version with a faster spatial development close to the inflow boundary but a longer adjustment fetch required to achieve an equilibrium state.
Alternatives to volumeforcing approaches are socalled filtering approaches, where spatially and temporally correlated perturbations are imposed only onto the velocity components at the lateral boundaries (e.g. Klein et al., 2003; Xie and Castro, 2008). Gronemeier et al. (2015) have originally implemented the synthetic turbulence generation method by Xie and Castro (2008) into PALM and found good agreement of the turbulent flow development in an urban environment compared to using the turbulence recycling method according to Kataoka and Mizuno (2002). The main challenge of this approach, however, is to adequately infer Reynoldsstress components, as well as turbulent length scales and timescales of the flow to generate appropriate inflow turbulence.
Beside the necessity to add perturbations at the boundaries, modellers should also be aware that numerical artefacts from the mesoscale model may propagate into the LES; e.g. Mazzaro et al. (2017) and MuñozEsparza et al. (2017) showed that underresolved flow structures propagate from a mesoscale WRF simulation into the LES. Honnert et al. (2011) found that in convectionpermitting mesoscale simulations, resolvedscale convection on the grid scale can develop when the boundarylayer depth approaches the horizontal grid resolution. When boundarylayer convection is explicitly resolved in mesoscale models, this is often referred to as the grey zone, or terra incognita (Wyngaard, 2004), where both mesoscale and LES model assumptions break; i.e. the grid spacing is already too small compared to the dominant length scales of the flow to justify usage of fully parameterized boundarylayer schemes but still too large to reliably resolve convective structures. Ching et al. (2014) and Zhou et al. (2014) showed that the strength and spatial scale of the resolvedscale convection strongly depends on the horizontal grid resolution, while Shin and Dudhia (2016) also confirmed a dependence on the applied boundarylayer scheme. With a grid nesting of a turbulenceresolving WRF simulation into a mesoscale WRF simulation, Mazzaro et al. (2017) showed that such underresolved flow structures propagate into the LES, delaying the spatial development of turbulence. For a strongly convective case, Mazzaro et al. (2017) further showed that firstorder statistics in the LES are not significantly affected by imposed underresolved convection from the parent mesoscale simulation when the flow has been fully adjusted, though variances, turbulent vertical fluxes, and length scales in the LES tend to become larger for stronger underresolved mesoscale convection. Further, they showed that the signals of the imposed up and downdrafts from the mesoscale model vanish after about 40 km downstream of the inflow boundary, even though they also noted that under less convective conditions the signals may even persist longer. However, this implies that in the event of underresolved convection in the mesoscale model, the turbulent transport in the LES as well as the location where up and downdrafts occur, depend on the mesoscale model setup, i.e. horizontal resolution, boundarylayer parameterization, etc. In our test scenario we also found underresolved rolllike convective structures that propagate into the LES domain.
Another issue that emerges when nesting LES in mesoscale models concerns the representation of the atmospheric boundary layer. Due to different treatment of turbulent transport, i.e. boundarylayer parameterizations in the mesoscale model vs. an explicit representation of turbulent eddies in the LES model, the vertical transport of energy, water, and momentum may differ considerably. In situations where this is the case, the mean state of the LES solution, which is generally more credible due to a wider range of explicitly resolved turbulent scales, will be pushed towards the mesoscale solution including any possible model biases.
In this paper, we present the mesoscale nesting interface of the PALM 6.0 model system. It provides timedependent spatially heterogeneous boundary conditions for PALM obtained from the mesoscale model COSMO (see for instance Baldauf et al., 2011) and includes a synthetic turbulence generation method to accelerate generation of turbulent fluctuations at the model boundaries. COSMO has been developed by the Consortium for Smallscale Modeling (http://cosmomodel.org, last access: 13 August 2021) and currently serves as the operational regional weather forecasting model at the German Meteorological Service (Deutscher Wetterdienst, DWD). PALM's mesoscale nesting interface consists of two major parts: (i) the preprocessor INIFOR, which derives initial and boundary conditions from mesoscale model output, and (ii) PALM's internal routines for reading these forcing data and superimposing synthetic turbulence. In particular, we impose synthetic turbulent structures at the lateral boundaries following Xie and Castro (2008), while the required turbulence statistics are parameterized based on mesoscale model output. At the moment, INIFOR is tailored towards the COSMO model, but extensions to WRF and ICON (Zängl et al., 2015; Reinert et al., 2020) are planned in the future. Note that the scope of this paper is on the description of our nesting approach and on the demonstration of its effectiveness. We defer indepth analyses and comparisons to other methods to further publications.
This approach provides a oneway nesting capability of PALM into a mesoscale simulation, where boundary conditions are only set for child model. At this point, we want to distinguish the mesoscale nesting interface from PALM's self nesting capabilities (Hellsten et al., 2021). While selfnesting allows a twoway coupling of a PALM child domain within a parent PALM domain, the mesoscale nesting interface presented in this paper realizes a oneway or offline nesting of PALM within a mesoscale model. That means, while PALM obtains timedependent boundary conditions from the mesoscale model, information produced by PALM is not fed back into the mesoscale model. Both nesting features may, however, be combined to carry out LES nested in COSMO with one or multiple twoway coupled child nests within PALM.
This paper is structured as follows. We describe the mesoscale nesting approach in Sect. 2, including the relevant mesoscale–microscale model differences, the resulting transformation and interpolation methodology implemented in INIFOR, as well as the synthetic inflow turbulence generation method with its underlying turbulence parameterizations. We demonstrate and verify our mesoscale nesting approach in a semiidealized benchmark simulation of a convective boundary layer under evolving synoptic conditions. We describe the setup in Sect. 3 and analyse the case thereafter in Sect. 4. We conclude this paper with a summary of our findings and an outlook on future developments in Sect. 5.
The PALM model is nested into the mesoscale model by prescribing initial conditions and timedependent Dirichlet boundary conditions derived from output of the parent mesoscale model. Boundary conditions for PALM are given for the top and lateral domain boundaries. The boundary conditions at the surface are provided by PALM's urban and landsurface model, which are initialized from the mesoscale data.
The boundary conditions enter PALM via the mesoscale nesting interface, which consists of two major components: (i) the preprocessor INIFOR and (ii) PALM's internal boundary condition routines. The workflow of a model run using the mesoscale nesting interface is illustrated in Fig. 1. First, the forcing data are interpolated in a preprocessing step using INIFOR and stored in a netCDF driver file. In analogy to the static driver (Maronga et al., 2020), which contains all static geospatial information such as topography, building and surface parameters, etc., we refer to this forcing file as the “dynamic driver”. During the simulation, PALM successively reads and processes the dynamic driver data. This includes temporal interpolation of the boundary data, removal of any residual divergence, as well as the superposition of synthetic turbulent fluctuations (see Sect. 2.3).
The required prognostic variables for which the dynamic driver provides initial and boundary conditions are listed in Table 1 next to their equivalents in the COSMO model output. Note that we use uppercase letters to denote COSMO's dependent variables and lowercase ones for PALM. In particular, INIFOR provides data for the state of the atmosphere ($u,v,w,\mathit{\theta}$, and q_{v}) at model start, which can be supplied either as onedimensional vertical profiles (level of detail, LOD = 1) or as threedimensional fields (LOD = 2). Since the initial atmospheric data are already interpolated onto the PALM's Cartesian grid by INIFOR (see Sect. 2.2), it can be directly copied onto the respective internal PALM arrays after it is read from the dynamic driver. Further, the dynamic driver contains the initial state of soil moisture (m_{soil}) and temperature (t_{soil}), again either as onedimensional profiles (LOD = 1) or as horizontally heterogeneous threedimensional data (LOD = 2). To allow for a different number of soil layers in the PALM domain depending on the local soil type, we decided to linearly interpolate the provided soil data during soilmodel initialization rather than in a preprocessing step done by INIFOR as it is done for the initial atmospheric quantities. At this point, we note that the provided initial soil data only contain values aggregated over a mesoscale grid cell, which in reality may feature surfaces with various soil types and different land use across which soil moisture and temperature can vary significantly.
Hence, we recommend to run the soilmodel spinup mechanism as described in Maronga et al. (2020) to obtain individual soil moisture and temperature profiles that are in equilibrium with the local conditions at each model surface. In the case of self nesting, where fineresolution child domains are nested within a coarserresolution outer parent domain, it is sufficient to provide initial mesoscale model data for the outermost parent domain only, while the respective initial data are propagated to the nested child domains as described in Hellsten et al. (2021). However, the user may also provide a separate dynamic driver for PALM to initialize atmosphere and soil quantities in the respective child domain, which is useful, for instance, if highresolution initial soil data for a limited area are available.
In addition to the initial state, the dynamic driver provides the timedependent boundary conditions for PALM's atmospheric prognostic variables $(u,v,w,\mathit{\theta},$ and q_{v}) at the top and the four lateral boundaries at certain points in time (hourly data are provided from COSMO output). These timedependent boundary data are read from the dynamic driver and are linearly interpolated onto the respective model time level, while the data are copied onto the respective model boundaries. In order to save memory, only the boundary data at LES time levels t_{i} and t_{i+1} are read, with ${t}_{i}\le {t}_{\mathrm{s}}<{t}_{i+\mathrm{1}}$, while t_{s} being the actual simulation time in the model. The boundary data can be provided either as onedimensional vertical profiles (one value for the top boundary) that are identical at each of the lateral boundaries (LOD = 1) or as individual twodimensional x–z (north and south lateral boundaries), y–z (east and west boundaries), and x–y (top boundary) crosssections.
The velocity boundary conditions and the associated massflux fields obtained from a compressible mesoscale model such as COSMO do not generally satisfy the divergencefree condition of incompressible models such as PALM. To overcome this, we correct the velocity ${w}_{\mathrm{b}}^{\mathrm{t}}$ at the top boundary similar to the approach described by Hellsten et al. (2021). The correction is calculated from the integrated mass flux through the lateral and top boundaries as
where u_{b} denotes the velocity vector at the boundary, n the boundary normal vector and Ω the surface area of the model boundaries. We obtain divergencefree boundary conditions by using the corrected vertical velocity:
instead of the preliminary boundary condition ${w}_{\mathrm{b}}^{\prime \mathrm{t}}(x,y)$ at the top boundary. With this correction, we satisfy the massflux continuity globally. Local continuity is continuously maintained by the pressure correction step (cf. Appendix A).
We do not currently use any damping zones near the lateral boundaries to relax the solution towards the boundary conditions as is done, for instance, in the WRF model. There, a relaxation term according to Davies and Turner (1977) is added in the momentum equations near the boundaries to suppress reflections of waves. For the present study, we tested the effect of such damping zones on the flow (not shown) but we could not observe any significant differences in the results nor any wave reflection. However, this may change in the future when PALM gains support for a compressible system of equations, which would add sound waves to the solution.
2.1 Model differences
In the following, we describe the relevant model properties and point out the relevant differences, which yield the conceptual steps that need to be carried out by INIFOR. Here, we omit indepth descriptions of both models and refer the reader to additional publications. In particular, more information about the formulation and numerics of COSMO can be found in the model documentation by Doms and Baldauf (2018). Details and verification studies of COSMODE – a particular model configuration used at DWD – have been published by Baldauf et al. (2011). For details about the PALM model system, please see the descriptions by Maronga et al. (2015, 2020) and the publications cited therein.
PALM and COSMO differ in a number of ways, between which INIFOR needs to translate in order to derive PALM forcing data. The first difference lies in the physical formulation of the models. COSMO is a nonhydrostatic limitedarea atmospheric model. It is based on fully compressible equations, which are formulated in terms of the three spherical wind components, absolute temperature and pressure, density, and multiple water phases. PALM, on the other hand, solves incompressible equations for the moist atmosphere, where either the Boussinesq or an anelastic limit of the Navier–Stokes equations may be used. The model is formulated in terms of the three Cartesian wind components, the potential temperature, and the water vapour mixing ratio. The continuity equation in the anelastic and Boussinesq approximations reduces to divergence constraint $\mathrm{\nabla}\cdot \left(\mathit{\rho}\mathit{v}\right)\equiv \mathrm{0}$. This restriction is not present in COSMO's compressible formulation and this difference is accounted for in PALM's side of the nesting interface by Eqs. (1) and (2). Furthermore, at horizontal grid spacings of several kilometres, turbulence in COSMO is fully parameterized such that its flow fields are essentially free of turbulent fluctuations. PALM, on the other hand, explicitly resolves the energetic part of the turbulent spectrum.
Secondly, due to their different domain extents, both models use different approximations of Earth's surface and, as a result, use different coordinate systems. COSMO represents the planet as a perfect sphere with radius R=6371.229 km and terrain layered on top of it. Consequently, it uses a spherical coordinate system, in particular a rotatedpole system as depicted in Fig. 2a. The origin of COSMO's coordinate system is rotated to the region of interest in order to minimize grid heterogeneity. The rotation is defined in terms of the location of the rotated North Pole with the restriction that the prime meridian continues to intersect with Earth's axis of rotation and, thus, with the geographical North Pole and South Pole. In contrast, PALM is designed as a tool for simulating the atmospheric boundary layer, which implies domain sizes that are small compared to Earth's radius. Hence, Earth's surface is represented as a tangential plane and the governing equations are formulated in a Cartesian frame of reference with the z coordinate aligned with the uniform gravitation vector field and the y coordinate facing north.
Lastly, COSMO and PALM use different numerical grids which require interpolation. Both models discretize their respective governing equations on structured grids aligned with their respective coordinate axes and equidistant horizontal spacings – in the case of COSMO equidistant in rotated latitudes and longitudes, and in the case of PALM equidistant in Euclidean length. Both are based on the ArakawaCtype grid structure, where scalars are defined at the mass points at the cell centre and velocity components are staggered one half grid cell. In the vertical, both models allow for grid stretching. COSMO uses a hybrid z coordinates, with levels in the lower region following the terrain and gradually approaching an upper region with horizontally homogeneous spacings. The grid is constructed starting with the staggered velocity points, the socalled “half layers”. The “full layers”, where mass points are located, are defined as the arithmetic mean of two neighbouring half layers. PALM, on the other hand, uses a horizontally uniform grid that may contain both, parts with vertically stretched as well as equidistant grid spacings. With PALM, typical grid spacings are on the order of 100 to 1 m, while COSMO is designed for horizontal grid spacings on the order of 10 to 1 km.
Currently, INIFOR is designed to process model output of DWD's current operational configuration COSMOD2 (Baldauf et al., 2018) and its predecessor COSMODE (Baldauf et al., 2014). Both configurations operate on rotatedpole grids with the rotated North Pole located at $({\mathit{\lambda}}^{\mathrm{N}},{\mathit{\phi}}^{\mathrm{N}})=(\mathrm{170}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\mathrm{W},\mathrm{40}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\mathrm{N})$, placing the origin at $({\mathit{\lambda}}^{\mathrm{O}},{\mathit{\phi}}^{\mathrm{O}})=(\mathrm{10}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\mathrm{E},\mathrm{50}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\mathrm{N})$, close to the centre of Germany (see Fig. 2b). COSMOD2 extends the prior COSMODE domain slightly towards the north, west, and south. With horizontal grid spacings of 2.2 and 2.8 km, respectively, both configurations run at convectionpermitting resolutions (cf. Baldauf et al., 2011). COSMOD2 uses 65 vertical levels, which is up from 50 levels in COSMODE. The vertical grid spacing of the lowest cell at sea level is 20 m for both configurations, which gets further compressed over elevated terrain. The particular rules governing the vertical grid generation can be found in DWD's database manuals (Baldauf et al., 2011, 2018) and in the COSMO model documentation (Doms and Baldauf, 2018), but in the context of data interpolation on that grid, knowledge of the underlying rules is not necessary. Rather, the vertical grid is completely defined by the threedimensional field of the halflayer heights, which is static in time and available in the model output.
2.2 INIFOR preprocessing
PALM and COSMO differ in their physical formulation, i.e. their prognostic and available output variables, the representation of Earth's surface, the coordinate systems, and the structure and resolution of the numerical grids used. To translate those differences, INIFOR needs to carry out the following conceptual steps:

convert between the sets of prognostic variables (see Table 1),

project PALM's Cartesian domain onto COSMO's spherical Earth,

transform PALM's projected Cartesian coordinates to COSMO's rotatedpole system, and

interpolate COSMO data onto PALM's grid in the COSMO rotatedpole system.
In the following sections, we describe how INIFOR addresses each of these steps in detail.
Note that the data interpolation could be carried out in different coordinate systems. With INIFOR, we decided to interpolate in COSMO's rotatedpole system where the required interpolation neighbours are located on a rectangular grid leading to simple and efficient interpolation rules. We obtain the COSMO coordinates for the PALM grid points using a twostep transformation as shown in Fig. 3. First, we project the PALM grid onto COSMO's geoid (see Sect. 2.2.2), resulting in a rotatedpole system similar to COSMO's but with a different rotated North Pole. Then, we transform from the rotated PALM system to the rotated COSMO system.
2.2.1 Conversion of prognostic variables
Differences in the model formulations require conversions between the sets of prognostic equations. In our case, this includes the computation of the potential temperature and the volumetric soil moisture. INIFOR converts both quantities before interpolating them onto the PALM grid.
As for the air temperature preprocessing, INIFOR replaces the absolute temperature T provided in the COSMO output by the potential temperature given by
Here, P is the corresponding mesoscale pressure, p_{ref}=10^{5} Pa is the reference pressure, and R_{d}=287 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$ and c_{p}=1005 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$ are the ideal gas constant and specific heat at constant pressure of dry air, respectively.
For soil data, preprocessing is slightly more involved because on sea or inland water cells, COSMO's soil data are not defined. Due to the coarser grid resolution, shorelines or inland lakes do not necessarily correspond to the highresolution surface input required by PALM. In order to provide soil data at each PALM grid point, the missing information is iteratively generated by horizontal averaging of soil data from neighbouring land cells. Every iteration of this procedure generates new virtual land cells. By repeating this procedure using both the original and newly generated virtual cells, the virtual shoreline moves one COSMO cell per iteration. This procedure is currently repeated five times, before the field is used for interpolation. After the data extrapolation on the COSMO grid, the units of COSMO's soil moisture are converted to PALM's formulation. COSMO provides soil moisture as vertically integrated water density, while PALM requires the volumetric water content. The conversion is given by
where WS_{k} and Δd_{k} are the columnintegrated soil moisture and thickness of the kth soil layer in COSMO, respectively, and ρ_{w}=1000 kg m^{−3} is the approximate density of water.
2.2.2 Inverse map projection
There are multiple ways as to how the differences in the representation of Earth's surface can be resolved. Two options are illustrated in Fig. 4: (i) by a direct spatial transformation between coordinates of the rotatedpole coordinates and the tangential Cartesian system or (ii) by projecting the Cartesian system onto COSMO's geoid surface. The first approach, however, implies a change in the gravitational field: while the isosurfaces of the gravitational potential are concentric spheres (gravitation vectors point towards the geoid centre), they are parallel planes in PALM's Cartesian system. As a result, a balanced stratified atmosphere in COSMO would produce baroclinic instabilities if it was directly transformed into PALM's Cartesian system. With INIFOR, we avoid this effect by choosing the second approach and project PALM's system onto COSMO's geoid. This corresponds to the inverse of a map projection.
In particular, we use an inverse plate carrée projection which linearly maps between the Cartesian coordinates $(x,y,z)$ and the spherical coordinates $({\mathit{\lambda}}^{\mathrm{p}},{\mathit{\phi}}^{\mathrm{p}},{z}^{\mathrm{p}})$ on a sphere of radius R according to
where the superscript “p” refers to PALM. This projection defines a rotatedpole system, the Equator and prime meridian of which pass through PALM's Cartesian origin (see Fig. 2c). By requiring the y axis to point towards geographical North, we obtain a rotatedpole system of the same kind as COSMO's rotatedpole system where the prime meridian also intersects with the Earth's North Pole.
2.2.3 Coordinate transformation
When transforming the PALM to the COSMO rotatedpole coordinates, we consider the PALM system a rotatedpole system relative to the COSMO rotatedpole system, the same way the COSMO system is a rotatedpole system relative to the geographical system. Because, as we discuss below, the definition and evaluation of the transformation from PALM's to COSMO's coordinates involves forward and backward transformations between rotated systems, we begin with the general forward and backward transformations. The forward transformation from a geographical (λ,φ) to a rotatedpole system (λ^{r},φ^{r}) is obtained from spherical geometry as (Baldauf et al., 2014)
where (λ_{N},φ_{N}) are the geographical coordinates of the rotated pole. The inverse transformation is given by
The definition of the PALM rotatedpole system starts with the specification of its origin in terms of its geographical coordinates (λ_{O},φ_{O}). In order to define the transformation to the rotated COSMO system, we need to translate the PALM origin into the corresponding rotated North Pole in the COSMO system. We do this by first computing the location of the rotated North Pole (λ_{N},φ_{N}) in the geographical system as
and then, using the forward transformation in Eqs. (6) and (7), we obtain the rotated North Pole coordinates ${\mathit{\lambda}}_{\mathrm{N}}^{\mathrm{r}}={\mathit{\lambda}}^{\mathrm{r}}({\mathit{\lambda}}_{\mathrm{N}},{\mathit{\phi}}_{\mathrm{N}})$ and ${\mathit{\phi}}_{\mathrm{N}}^{\mathrm{r}}={\mathit{\phi}}^{\mathrm{r}}({\mathit{\lambda}}_{\mathrm{N}},{\mathit{\phi}}_{\mathrm{N}})$ in COSMO's frame of reference. Now the horizontal transformation between the PALM and COSMO is fully defined and we can transform PALM rotatedpole coordinates to COSMO's rotatedpole system using the backward transformation in Eqs. (8) and (9) using the PALM coordinates (λ^{p},φ^{p}) as the rotatedpole coordinates (λ^{r},φ^{r}) and COSMO's rotatedpole coordinates (λ^{r},φ^{r}) as the geographical ones (λ,φ).
Finally, the PALM domain may generally be shifted a.s.l. by specifying a nonzero domain base z_{0} in order to vertically align COSMO and PALM orography. The COSMO heights (a.s.l.) of the vertical PALM levels at z^{p} are then given by
2.2.4 Spatial interpolation
Having the COSMO rotatedpole coordinates for each PALM grid point available, we can interpolate COSMO fields by locating the appropriate interpolation neighbours and by computing the corresponding interpolation weights. We use the coordinate symbols laid out in Table 3 to describe the interpolation methodology. In particular, the COSMO rotatedpole coordinates are denoted by $({\mathit{\lambda}}^{\mathrm{r}},{\mathit{\phi}}^{\mathrm{r}},h)$, while the Cartesian PALM coordinates are $\mathit{x}=(x,y,z)$. Grid points are referenced with the indices $i,j,k$ for the PALM grid, while points on COSMO's grid are denoted by an additional hat.
Using this convention, a general interpolation scheme for an arbitrary scalar s on the COSMO grid can be formulated in terms of the weighted sum of N_{l} neighbouring values S:
Here, the indices ${\widehat{i}}_{ijk}^{l},{\widehat{j}}_{ijk}^{l}$, and ${\widehat{k}}_{ijk}^{l}$ identify the lth neighbour on the COSMO grid for the PALM grid point at x_{ijk} and ${\mathcal{W}}_{ijk}^{l}$ are the corresponding interpolation weights, which satisfy ${\sum}_{l=\mathrm{1}}^{{N}_{l}}{\mathcal{W}}_{ijk}^{l}=\mathrm{1}$. Since the scalar's values on the mesoscale grid are known, the remaining task is to compute the values of those four fields. In INIFOR, we use bilinear and trilinear interpolation requiring only the four or eight closest neighbours, respectively, but the approach may be extended to higherorder schemes by including more points. INIFOR separates horizontal and vertical interpolation, which (i) simplifies the treatment of COSMO's terrainfollowing vertical grid and (ii) enables us to adapt the horizontal scheme to other grid structures in the future, such as the triangular horizontal grid of ICON, the Icosahedral Nonhydrostatic Model. (As of writing this paper, ICON is being used as the operational global weather prediction model at DWD, and ICONLAM – its limitedarea model variant – is designated to supersede COSMO as the regional model. For a description of ICON's grid, see, for instance, the paper by Wan et al., 2013.)
Twodimensional horizontal interpolation
In the case of bilinear interpolation, Eq. (12) reduces to two dimensions, and we can drop the vertical indices $\widehat{k}$ and k and N_{l}=4. The indices ${\widehat{i}}_{ij}^{l},{j}_{ij}^{l}$ of the four neighbours $l\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}\mathit{\}}$ are
The reference COSMO indices ${\widehat{I}}_{ij}$ and ${\widehat{J}}_{ij}$ mark the bottom left neighbour point (see Fig. 5) and are obtained from the rotatedpole coordinates of the PALM grid point according to
where ${\mathit{\lambda}}_{\mathrm{0}}^{\mathrm{r}}$ and ${\mathit{\phi}}_{\mathrm{0}}^{\mathrm{r}}$ mark the lowest longitude and latitude of the COSMO grid, and Δλ^{r} and Δφ^{r} are the equidistant grid spacings in the respective directions.
Using the location of the neighbour grid points, we can compute the corresponding bilinear interpolation weights based on the nondimensional coordinates:
along the COSMO cells faces. The bilinear interpolation weights at the four neighbour points are given by
which lets us interpolate scalars using to Eq. (12).
Threedimensional interpolation
The interpolation in three dimensions is split in two steps in INIFOR: (i) a bilinear horizontal interpolation onto an intermediate grid and (ii) a linear vertical interpolation in each of its columns. The intermediate grid, hereafter indicated by an overbar, shares PALM's fine horizontal grid but features COSMO's coarser vertical levels (see Fig. 5). Concretely, the vertical levels ${\stackrel{\mathrm{\u203e}}{h}}_{\widehat{i}\widehat{j}k}$ of the intermediate grid – as well as values of the corresponding interpolation quantity ${\stackrel{\mathrm{\u203e}}{s}}_{\widehat{i}\widehat{j}k}$ – are computed using the bilinear scheme above, i.e.
where h_{ijk} indicates the COSMO grid levels and $l\in [\mathrm{1},\mathrm{4}]$ iterates over the four neighbouring COSMO columns. In the second step, the interpolated values $\widehat{s}$ are interpolated vertically from the intermediate to the PALM target grid:
Below the lowest intermediate grid level ${\overline{h}}_{ij\mathrm{1}}$ of each column, all variables are extrapolated downwards as a constant. Beyond that, there is currently no further terrain adaptation made to blend the terrain from the coarse mesoscale to the fine microscale resolution.
Interpolation of velocities
The transformation in Eqs. (8) and (9) between the two rotatedpole systems (see Fig. 2b) involves a rotation as the meridians of the original and rotated system are generally not parallel. This deviation angle, the socalled “meridian convergence”, increases as we move away from the reference meridian, and its distribution is given by Baldauf et al. (2014) as
We obtain the Cartesian velocity components in the PALM system by rotating COSMO's spherical velocity components by the local meridian convergence according to
Since on the staggered ArakawaC grid U and V are not defined at the same location, INIFOR first interpolates horizontal velocities onto COSMO's mass points and then rotates the interpolated velocity vectors using Eq. (20). Apart from this preprocessing, velocities are interpolated the same way scalars are. The resulting interpolation neighbours and weights for velocities, however, do differ from those of scalars because u and v on PALM's staggered grid are in turn defined half a grid cell away from PALM's mass points.
Averaging of profile data
INIFOR provides the option to initialize and force PALM with threedimensional atmospheric data (LOD = 2) or with averaged profiles (LOD = 1). The latter has the advantage that, for large setups, INIFOR preprocessing is easier to handle in practice because less memory is required on the preprocessing machine and the resulting dynamic driver is greatly reduced in size because threedimensional arrays are omitted. INIFOR produces profile data by first averaging along COSMO levels and then interpolating in the vertical direction.
Concretely, this is done carrying out the following steps. We first define the averaging region as a horizontal box bounded by the minimum and maximum rotated longitude and latitude of the PALM domain. Once all COSMO columns in the region are identified, we compute the average vertical grid levels of the terrainfollowing COSMO grid and then compute the vertical neighbours and weights for every PALM level relative to the averaged COSMO levels. The average profile (denoted in the following by the double bar) is then formed by scanning through the N_{c} columns of the averaging region and adding the vertically interpolated values on every PALM level according to
with $({\widehat{i}}_{\mathrm{c}},{\widehat{j}}_{\mathrm{c}})$ being the indices of the N_{c} COSMO columns in the averaging region.
2.2.5 Program structure
INIFOR's program structure is organized around the set of netCDF variables that are to be computed for the dynamic driver. The dynamic driver contains individual netCDF variables for each combination of prognostic variable and model boundary, e.g. netCDF variables for the u velocity component at the east, south, west boundaries, and so forth. Internally, INIFOR maintains a list with representations of each of those netCDF variables. Each is associated with the netCDF metadata required to handle data input and output, the computational task – averaging profiles or interpolating fields in 2D or 3D – as well as with the appropriate interpolation grids. The latter contain both grid point coordinates and interpolation neighbours and weights. Generally, different variables that are defined at the same grid point type also share the interpolation parameters. For instance, the horizontal (intermediate) interpolation grid for scalars is shared among netCDF variables for the initial soil moisture and temperature fields as well as the top boundary conditions for w, θ, and q_{v}. Consequently, interpolation grids with their corresponding interpolation parameters are computed once and reused each time step and shared among multiple output variables.
This is reflected in INIFOR's program flow, which is depicted in Fig. 6. It is divided into two main sections: a setup phase and the main loop. During the setup phase, INIFOR constructs the required model and interpolation grids. This involves defining and transforming the coordinates of the PALM interpolation grids as well as precomputing interpolation neighbours and weights for every grid point. During the main loop, INIFOR then iterates through the output netCDF variables and time steps, reusing precomputed interpolation parameters that are associated with each variable. Each main loop iteration is structured into reading input data, preprocessing input data, interpolation, and data output. The preprocessing step is dependent on the kind of input and includes the extrapolation of soil data into water cells, conversions between model formulations (such as the computation of the potential temperature or the computation of volumetric soil moisture; see Sect. 2.2.1), and velocity vector rotation (see Sect. 19).
As input data, INIFOR reads hourly netCDF files containing COSMO analyses. Each hourly input is processed separately and translated into one instantaneous boundary condition in the dynamic driver. Input data are not interpolated temporally in INIFOR but rather in PALM during the simulation as described in Sect. 2.
2.3 Superposition of boundary conditions with synthetic turbulence
With the generation of timedependent boundary conditions from mesoscale model output in a preprocessing step and online processing of the boundary data, PALM is enabled to simulate more realistic scenarios considering timeevolving synoptic conditions. However, due to the nature of RANS models, turbulence is parameterized and thus the boundary values are free of any turbulent fluctuations. Mirocha et al. (2014) showed that without adding perturbations the turbulent flow needs several tens of kilometres to sufficiently develop. In order to accelerate the spatial development of turbulence in PALM in our mesoscale nesting approach, we employed the synthetic turbulence generator by Xie and Castro (2008) where perturbations are added onto $u,v,w$ – components imposed at the lateral boundaries. In the following, the preliminary boundary values without any perturbations added are indicated by an overbar.
To obtain turbulent flow components u_{i,b} on the lateral boundaries, spatially and temporally correlated disturbances ${u}_{i}^{\prime \prime}$ are imposed onto the preliminary velocity components ${\stackrel{\mathrm{\u203e}}{u}}_{i,\mathrm{b}}$:
a is the amplitude tensor that is calculated from the Reynoldsstress tensor r. To consider crosscorrelations between the velocity components, Lund et al. (1998) suggested a Cholesky decomposition to compute a recursively by
Depending on characteristic length scales L and timescales T of the flow, which are defined individually for each velocity component in each spatial direction, ${u}_{i}^{\prime \prime}$ computes as
with Δt being the actual LES time step and the twodimensional spatially correlated disturbances
The subscripts m and l indicate grid positions at the lateral boundary, ${N}_{i}=\mathrm{2}{L}_{i}/\mathrm{\Delta}{x}_{i}$, with Δx_{i} being the grid spacing. ζ_{i} indicates a set of equally distributed random velocities with zero mean and unit variance that are individually computed for each velocity component. Finally, the spatial filter function computes as
with ${b}_{i}^{*}=\mathrm{exp}\left(\frac{\mathit{\pi}\leftk\right\mathrm{\Delta}{x}_{i}}{{L}_{i}}\right)$. With this approach, the imposed ${u}_{i}^{\prime \prime}$ reflects the prescribed Reynolds stress as well as the spatial and temporal correlation according to L_{i} and T_{i}, respectively. At this point we want to note that Xie and Castro (2008) assumed an exponentially decaying autocorrelation function in Eq. (24), as well as for the formulation of the filter coefficients ${b}_{i}^{*}$. Mordant et al. (2001) have shown that this is a valid approach for sheardriven flows. To our knowledge, it has not been proven yet whether an exponentially decaying autocorrelation function is an appropriate choice for stratified flows. Due to the lack of universally valid alternatives, however, we employed the formulation described above also for stratified flows despite its possible limitations.
From a mathematical point of view, the imposed fluctuations should have zero mean. Due to a finite sample of random numbers and the finite number of discrete grid points, however, the fluctuations have mean values slightly different from zero in practice. In order to overcome this, Kim et al. (2013) proposed a correction for the boundary normal flow component in order to maintain constant mass flux through the boundary. In order to avoid the perturbations imposed onto the nonnormal components having nonzero mean too, e.g. nonzero mean w and v components at the western model boundary, we correct the imposed turbulent velocity components as
with S being the surface area of the respective lateral boundary.
For nonstationary flows, an inflow boundary can become an outflow boundary, and vice versa. Hence, the turbulence generator is applied at each lateral boundary simultaneously, while at opposite boundaries (west and east, as well as north and south) we use the same Ψ_{i} and thus the same set of random numbers (velocities). By doing so, we save computational resources because the same set of Ψ_{i} is already available on the west/east and north/south boundary according to our parallelization strategy. Further, perturbations are imposed at the end of each LES time step at the last Runge–Kutta substep right before the Poisson equation is solved to fulfil divergencefree flow.
From Eq. (25) it becomes obvious that the computational demand to calculate the spatially correlated disturbances is a function of the turbulent length scales – doubling the turbulent length scale leads to a quadrupling of the number of elements in the summation. For example, turbulent length scales vary significantly in time, reaching values >2000 m (see Fig. 13) and altering the computational demand of the turbulence generator within the course of the day. Especially for nonparallelized implementations of the method (e.g. Zhong et al., 2021), this may become a limiting factor with respect to the computational demand. In order to overcome this and to balance the computational load over various processes, we parallelized the synthetic turbulence generator using the Message Passing Interface (MPI). To achieve this, we made use of the 2D domain decomposition used in PALM (Maronga et al., 2015). The imposed disturbances are computed locally on each MPI process that belongs to a lateral boundary. In order to avoid that only processes residing at the domain boundaries execute the computationally heavy code of Eqs. (25) and (26), while other processes are on hold, the computation of the filtered random numbers is divided in vertical direction by the number of processes in normal direction to the inflow boundary (e.g. on the left boundary, computation is distributed over npex parts, where npex is the number of subdomains, or MPI processes, along x), while the filtered random numbers are gathered on the boundary process, subsequently. In the case of large ${L}_{i}/\mathrm{\Delta}{x}_{i}$, this significantly reduces the required computational time of the synthetic turbulence generator (see Sect. 4.6) as each MPI process needs to compute only a subset of perturbations. If ${L}_{i}/\mathrm{\Delta}{x}_{i}$ is of only low value, however, the additional parallelization can be omitted by choice as the additional MPI communication can consume the benefit of the distributed calculation. Results from a performance test for the turbulence generator are discussed in Sect. 4.6.
The random numbers ζ, which are defined on the discretized grid, are distributed over several MPI processes, and each process only knows its required set of random numbers. For example, suppose 𝚗𝚡𝚕_{𝙰} (𝚗𝚡𝚕_{𝙱}) and 𝚗𝚡𝚛_{𝙰} (𝚗𝚡𝚛_{𝙱}) are the left and right local domain boundary indices along the x direction on MPI process A (B), respectively, with ${\mathtt{\text{nxl}}}_{\mathtt{\text{B}}}={\mathtt{\text{nxr}}}_{\mathtt{\text{A}}}+\mathrm{1}$. On process A, random numbers need to be defined for the index range ${\mathtt{\text{nxl}}}_{\mathtt{\text{A}}}{\mathtt{\text{N}}}_{\mathtt{\text{x}}}:{\mathtt{\text{nxr}}}_{\mathtt{\text{A}}}+{\mathtt{\text{N}}}_{\mathtt{\text{x}}}$; equivalently, on process B, ${\mathtt{\text{nxl}}}_{\mathtt{\text{B}}}{\mathtt{\text{N}}}_{\mathtt{\text{x}}}:{\mathtt{\text{nxr}}}_{\mathtt{\text{B}}}+{\mathtt{\text{N}}}_{\mathtt{\text{x}}}$, meaning that on processes A and B random numbers partly overlap, e.g. within the index range ${\mathtt{\text{nxr}}}_{\mathtt{\text{A}}}{\mathtt{\text{N}}}_{\mathtt{\text{x}}}+\mathrm{1}$ and 𝚗𝚡𝚕_{𝙱}−𝙽_{𝚡}. In order to obtain the same ζ within the overlapping index range, we have to assure that the set of random numbers do not depend on the horizontal domain decomposition. This is achieved by linking the seed of the employed randomnumber generator to the grid index which is then independent of the domain decomposition. With respect to the computation of ζ, MPI communication can thus be reduced to only exchange data to compute its mean and variance. Note that according to Eq. (25), random numbers are also required for locations outside the model domain to allow for the computation of the spatial correlations, especially near the domain boundaries and for larger values of ${L}_{i}/\mathrm{\Delta}{x}_{i}$. Therefore, the required randomnumber arrays are allocated with an offset so that all required values fit into the respective arrays. In the event that ${L}_{i}/\mathrm{\Delta}{x}_{i}$ increase or decrease during the simulation, the respective arrays are resized.
In order to create time and heightdependent synthetic turbulence, respective information about the Reynolds stresses, as well as turbulent length scales and timescales for the velocity components are required. For stationary flows, this information can be deduced from observations or from cyclic precursor simulations (Xie and Castro, 2008). However, for nonstationary flows with pronounced diurnal cycles and/or changing synoptic conditions, running precursor simulations is practically not feasible. Also, to take this information from the mesoscale model output is also not possible since this detailed information is often neither available nor part of the operational output. Hence, to allow for an adjustment of the synthetic inflow turbulence to changing atmospheric conditions, we parametrize the Reynolds stresses based on the timedependent mesoscale inflow profiles. We follow the set of parameterizations presented by Rotach et al. (1996) which they employed in stochastic dispersion modelling. Please note that the following set of parameterizations refers to the stream and spanwise components of the Reynolds stress that are not necessarily parallel to the x or y axis, respectively. In order to emphasize this, we indicate stream and spanwise components with a tilde in the following.
Rotach et al.'s 1996 parameterizations are based on parameterizations Brost et al. (1982) derived from observations in stratocumulustopped marine boundary layers, which often differ in their vertical structure and turbulence production compared to boundary layers over land. However, since Rotach et al. (1996) have successfully validated the set of parameterizations against observations over land for a wide range of stability regimes, we are confident that the chosen set of parameterizations can be universally employed. Based on the original formulation by Brost et al. (1982), the variance of the streamwise flow component ${\stackrel{\mathrm{\u0303}}{r}}_{\mathrm{11}}$ is parameterized following Rotach et al. (1996):
who added a correction term (first term in Eq. 28) proposed by Gryning et al. (1987) to account for unstable nearsurface stratification. Here, u_{∗} is the friction velocity, κ=0.4 the von Kármán constant, L_{o} the Obukhov length, z_{i} the boundarylayer depth, and z the height above the ground. For neutral and stable situations, the first term is ignored. Similarly, we estimate the variance of the spanwise flow component by adding a correction term to the original formulation proposed by Brost et al. (1982):
The profile of vertical velocity variance is taken from Gryning et al. (1987) as
with w_{∗} being the convective velocity scale. The vertical transport of horizontal streamwise and spanwise momentum is estimated by Brost et al. (1982) as
and
respectively. To our knowledge, there exists no comparable formulation to estimate ${\stackrel{\mathrm{\u0303}}{r}}_{\mathrm{21}}$ in the literature. Hence, we decided to simply set ${\stackrel{\mathrm{\u0303}}{r}}_{\mathrm{21}}=\sqrt{{\stackrel{\mathrm{\u0303}}{r}}_{\mathrm{31}}^{\mathrm{2}}+{\stackrel{\mathrm{\u0303}}{r}}_{\mathrm{32}}^{\mathrm{2}}}$, assuming isotropy of horizontal and vertical transport of horizontal momentum. To estimate the boundarylayer depth for a wide range of stability regimes, including buoyancy and purely sheardriven boundary layers, we calculated z_{i} from a bulk Richardson number criterion according to Heinze et al. (2017) based on the bulk Richardson number:
Starting at the surface, z_{i} is defined as the height where Ri_{b} first exceeds the critical bulk Richardson number Ri_{b,c}=0.25, which revealed to be a robust criterion to estimate the depth of the layer with significant turbulent transports caused by the presence of the surface (Heinze et al., 2017). Here, u_{h} denotes the horizontal wind speed from mesoscale model input, θ_{v} the virtual potential temperature, θ_{v,s} the virtual potential surface temperature inferred from the second prognostic level above the surface, following Heinze et al. (2017), and g the acceleration of gravity. In the case of LOD = 1 input, z_{i} is determined based on the mean profiles prescribed at the lateral boundaries, while in the case of LOD = 2 input (xz and yz slices of boundary data), z_{i} is determined locally at each (x,y) boundary grid point and averaged horizontally afterwards.
In contrast to z_{i}, which can be well estimated from the mean boundarylayer profiles, the friction velocity u_{∗}, used in Eqs. (28) to (30), strongly depends on the surface roughness and local surface properties. To generate turbulence that roughly reflects the mean conditions within the LES domain, we decided to estimate u_{∗} by horizontally averaging the values as used in the surface parametrization according to the Monin–Obukhov similarity theory (see Maronga et al., 2020, Eq. 28) within the LES model domain. The same is done also to obtain L_{o}. By doing this, we are aware that u_{∗} and L_{o}, and thus also the parameterized Reynolds stress, are not entirely independent of each other, since adjustment effects of the turbulent flow may modify u_{∗} and L_{o} locally near the inflow boundaries, which in turn feeds back into the Reynoldsstress parametrization again modifying u_{∗} and L_{o} near the inflow boundaries. However, for sufficiently large model domains, this feedback loop is negligible, as we discuss in Sect. 4.3. The convective velocity is computed as ${w}_{\ast}={\left(g{H}_{\mathrm{0}}{z}_{i}/{\mathit{\theta}}_{\mathrm{s}}\right)}^{(\mathrm{1}/\mathrm{3})}$, which is only defined for positive values of the mean surface sensible heat flux H_{0}, else it is set to zero so that r_{3,3} remains defined. Even in neutral or stable boundary layers, a vertical velocity variance can be observed. Hence, in order to account the Reynoldsstress parametrization for a wide range of stability regimes, we followed Rotach et al. (1996), who replaced w_{∗} in Eq. (30) by the mixed velocity scale w_{m} (Troen and Mahrt, 1986) with
with β=0.6 according to Holtslag and Boville (1993).
Note that Eqs. (28), (29), (31), and (32) describe the flow characteristics in the streamwise and spanwise framework, which is indicated by the tilde. In a mesoscale nesting with changing wind directions, however, the streamwise and spanwise flow directions do not necessarily coincide with the Cartesian grid axis which the prognostic velocity components relate to. Hence, the individual components of $\stackrel{\mathrm{\u0303}}{r}$ are projected back onto the Cartesian grid by rotation about the vertical axis by the rotation angle defined by $\mathrm{arctan}\left(v/u\right)$.
Further, the synthetic turbulence generation requires information about the turbulent length scales and timescales. The turbulent timescale of the flow is estimated according to Brost et al. (1982) with
Parameterizations of turbulent length scales exist for the lower part of the boundary layer, (e.g. Flay and Stevenson, 1988; Salesky et al., 2013; Li et al., 2017; Emes et al., 2019); however, no parametrization of turbulent length scales that cover the entire depth of the boundary layer and all stability regimes exists to date to our knowledge. Hence, we calculate turbulent length scales of the flow components according to Tennekes and Lumley (1972) using the parameterized Reynolds stress and the timescale:
with $i\in \mathrm{1},\mathrm{2},\mathrm{3}$ indicating the streamwise, spanwise, and vertical directions, respectively.
Assuming that turbulence is only present within the boundary layer, r, T, and L are faded for z>z_{i} with
Here, the fading function is designed so that Φ(z) rapidly decreases above the boundary layer.
In the case of nonstationary flows, the turbulence parameters r,T, and L are adjusted hourly by default, but the frequency can be also modified by the user. We note that updating the turbulence parameters violates the temporal correlation expressed in Eq. (24). Hence, we performed a test simulation where we omitted the first term in Eq. (24) after updating the turbulence parameters and we found no difference with respect to the spatial development of the flow (not shown), indicating that occasional violations of the temporal correlation have no significant effect on the development of the flow.
In order to test the implemented mesoscale nesting approach, we selected a particular weather scenario with a developing daytime convective boundary layer (CBL) that features advective conditions with moderate wind speeds and changing meanwind direction. Moreover, the scenario is characterized by little to no cloud cover which is attributed to the fact that PALM cannot capture highaltitude clouds yet (due to missing icephase physics, planned) and thus cannot realistically reproduce the prevailing radiative forcing. We simulated one diurnal cycle of the evolving CBL starting at 00:00 UTC on 7 May 2016, for a domain located east of Berlin, Germany. The boundary layer on this day was characterized by clearsky conditions and moderate mean boundarylayer wind speeds of about 7–8 m s^{−1} from the east, later turning to the southeast during the morning hours. Figure 7 shows horizontal mean vertical profiles of horizontal wind u_{h}, meanwind direction, potential temperature θ, and specific humidity q_{v} at different points in time obtained from COSMO. During nighttime, the profiles indicate a stably stratified layer up to z=800 m, while at 08:00 UTC the stably stratified layer gets successively eroded by the beginning surface heating. Later in the day, a wellmixed CBL develops with maximum z_{i}≈2400 m. At about 16:00 UTC, the evening transition sets in and again a stably stratified layer develops near the surface. We assumed a horizontally homogeneous and flat surface instead of the particular terrain, land use, and buildings in the Berlin area. We made this idealization in order to be able to determine adjustment fetch lengths under timedependent inflow conditions. Surface heterogeneities in terms of land use or terrain would modify the turbulent flow field locally, making it impossible to disentangle changes in the turbulent flow due to adjustment effects and surface heterogeneity. We employed the embedded landsurface model (see Maronga et al., 2020; Gehrke et al., 2020) to obtain sensible and latent heat fluxes at the surface, which we assumed to be fully covered with short grass. We chose this setup as a tradeoff roughly reflecting the prevailing land use in this area with distributed farm and grassland, forest patches, and urbanized environments. The soil was initialized with horizontally homogeneous profiles of soil temperature and soil moisture taken from COSMO. Since the soil properties in COSMO are aggregated over various surfaces and soil types, the soil conditions are not necessarily in equilibrium with the assumed grass surfaces as well as selected atmospheric conditions. Hence, we ran a 2 d surface spinup as a precursor to the 3D simulation (Maronga et al., 2020) in order to bring the soil into equilibrium with the atmospheric conditions and to avoid spinup effects that may lead to varying heat fluxes at the beginning of the simulation. The incoming short and longwave radiation was modelled using the Rapid Radiative Transfer Model for Global Models (RRTMG, Clough et al., 2005).
The simulated domain is located at 52.5^{∘} N, 13.7^{∘} E (PALM origin), and extends over $\mathrm{43.2}\times \mathrm{43.2}\times \mathrm{4.7}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{3}}$ in the x, y, and z direction, respectively, with an isotropic grid spacing of 25 m. Above z=3 km – approximately 600 m above maximum boundarylayer depth – the vertical grid was successively stretched up to 50 m vertical grid spacing. This allows for proper resolution of the convective boundary layer, though we note that the nighttime stably stratified boundary layer is only poorly represented with the chosen grid spacing. The advection terms were discretized with a fifthorder upwind scheme (Wicker and Skamarock, 2002); for the time stepping, we applied a thirdorder Runge–Kutta method according to Williamson (1980).
We performed two simulations with different lateral boundary conditions. In the first simulation, hereafter referred to as REF, the boundary conditions were given as LOD = 1, i.e. horizontally averaged profiles. Unless it is not further noted, we refer to this simulation in the following analysis. In the second simulation, heterogeneous boundary values were prescribed (LOD = 2). This second simulation will be used to check whether the LES simulation results depend rollconvection artefacts in the convection grey zone.
Hence, we calculated resolvedscale variances of the velocity components by
while the angle brackets indicate a time average over half an hour and the prime indicates the turbulent fluctuation. The resolvedscale TKE was computed as $\text{TKE}=\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\cdot \phantom{\rule{0.125em}{0ex}}\sum \phantom{\rule{0.125em}{0ex}}\u2329{u}_{i}^{\prime}\phantom{\rule{0.125em}{0ex}}{u}_{i}^{\prime}\u232a$. For each grid point location, we determined its distance to the inflow boundary at the given wind direction. Therefore, we calculated virtual backward trajectories for each halfhour interval from the current meanwind direction; subsequently, we determined the distance d between the sampling location and the intersection point of the backward trajectory with the closest inflow boundary. Note that this analysis can be simplified in the case of LOD = 1 forcing, where the distance of each grid location to the next inflow boundary can be inferred directly from a linear equation without using backward trajectories. For LOD = 2 forcing, on the other hand, backward trajectories are still needed since the lateral inflow, and thus the wind direction can change significantly along the lateral boundaries. Finally, variances were averaged over similar distances to the inflow boundary, while we sorted similar distances into equally sized bins of 100 m to obtain sufficiently large sample size for each discrete distance.
Please note that in this study we will mainly focus on convective conditions, especially with respect to the spatial development of the flow. The nighttime stable flow is only poorly resolved at the given grid spacing, making it difficult to make reliable conclusions concerning the flow adjustment. Here, we will refer to followup studies.
In the following section, we show results from a mesoscale nested LES for a diurnal cycle. In order to better guide the reader through this section, we will first give a short outline of what to expect. First, we describe the boundarylayer structure and its development over the diurnal cycle in the LES as well as in COSMO. Subsequently, we discuss differences between the LES and COSMO with respect to the boundarylayer representation and its implications in a nested simulation. In the following, triggered by imposed timedependent synthetic turbulence, we focus on the spatial development of the turbulent flow within the LES domain and determine adjustment lengths where the turbulent flow is fully developed. In addition, we present results on how rolllike structures emerging in the COSMO simulation propagate into the LES. Moreover, we discuss implications near the LES domain inflow and outflow boundaries. Finally, we look at a more technical issue and demonstrate the computational efficiency of the synthetic turbulence generation.
4.1 Boundarylayer structure
Figure 8 shows vertical profiles of the velocity variances horizontally averaged over a 10 × 10 km^{2} at the centre of the domain where the turbulent flow has been already adjusted (see Sect. 4.3). The variances show a pronounced diurnal cycle, with only small values during nighttime and the morning hours. At 08:00 UTC, when the surface heating sets in, a doublepeaked profile can be observed. Later, the variances increase, with the horizontal variances exhibiting a maximum near the surface. With increasing height, the horizontal variances decrease, while at 13:00 UTC a secondary peak can be observed near the boundarylayer top. The vertical variances peak in the middle part of the boundary layer and approach zero at boundarylayer top. In the afternoon and evening hours, the value of the variances again decreases and the boundary becomes shallower. Overall, the variances show a typical diurnal cycle for clearsky conditions (André et al., 1978).
Figure 9 shows vertical profiles of the potential temperature at the inflow boundary and from the inner part of the domain, where the turbulent flow is spatially fully developed. The averaging region for the innerdomain profiles is indicated in Fig. 12a. At 09:00 UTC (dashed lines), the COSMO inflow profile indicates a warmer boundary layer within the lowest 500 m compared to the innerdomain profile, while further above (within the residual layer) the profiles from COSMO and the inner domain coincide. At 10:00 UTC, the imposed COSMO inflow profile already indicates an unstable stratification within the lower part of the boundary layer and a wellmixed layer up to z=2000 m, while the innerdomain profile indicates an unstable stratification only within the surface layer and a wellmixed layer above reaching only up to z=1400 m, where the potential temperature profile indicates similar values compared to the profile 1 h before. This means that between 09:00 and 10:00 UTC the boundary layer in COSMO develops more rapidly, where the stably stratified layer gets completely eroded, the residual layer becomes convective and the boundary layer grows significantly, while in PALM the boundary layer develops less rapidly and only the stably stratified layer gets eroded. At 13:00 UTC, the shapes of the COSMO inflow and innerdomain profiles are similar, indicating similar boundarylayer depth, though the inflow profile indicates a warmer boundary layer compared to the region further downwind of about 0.25 K. At 16:00 UTC, the COSMO inflow profile indicates already a weakly stable stratification below z=900 m, while the innerdomain profile still shows a vertically wellmixed boundary layer. At all points in time shown, the potential temperature profiles above the boundary layer do not change significantly, indicating only small horizontal temperature advection on the mesoscale. We calculated the advection tendency at 10:00 and 16:00 UTC to −0.1 and 0.05 K h^{−1} (not shown), respectively, meaning that largescale horizontal advection of temperature cannot explain the mismatch of the temporal boundarylayer development between COSMO and PALM. However, this does not necessarily exclude local advection on the mesoscale within the boundary layer.
Figure 10 shows the surface net radiation and surface heat fluxes for COSMO and simulation REF during the course of the day. During the night, the surface net radiation and sensible heat flux are significantly smaller in COSMO, indicating more cooling of the surface which results also in more negative surface sensible heat flux. During daytime, COSMO and PALM show a similar diurnal cycle of surface net radiation with comparable peak values at noon and only small differences during the course of the day. The available energy at the surface, however, is differently partitioned into the surface latent and sensible heat flux between PALM and COSMO. The sensible heat flux in COSMO shows slightly higher values compared to PALM between 10:00 and 12:00 UTC but significantly lower values during the afternoon hours where the sensible heat flux approaches zero at about 16:00 UTC corresponding to the stabilization of the COSMOsimulated boundary layer (see Figs. 7 and 9), whereas PALM still simulates a positive sensible heat flux of $>\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}$. In contrast, the surface latent heat flux in COSMO shows larger values in the afternoon compared to PALM, meaning that the bulk of the available energy is partitioned into the surface latent heat flux being not available for heating the boundary layer.
In addition, we performed a simulation where we prescribed the diurnal cycles of H_{0} and LE_{0} with values taken from COSMO as shown in Fig. 10 rather than computing them using the landsurface model. Hereafter, we refer to this simulation as PSF. This test was motivated to check whether the discrepancy between the inner domain and the COSMO inflow profile can be attributed to a possible misrepresentation of the surface energy balance attributed to the idealized setup with a homogeneous grass surface rather than a more realistic surface. However, except for minor differences with a slightly cooler boundary layer in the morning (see Fig. 9), a slightly shallower boundary layer at noon, and a slightly cooler boundary layer in the afternoon, we found no significantly different structure of the boundary layer between simulations PSF and REF. This shows that the discrepancy between the inner domain and the COSMO inflow profile can not be attributed to any misrepresentation in the surface energy balance compared to COSMO. Another possible explanation is the advection of boundarylayer characteristics that have already developed in COSMO locally further upstream and a timelagged representation thereof due to horizontal averaging and the only hourly resolution of the forcing data. However, we have not analysed this issue further and postpone investigations of this issue until future studies.
As the COSMO profiles are mapped onto the inflow boundary, the more rapid evolution or the earlier stabilization of the boundary layer at 10:00 and 16:00 UTC, respectively, also propagate into the PALM model domain. Figure 11 shows vertical crosssections of the potential temperature and corresponding boundarylayer height. At 10:00 and 13:00 UTC, according to the higher potential temperature at the inflow boundary compared to the inner part of the domain as shown in Fig. 9, the potential temperature within the boundary layer and the boundarylayer height decrease with increasing distance to the inflow boundary. This indicates that at these points in time, a deeper and warmer boundary layer is advected into the PALM domain. In contrast, at 16:00 UTC, when the inflow potential temperature profile (see Fig. 9) already indicates a stable inflow with lower values of potential temperature, the potential temperature and boundarylayer height increase with increasing distance to the inflow boundary (Fig. 11c, d). Especially the spatial gradient of the boundarylayer height close to the inflow boundary indicates that a significantly shallower boundary layer is advected into the PALM domain which further propagates downwind (Fig. 11d) later on. Also, with increasing distance to the inflow boundary, more and deeper convective updrafts, indicated by higher values of potential temperature, can be observed, while close to the inflow boundary only shallow convective updrafts occur. This suggests that the stable stratification near the inflow boundary suppresses convection in the later afternoon. In particular, the horizontal difference in the boundarylayer structure at 16:00 and 16:30 UTC shows that temporal changes on the inflow temperature (and humidity, not shown) reach the inner part of the model domain with a time lag. In this setup, it takes about 1 to 1.5 h until the signals imposed at the inflow boundary reach the outflow boundary, meaning that the boundary layer becomes horizontally heterogeneous during the transition phase. We note that this is in contrast to the largescale forcing approach by Heinze et al. (2017) where largescale advection and nudging terms were considered at each location at the same time so that the LES solution can approach the mesoscale solution as a whole, though also with this approach the transition of the LES towards the mesoscale mean state is time lagged according to the applied nudging timescale.
The temporal change in inflow conditions can also be observed visually in the vertical velocity shown in Fig. 12. At 10:00 UTC, where a deeper boundary layer accompanied with more energetic synthetic disturbances is advected into the model domain, the up and downdrafts close to inflow boundary show a larger amplitude compared to the up and downdrafts further downwind. In contrast, at 16:00 UTC, where a more shallow and stable boundary layer accompanied with only small synthetic disturbances (see Fig. 13) is advected into the model domain, the amplitude of the up and downdrafts is only small near the inflow boundary and increases farther downwind.
4.2 Characteristics of imposed turbulence
The parameterized Reynoldsstress components depend on the inflow profiles obtained from the mesoscale model input, i.e. z_{i}, as well as on u_{∗}, w_{∗}, and H_{0}, which were obtained from horizontal averaging over the entire model domain. While z_{i} determines the depth of the boundary layer, and thus the relevant window for the Reynolds stresses, the latter three determine the amplitudes of their components.
Figure 13 shows the characteristics of the synthetic turbulence imposed at the lateral boundaries. During nighttime, Reynolds stresses and turbulent length scales are only small and defined in a shallow layer up to about z=200 m, attributed to the shallow boundarylayer depth of the stably stratified layer. Later, when convection sets in and the boundarylayer depth increases, the values of the Reynoldsstress components and the length scales increase as well, with maximum length scales of about L=2500 m at 13:00 UTC. The shape of the Reynoldsstress and lengthscale profiles does only change slightly during the simulation due to the turning wind direction where the contribution from the streamwise and the spanwise parametrization to the components on the Cartesian grid slightly changes. Above the boundarylayer top, where the fading function in Eq. (37) becomes active (see the kink in the profiles), the Reynoldsstress components and length scales rapidly approach zero. For instance, at 13:00 UTC at about z=2700 m, the length scales approach zero, while the amplitudes of the imposed disturbances are already nearly zero, meaning that only small perturbations are added above the boundarylayer top, which will be quickly dissipated. Qualitatively, the parameterized Reynoldsstress components resemble the variance profiles shown in Fig. 8, though r_{11} (r_{22}) are underestimated (overestimated) compared to ${u}^{\prime}{u}^{\prime}$ (${v}^{\prime}{v}^{\prime}$) around noon, respectively, and do not account for the secondary peak near the boundarylayer top. Furthermore, at 10:00 UTC (and 16:00 UTC), it strikes that the Reynoldsstress components indicate a deeper (shallower) boundary layer compared to the velocity variances, respectively, which is in accordance with the horizontally heterogeneous structure of the boundary layer during transition periods (see Sect. 4.1), which again is due to the fact that the onset and offset of convection are shifted between COSMO and PALM.
In summary, the parameterized Reynolds stresses resemble the variances profiles created by the LES itself reasonably well during the course of the day. However, we emphasize that the imposed turbulence is only considered to be a rough description to resemble the secondorder statistics of a fully adjusted flow, while its spectral distribution or higherorder moments are not accounted for.
4.3 Spatial development of the flow
Figure 14 shows the spatial development of the resolvedscale TKE depending on the distance to the inflow boundary at different heights and points in time. Please note that the upper abscissa depicts the dimensionless distance to the inflow boundary using convective scaling, which was originally developed by Willis and Deardorff (1976) to scale Lagrangian dispersion experiments. Here, we use this to scale the travel time $d/{u}_{\mathrm{h}}$ of a signal imposed at the lateral boundary with the eddyturnover time ${z}_{i}/{w}_{\ast}$ in the CBL, where d indicates the distance to the inflow boundary in wind direction and u_{h} indicates the mean boundarylayer wind speed (averaged over the depth of the boundary layer). The use of this scaling is in accordance with MuñozEsparza and Kosović (2018), who argue that the flow development primarily scales with ${u}_{\mathrm{h}}/{w}_{\ast}$, describing well the dominant transition mechanism. At 10:00 UTC, the TKE peaks at about 3 and 5 km downstream of the inflow boundary within the lower and upper parts of the CBL, respectively. This can also be observed visually in Fig. 12, where the amplitude of the up and downdrafts close to the inflow boundary appears stronger. This TKE overshoot is a result of the boundarylayer adjustment process, where the peak location corresponds to the distance where the bulk of the initially uprising thermals reaches the inversion layer, starting to entrain warmer air from the free atmosphere into the boundary layer. This entrainment then slightly stabilizes the boundary layer, which in turn also decelerates the uprising thermals. Further downstream, the TKE gradually decreases, approaching a constant value at about 20 to 30 km downstream of the inflow boundary, with the upper parts of the boundary layer requiring the largest fetch. (Here, the required fetch length is quantified as the distance where the TKE does not deviate by more than 10 % from the target TKE, defined as the spatial average over $\mathrm{30}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\le d\le \mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{km}$.) At 13:00 UTC, the TKE peaks again close to the inflow boundary, though the amplitude of the peak value is smaller and closer to it equilibrium value where the flow has been spatially fully developed. Especially within the upper part of the CBL, the TKE reaches a nearly constant value at about 10 to 15 km downstream of the inflow boundary. The dimensionless distance until the flow has been developed is about $\mathrm{2.0}\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{h}}{z}_{i}/{w}_{\ast}$ at 10:00 and 13:00 UTC, meaning that the flow needs about two eddy turnovers to become fully developed.
At 16:00 UTC the situation becomes qualitatively different. Even though the TKE within the lower part of the CBL peaks again close to the inflow boundary and approaches a nearly constant value farther downstream, it gradually increases starting at about 20 km downstream. Within the middle and upper parts of the CBL, the TKE is close to zero near the inflow boundary, according to the only small amplitude of the imposed synthetic turbulence (see Fig. 13), and gradually increases towards the outflow boundary. This can also be observed in Fig. 12, where the up and downdrafts near the inflow boundary are only weak and become stronger further downstream, indicating that turbulence first needs to develop spatially. Due to the horizontally heterogeneous boundary layer with already weakly stable boundary conditions near the inflow boundary and still convective conditions farther downstream, the TKE does not reach an equilibrium value and a spatial adjustment length cannot be accurately determined for this point in time.
In order to investigate how the structure of the turbulent flow develops, Fig. 15 shows the corresponding horizontal profiles for the skewness of the vertical velocity component. As typically observed in a clearsky convective boundary layer, the skewness is positive and increases with height (Moeng and Rotunno, 1990). Near the inflow boundary, the skewness is close to zero or even slightly negative, indicating that the imposed up and downdrafts are equally distributed with respect to their area contribution. Similar to the TKE, the skewness peaks close to the inflow boundary and adjusts towards a constant positive value further downstream where strong/narrow thermal updrafts and weaker/wider downdrafts are in equilibrium. This equilibrium value is reached earlier for the skewness compared to the TKE (see Fig. 14), which is in accordance to the results shown in MuñozEsparza and Kosović (2018). In other words, the flow rapidly develops coherent turbulent structures, albeit these are still too energetic, as indicated by the TKE.
To investigate of how fast land–atmosphere interactions adjust, Fig. 16 shows horizontal profiles of LE_{0}, H_{0}, and u_{∗} depending on the distance to the inflow boundary. At 10:00 and 13:00 UTC, right behind the inflow boundary, H_{0} and u_{∗} increase with increasing distance and approach a nearly constant value after about 1.5–3 km, indicating almost homogeneous fluxes of sensible heat and horizontal momentum. Likewise, LE_{0} approaches a nearly constant value after about 1.5–3 km, though, in contrast to H_{0} and u_{∗}, it decreases slightly behind the inflow boundary. At 16:00 UTC, the fluxes behave similar though it takes a slightly longer distance of about 5–6 km to approach a nearly constant value. Compared to the TKE, the land–atmosphere exchange adjusts faster and does not show a significant dependence on the distance to the inflow boundary, which in turn is a necessary prerequisite that the turbulent flow can adjust. Especially the fact that the fluxes rapidly approach a homogeneous value is important for the parametrization of the Reynoldsstress components. The surface fluxes directly enter the parametrization – see Eqs. (28)–(32) and (34) – so that the imposed synthetic turbulence depends on the domainaveraged surface fluxes including the region near the inflow boundary. However, since the fluxes rapidly approach a homogeneous value, the error made by averaging is not significant for sufficiently large model domains.
Finally, we would like to note that we also simulated different scenarios with higher (2UH) and lower (05UH) wind speeds, as well as with increased (15RS) and reduced (075RS) shortwave solar radiation, in order to test the convective scaling. Even though the peak amplitudes in the horizontal TKE profile are different due to the modified forcing (see Fig. 17), the peak locations coincide with respect to $d{w}_{\ast}/\left({u}_{\mathrm{h}}{z}_{i}\right)$ at the respective height levels. This, in turn, indicates that the required distance needed to allow the flow to spatially develop is not an universal number but scales with ${u}_{\mathrm{h}}{z}_{i}/{w}_{\ast}$ in a convective boundary layer, meaning that with higher wind speeds or less surface heating the required adjustment fetch increases (see Table 4).
Summarized, under convective conditions, the turbulent flow is fully developed within the boundary layer after about $\mathrm{2.0}\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{h}}{z}_{i}/{w}_{\ast}$ and further adjustment effects further downstream are only small. However, we note that the absolute distance required to allow for fully spatially developed turbulence is still on the order of kilometres, meaning that the model domain should be sufficiently large to place the region of interest sufficiently apart from the inflow boundaries.
4.4 Effect of heterogeneous inflow conditions
Figure 18c shows a horizontal crosssection of the vertical velocity component from COSMO within the middle part of the boundary layer at 12:00 UTC. The COSMO simulation shows elongated structures that are mainly orientated along the meanwind direction with up and downdrafts on the order of m s^{−1}. Visually estimated, the wavelength of these structures is ≈0.1–0.15^{∘}, which is on the order of COSMO's horizontal grid spacing of 0.025^{∘} (2.8 km). Previous studies with the WRF model revealed that these kind of up and downdrafts are a numerical artefact rather than a realistic feature of the boundary layer when the boundarylayer depth is within the range of the horizontal grid spacing (Ching et al., 2014; Zhou et al., 2014; Shin and Dudhia, 2016). Even though results from mesoscale WRF simulations are not necessarily transferable one to one to COSMO, we nevertheless assume a similar behaviour here. With a boundarylayer depth of ≈2.4 km at 12:00 UTC (see Fig. 7), the dominant length scales of the flow approach the horizontal grid spacing, meaning that convection can partly be resolved on the COSMO grid. Figure 18a, b show corresponding horizontal crosssections of the vertical velocity component from a PALM simulation with heterogeneous and homogeneous boundary values prescribed, respectively. After some adjustment behind the inflow boundaries (east and south boundaries) where turbulent structures are weaker and appear on smaller scales (see Sect. 4.3), elongated structures orientated along the meanwind direction form in both simulations, with the typical strength of up and downdrafts for a convective boundary layer. In the heterogeneous case, however, the elongated turbulent up and downdrafts appear more clustered with similar wavelength as in the COSMO simulation, while in the homogeneous case the turbulent up and downdrafts are more homogeneously distributed. These clustering of up and downdrafts in the PALM simulation indicates that griddependent flow structures resolved by COSMO propagate into the LES domain and trigger the development of elongated structures in PALM with a similar wavelength that persists throughout the entire model domain. This is in contrast to Mazzaro et al. (2017), who found that with the cellperturbation method griddependent structures do not significantly bias the turbulence behind the turbulenceadjustment region, though they might affect the rate of evolution of turbulence near the inflow boundaries. This might be attributed to the different coupling time steps. In this study, the coupling time step was 1 h, while Mazzaro et al. (2017) used 1 min. Since the underresolved rolllike structures are not necessarily stationary, the horizontal movement can be well considered with a 1 min coupling, while with only 1hourly coupling these imposed temporary structures are present over a longer time interval at the lateral boundaries so that the underresolved structures can effectively propagate into the model domain. As Fig. 18a indicates, this may introduce a location bias to the turbulent flow, so that the PALM solution becomes dependent on the presence of unrealistic flow structures in the mesoscale model. In such a case, modellers may consider using homogeneous boundary conditions which, as Fig. 18b shows, avoid the artificial generation of persisting structures. However, one should be aware that, especially for large LES domains approaching the size of mesoscale structures, largescale gradients vanish and the mean mesoscale boundarylayer development may not necessarily represent local conditions any more.
4.5 Implications near the inflow and outflow boundaries
In this section, we focus on the flow near the inflow and outflow boundaries. Due to different model representations of surface processes and different surface input data, the mesoscale nearsurface wind and temperature profiles can deviate from the one the LES would simulate. As an example, Fig. 19b shows the nearsurface wind profiles taken at the inflow boundary (COSMO solution) and taken from the inner part of the LES model domain at about 20 km downstream of the inflow boundary. Within the lower 200 m, the PALM profile shows a higher wind speed compared COSMO, indicating that the imposed inflow profile is not in equilibrium with the surface in PALM. As a consequence, the horizontal nearsurface flow behind the east and south inflow boundaries is accelerated, causing a mean downdraft to maintain continuity, which can be observed in Fig 19a. Likewise, to maintain continuity, the flow needs to adjust the mesoscale profile at the outflow boundary (which resembles the inflow profile). Here, the horizontal flow needs to be decelerated which causes a mean updraft close to the outflow boundary. We note that these mean up and downdrafts are most pronounced at nighttime, ranging between 0.1–0.4 m s^{−1} in this setup. Even though during daytime these mean up and downdrafts cannot be detected visually in the instantaneous vertical wind speed (see, e.g. Fig. 12), we also found mean up and downdrafts near the in and outflow boundaries but with a lower amplitude on the order of 0.05–0.1 m s^{−1}. Note that these nearboundary adjustment effects cannot be avoided by the massflux correction in Eqs. (1) and (2), which only acts on the global scale and ensures divergencefree boundaries, but are a result of the interconnection between the surface friction and the pressure solver to maintain incompressibility on all considered scales.
4.6 Computational efficiency of the synthetic turbulence generator
With respect to the mesoscale nesting, the generation of synthetic turbulence represents the computationally most expensive part, while setting the boundary conditions and enforcing a divergence flow field is computationally much less expensive. In order to examine the efficiency of the turbulence generator implementation and estimate its computational cost, we have carried out a performance test.
The most expensive part of the turbulence generation is the computation of the filtered random numbers – see Eqs. (25) and (26) – which depend on the turbulent length scales. Therefore, based on the described setup in Sect. 3, we ran a set of idealized simulations where we varied the turbulent length scales. The length scales were set to a vertically constant value up to z=2500 m and to zero above. We performed 100 time integrations with a time step that was held constant at 0.5 s for all simulations. The simulations with different length scales were performed twice: one set where the computation of the filtered random numbers is distributed over all available MPI processes (hereafter referred to as distributed) and a second set where the filtered random numbers are only computed on boundary MPI processes (hereafter referred to as nondistributed), and the rest of the MPI processes were on hold. The number of MPI processes used for this scaling test was n=1296. Figure 20 shows the consumed CPU time by the synthetic turbulence generator for different length scales. As expected, the consumed CPU time increases with increasing length scale. For small ${L}_{i}/\mathrm{\Delta}{x}_{i}$ the consumed CPU time for both, nondistributed and distributed computation, is below 1 % of the CPU time consumed in the time stepping (see red lines), and no significant difference between both cases can be observed (black lines); this also shows that the additional MPI communication required in the distributed case does not deteriorate the computational performance of the turbulence generator for small ${L}_{i}/\mathrm{\Delta}{x}_{i}$. For larger ${L}_{i}/\mathrm{\Delta}{x}_{i}$, the turbulence generator consumes significant portions of the available resources in the nondistributed case with relative contributions of >60 % for length scales of about 2000 m (${L}_{i}/\mathrm{\Delta}{x}_{i}=\mathrm{100}$). In contrast, in the distributed case, the CPU consumption only increases moderately, reaching only up to 10 % of the total CPU time consumed in the time stepping for length scales of about 2000 m. This shows that parallelizing the tasks needed for the synthetic turbulence generation saves significant computational resources.
Finally, we note that for the simulation covering an entire diurnal cycle, the length scales vary significantly with lower values during nighttime and larger values around noon. For these simulations, the turbulence generation consumed about 2.5 % of the CPU time.
In this paper, we presented a mesoscale nesting interface for the PALM 6.0 model system that extends PALM's capabilities to simulate atmospheric boundary layers under evolving synoptic conditions. The mesoscale nesting interface, which currently relies on output of the COSMO model, consists of two components: (i) the preprocessing interpolation tool INIFOR which provides initial and boundary conditions as a netCDF file, and (ii) PALM's internal boundary condition routines which read and process the initial and boundary conditions as well as imposed additional synthetically turbulent fluctuations. We described INIFOR's interpolation methodology in detail, beginning with the relevant model differences between PALM and COSMO, leading to the conceptual steps needed to interpolate COSMO model output onto the PALM grid. Since the interpolated mesoscale boundary conditions are essentially free of turbulent fluctuations, the flow first needs to spatially develop before the turbulent transport of momentum, energy, and water can be analysed. In order to minimize the extent of development zones near the lateral inflow boundaries which the LES model would otherwise require to generate turbulence via shear and convective instabilities by itself (Mirocha et al., 2014), we employed a synthetic turbulence generation method according to Xie and Castro (2008). Using this approach, spatially and temporally correlated fluctuations of all three velocity components are generated based on parameterizations of the Reynolds stresses as well as turbulent length scales and timescales.
We demonstrated the nesting interface and the effectiveness of the synthetic turbulence generation using a semiidealized benchmark case: we simulated a convective boundary layer developing near Berlin, Germany, on a clearsky spring day using initial and boundary conditions derived from DWD's operational COSMODE analysis. For the sake of analysing the spatial development of the flow, the case was idealized in that we assumed flat terrain with homogeneous grassland instead of using more realistic landsurface heterogeneity, in order to disentangle turbulence built up due to the synthetic turbulence generation and convective and shear instabilities from effects of the particular surface heterogeneities of the Berlin area. We found that the flow rapidly develops up and downdrafts, whereas the adjustment of the TKE takes a longer distance of about 2–$\mathrm{3}\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{h}}{z}_{i}/{w}_{\ast}$, meaning that the turbulent flow needs a fetch length that corresponds to at least two eddyturnover times to be fully adjusted. Even though the adjustment distance could be significantly reduced, it is still on the order of several kilometres, which means that significant parts of the computational resources are still required only for the spatial development of the flow. To further reduce the computational effort, an alternative could be to combine the mesoscale nesting together with PALM's selfnesting (Hellsten et al., 2021), i.e. a relatively coarse grid resolution in the outermost parent domain and finer grid resolutions within the nested child domains. Another worthwhile branch of development would be the implementation of the cellperturbation method by MuñozEsparza and Kosović (2018) in PALM as an alternative to the synthetic turbulence generation to profit from recent promising developments in this regard (Mazzaro et al., 2019; MuñozEsparza and Kosović, 2018). For a stationary boundary layer, MuñozEsparza et al. (2015) showed that the cellperturbation method requires shorter fetches compared to synthetically generated turbulence according to Xie and Castro (2008).
In our benchmark case, the boundary layer in the mesoscale COSMO model does not develop synchronously with the boundary layer in the LES. For example, in COSMO, the boundary layer develops more rapidly before noon and the evening transition starts earlier compared to the LES simulation. As the signals due to nonsynchronous boundarylayer development are imposed to the inflow boundary, these propagate through the LES domain creating a horizontally heterogeneous boundary layer during the morning and evening transition phase. Furthermore, we observed underresolved convective rolls emerging in the mesoscale model that, similar to Mazzaro et al.'s (2017) findings, propagate into the LES domain. In the present study, we eliminated these rolllike structures by averaging over the inflow boundary, being aware that especially for larger domains synopticscale horizontal gradients or effects of mesoscale topography cannot be considered (Mazzaro et al., 2019). To eliminate spurious underresolved convection already in the mesoscale WRF simulations, MuñozEsparza et al. (2017) increased the vertical diffusion by increasing the Smagorinsky constant in the turbulence closure and, thus, damping the vertical up and downdrafts without changing the general boundarylayer structure significantly. However, the choice of the Smagorinsky constant may be case dependent and the general applicability of this approach yet needs to be investigated in detail. Another alternative could lie in the filtering the of the boundary conditions using a filter width corresponding to the horizontal grid spacing of the mesoscale model may help to eliminate such spurious flow structures.
Overall, especially the nonsynchronous boundarylayer development and the imposed rolllike convection emphasize that the representation of the boundary layer in the LES and accompanied vertical gradients of wind velocity, potential temperature, etc. depend on the boundarylayer representation in the mesoscale model. Suppose the boundary layer is not well captured in the mesoscale model, e.g. due to misrepresented convection and turbulent mixing, cloud cover, or atmosphere–surface exchange; the boundary layer resolved in the LES will also be affected by this. In such cases, the physically more credible LES solution (with respect to the boundarylayer representation) will be continuously pushed towards the mesoscale solution. Here, further research is required to better understand the causes for such model discrepancies, under which circumstances they arise, and what the implications are for the representation of the turbulent boundarylayer flow.
Further branches of future development will be to enable INIFOR to also process WRF (Skamarock et al., 2008) and ICON (Zängl et al., 2015; Reinert et al., 2020) output, as well as to add further prognostic quantities to the mesoscale nesting interface, e.g. chemical compounds, aerosols, liquid and frozen water. This is especially important to properly consider clouds and precipitation in the LES, which in turn also affect the surface net radiation and thus the entire boundarylayer development. However, we expect that in many future applications with mesoscale nesting the outermost parent domain will only run with relatively coarse grid resolution, so that cloud physics will not be captured well in the LES, especially for highaltitude clouds. Hence, we also plan to enable INIFOR to also provide incoming short and longwave radiation fluxes, which, to date, PALM is already enabled to consider either from observations or manually from observations or mesoscale model output.
The main focus of this study was to demonstrate the capability of the mesoscale nesting approach and to confirm the effectiveness of the synthetic turbulence generation to reduce the fetch length needed for the model to develop balanced turbulence characteristics. Dedicated evaluation runs of the PALM 6.0 model system including the mesoscale nesting interface are currently on their way within the project [UC]^{2} (Scherer et al., 2019).
As opposed to our comment in the PALM 6.0 overview paper (Maronga et al., 2020), the geostrophic wind forcing is not required in the present mesoscale nesting interface, which we think warrants further explanation. In idealized model setups, the assumption of periodicity is often used. However, using periodic boundary conditions prevents the model from developing any mean horizontal pressure gradient. Thus, if largescale pressure gradients are important for a given problem, they need to be externally prescribed. This is often done by using an equivalent geostrophic wind profile that is obtained from the mesoscale pressure and density fields P and ρ, respectively:
which enters the model in the form of the additional forcing tendency
in the horizontal momentum equations. With this external mesoscale forcing, even an atmosphere initially at rest will eventually develop a mean horizontal flow representative of the real conditions. In the case of inflow Dirichlet boundary conditions, the situation is reversed: the dynamic pressure develops a mean horizontal gradient as a result of internal forces in order to maintain continuity under the prescribed inflow boundary conditions.
In incompressible formulations, the pressure solution is obtained by constructing and solving a Poissontype equation which can be obtained by applying the divergence operator to the momentum equation. The equation is simplified by exploiting the incompressible continuity equation which represents a divergence constraint on the mass flux ρv. As a result, the pressure solution of the Poisson equation acts as to enforce this divergence constraint onto the flow field. In the case where Dirichlet boundary conditions are used for the velocity, any divergence resulting from the tendencies in the momentum equation will be compensated by a corresponding gradient in the dynamic pressure solution. For instance, mean friction and mean Coriolis forces will result in pressure gradients opposing those effects.
The PALM model system 6.0 is freely available at http://palmmodel.org (PALM, 2021) and distributed under the GNU General Public License v3 (http://www.gnu.org/copyleft/gpl.html, last access: 13 August 2021). The preprocessor INIFOR is included in the PALM software repository as a utility and is currently available at https://palm.muk.unihannover.de/trac/browser/palm/trunk/UTIL/inifor (Kadasch, 2021). The simulations documented in the present article were performed using revision 4564 of the PALM model system, which includes INIFOR version 1.4.15. A complete archive of the software used for this publication, including the input data used, analysis, and plotting scripts, as well as a stepbystep reproduction guide, is available at https://doi.org/10.25835/0084787 (Kadasch and Sühring, 2020).
EK was responsible for conceptualization and development of INIFOR, writing parts of the manuscript, conceptualization of the study, analysis of the results; MS was responsible for conceptualization and implementation of mesoscale nesting into PALM, adaption of the synthetic turbulence generator, writing parts of the manuscript, conceptualization of the study, analysis of the results; TG was responsible for adaptation of synthetic turbulence generator and writing parts of the manuscript; SR was responsible for conceptualization of INIFOR and mesoscale nesting; all authors were responsible for revision of the manuscript.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work is part of the [UC]^{2} project. Financial support was provided by the German Federal Ministry of Education and Research (BMBF) within the framework of Research for Sustainable Development (FONA; http://www.fona.de, last access: 13 August 2021), which is gratefully acknowledged. All simulations with PALM have been performed on the supercomputers of the North German Supercomputing Alliance (HLRN). We thank Heike SchauNoppel and the three anonymous reviewers for their valuable comments on the manuscript.
This research has been supported by the German Federal Ministry of Education and Research (Bundesministerium für Bildung und Forschung, BMBF) under grant no. 01LP1601.
This paper was edited by Xiaomeng Huang and reviewed by three anonymous referees.
André, J. C., De Moor, G., Lacarrère, P., and du Vachat, R.: Modeling the 24Hour Evolution of the Mean and Turbulent Structures of the Planetary Boundary Layer, J. Atmos. Sci., 35, 1861–1883, https://doi.org/10.1175/15200469(1978)035<1861:MTHEOT>2.0.CO;2, 1978. a
Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M., and Reinhardt, T.: Operational ConvectiveScale Numerical Weather Prediction with the COSMO Model: Description and Sensitivities, Mon. Weather Rev., 139, 3887–3905, https://doi.org/10.1175/MWRD1005013.1, 2011. a, b, c, d
Baldauf, M., Förstner, J., Klink, S., Reinhardt, T., Schraff, C., Seifert, A., and Stephan, K.: Kurze Beschreibung des LokalModells Kürzestfrist COSMODE (LMK) und seiner Datenbanken auf dem Datenserver des DWD, Tech. rep., Deutscher Wetterdienst, available at: https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_de/cosmo_de_dbbeschr_version_2_3_201406.pdf?__blob=publicationFile&v=5 (last access: 13 August 2021), version 2.3, 2014. a, b, c
Baldauf, M., Gebhardt, C., Theis, S., Ritter, B., and Schraff, C.: Beschreibung des operationellen Kürzesfristvorhersagemodells COSMOD2 und COSMOD2EPS und seiner Ausgabe in die Datenbanken des DWD, Tech. rep., Deutscher Wetterdienst, available at: https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_version_1_0_201805.pdf?__blob=publicationFile&v=3 (last access: 13 August 2021), version 1.0, 2018. a, b
Brost, R. A., Wyngaard, J. C., and Lenschow, D. H.: Marine Stratocumulus Layers. Part II: Turbulence Budgets, J. Atmos. Sci., 39, 818–836, https://doi.org/10.1175/15200469(1982)039<0818:MSLPIT>2.0.CO;2, 1982. a, b, c, d, e
Ching, J., Rotunno, R., LeMone, M., Martilli, A., Kosovic, B., Jimenez, P. A., and Dudhia, J.: Convectively Induced Secondary Circulations in FineGrid Mesoscale Numerical Weather Prediction Models, Mon. Weather Rev., 142, 3284–3302, https://doi.org/10.1175/MWRD1300318.1, 2014. a, b
Clough, S. A., Shephard, M. W., Mlawer, E. J., Delamere, J. S., Iacono, M. J., CadyPereira, K., Boukabara, S., and Brown, P. D.: Atmospheric radiative transfer modeling: A summary of the AER codes, Short Communication, J. Quant. Spectrosc. Ra., 91, 233–244, 2005. a
Davies, H. C. and Turner, R. E.: Updating prediction models by dynamical relaxation: an examination of the technique, Q. J. Roy. Meteor. Soc., 103, 225–245, https://doi.org/10.1002/qj.49710343602, 1977. a
Doms, G. and Baldauf, M.: A Description of the Nonhydrostatic Regional COSMOModel, Tech. rep., Deutscher Wetterdienst, available at: http://www.cosmomodel.org/content/model/documentation/core/cosmoDyncsNumcs.pdf (last access: 13 August 2021), cOSMOModel 5.5, 2018. a, b
Emes, M. J., Arjomandi, M., Kelso, R. M., and Ghanadi, F.: Turbulence length scales in a lowroughness nearneutral atmospheric surface layer, J. Turbul., 20, 545–562, https://doi.org/10.1080/14685248.2019.1677908, 2019. a
Flay, R. and Stevenson, D.: Integral length scales in strong winds below 20 m, J. Wind Eng. Ind. Aerod., 28, 21–30, https://doi.org/10.1016/01676105(88)900980, 1988. a
Gehrke, K. F., Sühring, M., and Maronga, B.: Modeling of landsurface interactions in the PALM model system 6.0: Land surface model description, first evaluation, and sensitivity to model parameters, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd2020197, in review, 2020. a
Gronemeier, T., Inagaki, A., Gryschka, M., and Kanda, M.: Largeeddy simulation of an urban canopy using a synthetic turbulence inflow generation method, Journal of Japan Society of Civil Engineers, Ser. B1 (Hydraulic Engineering), 71, I_43–I_48, https://doi.org/10.2208/jscejhe.71.I_43, 2015. a
Gronemeier, T., Raasch, S., and Ng, E.: Effects of Unstable Stratification on Ventilation in Hong Kong, Atmosphere, 8, 168, https://doi.org/10.3390/atmos8090168, 2017. a
Gryning, S., Holtslag, A., Irwin, J., and Sivertsen, B.: Applied dispersion modelling based on meteorological scaling parameters, Atmos. Environ., 21, 79–89, https://doi.org/10.1016/00046981(87)902733, 1987. a, b
Heinze, R., Moseley, C., Böske, L. N., Muppa, S. K., Maurer, V., Raasch, S., and Stevens, B.: Evaluation of largeeddy simulations forced with mesoscale model output for a multiweek period during a measurement campaign, Atmos. Chem. Phys., 17, 7083–7109, https://doi.org/10.5194/acp1770832017, 2017. a, b, c, d, e
Hellsten, A., Ketelsen, K., Sühring, M., Auvinen, M., Maronga, B., Knigge, C., Barmpas, F., Tsegas, G., Moussiopoulos, N., and Raasch, S.: A nested multiscale system implemented in the largeeddy simulation model PALM model system 6.0, Geosci. Model Dev., 14, 3185–3214, https://doi.org/10.5194/gmd1431852021, 2021. a, b, c, d
Holtslag, A. A. M. and Boville, B. A.: Local Versus Nonlocal BoundaryLayer Diffusion in a Global Climate Model, J. Climate, 6, 1825–1842, https://doi.org/10.1175/15200442(1993)006<1825:LVNBLD>2.0.CO;2, 1993. a
Honnert, R., Masson, V., and Couvreux, F.: A Diagnostic for Evaluating the Representation of Turbulence in Atmospheric Models at the Kilometric Scale, J. Atmos. Sci., 68, 3112–3131, https://doi.org/10.1175/JASD11061.1, 2011. a
Jähn, M., MuñozEsparza, D., Chouza, F., Reitebuch, O., Knoth, O., Haarig, M., and Ansmann, A.: Investigations of boundary layer structure, cloud characteristics and vertical mixing of aerosols at Barbados with large eddy simulations, Atmos. Chem. Phys., 16, 651–674, https://doi.org/10.5194/acp166512016, 2016. a
Jiang, P., Wen, Z., Sha, W., and Chen, G.: Interaction between turbulent flow and sea breeze front over urbanlike coast in largeeddy simulation, J. Geophys. Res.Atmos., 122, 5298–5315, https://doi.org/10.1002/2016JD026247, 2017. a
Kadasch, E.: INIFOR [code], available at: https://palm.muk.unihannover.de/trac/browser/palm/trunk/ UTIL/inifor, last access: 13 August 2021. a
Kadasch, E. and Sühring, M.: Supplementary material to “Mesoscale nesting interface of the PALM model system 6.0”, Leibniz Universität Hannover [data set], https://doi.org/10.25835/0084787, 2020. a
Kataoka, H. and Mizuno, M.: Numerical flow computation around aeroelastic 3D square cylinder using inflow turbulence, Wind Struct., 5, 379–392, https://doi.org/10.12989/WAS.2002.5.2_3_4.379, 2002. a, b
Kim, Y., Castro, I. P., and Xie, Z.T.: Divergencefree turbulence inflow conditions for largeeddy simulations with incompressible flow solvers, Comput. Fluids, 84, 56–68, https://doi.org/10.1016/j.compfluid.2013.06.001, 2013. a
Klein, M., Sadiki, A., and Janicka, J.: A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations, J. Comput. Phys., 186, 652–665, https://doi.org/10.1016/S00219991(03)000901, 2003. a
Lee, G.J., MuñozEsparza, D., Yi, C., and Choe, H. J.: Application of the Cell Perturbation Method to LargeEddy Simulations of a Real Urban Area, J. Appl. Meteorol. Clim., 58, 1125–1139, https://doi.org/10.1175/JAMCD180185.1, 2019. a, b
Letzel, M. O., Helmke, C., Ng, E., An, X., Lai, A., and Raasch, S.: LES case study on pedestrian level ventilation in two neighbourhoods in Hong Kong, Meteorol. Z. 21, 575–589, https://doi.org/10.1127/09412948/2012/0356, 2012. a
Li, S., Hu, Z., Chan, P., and Hu, G.: A study on the profile of the turbulence length scale in the nearneutral atmospheric boundary for sea (homogeneous) and hilly land (inhomogeneous) fetches, J. Wind Eng. Ind. Aerod., 168, 200–210, https://doi.org/10.1016/j.jweia.2017.06.008, 2017. a
Lund, T. S., Wu, X., and Squires, K. D.: Generation of Turbulent Inflow Data for SpatiallyDeveloping Boundary Layer Simulations, J. Comput. Phys., 140, 233–258, https://doi.org/10.1006/jcph.1998.5882, 1998. a
Maronga, B. and Raasch, S.: LargeEddy Simulations of Surface Heterogeneity Effects on the Convective Boundary Layer During the LITFASS2003 Experiment, Bound.Lay. Meteorol., 146, 17–44, https://doi.org/10.1007/s105460129748z, 2013. a
Maronga, B., Gryschka, M., Heinze, R., Hoffmann, F., KananiSühring, F., Keck, M., Ketelsen, K., Letzel, M. O., Sühring, M., and Raasch, S.: The Parallelized LargeEddy Simulation Model (PALM) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments, and future perspectives, Geosci. Model Dev., 8, 2515–2551, https://doi.org/10.5194/gmd825152015, 2015. a, b
Maronga, B., Banzhaf, S., Burmeister, C., Esch, T., Forkel, R., Fröhlich, D., Fuka, V., Gehrke, K. F., Geletič, J., Giersch, S., Gronemeier, T., Groß, G., Heldens, W., Hellsten, A., Hoffmann, F., Inagaki, A., Kadasch, E., KananiSühring, F., Ketelsen, K., Khan, B. A., Knigge, C., Knoop, H., Krč, P., Kurppa, M., Maamari, H., Matzarakis, A., Mauder, M., Pallasch, M., Pavlik, D., Pfafferott, J., Resler, J., Rissmann, S., Russo, E., Salim, M., Schrempf, M., Schwenkel, J., Seckmeyer, G., Schubert, S., Sühring, M., von Tils, R., Vollmer, L., Ward, S., Witha, B., Wurps, H., Zeidler, J., and Raasch, S.: Overview of the PALM model system 6.0, Geosci. Model Dev., 13, 1335–1372, https://doi.org/10.5194/gmd1313352020, 2020. a, b, c, d, e, f, g
Mazzaro, L. J., MuñozEsparza, D., Lundquist, J. K., and Linn, R. R.: Nested mesoscaletoLES modeling of the atmospheric boundary layer in the presence of underresolved convective structures, J. Adv. Model. Earth Sy., 9, 1795–1810, https://doi.org/10.1002/2017MS000912, 2017. a, b, c, d, e, f, g
Mazzaro, L. J., Koo, E., MuñozEsparza, D., Lundquist, J. K., and Linn, R. R.: Random Force Perturbations: A New Extension of the Cell Perturbation Method for Turbulence Generation in Multiscale Atmospheric Boundary Layer Simulations, J. Adv. Model. Earth Sy., 11, 2311–2329, https://doi.org/10.1029/2019MS001608, 2019. a, b, c
Mirocha, J., Kosović, B., and Kirkil, G.: Resolved Turbulence Characteristics in LargeEddy Simulations Nested within Mesoscale Simulations Using the Weather Research and Forecasting Model, Mon. Weather Rev., 142, 806–831, https://doi.org/10.1175/MWRD1300064.1, 2014. a, b, c, d
Moeng, C.H. and Rotunno, R.: VerticalVelocity Skewness in the BuoyancyDriven Boundary Layer, J. Atmos. Sci., 47, 1149–1162, https://doi.org/10.1175/15200469(1990)047<1149:VVSITB>2.0.CO;2, 1990. a
Mordant, N., Metz, P., Michel, O., and Pinton, J.F.: Measurement of Lagrangian Velocity in Fully Developed Turbulence, Phys. Rev. Lett., 87, 214501, https://doi.org/10.1103/PhysRevLett.87.214501, 2001. a
Munters, W., Meneveau, C., and Meyers, J.: Shifted periodic boundary conditions for simulations of wallbounded turbulent flows, Phys. Fluids, 28, 025112, https://doi.org/10.1063/1.4941912, 2016. a, b
MuñozEsparza, D. and Kosović, B.: Generation of Inflow Turbulence in LargeEddy Simulations of Nonneutral Atmospheric Boundary Layers with the Cell Perturbation Method, Mon. Weather Rev., 146, 1889–1909, https://doi.org/10.1175/MWRD180077.1, 2018. a, b, c, d, e, f
MuñozEsparza, D., Kosović, B., Mirocha, J., and van Beeck, J.: Bridging the Transition from Mesoscale to Microscale Turbulence in Numerical Weather Prediction Models, Bound.Lay. Meteorol., 153, 409–440, https://doi.org/10.1007/s1054601499569, 2014. a
MuñozEsparza, D., Kosović, B., van Beeck, J., and Mirocha, J.: A stochastic perturbation method to generate inflow turbulence in largeeddy simulation models: Application to neutrally stratified atmospheric boundary layers, Phys. Fluids, 27, 035102, https://doi.org/10.1063/1.4913572, 2015. a, b
MuñozEsparza, D., Lundquist, J. K., Sauer, J. A., Kosović, B., and Linn, R. R.: Coupled mesoscaleLES modeling of a diurnal cycle during the CWEX13 field campaign: From weather to boundarylayer eddies, J. Adv. Model. Earth Sy., 9, 1572–1594, https://doi.org/10.1002/2017MS000960, 2017. a, b, c
PALM: The PALM model system web pages [code], available at: http://palmmodel.org, last access: 13 August 2021. a
Park, S.B., Baik, J.J., and Lee, S.H.: Impacts of Mesoscale Wind on Turbulent Flow and Ventilation in a Densely Builtup Urban Area, J. Appl. Meteorol. Clim., 54, 811–824, https://doi.org/10.1175/JAMCD140044.1, 2015. a
Reinert, D., Prill, F., Frank, H., Denhard, M., Baldauf, M., Schraff, C., Gebhardt, C., Marsigli, C., and Zängl, G.: DWD Database Reference for the Global and Regional ICON and ICONEPS Forecasting System, Tech. rep., Deutscher Wetterdienst, available at: https://www.dwd.de/DWD/forschung/nwv/fepub/icon_database_main.pdf, version 2.1.1, last access: 6 June 2020. a, b
Rotach, M. W., Gryning, S.E., and Tassone, C.: A twodimensional Lagrangian stochastic dispersion model for daytime conditions, Q. J. Roy. Meteor. Soc., 122, 367–389, https://doi.org/10.1002/qj.49712253004, 1996. a, b, c, d, e, f
Salesky, S. T., Katul, G. G., and Chamecki, M.: Buoyancy effects on the integral lengthscales and mean velocity profile in atmospheric surface layer flows, Phys. Fluids, 25, 105101, https://doi.org/10.1063/1.4823747, 2013. a
Schalkwijk, J., Jonker, H. J. J., Siebesma, A. P., and Van Meijgaard, E.: Weather Forecasting Using GPUBased LargeEddy Simulations, B. Am. Meteorol. Soc., 96, 715–723, https://doi.org/10.1175/BAMSD1400114.1, 2015. a, b
Scherer, D., Antretter, F., Bender, S., Cortekar, J., Emeis, S., Fehrenbach, U., Gross, G., Halbig, G., Hasse, J., Maronga, B., Raasch, S., and Scherber, K.: Urban Climate Under Change [UC]2 – A National Research Programme for Developing a BuildingResolving Atmospheric Model for Entire City Regions, Meteorol. Z. 28, 95–104, https://doi.org/10.1127/metz/2019/0913, 2019. a
Shin, H. H. and Dudhia, J.: Evaluation of PBL Parameterizations in WRF at Subkilometer Grid Spacings: Turbulence Statistics in the Dry Convective Boundary Layer, Mon. Weather Rev., 144, 1161–1177, https://doi.org/10.1175/MWRD150208.1, 2016. a, b
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D., Duda, M. G., Huang, X., Wang, W., and Powers, J. G.: A Description of the Advanced Research WRF Version 3 (No. NCAR/TN475+STR), Tech. rep., University Corporation for Atmospheric Research, https://doi.org/10.5065/D68S4MVH, 2008. a, b
Tennekes, H. and Lumley, J.: A first course in turbulence, The MIT Press, Cambridge, Mass., 1972. a
Troen, I. B. and Mahrt, L.: A simple model of the atmospheric boundary layer; sensitivity to surface evaporation, Bound.Lay. Meteorol., 37, 129–148, https://doi.org/10.1007/BF00122760, 1986. a
Wan, H., Giorgetta, M. A., Zängl, G., Restelli, M., Majewski, D., Bonaventura, L., Fröhlich, K., Reinert, D., Rípodas, P., Kornblueh, L., and Förstner, J.: The ICON1.2 hydrostatic atmospheric dynamical core on triangular grids Part 1: Formulation and performance of the baseline version, Geosci. Model Dev., 6, 735–763, https://doi.org/10.5194/gmd67352013, 2013. a
Wicker, L. J. and Skamarock, W. C.: TimeSplitting Methods for Elastic Models Using Forward Time Schemes, Mon. Weather Rev., 130, 2088–2097, https://doi.org/10.1175/15200493(2002)130<2088:TSMFEM>2.0.CO;2, 2002. a
Williamson, J. H.: Lowstorage RungeKutta schemes, J. Comput. Phys., 35, 48–56, https://doi.org/10.1016/00219991(80)900339, 1980. a
Willis, G. E. and Deardorff, J. W.: A laboratory model of diffusion into the convective planetary boundary layer, Q. J. Roy. Meteor. Soc., 102, 427–445, https://doi.org/10.1002/qj.49710243212, 1976. a
Wu, X.: Inflow Turbulence Generation Methods, Annu. Rev. Fluid Mech., 49, 23–49, https://doi.org/10.1146/annurevfluid010816060322, 2017. a
Wyngaard, J. C.: Toward Numerical Modeling in the “Terra Incognita”, J. Atmos. Sci., 61, 1816–1826, https://doi.org/10.1175/15200469(2004)061<1816:TNMITT>2.0.CO;2, 2004. a
Xie, Z. and Castro, I.: Efficient Generation of Inflow Conditions for Large Eddy Simulation of StreetScale Flows, Flow Turbul. Combust., 81, 449–470, https://doi.org/10.1007/s1049400891515, 2008. a, b, c, d, e, f, g, h
Zhong, J., Cai, X., and Xie, Z.T.: Implementation of a synthetic inflow turbulence generator in idealised WRF v3.6.1 large eddy simulations under neutral atmospheric conditions, Geosci. Model Dev., 14, 323–336, https://doi.org/10.5194/gmd143232021, 2021. a
Zhou, B., Simon, J. S., and Chow, F. K.: The Convective Boundary Layer in the Terra Incognita, J. Atmos. Sci., 71, 2545–2563, https://doi.org/10.1175/JASD130356.1, 2014. a, b
Zängl, G., Reinert, D., Rípodas, P., and Baldauf, M.: The ICON (ICOsahedral Nonhydrostatic) modelling framework of DWD and MPIM: Description of the nonhydrostatic dynamical core, Q. J. Roy. Meteor. Soc., 141, 563–579, https://doi.org/10.1002/qj.2378, 2015. a, b
 Abstract
 Introduction
 Mesoscale nesting interface
 Simulation setup and statistical analysis
 Results
 Summary and conclusions
 Appendix A: A note on largescale pressure forcing
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Mesoscale nesting interface
 Simulation setup and statistical analysis
 Results
 Summary and conclusions
 Appendix A: A note on largescale pressure forcing
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References