Steady-State Large-Eddy Simulations of Convective and Stable Urban Boundary Layers

A comprehensive investigation is carried out to establish best practice guidelines for the modelling of statistically steady-state non-neutral urban boundary layers (UBL) using large-eddy simulation (LES). These steady-state simulations enable targeted studies under realistic non-neutral conditions without the complications associated with the inherently transient nature of the UBL. An extensive set of simulations of convective and stable conditions is carried out to determine which simplifications, volumetric forcings, and boundary conditions can be applied to replicate the mean and turbulent (variance and covariance) statistics of this intrinsically transient problem most faithfully. In addition, a new method is introduced in which a transient simulation can be ‘frozen’ into a steady state. It is found that non-neutral simulations have different requirements to their neutral counterparts. In convective conditions, capping the boundary-layer height h with the top of the modelled domain to h/5 and h/10 (which is common practice in neutral simulations) reduces the turbulent kinetic energy by as much as 61% and 44%, respectively. Consistent with the literature, we find that domain heights lz≥5|L|\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l_z \ge 5 |L|$$\end{document} are necessary to reproduce the convective-boundary-layer dynamics, where L is the Obukhov length. In stably stratified situations, the use of a uniform momentum forcing systematically underestimates the mechanical generation of turbulence over the urban canopy layer, and therefore leads to misrepresentations of both the inner- and outer-layer dynamics. The new ‘frozen-transient’ method that is able to maintain a prescribed flow state (including entrainment at the boundary-layer top) is shown to work well in both stable and convective conditions. Guidelines are provided for future studies of the capped and uncapped convective and stable UBL.


Introduction
As cities continue to grow and become more dense, it becomes increasingly important to understand the urban climate, particularly in regard to the sustainability of urban life and global challenges such as urban air quality and the urban-heat-island effect. The ability to model urban environments is integral to developing this understanding. The urban fabric and the anthropogenic processes occurring therein interact with the planetary boundary layer to produce a complex modelling problem over a large range of spatial and temporal scales. High-resolution modelling of the urban canopy layer (UCL) is imperative to resolving the unsteady and heterogeneous urban flow field, rendering large-eddy simulation (LES) a key modelling tool (Tominaga and Stathopoulos 2013;Barlow 2014;Lateb et al. 2016). The parameter space for a realistic urban boundary layer (UBL) with a given urban form (morphology, trees, etc.) is very large, resulting in high computational costs with limited generality. Urban studies must therefore be conducted with the aims of, (1) isolating particular phenomena or processes of interest and (2) accounting for horizontal, vertical, and temporal variability (targeted studies; see Oke et al. 2017). The production of statistical steady states is integral to satisfying both of these aims and allows for converged statistics to be obtained (directly with periodic boundary conditions and when used as driver simulations; see Tomas et al. 2015).
Extensive research has been carried out in order to understand the fundamental physics and processes of the urban flow field under neutral conditions in a statistical steady state (Walton and Cheng 2002;Belcher 2005;Letzel et al. 2008;Xie and Castro 2009). As our understanding progresses, additional levels of complexity have been incorporated into urban LES investigations under the same framework. Much of this work has involved modelling thermal effects, such as atmospheric stability (Xie et al. 2013;Boppana et al. 2014), the use of radiative models (Resler et al. 2017), and the incorporation of trees (Li and Wang 2018;Matsuda et al. 2018). Recent studies have moved to combine all of these additional effects to utilize the LES approach to study realistic integrated urban environments (Kurppa et al. 2018). Hanna (1989) highlighted the likelihood of degrading modelling performance with additional complexity. Blocken (2015) used this assertion to indicate the need for more extensive best practice guidelines for urban LES and this is particularly prevalent in the case of non-neutral conditions.
The fundamental physics of the UBL is altered by the role of buoyancy and an additional degree of freedom is introduced into the modelling problem as a result of the slow transients associated with convective (CBL) and stable (SBL) boundary-layer dynamics. Mesoscale meteorological conditions, radiative transfer, and anthropogenic heat sources combine to produce a complex urban surface energy balance, which drives a diurnal cycle in the UBL. The development of the CBL and SBL during the day and night, respectively, and the transitions between them result in time-dependent boundary-layer heights and degrees of atmospheric stability driven by turbulent entrainment at the UBL top and temperature changes at the surface (see Fig. 1). The UBL tends to be less stable and deeper than its rural counterpart (Roth 2000;Barlow 2014).
Field, experimental and numerical investigations have indicated the importance of atmospheric stability on the urban climate. Numerical studies have illustrated the effects upon, for example, the urban heat island effect (Schrijvers et al. 2015), the form of the roughness sublayer (Barlow 2014) and the dispersion of pollutants (Caton et al. 2003;Salizzoni et al. 2009;Li et al. 2016). Wind-tunnel experiments have shown the effect of local stability on the turbulent structure within the UCL (Uehara et al. 2000;Kanda and Yamao 2016). Numerous field studies have identified atmospheric stability as the prevailing meteorological factor in determining urban air quality (Sánchez-Ccoyllo and De Fátima Andrade 2002;Mayer 1999;Hien et al. 2002;Perrino et al. 2008). The form of the UBL can be characterised by three key length scales: the boundary-layer height h, the Obukhov length L and the average building height H . A range of stability regimes exists under both the convective and stable classifications. The UBL is generally non-neutral and sheared convective conditions typically represent the most prevalent regime (Kotthaus and Grimmond 2018). In terms of the SBL, weak stability with continuous turbulence is the dominant regime over urban areas (Uno et al. 1988;Barlow 2014). Developing the ability to model non-neutral conditions is thus fundamental to understanding the urban environment and facing the challenges it presents.
A key challenge in this regard is to develop a model of a non-neutral UBL that is statistically steady. Previous LES investigations have taken different levels of abstraction in modelling the heat fluxes at the urban surface, from idealized isothermal boundary conditions (Boppana et al. 2014) to full three-dimensional radiative transfer models (Resler et al. 2017). A wide range of forcings and boundary conditions have similarly been applied to enforce steady states on the resulting unstable and stable conditions (summarized in Table 1 and discussed in detail in Sect. 2). Many of the simplifications that are used for neutral conditions have been applied to the non-neutral case with little consideration of their effects on the flow statistics. The use of these different forcings, boundary conditions and domain heights fundamentally alters the resulting mean and turbulent characteristics of the boundary layer. Studies citing the same bulk Richardson numbers or Obukhov lengths can exhibit different turbulent structures as a result of the applied steady-state methodologies. The aim of the present study is to outline the best practice guidelines (specifically the applied forcing terms, boundary conditions and computational domain sizes) for producing a statistically steady-state and physically realistic CBL or SBL using the urban LES approach (using periodic boundary conditions). Much of the presented results and analyses are equally applicable to generic (non-urban) boundary layers and the findings can be applied as such. The urban context is implicit within the approach taken here and the formation of the devised guidelines (e.g. the aim of simplifying the problem and the requirement to have sufficient resolution to capture flows within the UCL; Blocken 2015).
The paper is structured as follows: Sect. 2 introduces the prognostic equations applicable to the UBL and provides a comprehensive review on what forcings and boundary conditions can be applied to create a statistically steady state in non-neutral conditions. A description of the LES model and the suite of simulations carried out as part of this study is provided in Sect. 3. Two transient simulations of the developing CBL and SBL are discussed in Sect. 4.1. Sets of steady-state simulations that aim to reproduce the conditions of the transient CBL and SBL after 10 h using different steady-state methodologies are displayed in Sects. 4.2 and 4.3. The shortcomings of particular methodologies and best practice for conducting non-neutral urban LES investigations are outlined from this analysis. Section 5 discusses a novel methodology in which transient simulations can be 'frozen' in a specific state without modifying their statistical properties. Finally, the best practice guidelines are summarized and conclusions are drawn in Sect. 6.

Requirements for Steady-State Urban Boundary Layers
The real UBL exhibits horizontal, vertical and temporal variation over a range of scales, driven by, for example, meteorological and morphological changes. A steady-state UBL is defined hereafter as a flow that satisfies statistical stationarity and horizontal homogeneity [neglecting local-scale heterogeneity in the roughness sublayer (RSL)]. By definition, the mean and turbulent properties of the boundary layer must therefore be temporally and horizontally invariant given sufficiently long sampling (Pope 2000). The criterion for horizontal homogeneity is satisfied by the use of periodic boundary conditions.

Integral Description of the Urban Boundary Layer
The prognostic boundary-layer equations illustrate the inherent transiency of the problem and allow us to understand how statistical steady states can be achieved in urban LES. Defining the horizontal-averaging operator · and assuming horizontally homogeneous conditions, the equations for the horizontal velocity components and the potential temperature in the UBL are given by where F u and F v are forcing terms applied to model the mesoscale momentum forcing, τ x and τ y are the total shear stresses in the two horizontal directions, F θ denotes a thermal forcing term and φ is the total vertical heat flux. Mesoscale forcings such as subsidence and (m) Fig. 1 Diagrammatic view of a the diurnal cycle of the planetary boundary layer, b a canonical CBL potential temperature profile (and the free atmosphere (FA) aloft) and c a canonical range of SBL potential temperature profiles (with the residual layer (RL) aloft). Black solid arrows indicate temporal development of the potential temperature profiles, dotted grey lines denote fluxes of potential temperature into or out of the boundary layer and grey arrows are possible volumetric thermal-forcing terms (discussed in more detail in Table 1 and Sect. 5) radiative divergence are neglected here and F θ can thus be expected to be equal to zero in Eq. 1c under normal conditions. The role of the heterogeneous urban canopy (form drag) and other meteorological forcings have not been included here but can be added without loss of generality (Raupach and Thom 1981;Xie and Fuka 2018). In Eqs. 1a and 1b, it is evident that the time derivatives/storage terms reduce to zero under the conditions that the magnitude of the local momentum forcing balances the divergence of the local shear stress. However, there is no equivalent forcing term for potential temperature (since F θ ≈ 0 under normal circumstances) meaning that the storage term is directly determined by the local divergence of the vertical heat flux, which can therefore only reach a steady-state in the limit of a uniform heat-flux profile ∂ φ /∂z = 0. This assertion indicates that the mean potential temperature is inherently transient in the presence of a non-zero local divergence in the vertical heat flux, which is expected across a typical CBL and SBL.
Further insight can be obtained by integrating Eqs. 1a-1c across the boundary layer, which results in where the Leibniz integration rule has been used. The hat symbol (·) denotes an integral property of the boundary layer and subscripts h and 0 denote the value at the boundary-layer top (hereafter defined as the top of the inversion) and urban surface, respectively. The terms on the left-hand side of Eqs. 2a-2c demonstrate two distinct ways in which the UBL can be transient: (1) via the bulk change in the boundary layer's integral properties (e.g. bulk cooling or heating), and (2) via the vertical growth of the boundary layer due to entrainment/encroachment (also diagrammatically shown for both the CBL and SBL in Fig. 1). The presence of two time derivatives illustrates the challenges associated with the variabledepth boundary-layer equations. All steady-state methodologies must therefore tackle both mechanisms as well as the local balance determined by the divergence of vertical flux terms.

Steady-State Methods
Investigations into the non-neutral UBL face the same challenges as the neutral case but are complicated by the inherent imbalance exhibited by Eqs. 1c and 2c and the fundamental differences in the turbulent structure and boundary-layer heights of the modelled flows. A range of simplifications, forcings and boundary conditions can be applied to form steady states, but these methods affect the mean and turbulent statistics of the modelled flow. They must be formulated to produce physically realistic UBL dynamics whilst also satisfying the aims and computational demands of a given study. Steady-state methods in this context refer to the volumetric terms and top boundary conditions that can be applied such that the modelled boundary layer reaches a statistical steady state (periodic lateral boundary conditions are assumed in the lateral directions). Both the imposed thermal and momentum forcings are important. Thermal forcings, both as defined within the literature and those presented here, are summarized in Table 1 and discussed in detail in Sects. 2.2.2 and 2.2.3. The momentum forcings are presented in Sect. 2.2.1 and in Table 2. In the context of running a steady-state urban LES model, it is important to note herein the use of both the spatial and temporal (·) averaging operators (following the ergodic principle; Monin and Yaglom 2013;Ghannam et al. 2015). Averaging over the horizontal extent of the domain alone is often insufficient to fully resolve the scales of turbulent motion.

Momentum Forcing
The momentum forcings F u and F v require consideration as they directly affect the shear stresses τ x and τ y . Four momentum forcings are used as part of this study, labeled a to d, which are discussed below and are detailed in Table 2. Real UBL momentum forcings arise from mean mesoscale pressure gradients and the Coriolis effect, ∂ y where f is the Coriolis parameter (equivalent to forcing a). This results in an Ekman spiral through the UBL (under stable and neutral conditions), meaning that the direction and magnitude of the momentum forcing changes with height. In neutral urban LES studies, the Ekman layer and non-uniformity is typically circumvented by neglecting the Coriolis force and therefore simplifying the forcing term to a uniform pressure gradient F u = F θ = 0, φ l z = 0 (Unforced) Flow is capped by the zero-flux condition at the top of the domain. The mean temperature profile develops in time but the turbulent heat flux reaches a steady state assuming self-similar behaviour (Saiki et al. 2000;Li and Wang 2018) (i)

Boundary conditions
Boundary layer capped by a heat flux at the top of the domain (equal to flux from urban surface). Results in a uniform heat-flux profile over domain height (Boppana et al. 2014) (ii) θ l z = const. (Isothermal) Boundary layer capped by an isothermal top boundary condition. Equivalent to steady-state constant-flux condition (Li et al. 2010(Li et al. , 2016Chan and Liu 2013) (ii)

Volumetric source/sink
Volumetric term applied to balance the bulk heating/cooling of the boundary layer. Results in a linear turbulent heat-flux profile up to top of the domain (Boppana et al. 2014;Tomas et al. 2016;Shen et al. 2017) As above, except the term R is redefined at each timestep by nudging towards a predefined bulk temperature. Consideration required on the relaxation time scale t N (iii) Volumetric term that permits non-linear turbulent heat-flux profiles. Must satisfyF θ = φ l z − φ 0 . Can be calculated from a transient simulation (e.g., Δ θ /Δt) or set a priori to produce a predefined heat-flux profile (iv) Volumetric term defined by nudging towards predefined temperature profile, θ i . Consideration required for the magnitude of t N and translation of temperature profile that arises. If applied in top cells, uniform heat flux below the nudged layer (Matsuda et al. 2018) (v)

Volumetric source/sink and stabilizing velocity
van Driel and Jonker (2011) utilized the self-similarity of the CBL and the flux-ratio method to assign uniform values of w s and R that balance the bulk heating and vertical growth of the CBL and enforce an assigned lapse rate Γ in the free atmosphere. Developed to work in the limit of free convection, no similar closure method or self-similarity for application to the SBL Indices are repeated for methods with equivalent time-averaged effects. Periodic lateral boundary conditions are assumed and any bottom boundary condition is possible given that some statistically steady heat flux φ 0 exists The volumetric source/sink term R and the stabilizing velocity w s are shown diagrammatically in Fig. 1 Table 2 A summary of the momentum forcings (F u and F v ) Momentum forcing defined as a function of the mesoscale pressure gradient and Coriolis force via geostrophic velocity components (u g and v g ) and the Coriolis parameter, f (a) Uniform and unidirectional forcing defined by idealized mesoscale pressure gradients (using a reference friction velocity u * ; neglects Coriolis force) Momentum forcing is defined using the absolute geostrophic deficit of a transient simulation. This acts as an approximation of the vertical variability of the momentum forcing whilst neglecting its rotationality (c) Momentum forcing is predefined following the stationary solution of Nieuwstadt (1984) (using a reference u * ). Results in non-linear total-vertical-momentum-flux profile Coceal et al. 2006;Xie et al. 2008). This simplification allows for the wind direction and friction velocity/volume flow rate (depending on its formulation) to be user-defined. The effect of changing from an Ekman layer to a unidirectional flow under neutral conditions has been shown to have a minimal effect within the inner layer of the boundary layer but to cause significant changes to the turbulent structure in the outer layer due to the role of vortex tilting (Ansorge and Mellado 2014). The uniform forcing therefore only holds for analysis of the inner layer of the neutral UBL. In fact, it is not possible to model the free atmosphere/residual layer above the UBL in a steady state with a uniform forcing, as the flow accelerates uniformly in the absence of the Coriolis force. The modelled boundary layer is, therefore, capped by the domain top in a steady state by definition when using a uniform pressure gradient. Enforcing a constant pressure gradient also results in a linear shear-stress profile. This forcing is therefore consistent with for example the mixed layer of the CBL. However, the shear-stress profiles of the SBL have been shown to be far more challenging to simplify or parametrize in such a way (Mahrt 1999).
Forcings c and d similarly neglect the rotationality of the flow through the UBL and require the boundary layer to be capped, but they provide the capability to model non-linear shear-stress profiles. Forcing c defines F u and F v using the magnitude of the geostrophic deficit of a transient simulation. Forcing d utilizes the solution of Nieuwstadt (1984) for the homogeneous and stationary SBL. It defines the momentum forcing such that the expected 3/2 power-law shear-stress profile is enforced (Glazunov 2014).

Capped Urban Boundary Layers
The majority of neutral urban LES studies use a free-slip boundary condition at the top of the domain (τ l z = 0 where l z is the domain height ) to truncate the modelled boundary layer such that dh/dt = 0 by construction and the boundary-layer height is equal to the vertical extent of the domain (l z = h; hereafter referred to as a capped UBL). The conditions for the integral properties of the boundary layer to reach a steady state subsequently simplify toF u = τ x 0 andF v = τ y 0 , which has the additional advantage of being computationally efficient: the boundary-layer height can be predefined and set specifically for the problem of interest. The quiescent layer above the boundary layer, which is often not of interest in targeted studies, is neglected. Domain heights are typically ≈100-200 m, assuming that it is sufficient to model the surface layer alone (≈ h/10). Capping the boundary layer in this manner is similarly used in methods i-v in Table 1, although it is unclear whether these assumptions hold under non-neutral conditions.
Capping the boundary layer and applying a uniform forcing have been routinely used for both convective and stable conditions. Li and Wang (2018) applied this condition with a zeroflux top boundary condition for heat and no thermal forcing (φ l z = 0 and F θ = 0; method i). These conditions do not equilibrate the inherent transiency of Eqs. 1c and 2c, allowing the mean temperature profile to change in time. Self-similar behaviour was assumed for the turbulence profiles once the temperature profile was increasing uniformly in time. By definition, the steady-state criteria are not satisfied by this method.
Following Eq. 2c, if F θ = 0, the only way to produce a steady state is to ensure that φ 0 = φ h . Constant-flux or isothermal boundary conditions at the top of the domain can be applied to satisfy this condition (method ii; Li et al. 2010;Boppana et al. 2014;Chan and Liu 2013). These conditions result in a uniform vertical heat-flux profile over the domain, which is consistent with the established surface-layer assumption (h/10). This description holds for the CBL but, by definition, the height of the SBL surface layer of is of similar order of magnitude to the building height. The changes in the scaling of turbulence under nonneutral conditions also challenge the justification for modelling the surface layer alone (also discussed in Sect. 2.2.1). For example, Inagaki et al. (2012) highlight that coherent structures in the convective mixed layer significantly impact the nature of flow within the UCL and Fodor et al. (2019) illustrate that large-scale updrafts and downdrafts cause deviations from Monin-Obukhov similarity theory within the freely convective surface layer and quantify this effect for different top boundary conditions. The applicability of this assumption is analyzed in detail in Sects 4.2 and 4.3.
Introducing a thermal-forcing term F θ = 0 into Eqs. 1c and 2c provides the capability to produce steady states that feature non-uniform heat-flux profiles. The term F θ must be defined such that it does not affect the turbulent dynamics of the modelled flow. However, this criterion is generally satisfied due to the slow transients associated with the CBL and SBL dynamics (consideration is required for time-dependant terms). Continuing under the assumption of a capped boundary layer, the condition for a steady state in Eq. 2c becomeŝ F θ = φ l z − φ 0 and the local balance in Eq. 1c becomes F θ = ∂ φ ∂z (equivalent to the relationship between momentum forcing and shear stress). Under these conditions, a volumetric heat source/sink term, F θ = constant, can be defined to balance the bulk and local cooling/heating of the boundary layer. The magnitude and vertical distribution of this forcing term directly impacts the resulting heat-flux profile. Three methods are defined in Table 1 that follow this methodology (methods iii-v). Boppana et al. (2014) and Shen et al. (2017) imposed a constant, uniform profile of the forcing term F θ , which results in a linear heat-flux profile (method iii). This profile is consistent with the expected linear heat flux within the convective mixed layer and also with the stationary solution defined by Nieuwstadt (1984) for the SBL. In both of these studies, a zero-flux top boundary condition φ l z = 0 was applied for heat, analogous to modelling up to the top of the UBL. However, the domain heights were considerably smaller than the typical height of the CBL. It is possible to apply this method with non-zero values of φ l z in order to simulate, for example, the role of entrainment at the top of the CBL or to enforce profiles of the turbulent heat flux that are consistent with the CBL height but modelled in a smaller domain. If the surface heat flux is not known a priori, it is possible to define an equivalent thermal forcing by applying a time-dependent bulk-nudging term F θ = −(θ − θ i )/t N (where θ i is a predefined temperature profile and t N is the relaxation time scale).
For typical meteorological conditions over cities, particularly for the SBL, heat-flux profiles are non-linear. The ability to produce non-linear turbulent flux profiles using a non-uniform forcing can therefore be a useful tool (methods iv and v). Following Eq. 1c, the turbulent heat flux can be defined as φ = F θ dz + C, with applied top and bottom boundary conditions (method iv; following a similar methodology to forcing d and still The term F θ can also be defined such that it is a function of both height and time. The mean temperature profile can be nudged towards a predefined profile by using Matsuda et al. (2018) applied nudging of this form over the top cells in the domain to produce steady-state unstable conditions. If applied over a limited number of cells, a uniform heat flux is present up to the lowest nudged cell and the nudged cells act to form a buffer layer at the top of the domain (therefore acting equivalently to method ii but defining a buffer region between the bulk boundary-layer flow and the top of the domain). For the CBL, a steady-state mean temperature profile can be obtained, but this forcing would behave fundamentally differently under stable conditions. Time-dependent thermal forcings are useful for application in simulations where the surface heat flux is not known a priori (e.g. with isothermal boundary conditions or for full three-dimensional radiative-transfer models) or when the aim is to reproduce, for example, a specific mean temperature profile. These forcings follow the form of 'nudging' or data assimilation that has been applied extensively in transient studies (Pleim 2007). In the context of the present work, the forcing serves to nudge the potential temperature towards a steady state as opposed to towards some transient profile as is typically done. In a steady state, the nudging term, as defined in Table 1, results in a translation of the mean temperature profile as a function of the heat-flux divergence and the relaxation time scale (Cai et al. 1997;Lakshmivarahan and Lewis 2013). The relaxation time scale must also be sufficiently large to ensure that the nudging term does not directly affect the turbulent state of the flow.

Uncapped Urban Boundary Layers
Uncapped UBL simulations extend vertically into the free troposphere, well away from the top of the UBL, which implies that both the boundary layer and its interaction with the free atmosphere/residual layer aloft are included. Sorbjan (1996) studied the impact of modelling a 'solid-lid' non-penetrative CBL, and identified a significant difference between the turbulent structure of modelled boundary layers where entrainment through the inversion layer is captured and where the boundary layer is capped by a solid and adiabatic wall. If entrainment processes are to be incorporated faithfully in a steady-state formalism, the UBL must not be capped by the top of the domain.
van Driel and Jonker (2011) derived a convective fixed-point solution to create steady states in uncapped simulations of the CBL. The vertical growth of the boundary layer was balanced by the addition of a stabilizing velocity, w s , in the thermal forcing term (F θ = R − w s ∂ θ ∂z ; method vi). The stabilizing velocity is analogous to a subsidence velocity and provides a forcing component that is dependent on the local thermal stratification. The magnitude of the uniform profiles of the terms R and w s can be defined through the convective fixedpoint solution that applies the flux-ratio method and mixed-layer model to Eq. 2c (see Table  1; Ball 1960;Tennekes 1973). The thermal forcings are also applied above the boundary layer where they act to further maintain the defined lapse rate Γ following Γ = −R/w s . This method is designed for use with free, as opposed to mixed or forced, convection (the entrainment ratio A will change for a sheared CBL; Liu et al. 2018;Haghshenas and Mellado 2019) and there is no similar closure method for use with the SBL. A novel 'frozen-transient' method is developed in Sect. 5 where non-uniform profiles of the radiative transfer R and stabilizing velocity w s are defined such that the method is suitable for both convective and stable stratifications (method vii).

Large-Eddy Simulation Set-Up
Full-scale simulations are performed with an urban version of the Dutch Large-Eddy Simulation model (Heus et al. 2010) called the uDALES model (Tomas et al. 2015(Tomas et al. , 2016Grylls et al. 2019), whose governing equations are where ... denotes low-pass filtering at the scale of the grid, Π is the modified pressure, g is the acceleration due to gravity, δ i3 is the Kroenecker delta, T is the deviatoric component of the subgrid-scale (SGS) momentum flux, and R is the SGS flux vector (the filtered notation is omitted for the remainder of the paper for clarity; Heus et al. 2010). These filtered Navier-Stokes equations are solved using an explicit third-order Runge-Kutta time integration scheme and a second-order central-differencing advection scheme. The SGS model used follows Vreman (2004). The immersed boundary method has been incorporated into the uDALES model in order to provide the capability to model buildings within the domain (Tomas et al. 2015). Logarithmic wall functions are used to resolve the SGS dynamics close to the walls (negating the necessity to refine the grid structure close to the modelled buildings but requiring a parametrization to approximate the near-wall dynamics; Uno et al. 1995;Suter 2018). The numerical scheme consists of a spatial discretization on an Arakawa C-grid, utilizing a second-order central-differencing scheme for all field variables and a third-order Runge-Kutta time-integration scheme.
Simulations of both convective (-C) and stable (-S) conditions are performed to evaluate the different steady-state methodologies (simulation details are provided in Table 3). Two transient simulations (tC and tS), applying a geostrophic momentum forcing and no volumetric thermal forcing, are used to model the development of a typical CBL and SBL over a 10-h time period. A suite of steady-state simulations (s-) are set up with the aim of reproducing the state of their corresponding transient simulation after 10 h of run time (hereafter referred to as the reference state). Simulations: • sC1-3 apply established simplifications from the literature to model the assumed surface layer of the CBL. The boundary layers are capped at the heights h/10 and h/5, a uniform momentum forcing is applied, and methods ii and iii are used to enforce uniform or linear vertical heat-flux profiles, respectively. • sC4-6 similarly cap the boundary layer and apply a uniform momentum forcing. In these simulations, the top heat flux is defined such that the vertical heat-flux profile matches that of the reference state and the domain height is varied to investigate its effect on the modelled flows (l z = h/5, h/2 and h).
Indices for the momentum forcings are outlined in Table 2 and indices for the thermal forcing refer to those in Table 1 a Using rotation rate of the Earth Ω = 7.292 × 10 −5 rad s −1 and latitude ϕ = 51 • N b The surface heat flux was scaled from 5φ 0 to φ 0 over the first 1800 s of simulation time in order to avoid enhanced boundary-layer development resulting from the initial conditions • sC7-8 do not cap the boundary layer and correspondingly use a geostrophic momentum forcing. Methods vi and vii are used to enforce a steady state. • sS1-2 apply methodologies as found in the literature to model the SBL. The boundary layers are capped, a uniform momentum forcing is applied and methods ii and iii are used to enforce a uniform or linear vertical heat-flux profile, respectively. • sS3-5 cap the boundary layer and apply different combinations of non-uniform flow and thermal forcings. The simulation sS3 satisfies the solution of Nieuwstadt (1984) by applying a uniform thermal forcing (method iii) and a 3/2 power-law momentum forcing (forcing d). Simulation sS4 uses the statistics of tS to obtain non-uniform profiles of both the thermal and momentum forcing (method vi and forcing c). Simulation sS5 applies the 3/2 power-law momentum forcing and method v to nudge the mean temperature profile towards the temperature profile of the reference state. • sS6 does not cap the boundary layer and correspondingly uses a geostrophic momentum forcing. Method vii is applied to enforce a steady-state.
The surface heat flux (φ 0,C = 0.04 and φ 0,S = − 0.005 K m s −1 ), friction velocity (u * ,C = 0.45 and u * ,S = 0.16 m s −1 ) and therefore the Obukhov length (L C = −163 m and L S = 63 m) are matched to that of the corresponding reference state in all convective and stable steady-state simulations. In simulations where a unidirectional forcing is applied, the wind direction is defined to match that of the reference state at building height (α C | z=H = 28 • and α S | z=H = 51 • ). The boundary-layer heights h of the reference states, which are used to determine the domain heights l z of the capped simulations, are h C = 1234 m and h S = 120 A base-simulation set-up is defined such that the other simulation conditions remain constant (including the urban morphology) whilst the forcings, top thermal boundary condition, and domain configurations are varied as discussed. The base set-up applies periodic horizontal boundary conditions for both momentum and heat and a free-slip condition at the top of the domain. A staggered cubic array is used across all simulations, diagrammatically shown in Fig. 2, with a building height H = 16 m (therefore also matching the building-height length scale H ). Constant heat fluxes are defined at the road and roof-top levels with no heat flux defined through the vertical building walls. The horizontal sizes of the computational domains are defined to be sufficiently large with respect to their respective vertical domain heights such that one can average over various large coherent structures. It is important to note that best practice in the urban LES approach is to require a minimum of ≈10 grid cells per building height H to resolve flows within the UCL (Blocken 2015). This criterion is not satisfied in some simulations, which highlights the advantage of capping the boundary layer at lower heights in urban LES studies (this criterion is discussed in more detail in Sect. 4.2).
The initial profiles of horizontal velocity components and potential temperature also vary by simulation. Simulation tC applies an initial lapse rate, Γ = − 0.003 K m −1 , in order to represent the free atmosphere, whereas simulation tS uses a uniform initial temperature profile of 288 K to represent the residual layer. Both transient simulations utilize uniform profiles of the specified geostrophic velocity components (u g = 4 m s −1 , v g = 0) as initial conditions. Steady-state simulations that apply method vii (sC8 and sS6; see Sect. 5 for details) are run directly from the state of the corresponding transient simulation after 10 h run time. Following van Driel and Jonker (2011), simulation sC7 sets a uniform initial potential temperature profile up to the defined boundary-layer height h C , an increase in potential temperature of the inversion strength, Δθ = 0.78 K, at h C , and an increase in potential temperature at the defined lapse rate, −Γ , for z > h C . The other steady-state convective simulations (sC1-6) utilize a uniform initial potential-temperature profile of 288 K. Simulations sC1-7 apply uniform initial velocity profiles that represent the geostrophic velocity adjusted to the wind direction at building height as discussed. The remaining steadystate stable simulations (sS1-5) use profiles of the mean potential temperature and the wind speed (adjusting the wind angle) from the reference state. Using the values from the reference state helped to reduce the required simulation times in the presence of large time scales in the outer layer.
Statistics are obtained by spatially averaging in the horizontal directions (a mask is applied in the positions of the buildings; c.f. intrinsic spatial averaging; Xie and Fuka 2018) and temporally averaging with a sampling time of 5 s. The transient simulations are averaged over a time period of 1 h and the steady-state simulations are averaged over a time period sufficient to produce converged statistics. For the convective simulations, a maximum averaging period of ≈ 20T C is used (T C = h C / gh C φ 0 /θ C 1/3 = 1041 s) and for the SBL, an averaging period of ≈ 10T S (T S = h S /u * ,S = 750 s) is used across all simulations. The simulations are run until a negligible change is observed between the resulting mean and turbulent statistics over consecutive averaging periods.

Results
The two transient simulations are discussed first to provide a reference state for a typical urban CBL and SBL. Mean and turbulent profiles of the steady-state simulations are then plotted alongside those of the reference state in order to evaluate the range of applied steady-state methodologies. Deviations from the reference state in these profiles indicate a discrepancy in the representation of the modelled flow. The boundary conditions, forcings and simplifications applied to produce these steady states are reviewed and guidelines are outlined for best practice in modelling the non-neutral steady-state UBL in future studies. Much of the following analysis is equally applicable to generic boundary layers in a non-urban context, however the discussion and guidelines are based around modelling approaches for urban LES models since this is where most applications originate.

Transient Non-Neutral Urban Boundary Layers
Simulations tC and tS capture the development of a typical CBL and SBL over an urban area throughout the day and night, respectively. The mean and turbulent statistics of the two developing non-neutral UBL simulations are shown in Fig. 3. The CBL height grows in time, achieving a height h C = 1234 m after 10 h simulation time. This growth is driven by the formation of an inversion through which warm air is entrained down from the stratified free atmosphere aloft and encroachment due to the surface heating. The mean temperature and turbulent heat-flux profiles exhibit self-similar behaviour consistent with classical mixedlayer theory, as has been found in urban field studies (Wood et al. 2010). The development of the velocity and turbulent momentum-flux profiles is consistent with a weakly sheared CBL flow (Fedorovich and Conzemius 2008). The flow veers within the RSL (26.6 • ) and entrainment zone (27.8 • ) and is relatively uniform within the mixed layer. Park and Baik (2014) provide detailed analysis of a similar urban CBL.
The SBL grows to a quasi-equilibrium height, h S ≈ 120 m, within the first 4 h of simulation time, which is an order of magnitude smaller than the convective case. An inversion forms and strengthens, limiting turbulent entrainment through the boundary-layer top. As the inversion strengthens, the cooling is redistributed across the boundary layer, and a nocturnal jet forms towards the top of the velocity profile. The wind direction rotates down through the boundary layer, with a maximum of α S = 77 • within the UCL after 10 h simulation time. The temperature and heat-flux profiles do not exhibit self-similar behaviour. Classifying or applying canonical results and similarity theory to the SBL has been shown to be challenging  (Mahrt 1999;Barlow et al. 2011). The development of simulation tS is typical of a weak SBL characterised by continuous turbulence and an elevated inversion. Approaching 10 h simulation time, the turbulent profiles are relatively consistent with the stationary solution of Nieuwstadt (1984).
Three key length scales, the boundary-layer height h, Obukhov length L, and building height H , are prevalent in UBL dynamics. These two transient simulations provide typical reference cases of sheared convective and weakly stable conditions over urban areas (Kotthaus and Grimmond 2018). In future investigations, it would be useful to investigate more combinations of the length scales h, L and H ; however applying the steady-state methods to these reference states covers the predominant cases of interest and also provides a platform for future works.

Steady-State Convective Urban Boundary Layers
The development of simulation tC highlights several points of interest that require consideration when forming a steady-state urban CBL: (1) the large boundary-layer height and corresponding turbulent length scale, h C H , (2) the significant vertical growth of the boundary layer via encroachment and turbulent entrainment at the boundary-layer top, dh C /dt, and (3) the applicability of mixed-layer theory and therefore the expected uniform mean profiles and linear turbulent-flux profiles.
The behaviour of capped CBL simulations, in which only the surface layer is considered, is explored first. Simulations sC1-3 apply methodologies as seen in the literature to model the CBL up to heights of h C /10 or h C /5. The resulting mean and turbulent profiles are shown in Fig. 4 (where · denotes the fluctuating component of a variable and therefore Fig. 4b, e display the resolved turbulent components of the vertical temperature and momentum fluxes, respectively). The application of a uniform momentum forcing and enforcement of either uniform (method ii) or linear (method iii) vertical heat-flux profiles are self-evident from the turbulent-flux profiles. By definition, the limited domain heights and these forcings lead to deviations from the reference state in the turbulent-flux profiles. The use of a non-zero top thermal boundary condition in simulations sC1 and sC2 results in a thin layer at the top of the domain with an unstable stratification and a corresponding local increase in the potential temperature variance. In simulation sC1, this deviation is significant from as low as 60 m. Simulation sC3 does not exhibit the same thermal interaction with the domain top. However, the reduced turbulent heat flux leads to a corresponding reduction in the potentialtemperature variance (bulk reduction of 35% in comparison with the reference state). The increased potential-temperature variances at the surface are a result of the differences in the modelled resolutions (Sullivan and Patton 2011). The mean wind-speed profiles are similar to the reference state although the steady-state simulations exhibit a reduction in the mean wind speed above the RSL (maximum difference in wind speed of 24% at 114 m in simulation sC1; the RSL is the layer up to which the dispersive component of the total momentum flux is non-zero).
The most pronounced deviation from the reference state is the reduced magnitude of turbulence kinetic energy (TKE; e) at building height and above in all three simulations in Fig. 4f. A maximum reduction of 61%, 44% and 72% is exhibited in simulations sC1-3, respectively. The magnitude of TKE is indicative of the intensity of the local turbulent motion and therefore, for example, the turbulent mixing (it is the defining parameter in terms of the vertical exchange of pollutants out of the UCL; Salizzoni et al. 2011). A reduction in TKE thus signifies a fundamental change in the turbulent structure of these flows in comparison with the reference state, and is therefore a significant potential source of error in urban LES studies under unstable conditions. The TKE reductions are a result of the limited domain heights. The reduced turbulent-flux profiles directly affect the mechanical (sC1-3) and buoyant (sC3) production of turbulence above the RSL, but these profiles alone cannot account for the nature of the TKE reductions shown in Fig. 4 (e.g. at building height). A more fundamental difference in the turbulent structure of the modelled flows exists that affects both the dissipation and transport of TKE within these boundary layers. Typically, within a mixed convective surface layer, large-scale eddy structures (e.g. temperature-ramp structures; Wilczak 1984) form and, within and just above this layer, these structures merge to produce the convective thermals that drive the formation of the mixed layer aloft. The relationship between the boundary-layer height and the Obukhov length and their effect on TKE is discussed in Lenschow (1974). Within the convective mixed layer, there is near uniform spectral behaviour of the velocity components with height. However, underneath the mixed layer, low-frequency spectral components do not show such uniformity and exhibit a scaling with the dimensionless height z/L (Kaimal et al. 1976). The limited domain heights in simulations sC1-3 inhibit the formation of the convective mixed layer, which evidently affects the turbulent structure of the modelled flow.
The domain heights of these simulations are equal to −0.8L C and −1.6L C , indicating the prominence of mechanical over buoyant effects in these boundary layers. The scaling regions for the CBL dictate that the parameter h/|L| ≥ 5 is the criterion to drive the CBL in a convective state (with the surface layer then defined at −L/2; Deardorff 1974; Holtslag and Nieuwstadt 1986). The importance of capturing the large-scale structures of the convective mixed layers and their interaction with the surface layer has previously been highlighted in other LES studies (Mason 1989;Inagaki et al. 2012). These results clearly indicate that modelling unstable conditions by significantly inhibiting the boundary-layer height results in flows that do not capture the turbulent dynamics of a realistic CBL.
Having established that it is not sufficient to resolve the surface layer only for a CBL, the following set of simulations use larger domain heights and more realistic thermal forcings. Figure 5 shows the results for simulations sC4-8. Simulations sC4-6 were set up similar to simulation sC3, however method iii is applied such that the turbulent heat-flux profile follows that of the reference state (φ l z = φ C | z=l z ) and the domain height is varied l z = h C /5, l z = h C /2 and l z = h C . Simulations sC7-8 do not cap the boundary layer and these results are discussed in Sect. 5.
The mean, turbulent-flux and potential-temperature-variance profiles of simulations sC4-6 follow the trends identified in the analysis of simulations sC1-3. Simulation sC4 matches the turbulent heat-flux profile up to the height h C /5 but exhibits the same systematic underestimation of the TKE in Fig. 5f due to the impeded size of the boundary layer. Although simulation sC5 uses a limited domain height, l z = h C /2, it does not exhibit the same reduction over the surface layer of the modelled flow. The TKE shows a good agreement up to a height of 270 m, above which a small reduction is evident. This profile indicates that modelling up to a height of h C /2 and 4.5L C is sufficient to capture the convective nature of the boundary layer in this case, and sufficient to realistically model the turbulent dynamics over the UCL. Following this finding and the assertion in Holtslag and Nieuwstadt (1986) that the criterion h/L ≥ 5 is required to drive the convective state of the boundary layer, the scaling l z ≥ 5|L| provides a guideline through which to avoid systematic underestimations in the TKE in future urban LES studies. The applicability of this scaling for other values of the stability parameter h/L is investigated in the Appendix where the same criterion is found to work in the limit of free convection (provided the minimum Obukhov length is used; Businger 1973). A comprehensive study is needed to further investigate this scaling for different h/L and to determine the influence of the building height H (Lenschow 1974).
Simulation sC6 applied the same simplifications with a domain height equal to that of the reference state (l z = h C ). The agreement shown between simulation sC6 (as well as sC5) and the reference state in Fig. 5 illustrates that simplifying the CBL by neglecting the Coriolis force and by capping the boundary layer is reasonable (given that the domain height is sufficiently large). These simplifications are possible due to the well-mixed nature of the CBL. The changes in turbulent structure resulting from 'solid-lid' non-penetrative convection (Sorbjan 1996) are not evident within the potential-temperature variance or TKE profiles of the simulations sC5 and sC6. The application of a negative heat flux at the top of the domain in simulation sC6 [which was not imposed by Sorbjan (1996)] acts to minimize these differences by modelling the entrainment of heat down from the free atmosphere.
The challenge for urban LES investigations is that a trade-off exists between the physical size of the domain and the computational feasibility of a given simulation. High resolutions are required to effectively resolve the flow within the UCL (1-2 m; Tseng et al. 2006;Xie and Castro 2009). However, as has been highlighted, relatively large domain heights are necessary to effectively model the CBL. This criterion for the resolution was not satisfied in the simulations tC and sC6-8, but was obtained in the vertical direction in simulation sC5. The discussed set of steady-state convective simulations encompasses the range of setups that could be applied to model convective conditions in targeted urban LES studies. It is shown that capping the boundary layer and using a one-dimensional momentum forcing are reasonable simplifications, particularly if a study is focused on the inner layer of the UBL. However, a key finding is that the vertical extent of the domain must be equal or greater than five times the Obukhov length (l z ≥ 5|L|) to avoid a systematic underestimation of the TKE. In studies of, for example, a single infinite canyon, it may seem superfluous to model up to heights of 600 m and above; however, capturing the large-scale structures of the convective mixed layer is fundamental to understanding the urban climate under convective conditions.

Steady-State Stable Urban Boundary Layers
The SBL provides a different set of challenges in terms of modelling realistic steady-state conditions in an urban LES model. The difficulty in defining any generalized state of the SBL is alleviated by the use of a reference state that represents a typical weak SBL (Mahrt 1999). The modelled conditions also avoid the simulation difficulties associated with runaway cooling and strongly stable conditions (Van de Wiel et al. 2012). The challenges arising from the form of the reference state are, (1) the relatively small boundary-layer height (h S ≈ 8H ), (2) the non-uniform mean and non-linear turbulent-flux profiles, (3) the rotation of the flow down through the boundary layer, and (4) the long time scales associated with the outer-layer dynamics. Figure 6 shows time series of the bulk wind speeds from all of the steady-state simulations. Simulation sS2 exhibits global intermittency over a time period of ≈ 40T S which exceeds the defined sampling period (10T S ) and therefore does not satisfy the steady-state criterion (Mahrt 1999;Ansorge and Mellado 2014). The intermittency in simulation sS2 arises due to the use of a uniform momentum forcing. The momentum forcing of the reference state is a function of the geostrophic deficit, which is largest at the urban surface and reduces to zero at the boundary-layer top. Despite providing the same friction velocity, the uniform momentum forcing systematically underestimates the forcing within the UCL and therefore directly affects the local mechanical generation of turbulence. In this case, the underestimation permits the cooling induced by the surface heat flux to periodically dominate the flow dynamics. The use of this well-established simplification is therefore fundamentally limited because it has resulted in dynamics that were not present in the reference state. It should be noted that for other SBL states (e.g. different Obukhov lengths), the same change in forcing may not instigate global intermittency; however, it similarly underestimates the mechanical generation of turbulence and, therefore, leads to a different turbulent state of the boundary layer.
Simulation sS1 also applied a uniform momentum forcing but alongside method ii, which acts to enforce a uniform heat-flux profile. The bulk wind speed exhibited in Fig. 6 is 40% higher than that of the reference state's bulk wind speed, indicating a fundamental problem with this steady-state methodology. The effect of this thermal forcing is shown more clearly in Fig. 7, which exhibits the mean and turbulent profiles of the steady-state stable simulations (excluding simulation sS2). The enforced vertical heat flux in simulation sS1 leads to a lower and stronger inversion over which high wind speeds and potential-temperature variances arise. The fundamental problem of enforcing a uniform heat flux in an SBL is evident from the turbulent heat-flux profile in Fig. 7b as the flux in the reference state reduces to zero at the boundary-layer top. A further problem in simulation sS1 is indicated by the fact that the normalized turbulent heat flux is significantly less than unity above the RSL, which is due to the prevalence of the subgrid stress term.
Simulations sS3-6 exhibit only negligible changes in bulk wind speed over time in Fig. 6 and also show far better agreement with the reference state in Fig. 7. The increased magnitude of the bulk wind speed in simulations sS3-5 is due to the neglect of the wind veer (absence of the Coriolis force). Ansorge and Mellado (2014) highlight key differences in the outer-layer dynamics of neutral and stable boundary layers with and without the Ekman spiral (due to the role of vortex tilting; Ansorge and Mellado 2014, also discussed in Sect. 2.2.1). The increase in the wind speed in the outer layer of the modelled flows can be clearly seen in Fig. 7, where there is also a small overestimation in the turbulent momentum flux within the RSL that correspondingly leads to greater TKE in this region. However, the momentum forcings applied in these simulations (forcings c and d) generally show a good agreement with the turbulent momentum flux and TKE profiles of the reference state. In terms of conducting Fig. 6 Time series of the bulk wind speed of the steady-state stable simulations sS1-6 in comparison with the reference state t S| t=10h . Simulations sS1-2 apply simplifications as seen in the literature. Simulations sS3-5 use different steady-state methods to capture non-linear fluxes. Simulation sS6 models the uncapped steady-state SBL targeted urban LES investigations they present a significant improvement from the uniform momentum forcing applied in simulations sS1 and sS2 and in particular the agreement with the 3/2 power law is an important finding as it can be defined a priori and provides the numerical control to specify the wind direction and friction velocity of a given simulation.
Simulations sS3-6 all cap the boundary layer with a zero-heat-flux top boundary condition and therefore do not capture the boundary layer's interaction with the residual layer aloft. As a result of this, these simulations do not capture the local increase in the potentialtemperature variance over the inversion in Fig. 7c. Apart from this deviation, the applied thermal forcings (methods iii-v) all reproduce the mean temperature and turbulent heat-flux profiles of the reference state well. Simulation sS3 simply uses a uniform thermal forcing, which alongside the 3/2 power law, enforces flux profiles that follow the solution of Nieuwstadt (1984). This agreement may not be surprising due to the quasi-stationary state of the boundary layer at the reference state, but the ability to apply these simplified forcings to reproduce the turbulent statistics of a realistic reference SBL flow in a steady state is significant. The difference between the performance of simulations sS2 and sS3, which only differ by the applied momentum forcing, further indicates the importance of applying the 3/2 power law forcing here. Simulations sS4 and sS5 illustrate successful applications of nonuniform thermal forcings (derived from simulation tS and through nudging, respectively). The resulting non-uniform thermal forcings capture the non-linear details of the vertical heat flux, providing a small improvement in terms of the mean temperature profile in comparison with simulation sS3. The application of the nudging term alongside the 3/2 power law illustrates a methodology that can be useful in studies attempting to reproduce a specific mean temperature profile, for example, from observation.
The typical heights of the SBL (≈ 100-200 m) indicate that neglecting the interaction with the residual layer and the deviations in the outer layer dynamics is significant in some urban . 7 Mean and turbulent profiles of the steady-state stable simulations sS1 and sS3-6 in comparison with the reference state t S| t=10h . Simulation sS1 applies simplifications as seen in the literature. Simulations sS3-5 use different steady-state methods to capture non-linear fluxes. Simulation sS6 models an uncapped steady-state SBL studies. Prevalent current examples of this are studies of the role of tall buildings in the urban climate (building heights may even exceed the SBL height). Simulation sS6 applies method vii to capture the SBL flow in its entirety as well as the residual layer aloft and therefore presents a method that can be beneficial to these studies. This method and these results are discussed in detail in Sect. 5. For targeted LES studies where the dynamics at the top of the SBL are not of interest and a simplified representation of the SBL flow is advantageous, the most significant finding in this section is the limitation of using a uniform momentum forcing. The solution to this problem is the application of the 3/2 power-law momentum forcing. In terms of the thermal forcing, three methods were presented that achieved reasonable results in this context by applying volumetric source terms. The uniform thermal forcing of method iii is the simplest and easiest to implement in an urban LES model.

Steady-State Uncapped Non-Neutral Urban Boundary Layers
The discussion above focused on simulations where the boundary layer is capped by the top of the domain. It is of interest to form steady-states where the interaction with the quiescent layer above the UBL is captured. In terms of the CBL, it was shown by Sorbjan (1996) that the turbulent structure of the boundary layer is different when modelling penetrative (uncapped) as opposed to 'solid-lid' non-penetrative (capped) convection. The main differences are found to be in the upper half of the flow where the warmer air entrained down into the boundary layer alters the form of the updrafts and downdrafts [see Sorbjan (1996) for details]. A recent study by Fodor et al. (2019) also investigates the impact of different top boundary conditions (capped and uncapped) on the CBL (in the limit of free convection) and similarly identifies differences in the upper half of the flow (the surface layer is shown to be unaffected unless downdrafts are sufficiently cold and strong). It was shown in Sect. 4.3 that the outer-layer flow is significantly misrepresented when capping the SBL and that typical SBL heights are of similar magnitude to tall buildings. Modelling both the uncapped CBL and SBL in steady states can also be useful to study developing/transitional flows where the steady profiles can be used at the inlet and the development of the boundary layer including their interaction with the free atmosphere/residual layer aloft can be studied.

Fixed-Point Solution Method
The fixed-point solution of van Driel and Jonker (2011), discussed in Sect. 2.2 and in Table 1 (method vi), outlines a method by which an uncapped freely convective UBL can be modelled in a steady state. The advantage of this method is that the magnitude of the stabilizing velocity and volumetric source term (F θ = R − w s ∂ θ ∂z ; shown diagrammatically in Fig. 1) can be defined to produce a mean temperature profile with a predefined boundary-layer height and inversion strength a priori. However, this method is limited due to its employment of the entrainment ratio as this term changes under sheared-convective conditions (Liu et al. 2018;Haghshenas and Mellado 2019) and because no such simplification and closure scheme exists in the stable case. Method vi was applied to reproduce the mean temperature profile of the convective reference state. The resulting profiles of the applied thermal forcings are shown in Fig. 8a, b. The steady-state mean and turbulent profiles from the corresponding simulation (sC7) are shown in Fig. 5. A steady-state CBL flow is achieved with good agreement to the reference state despite the sheared nature of the reference CBL flow. A small increase in the boundary-layer height is shown and the velocity profile exhibits a maximum value within the entrainment zone.

Frozen-Transient Method
A novel method to freeze the state of the atmosphere in a specific statistical state is presented below (method vii), which works in both stable and convective situations, unlike method vi, and only requires information about the average potential-temperature profile of a transient simulation and its time derivative. Following van Driel and Jonker (2011), this method uses a thermal forcing consisting of both a stabilizing velocity w s , and a volumetric source/sink term R. However, instead of relying on a parametrization of the UBL to obtain estimates for the parameters w s and R, the frozen-transient method relies on the fact that the normal velocity component v n of an isosurface of a scalar χ is given by (Dopazo et al. 2007;Holzner and Lüthi 2011;van Reeuwijk and Holzner 2014) where v n = v · n is the normal velocity component of the isosurface and n = ∇χ/|∇χ| is the normal vector to the isosurface. For the case of a transient UBL flow, the aim is to create a steady-state case using information of the potential temperature θ (z, t) only. Setting χ = θ and using the homogeneity of the problem, Eq. 5 simplifies to where the relation w n = v n has been introduced to emphasize that the isolines of potential temperature move in the vertical direction only. If a statistically stationary steady state is required, then in principle it suffices to set the stabilizing velocity to w s = −w n such that the resulting volumetric thermal forcing term, w s ∂ θ ∂z , balances the temporal development of the mean potential temperature profile. However, care must be taken in the values for w n that are obtained: (1) in growing boundary layers, we require that w s < 0 (it will only be negative in the convective mixed layer), and (2) large values of w s must be avoided to ensure that it does not modify the turbulent statistics. We thus define a maximum allowable stabilizing velocity w n,max = w * , where w * is a turbulent velocity scale (defined as the Deardorff velocity scale for the convective case and equal to the friction velocity for the stable case), and = 0.1 is used here. The stabilizing velocity over the boundary-layer depth 0 < z < h can then be defined according to In the absence of radiative forcing or subsidence, the vertical velocity component of the potential-temperature isosurface is zero above the boundary layer by definition. However, maintaining a negative stabilizing velocity above the boundary layer (z > h) acts as a secondary system to impede the vertical growth of the boundary layer. Applying the stabilizing velocity above the boundary layer is particularly important in sheared flows where the effects of the thermal forcing must also act to equilibrate the development of the velocity profile. The most robust steady states are achieved by linearly increasing the stabilizing velocity from the boundary-layer top to the domain top. The conditions enforced in Eq. 7 and the possibility of non-zero lapse rates Γ above the boundary layer (as in simulation tC) necessitate an additional thermal-forcing term. A volumetric term R(z) can be defined that, (1) absorbs the remainder of the unsteadiness (arising due to the conditions in Eq. 7; w s = −w n ), and (2) enforces the lapse rate above the boundary layer (Γ = − R w s ; also used in van Driel and Jonker (2011) and equivalent to the radiative equilibrium with a constant large-scale divergence). This term can be obtained by substituting ∂ θ ∂t = −F θ in Eq. 6 giving Practically, there is a challenge associated with obtaining the vertical velocity component from the potential-temperature isosurface: the averaging period must be large enough that the effects of the smallest scale eddies are averaged out but small enough that the vertical development of the boundary layer does not significantly alter the time-averaged statistics. Maximizing the horizontal extent of the domain of the transient simulation is advantageous in minimizing the associated averaging errors. An averaging period of 1 h worked for both the stable and convective transient simulations reported here.
The mean temperature profiles of simulations tC and tS were used to 'freeze' the state of these simulations in time at the reference state following this methodology. Figure 8 shows the calculated velocity component of the potential-temperature isosurface w n , the stabilizing velocity w s and the volumetric source/sink profiles R which are applied in simulations sC8 (Fig. 8a, b) and sS6 (Fig. 8c, d). In the convective case, the defined criteria (Eqs. 7 and 8) partition the thermal forcing between • a relatively uniform volumetric cooling (that shows agreement with the R profile obtained using method vi; shown in Fig. 8b) • a stabilizing velocity w s that sharply decreases over the entrainment zone with a lower peak than that obtained from method vi, which is partly due to the use of an entrainment ratio of −0.25 in method vi (varies under sheared-convective conditions; Liu et al. 2018;Haghshenas and Mellado 2019) and partly due to the discussed limitations of the applied temporal averaging) • a linear decrease in both the parameters w s and R for z > h that forces the flow in the free atmosphere to the predefined lapse rate, Γ = − 0.003 K m −1 .
Note that the parameter w n reaches large negative values within the convective mixed layer in Fig. 8a, which is due to ∂ θ /∂z tending to zero within the convective mixed layer; this highlights the importance of using a maximum value w n,max (Eq. 7). In the convective case, the parameters w s and R therefore act to balance the two growth mechanisms identified in Eq. 2c independently. The parameters w s and R have fundamentally different profiles in the SBL. Figure 8c,d illustrates that the stabilizing velocity is used to balance the development of the potentialtemperature profile in this case. The fluid is stratified throughout the domain in simulation tS, and in contrast with the CBL, the criterion for w n,max does not need to be invoked (w s = −w n ); this in turn implies that R = 0 (Fig. 8d). The stabilizing velocity w s is largest towards the centre of the SBL where the stratification is weakest, and decreases linearly up through the residual layer aloft.
The resulting steady-state mean and turbulent profiles of simulations sC8 and sS6 are shown in Figs. 5 and 7, with both simulations displaying good agreement with the reference state; thus, the method has successfully 'frozen' a developing CBL and SBL. The set-up of these simulations means that the vertical distribution and directionality of the momentum forcing (inclusive of the Coriolis force) and the process of entrainment have therefore been captured in a steady state. Whilst this method demands a precursor transient simulation, a large domain size for the CBL and does not provide control over, for example, the wind direction, it is the only method through which the convective and stable UBL can be modelled faithfully in a steady state. The use of a precursor simulation can be advantageous in the stable case where predefining, for example, the boundary-layer height and turbulent-flux profiles is not trivial. For the CBL, this method enables the modelling of penetrative convection in a steady state and accounts for the effect of shear in the entrainment zone (which may limit the application of the van Driel and Jonker (2011) solution under more strongly sheared conditions). A growing area of interest in urban LES studies is the effect of tall buildings and this method provides the capability to study their interaction with the SBL in the specific case where the building height exceeds that of the domain height.

Conclusion
While thermal effects play an integral role in the urban environment, their incorporation into urban LES studies is complicated by corresponding changes in the state of the overlying UBL (fundamental differences in its mean and turbulent structure and the introduction of slow transients). Modelling a non-neutral UBL in a statistical steady state is key to obtaining converged statistics and to reducing the degrees of freedom of a given problem. Steady-state methodologies are therefore necessary to balance this inherent transiency. Previous investigations have typically applied the same simplifications, forcings and boundary conditions as for the well-established neutral case (e.g. the surface-layer assumption). Few studies have considered the altered demands arising from the fundamentally different forms of the convective and stable UBL and subsequently there is a lack of coherence on how best to produce realistic steady states. A comprehensive investigation was therefore conducted to determine best practice guidelines (specifically the applied forcing terms, boundary conditions and computational domain sizes; Blocken 2015) for the modelling of realistic steady-state convective and stable conditions in urban LES investigations. The results and findings are equally applicable to generic (non-urban) boundary layers but the guidelines are formulated in the context of urban-modelling approaches and best practices (e.g., sufficient resolution within the UCL) since this is where most applications are.
The boundary conditions and forcings that can be applied to enforce a steady state were reviewed in Sect. 2 and Table 1. The resulting steady-state methodologies were contextualized via an integral description of the UBL and the formulation of a generalized thermal forcing term F θ . A defining characteristic of these methodologies was the difference between those that cap the top of the boundary layer (ensuring dh/dt = 0 by definition) and those that model the interaction of the UBL with the quiescent layer aloft. A novel frozen-transient method was developed to model the uncapped steady-state UBL under both convective and stable conditions. Guidelines: apply either method vi (sets mean potential-temperature profile a priori but consideration is required on the entrainment ratio under sheared conditions) or method vii (use driver simulation and ensure that the averaging period does not significantly affect the parameter w n over the entrainment zone)

Capped
Challenges: difficult to categorize or define a typical SBL state over urban areas and the non-uniformity of the momentum forcing and wind direction Results: (1) Neglecting the veering of the flow results in a significant overprediction of wind speed in the outer layer, (2) the application of a uniform momentum forcing leads to an underestimation of the mechanical generation of turbulence over urban roughness and resulted in a fundamental misrepresentation of the modelled flow in simulation sS2 (global intermittency) and (3) several methods of varying complexity were shown to successfully reproduce non-linearities in the turbulent-flux profiles, providing good agreement within the inner layer (methods iii-v) and forcings c and d Guidelines: simplest method that also provides numerical control is application of a uniform thermal-forcing term (method iii) and a 3/2 power-law forcing (both follow the solution of Nieuwstadt (1984) for the stationary weak SBL) Uncapped Challenges: defining any generalized state of the SBL and parametrizing the growth of the SBL (e.g. determining boundary-layer height/entrainment flux a priori) Results: only the frozen-transient method, outlined as part of this study, was capable of capturing a steady-state uncapped SBL Guidelines: apply method vii to a transient developing SBL, consideration is required on the surface heat flux and geostrophic velocity to avoid intermittency/runaway cooling at the surface Two transient simulations of the developing CBL and SBL over an idealized urban morphology provided a realistic reference state. The mean and turbulent profiles of a series of steady-state simulations (set up to reproduce the conditions of the reference state) were then analyzed to assess the advantages and shortcomings of a range of steady-state methodologies. Table 4 outlines the main results from this analysis and the guidelines for best practice in future studies.
The results indicated that capping the UBL is suitable in studies where the dynamics within the inner layer or RSL are primarily of interest. This methodology enables user control over, for example, the boundary-layer height, wind direction and friction velocity, and can be set up to minimize the computational demands of a given simulation. However, care needs to be taken in terms of the domain height and the applied flow and thermal forcing in order to obtain the correct mean and turbulent statistics (simulations citing the same Obukhov length or Richardson number could produce different flow structures; see Table 4 for details). In studies where capturing the UBL in its entirety or the its interaction with the quiescent layer aloft is of importance, uncapped steady states are necessary.
The methods reviewed and developed here are based on statistical steady states in domains with periodic boundary conditions. However, other applications may demand an inflowoutflow configuration. In these cases, a steady-state precursor simulation may be executed to obtain a fully realistic turbulent flow field, of which one plane is then used as the inflow boundary condition for the actual simulation (thereby avoiding having to resort to synthetic turbulence at the inflow boundary; Tomas et al. 2015). A pertinent example of this would be the application of the method outlined in Sect. 5 to study the interaction between tall buildings and the SBL e.g. low-level inversions. A steady-state driver simulation using this method would enable this problem to be studied with a realistic SBL and even if the building height exceeded the UBL depth.
The methods developed here are also applicable to achieve steady states in other situations. Indeed, the issue of transient statistics is prevalent for any scalar in a periodic domain, including moisture and pollutants. As has been shown, consideration is necessary to define a meaningful steady state (e.g. realistic vertical profiles). However, with this in mind, these methods outline pathways with which to achieve the required steady state in a straightforward manner. Table 5 Properties of the transient (t-) and steady-state (s-), freely convective (-FC) simulations. Indices for the momentum forcing are presented in Table 2 and for indices of the thermal forcing refer to those in Table 1 Simulation Domain size Grid points φ 0 φ l z , (θ l z )

Appendix: Dependence on h/L
In Sect. 4.2, it was shown that capping the CBL height to 0.8|L| and 1.6|L| results in a systematic underestimation in the TKE profile and that a domain height, l z ≥ 5|L| is required to model the CBL in a convective state and allow the mixed layer to form. The reference simulation only provides an example of this condition for a single stability parameter h/L. It is important to consider whether this scaling still applies for different values of L, in particular shifts towards the near-neutral and free-convection regimes.
In the limit of neutral conditions, it is well established that the turbulence returns to a dominant scaling with the building height, H , as the heat flux φ 0 tends towards zero. The TKE can be derived as a function of the friction velocity u * , irrespective of the boundary layer height h (given that the value of h is sufficiently large to avoid blockage effects, i.e. h H ; Blocken 2015). Further simulations were performed to explore the sensitivity of the TKE reduction found in Sect. 4.2 for lower values of |L|. A transient convective simulation was performed (tFC) that is similar to simulation tC but with a higher surface heat flux, φ 0 = 0.12 km s −1 , and lower geostrophic wind speed u g = 1 m s −1 (details provided in Table 5). These conditions result in a smaller value of |L| in comparison with simulation tC, and, therefore, in a regime shift from sheared towards free convection that fundamentally affects the dynamics of the inner layer of the CBL. Simulation tFC develops more quickly than simulation tC and as such the reference profile is taken after 5 h simulation time. The reference boundary-layer height h FC is 1582 m, the friction velocity u * ,FC is 0.21 m s −1 , the Obukhov length L FC is −5.4 m, and the wind angle at the building height α FC | z=H is 23.8 • .
For the free CBL, convection-induced stress leads to an additional shear production of TKE that is not related to the mean flow (Businger 1973). Under these conditions, a minimum friction velocity, u * ,min = w * f (h/z 0 ), can be defined that captures this effect where w * is the Deardorff velocity scale and z 0 is the roughness length. Following the definition of u * ,min in Schumann (1988) for rough surface layers, we obtain u * ,min,FC = 0.24 m s −1 . This term is larger than the value of u * ,FC and therefore can be substituted to define the minimum Obukhov length L min,FC = −55.7 m (Businger 1973).
Four steady-state simulations (sFC1-4) were performed following the same methodology as outlined in Sect. 3 (matching the parameters H , φ 0 , u * and α| z=H to simulation tFC; details provided in Table 5). The uniform pressure gradient momentum forcing (b) is used across all simulations and the steady-state methods ii and iii are applied to set linear vertical heat-flux profiles in a steady state. Simulations sFC1-4 cap the height of the modelled boundary layer The vertical profiles of the mean potential temperature, turbulent heat flux and TKE are shown in Fig. 9. Figure 9c illustrates that simulations sFC2-4 do not systematically underestimate the TKE but that simulation sFC1 underestimates the mean TKE by 45% above the UCL. This finding shows agreement with the results in Sect. 4.2 that a domain height l z ≥ 5|L| is required to realistically model the CBL without underestimating the TKE.
Further work is required to fully explore these conditions in the context of the three key length scales h, L and H ; however, the results obtained from both a sheared and free CBL suggest that the inequaity l z ≥ 5|L| where L = L min for the freely convective regime is a robust criterion in the context of typical urban LES conditions.