Open boundary conditions for atmospheric large-eddy simulations and their implementation in DALES4.4

Open boundary conditions were developed for atmospheric large-eddy simulation (LES) models and implemented into the Dutch Atmospheric Large-Eddy Simulation model. The implementation was tested in a “Big Brother”-like setup, in which the simulation with open boundary conditions was forced by an identical control simulation with periodic boundary conditions. The results show that the open boundary implementation has minimal influence on the solution. Both the mean state and the turbulent structures are close to the control simulation, and disturbances at the in- and outflow boundaries are negligible. To emulate a setup in which the LES is coupled to a coarser model, the influence of coarse boundary input was tested by smoothing the output of the periodic control simulation both temporally and spatially before feeding it as input to the simulation with open boundary conditions. When smoothing is applied over larger spatial and longer temporal scales, disturbances start to form at the inflow boundary and an area exists where turbulence needs to develop. Adding synthetic turbulence to the smoothed input reduces the size of this area and the magnitude of the disturbances.

1 Introduction

Large-eddy simulation (LES) is a numerical simulation tool used to study turbulent motions in the atmospheric boundary layer (ABL). Employing resolutions ranging from 1–100 m, the largest turbulent eddies containing most turbulent kinetic energy (TKE) are resolved, whereas the effects of smaller unresolved eddies are parameterized. With most of the TKE being resolved, LES has the advantage over coarser limited-area models (LAMs) when it comes to representing the effects of boundary layer turbulence. This advantage comes at the cost of domain size and/or simulation time. Idealized ABL studies using LES started in the late 1960s–early 1970s (e.g. Lilly, 1966; Deardorff, 1972; Sommeria, 1976). Traditionally, LES was mainly used to study ABLs with idealized homogeneous forcings, employing periodic lateral boundary conditions (LBCs).

With the increase in computational power, the use of LES has shifted from idealized cases to more complex and realistic scenarios. Some examples are the simulation of urban areas (e.g. Giometto et al., 2016; Kurppa et al., 2018), wind farms (e.g. Mehta et al., 2014) and very large case studies pushing towards domain sizes of 𝒪(1000 km) (e.g. Schalkwijk et al., 2015; Heinze et al., 2017). For the latter it is especially important to capture the heterogeneity present in the domain. Periodic LBCs are by definition not suited for this (Moeng et al., 2007). Furthermore, it is often desired to couple LES to a regional weather model to transfer large-scale atmospheric structures. For these reasons, it is desirable to have open LBCs in place. Ideally, open LBCs allow the prescription of variables at inflow boundaries and propagate variables unperturbed out of the domain at outflow boundaries. Having the ability to use open BCs makes an LES model much more versatile in simulating a range of phenomena, especially over heterogeneous terrain. While periodic boundary conditions (BCs) can sometimes be used to study large-scale phenomena over such terrain, the large domains required to do so quickly become computationally prohibitive.

There is no consensus on the “best” implementation of open boundary conditions for anelastic turbulent flow. In 1991 two mini symposia were unsuccessfully dedicated to this topic, and the effort was summarized as a frustrating one (Sani and Gresho, 1994). A popular choice is an outflow condition based on the radiation condition of Sommerfeld (1949). The radiation BC states that waves generated in the interior of the domain should propagate outwards with no reflections at the boundaries. It takes the form of a propagating wave and replaces the Navier–Stokes equations at the boundaries. The difficulty lies in determining the phase speed of the wave, which is required for applying the radiation boundary condition. Different implementations for the phase speed have been defined (e.g. Orlanski, 1976; Klemp and Wilhelmson, 1978; Hedley and Yau, 1988). Orlanski (1976) uses a variable phase speed that is defined upwind of the boundary and propagated to the boundary. The results of a 2D test case show that the implementation works well and results in minimal reflections. Klemp and Wilhelmson (1978) use radiation LBCs in their 3D storm model and evaluate their influence in a 2D version of the model. They define their phase speed as a constant plus the local boundary-normal velocity component. Using a similar test setup to Orlanski (1976), Klemp and Wilhelmson (1978) show that their implementation is capable of producing realistic results. They do note that the results are sensitive to the choice of the fixed part of the phase speed. Hedley and Yau (1988) compare the implementations of Orlanski (1976) and Klemp and Wilhelmson (1978) with their new implementation, which is a hybrid version of the implementation of Orlanski (1976). They conclude that their hybrid implementation is superior to both. Craske and Van Reeuwijk (2013) give a summary of open BCs for incompressible turbulent flows and state that a radiation outflow condition results in the least amount of distortion for convection-dominated flows. Incompressible LES models such as PALM (Maronga et al., 2015, 2020) and Meso-NH (Lac et al., 2018) and the fully compressible WRF-LES (Skamarock et al., 2021) have the option to use radiation boundary conditions for the boundary-normal velocity components based on one of the previously mentioned implementations. For the other variables, (homogeneous) Neumann BCs, which specify the boundary-normal derivative, are often used at outflow boundaries. For inflow conditions, Dirichlet(-like) boundary conditions, which specify the fields at the boundary, are common. The implementation of open boundary conditions in these LES models is summarized in Table 1.

Orlanski (1976)Maronga et al. (2015)Maronga et al. (2020)Kadasch et al. (2021)Carpenter (1982)Lafore et al. (1998)Lac et al. (2018)Klemp and Wilhelmson (1978)Skamarock et al. (2021)Heinze et al. (2017)

Table 1 Summary of the open boundary implementations in the mentioned LES models.

The implementation of open LBCs make it possible to nest LES within both itself and mesoscale models (e.g. Moeng et al., 2007; Zhu et al., 2010; Talbot et al., 2012; Mazzaro et al., 2017; Heinze et al., 2017; Kadasch et al., 2021; Mirocha et al., 2014). These studies use prescribed boundary conditions instead of radiation BCs to nest their LES. Prescribed boundary conditions are a Dirichlet boundary condition, where the LBCs of the child simulation are directly prescribed by the parent simulation. These types of prescribed LBCs are similar to what is used in the mesoscale modelling community and are intuitive to implement. Dirichlet boundary conditions are, however, known to create reflections and perturbations at outflow boundaries for turbulent flows (e.g. Wesseling, 2009; Ol'shanskii and Staroverov, 2000). For this reason Moeng et al. (2007), Zhu et al. (2010) and Heinze et al. (2017) use a relaxation zone in combination with a prescribed boundary condition, in which the fields near the boundary are nudged towards the boundary values to dampen any numerical noise due to the LBCs. Moeng et al. (2007) use WRF-LES to conduct two-way nested simulations of LES nested within LES. They conclude that the nesting works well for LES within LES but state the challenges that will arise for both one-way and a two-way nesting of LES within a mesoscale model. Mesoscale models will have different vertical profiles due to their turbulent transport parameterizations in the planetary boundary layer (PBL) as opposed to the 3D resolved turbulence of LES. Furthermore, since mesoscale models are non-turbulence resolving, the lack of turbulence at inflow boundaries will result in a spinup area required for turbulence to develop. The spinup area is further increased by the implementation of a relaxation zone, as it dampens not only numerical artefacts but also turbulence. Zhu et al. (2010) tested both a one-way and a two-way nested setup with WRF-LES being forced by National Centers for Environmental Prediction (NCEP) reanalysis data. They found that the relaxation zone in the outermost model is able to mitigate potential problems with the large resolution jump between the coarsest WRF-LES domain and the NCEP data set. However, they found that the cloud fields can be strongly modulated by mesoscale organization, especially in high-wind conditions where the clouds align with the mean wind direction. They found little benefit of two-way nesting over one-way nesting. Talbot et al. (2012) coupled WRF-LES within WRF in a realistic one-way nesting setup using three LES domains with increasing resolution in three mesoscale domains. They found that the use of a nested LES setup mainly improves the surface fluxes and near-surface fields, but the bulk ABL dynamics such as the boundary layer height of the mesoscale models agreed better with observations. They also found that the initial and boundary forcings were most important for the results and had a much bigger influence than the choice of subgrid scheme. Mazzaro et al. (2017) studied the effect of unresolved mesoscale flows on LES. They forced three similar LES domains with different-resolution mesoscale simulations. They test their results both with and without the addition of the cell perturbation method of Muñoz-Esparza et al. (2014, 2015) and find that the LES is capable of overcoming erroneous features in the mesoscale output. The cell perturbation scheme helps to greatly reduce the distance required for turbulence to develop, especially for the coarser mesoscale forcings. The best results are obtained with the highest-resolution mesoscale model. Heinze et al. (2017) used a one-way nesting approach to employ realistic LES over Germany. Three domains were used to step down from 625 m horizontal resolution to 156 m with a constant grid refinement factor of 2. They compared their results to the observations of the HD(CP)² campaign and conclude that when it comes to small-scale to mesoscale variability, the use of LES drastically improves the results compared to their reference mesoscale model COSMO. PALM has also recently implemented an option for offline nesting within COSMO (Kadasch et al., 2021). They employ prescribed boundary conditions and impose synthetic turbulence in addition to the boundary fields. At the moment there is no relaxation zone implemented, but they do note that this might change in the future. In their test cases, they find that the boundary input has the largest impact on the main flow structures. Flow and updraughts rapidly develop with the help of the synthetic turbulence routine. Fully developed turbulence was found after 2 to 3 times the distance corresponding to the eddy turnover time.

In this research we develop a set of open LBCs for anelastic LES and implement them in the Dutch Atmospheric Large-Eddy Simulation (DALES) model. The goal of the paper is threefold. First, we will give a clear and extensive description of the open LBCs developed in this research. Second, we will show the influence of the LBCs on the mean fields and turbulent characteristics. Third, we will see how, in an idealized setup, the results depend on the temporal and spatial resolution of the input data, as one would encounter when embedding the LES in a coarser, non-turbulence-resolving LAM. The LBCs are developed to minimize reflections and the area needed for turbulence to develop and to allow for potential future one-way nesting with coarser LAMs. To minimize reflections, the outflow boundary conditions will be based on the radiation boundary condition of Sommerfeld (1949), and for the inflow boundary, a new set of Robin boundary conditions will be derived. To allow for one-way nesting with coarser LAMs, the open LBCs will be developed such that they allow time-varying input. The LBCs are tested with a simplistic dry convective case in a “Big Brother”-like setup (Denis et al., 2002). This allows us to single out the influence that the LBCs have on the fields in the interior of the domain. To study the influence of the spatial and temporal resolution of the boundary input data, the turbulence in the input data is filtered in both space and time simultaneously. This allows us to study the influence of the open LBCs in a setup where the LES is coupled to a non-turbulence-resolving model and to quantify the influence of the spatial and temporal resolution ratios between the parent and child model. We will investigate how long it takes for turbulence to fully develop. Furthermore, the influence of synthetic turbulence on generating inflow turbulence is explored.

2 Boundary condition implementation

This section will describe the implementation of the open boundary conditions in the Dutch Atmospheric Large-Eddy Simulation (DALES) model (Heus et al., 2010). The presented open boundary implementation is applicable to any incompressible atmospheric LES and, except for the discussion about mass conservation, could also be used for fully compressible LES. DALES solves the anelastic Navier–Stokes equations on a staggered Arakawa-C grid. The prognostic variables are the three velocity components ($u,v,w)$, liquid potential temperature (θ_l), total water specific humidity (q_t), rainwater specific humidity (q_r), rain droplet number concentration (N_r), subfilter-scale turbulence kinetic energy (e), and up to 100 active or passive scalars. Appropriate boundary conditions are required for all the prognostic variables at the resolution of the simulation. The velocity components are located at their respective cell faces and the rest of the variables at the cell centres. The boundary is defined as the cell faces of the outermost grid cells. Therefore, the boundary-normal velocity components are located at the boundary, whereas the other variables are located offset from the boundary. If the boundary input is not at the same time intervals as the simulation, the input data are linearly interpolated in time to the model time. First, the implementation for the boundary-normal velocity components will be given and conservation of mass will be discussed. Second, the implementation for the other variables is described. Third, the algorithm used to add synthetic turbulence at the boundaries will be discussed.

2.1 Boundary-normal velocity components

The boundary condition for the boundary-normal velocity components depends on whether the cell is an in- or outflow cell. An inflow cell for the boundary-normal velocity component is defined as ${\mathit{u}}^{\mathrm{B}}\cdot \widehat{\mathit{n}}<\mathrm{0}$, where u^B is the input velocity vector specified at the boundary, given by external data, and $\widehat{\mathit{n}}$ the outward-pointing boundary-normal unit vector. An outflow cell is defined by ${\mathit{u}}^{\mathrm{B}}\cdot \widehat{\mathit{n}}\ge \mathrm{0}$.

2.1.1 Outflow

The outflow boundary condition is based on the Sommerfeld radiation boundary condition (Sommerfeld, 1949), which states that disturbances should only be advected out of the domain with no reflections. The radiation boundary condition takes the form of a single propagating wave.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u_{\mathrm{n}}}{\partial t}}=\left\{\begin{array}{ll}-\frac{U}{\mathit{\rho }}\frac{\partial \mathit{\rho }{u}_{\mathrm{n}}}{\partial n}+\mathit{ϵ}& \text{for lateral boundaries}\\ -\frac{U}{\mathit{\rho }}\frac{\partial \mathit{\rho }{u}_{\mathrm{n}}}{\partial n}+g\frac{\mathit{\theta }-\langle \mathit{\theta }\rangle }{\langle \mathit{\theta }\rangle }+\mathit{ϵ}& \text{for top boundary}\end{array}\right.\end{array}$$

In Eq. (1) u_n is the boundary-normal velocity component; U the phase speed of the disturbances; $\partial /\partial n$ the boundary-normal derivative; ρ the reference density profile used by DALES; θ the potential temperature; g the gravitational acceleration; ϵ a correction factor required to conserve mass, which will be explained in more detail in Sect. 2.1.3; and 〈〉 a horizontal slab average. For the vertical component at the top boundary, the buoyancy force, which works as a damping factor for the top boundary in stably stratified flows, is added. The time derivative is discretized using DALES' third-order Runge–Kutta scheme (Heus et al., 2010). The spatial derivative is discretized using a first-order upwind scheme.

$$\begin{array}{}\text{(2)}& {\left.{\displaystyle \frac{\partial u_{\mathrm{n}}}{\partial n}}\right|}_i\approx \left\{\begin{array}{ll}\frac{u_i-u_{i-\mathrm{1}}}{\mathrm{\Delta }{x}_{\mathrm{n}}}& \text{for }{u}_B\ge \mathrm{0}\\ \frac{u_{i+\mathrm{1}}-u_i}{\mathrm{\Delta }{x}_{\mathrm{n}}}& \text{for }{u}_B<\mathrm{0}\end{array}\right.\end{array}$$

For non-dispersive waves with a phase speed equal to U, the 1D case of Eq. (1) without the correction factor ϵ will not generate any reflections. In the case of atmospheric simulations, which comprise a dispersive system, the phase speed U needs to be chosen carefully such that reflections are minimized. Popular implementations for the phase speed are given by Orlanski (1976) for those used in PALM (Maronga et al., 2015, 2020) and by Klemp and Wilhelmson (1978) for those used in Meso-NH (Lafore et al., 1998; Lac et al., 2018) and WRF-LES (Skamarock et al., 2021). Here we will use a slightly adjusted version of the implementation given by Hedley and Yau (1988). The implementation of Hedley and Yau (1988) is a hybrid version of the implementation given by Orlanski (1976) and is shown to work better. Similarly to Orlanski (1976), the velocity field and tendencies upstream of the boundary at the previous time step are used to define the local phase speed, which is then propagated to the boundary for the next time step. Additionally, Hedley and Yau (1988) set a fixed lower limit for the phase speed. We will set the lower limit to the boundary input normal velocity component, $u_{\mathrm{n}}^{\mathrm{B}}$.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& U^{*}={\left.U\right|}_{x_{\mathrm{n}}-\widehat{\mathit{x}}\cdot \widehat{\mathit{n}}\mathrm{\Delta }{x}_{\mathrm{n}}}^{t-\mathrm{\Delta }t}={〈-\mathit{\rho }{\displaystyle \frac{\partial u_{\mathrm{n}}}{\partial t}}{\left({\displaystyle \frac{\partial \mathit{\rho }{u}_{\mathrm{n}}}{\partial n}}\right)}^{-\mathrm{1}}〉}^{\text{int}},\\ & U=\left\{\begin{array}{ll}{u}_{\mathrm{n}}^{\mathrm{B}}& \text{if }\left|U^{*}\right|\le \left|u_{\mathrm{n}}^{\mathrm{B}}\right|\\ U^{*}& \text{if }\left|u_{\mathrm{n}}^{\mathrm{B}}\right|<\left|U^{*}\right|<\frac{\mathrm{\Delta }{x}_{\mathrm{n}}}{\mathrm{\Delta }t}\\ \text{sign}\left(U^{*}\right)\frac{\mathrm{\Delta }{x}_{\mathrm{n}}}{\mathrm{\Delta }t}& \text{if }\left|U^{*}\right|\ge \frac{\mathrm{\Delta }{x}_{\mathrm{n}}}{\mathrm{\Delta }t}\end{array}\right.\end{array}\end{array}$$

In Eq. (3) t−Δt denotes the previous time step and $x_{\mathrm{n}}-\widehat{\mathit{x}}\cdot \widehat{\mathit{n}}\mathrm{\Delta }{x}_{\mathrm{n}}$ the location one grid size upstream of the boundary. To avoid large fluctuations in the phase speed due to local gradients, the phase speed is averaged over the horizontal dimension perpendicular to the boundary vector over a distance of Δx^int (north and south boundaries) or Δy^int (west and east boundaries), denoted by 〈〉^int. This is similar to PALM, which averages laterally over the entire boundary (Maronga et al., 2015). For stability reasons the upper bound of the phase speed is set to the Courant–Friedrichs–Lewy (CFL) condition. Equation (3) is discretized using a first-order upwind scheme (Eq. 2).

2.1.2 Inflow

For inflow cells, the boundary-normal velocity at the boundary u_n is nudged towards the input value $u_{\mathrm{n}}^{\mathrm{B}}$ with a relaxation timescale equal to the integration timescale used by DALES (Δt). The discretization of the time derivative is given by the third-order Runge–Kutta scheme used by DALES (Heus et al., 2010).

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial u_{\mathrm{n}}}{\partial t}}={\displaystyle \frac{u_{\mathrm{n}}^{\mathrm{B}}-u_{\mathrm{n}}}{\mathrm{\Delta }t}}+\mathit{ϵ}\end{array}$$

2.1.3 Conservation of mass

The use of radiation boundary conditions means that continuity is not guaranteed and a correction factor, ϵ, needs to be added. Hedley and Yau (1988) enforce the rule that the height-integrated mass flux through each boundary does not change in time. This limits, however, the functionality for time-varying wind fields, in which inflow boundaries can become outflow boundaries and vice versa. Here we derive a correction term that forces the mass flux through the boundary to the boundary input on a defined length scale. This allows the wind field to change in magnitude and direction over time. To conserve mass, the following constraints are imposed.

1. The input boundary-normal velocity components integrated over the lateral and top boundaries S(B) satisfy the continuity equation that is calculated with the reference density profile used by DALES.

   $$\begin{array}{}\text{(5)}& \underset{S\left(B\right)}{\iint }\mathit{\rho }{\mathit{u}}^{\mathrm{B}}\cdot \widehat{\mathit{n}}\mathrm{d}S=\mathrm{0}\end{array}$$

2. The lateral and top boundaries are subdivided into patches S^int defined by Δy^int and Δz for the west and east boundaries, Δx^int and Δz for the north and south boundaries, and Δx^int and Δy^int for the top boundary. We enforce that the mass flux integrated over each patch equals the mass flux given by the input velocities integrated over the same patch.

   $$\begin{array}{}\text{(6)}& \underset{S^{\text{int}}}{\iint }\mathit{\rho }\mathit{u}\cdot \widehat{\mathit{n}}\mathrm{d}S=\underset{S^{\text{int}}}{\iint }\mathit{\rho }{\mathit{u}}^{\mathrm{B}}\cdot \widehat{\mathit{n}}\mathrm{d}S\end{array}$$

To obtain the correction factor ϵ, we define ϵ to be constant (in space) within a single integration patch S^int, but it can differ between patches. To obtain an expression for the correction term on a particular integration patch ϵ(S^int), we take the time derivative of Eq. (6). Further, we define $\frac{\partial {\stackrel{\mathrm{̃}}{u}}_{\mathrm{n}}}{\partial t}=\frac{\partial u}{\partial t}-\mathit{ϵ}$ as the tendency from either Eq. (1) or Eq. (4) minus the correction term. Within DALES, the tendencies for the boundary-normal velocities are first calculated without the correction term. These tendencies are then used to calculate the correction term ϵ for each integration patch using Eq. (7). The correction factor is then added to the tendencies before applying them to make sure mass is conserved.

$$\begin{array}{}\text{(7)}& \mathit{ϵ}\left(S^{\text{int}}\right)={\displaystyle \frac{{\iint }_{S^{\text{int}}}\mathit{\rho }\left(\frac{\partial u_{\mathrm{n}}^{\mathrm{B}}}{\partial t}-\frac{\partial {\stackrel{\mathrm{̃}}{u}}_{\mathrm{n}}}{\partial t}\right)\mathrm{d}S}{{\iint }_{S^{\text{int}}}\mathit{\rho }\mathrm{d}S}}\end{array}$$

The correction factor ϵ can be physically interpreted as the correction required to force the mass flux through the integration patch S^int to the mass flux integrated over the patch as given by the input. Since the constraint is set on the integrated quantity, fluctuations smaller than the set integration patch are conserved. Smaller values for Δx^int and Δy^int impose more strict boundary conditions, with Dirichlet conditions in the limit where Δx^int=Δx and Δy^int=Δy. When used in a nested simulation, Δx^int and Δy^int could be set to the grid size used by the parent model. In this setup, the total mass flux through a parent cell at the boundary of the child model (DALES) is conserved, while the child model is free to generate turbulence on smaller scales. This is illustrated in 2D in Fig. 1, in which the blue cells correspond to the parent model and have a resolution of Δx^parent and the brown cells correspond to the child model (DALES).

Figure 1 A 2D illustration of a nested setup in which the integration length scales are set to the grid size of the parent model. In this setup, the mass flux through a parent cell (blue) at the boundary of the child model (brown) is conserved, while the child model is free to generate turbulence on smaller scales.

The role of the correction term is to conserve mass integrated over the domain such that the pressure solver, which needs to find a solution that conserves mass locally, can find a solution. It is possible to implement the tendency from the correction factor as a non-homogeneous Neumann boundary condition for the modified pressure (defined in Heus et al., 2010) $\frac{\partial \mathit{\pi }}{\partial n}=-\mathit{ϵ}$ such that all the tendencies as a result of the continuity requirement are together. We chose, however, to add the term in the equations for the boundary-normal velocity components and use homogeneous Neumann boundary conditions for the modified pressure $\frac{\partial \mathit{\pi }}{\partial n}=\mathrm{0}$ because this allows us to keep using the Fourier pressure solver present in DALES (Heus et al., 2010) by using cosine basis functions only.

At the moment, the vertical length scale of the integration patch is fixed to the vertical grid resolution. This allows for a straightforward implementation when using stretched vertical grids. We have also experimented with setting the vertical length scale of the integration patch to the domain height. This couples the boundary layer with the column above the inversion layer and gave unwanted results. In the future, the implementation can be extended to allow for a variable vertical integration length scale as well.

2.2 Boundary-tangential velocity components and cell-centred variables

This section will discuss the boundary conditions for the cell-centred variables and the tangential velocity components. These variables are not computed at the boundary. Instead, ghost cells are used together with a second-order central discretization to determine the behaviour of the variable at the boundary. The implementation is different for in- and outflow boundaries. For the cell-centred variables and tangential velocity components, a boundary is defined as inflow if $\mathit{u}\cdot \widehat{\mathit{n}}<\mathrm{0}$ and as outflow otherwise. Note that this is different from the definition for the boundary-normal velocity components, where the nature of the boundary is determined by the input velocity $u_{\mathrm{n}}^{\mathrm{B}}$. These two can differ for outflow boundaries when the advection velocity is low and turbulence strong enough to reverse the local flow direction, as the radiation boundary condition does not enforce outflow on the local scale.

2.2.1 Outflow

For outflow cells, homogeneous Neumann conditions, Eq. (8), are specified at the lateral boundaries.

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}=\mathrm{0}\end{array}$$

In Eq. (8) ψ is any of the cell-centred variables (θ_l, q_t, q_r, N_r, e) or tangential velocity components. At the top of the domain, Neumann boundary conditions are set, which take the slab-averaged vertical derivative into account,

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}={\displaystyle \frac{\partial \langle \mathit{\psi }\rangle }{\partial z}},\end{array}$$

in which 〈〉 denotes a slab average. The decision to use homogeneous Neumann boundary conditions for all but the boundary-normal velocity components is based on the results of Sani and Gresho (1994) and Craske and Van Reeuwijk (2013). Sani and Gresho (1994) state that Neumann boundary conditions tend to produce fewer perturbations in comparison to a boundary condition placed on the variable itself (Dirichlet). Setting homogeneous Neumann conditions for the boundary-normal velocity components results in an ill-posed system with fluctuations in the pressure field and is not suited for turbulent flows (Sani and Gresho, 1994; Craske and Van Reeuwijk, 2013).

2.2.2 Inflow

For inflow boundaries, Dirichlet boundary conditions are a common choice (e.g. Maronga et al., 2015; Lac et al., 2018). However, for flows in which boundary cells change from in- to outflow boundaries and in which the outflow boundary is free to diverge from the boundary input, Dirichlet boundary conditions can result in large gradients over the boundary when they instantaneously set the value at the boundary to the boundary input value. For models that use radiation boundary conditions, this can result in unrealistically large tendencies at the boundary. Meso-NH poses a less strict Dirichlet inflow boundary condition by setting the boundary value to a weighted average between the input value and the nearest LES domain value, with a weight of 0.8 for the interior values (Lac et al., 2018). In this research we take a different approach and implement a Robin boundary condition, which will be derived in this section. The Robin boundary condition is a weighted average between a Dirichlet and Neumann boundary condition.

To derive the inflow boundary condition, we assume that advection is the only process taking place at the boundary.

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\psi }}{\partial t}}+u_{\mathrm{n}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}=\mathrm{0}\end{array}$$

We also impose the condition that the boundary value is nudged towards a given input value ψ^B over a timescale τ.

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathit{\psi }}{\partial t}}={\displaystyle \frac{{\mathit{\psi }}^{\mathrm{B}}-\mathit{\psi }}{\mathit{\tau }}}\end{array}$$

Combining these two constraints gives

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{{\mathit{\psi }}^{\mathrm{B}}-\mathit{\psi }}{\mathit{\tau }}}+u_{\mathrm{n}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}=\mathrm{0},\end{array}$$

which can be rewritten in the form of a Robin boundary condition:

$$\begin{array}{}\text{(13)}& \mathit{\psi }-u_{\mathrm{n}}\mathit{\tau }{\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}={\mathit{\psi }}^{\mathrm{B}}.\end{array}$$

The behaviour of Eq. (13) is determined by the value of u_nτ. Dirichlet and homogeneous Neumann conditions correspond to different limits.

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \underset{u_{\mathrm{n}}\mathit{\tau }\to \mathrm{0}}{lim}\mathit{\psi }={\mathit{\psi }}^{\mathrm{B}}\text{(Dirichlet)}\\ & \underset{u_{\mathrm{n}}\mathit{\tau }\to \pm \mathrm{\infty }}{lim}{\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}=\mathrm{0}\text{(homogeneous Neumann)}\end{array}\end{array}$$

The classical Dirichlet inflow conditions can thus be obtained by setting τ=0. When τ≠0, the boundary condition transitions from Dirichlet to homogeneous Neumann conditions as the velocity increases, avoiding large fluxes into the domain. At u_n=0, the transition point between in- and outflow conditions, the boundary condition changes from Dirichlet (inflow), u_nτ=0, to homogeneous Neumann (outflow). This transition can be smoothed by introducing a variable timescale for the inflow conditions. The inflow conditions were derived with the proposition that advection nudges the boundary over a fixed timescale. At very low velocities, advection plays a minor role and this assumption breaks down. To overcome this, the timescale needs to increase as the velocity approaches 0. The following requirements are set for τ:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& \underset{u_{\mathrm{n}}\to \mathrm{0}}{lim}{u}_{\mathrm{n}}\mathit{\tau }=\mathrm{\infty },\\ & \underset{u_{\mathrm{n}}\to \mathrm{\infty }}{lim}\mathit{\tau }={\mathit{\tau }}_{\mathrm{0}},\\ & {\mathit{\tau }}_{\mathrm{0}}=\mathrm{0}\Rightarrow \mathit{\tau }=\mathrm{0}.\end{array}\end{array}$$

The first condition is set such that Eq. (13) approaches homogeneous Neumann conditions for u_n=0, which removes the discontinuity. The second condition specifies that for large advection velocities we would like to have a constant nudging timescale τ₀. The third condition allows us to set the Robin inflow condition to Dirichlet inflow conditions when the nudging timescale is set to τ₀=0. A definition where τ is proportional to ${\left(\mathrm{1}/u_{\mathrm{n}}\right)}^p$ and p≥2 satisfies the conditions. The relation used is given by

$$\begin{array}{}\text{(16)}& \mathit{\tau }={\mathit{\tau }}_{\mathrm{0}}\left[\mathrm{1}+{\left|{\displaystyle \frac{u_{\mathrm{s}}}{u_{\mathrm{n}}}}\right|}^p\right],\end{array}$$

in which u_s is a subgrid velocity scale at the boundary. Here we used the square root of the subgrid turbulent kinetic energy taken from the subfilter-scale turbulence kinetic energy (SFS-TKE) scheme used by DALES (Heus et al., 2010). A different estimate can be used as well. When the resolved velocity is larger than the subgrid velocity, u_n≫u_s, the timescale reduces to τ=τ₀. When the resolved velocity drops below the subgrid velocity, the timescale will increase, providing a transition from the Robin boundary condition to the homogeneous Neumann condition at u_n=0. The final form of the Robin inflow boundary conditions is given by Eq. (17).

$$\begin{array}{}\text{(17)}& \mathit{\psi }-u_{\mathrm{n}}{\mathit{\tau }}_{\mathrm{0}}\left[\mathrm{1}+{\left|{\displaystyle \frac{u_s}{u_{\mathrm{n}}}}\right|}^p\right]{\displaystyle \frac{\partial \mathit{\psi }}{\partial n}}={\mathit{\psi }}^{\mathrm{B}}\end{array}$$

At the top of the domain the slab-averaged vertical gradient is taken into account. The Robin boundary condition at the top of the domain is given by Eq. (18).

$$\begin{array}{}\text{(18)}& \mathit{\psi }-w{\mathit{\tau }}_{\mathrm{0}}\left[\mathrm{1}+{\left|{\displaystyle \frac{u_s}{w}}\right|}^p\right]\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}-{\displaystyle \frac{\partial 〈\mathit{\psi }〉}{\partial z}}\right)={\mathit{\psi }}^{\mathrm{B}}\end{array}$$

2.3 Synthetic turbulence routine

To investigate the potential of synthetic turbulence in reducing the turbulence spinup area, the random flow generation (RFG) algorithm of Smirnov et al. (2001) is implemented. When used, ψ^B in Eqs. (18) and (17) and $u_{\mathrm{n}}^{\mathrm{B}}$ in Eq. (4) are replaced by ψ^B+ψ^R and $u_{\mathrm{n}}^{\mathrm{B}}+u_{\mathrm{n}}^{\mathrm{R}}$ respectively, where the superscript R denotes the perturbation given by the RFG algorithm. In the calculations for the mass conservation correction factor, ϵ (Eq. 7), $u_{\mathrm{n}}^{\mathrm{B}}$ is still used to satisfy condition Eq. (5). The RFG algorithm involves scaling and orthogonal transformation to create non-homogeneous anisotropic (near-)divergence-free velocity perturbations for a given covariance matrix $\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}$, turbulent length scale λ and turbulent timescale τ^R from the summation of N harmonic functions. The RFG routine is extended to give correlated potential temperature perturbations as well. From personal experience it is known that potential temperature perturbations are more effective in initiating turbulence than momentum perturbations. To create the potential temperature perturbations, a perturbation field is created from the summation of N harmonics,

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& \mathit{\alpha }=\sqrt{{\displaystyle \frac{\mathrm{2}}{N}}}\sum _{i=\mathrm{1}}^N{p}_i\cos \left({\mathbf{k}}_{\mathbf{i}}\cdot {\displaystyle \frac{\mathbf{x}}{\mathit{\lambda }}}+{\mathit{\omega }}_i{\displaystyle \frac{t}{{\mathit{\tau }}^{\mathrm{R}}}}\right)\\ & +q_i\sin \left({\mathbf{k}}_{\mathbf{i}}\cdot {\displaystyle \frac{\mathbf{x}}{\mathit{\lambda }}}+{\mathit{\omega }}_i{\displaystyle \frac{t}{{\mathit{\tau }}^{\mathrm{R}}}}\right),\\ & p,q,\mathit{\omega }\in N\left(\mathrm{0},\mathrm{1}\right),\\ & \mathbf{k}\in N\left(\mathrm{0},\mathrm{0.5}\right),\end{array}\end{array}$$

where x is the position vector, t the time, and N(μ,σ) samples from a normal distribution with mean μ and standard deviation σ. Next, the perturbation field is scaled for a given $\stackrel{\mathrm{‾}}{{{\mathit{\theta }}^{\prime }}^{\mathrm{2}}}$ and correlated to w^R for a given $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}$.

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathit{\theta }}^{\mathrm{R}}=\left(\mathit{\rho }{\displaystyle \frac{w^{\mathrm{R}}}{\sqrt{\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}}}}+\mathit{\alpha }\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}\right)\sqrt{\stackrel{\mathrm{‾}}{{{\mathit{\theta }}^{\prime }}^{\mathrm{2}}}}\\ & \mathit{\rho }={\displaystyle \frac{\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}}{\sqrt{\stackrel{\mathrm{‾}}{{{\mathit{\theta }}^{\prime }}^{\mathrm{2}}}\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}}}}\end{array}\end{array}$$

The RFG algorithm is easy to implement and is computationally inexpensive. A downside of the RFG algorithm is that it produces a Gaussian-model-like power spectrum. Huang et al. (2010) developed an improved algorithm that allows for any power spectra, but it comes with an additional computational cost. There are many other techniques and routines that are being used to help generate turbulence at inflow boundaries. However, the aim of this paper is not to study the performance of different inflow turbulence routines but rather to show the potential of adding perturbations to the boundary input fields in general.

3 Simulation setup and methodology

The test case setup used in this research is summarized in Fig. 2 and consists of a Big Brother-like setup (Denis et al., 2002) to test the performance of the open LBC implementation, simulations with spatially and temporally smoothed input to test the influence of turbulence present in the input data and simulations with added synthetic turbulence in addition to the smoothed input to see how these algorithms can help generate turbulence.

Figure 2 Illustration of the simulation setup. The solid blue rectangles show the different simulations and the sections in which their results are analysed.

The simulation case used in the test setup is the development of a dry convective boundary layer. This case is well understood, and DALES is known to produce realistic results (Heus et al., 2010). The dry convective boundary layer is forced with a constant surface heat flux of ${\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}}_{\mathrm{s}}=\mathrm{0.115}\mathrm{K}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, a zero surface momentum flux of $u^{*}=\mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and a geostrophic forcing in the east–west direction corresponding to $u_g=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The simulation is initialized with an east–west velocity of $U=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$; a north–south velocity of $V=\mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$; and an initial potential temperature profile that consists of a boundary layer with a temperature of 300 K, an inversion layer at 950 m and an inversion jump of Δθ=8 K over 120 m (linear interpolation between 300 and 308 K over 120 m) with a constant temperature gradient of $\frac{\partial \mathit{\theta }}{\partial z}=\mathrm{0.003}\mathrm{K}{\mathrm{m}}^{-\mathrm{1}}$ above. This corresponds to a convective velocity scale of $w^{*}=\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The domain size is $L_x\times L_y\times L_z=\mathrm{15.36}\times \mathrm{3.84}\times \mathrm{1.92}\mathrm{km}$ with a horizontal resolution of $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{60}\mathrm{m}$ and a vertical resolution of Δz=20 m. The simulations last 6 h and have an integration time step of Δt=5 s. The subgrid scheme used is the SFS-TKE scheme described in Heus et al. (2010). For the advection of all variables, DALES' second-order central scheme was used (Heus et al., 2010). This setup is very close to the dry (strong) convective boundary layer shown in Heus et al. (2010), which has already been studied by Sullivan et al. (1998). The differences are the addition of a mean background wind, a weaker surface heat flux, a higher horizontal resolution, the use of second-order advection schemes and a fixed integration time step. The initial profiles and the evolution over time of the potential temperature, east–west wind velocity, potential temperature flux and east–west wind variance are shown in Fig. 3.

Table 2 Setup parameters for the reference case. From left to right; grid spacing, domain size, integration time step, surface heat flux, surface momentum flux and geostrophic wind forcing.

Figure 3 Evolution of the periodic reference case from initial profiles to the end of simulation (6 h). Left to right; slab-averaged potential temperature, slab-averaged east–west velocity, and slab-averaged resolved and total heat flux, slab-averaged east–west velocity variance.

The Big Brother-like experiment, as was first proposed by Denis et al. (2002), consists of a simulation with open boundary conditions that is directly coupled to an identical reference simulation with periodic boundary conditions. The boundary fields of the periodic simulation are communicated at every time step to the simulation with open boundaries. This allows us to directly study the influence of the open boundary implementation, since both the periodic and the open boundary simulation are now identically forced and only differ in the implementation of their boundary conditions. The coupling is done offline, which means that the periodic simulation is done first and the boundary output is saved for every time step. This output is then used to force the simulation with open boundary conditions. In this setup, the west boundary is (mainly) an inflow boundary, the east boundary is (mainly) an outflow boundary, and the north and south boundaries will be in- and outflow boundaries changing for each grid cell and with time. The periodic simulation uses periodicity for the lateral boundaries and a no-stress boundary condition at the top (Heus et al., 2010). The simulation with open boundary conditions uses open boundary conditions for the lateral and top boundaries. First, we carry out a sensitivity analysis to study the dependence of the solution on the parameters introduced in the open boundary implementation. The parameters will be individually perturbed around a reference set. Next, a more in-depth analysis is conducted on the results of the simulation with the reference parameters. The parameters for the sensitivity analysis are listed in Table 3 with the default parameters in bold.

Table 3 Settings of the open boundary implementation for the sensitivity runs. The default settings are in bold.

In practice, the open boundary conditions will often be used to couple the LES to a coarser-resolution model, such as a mesoscale weather model. To study the impact of coarse-resolution (in space and time) boundary data, the periodic output is smoothed with a Gaussian filter before it is used to force the open boundary simulation. The simulation with open boundary conditions is repeated for different degrees of spatial and temporal smoothing. This setup emulates a one-way nesting setup and moves from the LES being nested in a turbulence-resolving model to a non-turbulence-resolving model. It also allows us to study the influence of resolution ratios between the parent and child model in a nested setup for both the spatial and the temporal resolutions. Since the smoothed fields come from the same model with the same model physics, resolution and subgrid parameterizations, any differences between the results of the simulation with the smoothed input and the reference (periodic) simulation must be caused by the boundary implementation and the smoothing. Comparison to the case without smoothing allows us to see the influence of smoothing, which relates to the resolution of and the turbulent scales present in the emulated parent model.

Different techniques exist to artificially add turbulence or increase the turbulent scales present in coarse data. To demonstrate the potential of one such technique, the synthetic turbulence algorithm of Smirnov et al. (2001) is implemented and expanded to give perturbations for the potential temperature as well (Sect. 2.3). The smoothed-input open boundary simulations are repeated with the addition of synthetic turbulence. The perturbations are created using height-dependent covariance matrices for u and θ obtained from the differences between the smoothed and non-smoothed-input fields. The turbulent length scale is set to the boundary layer height, which represents the largest turbulent eddies. The turbulent timescale is calculated as the turbulent length scale over the mean advection velocity. The covariance information would not be available in a real-case setup, but it allows us to see how the algorithm would perform in a best-case scenario. The purpose of these simulations is not to find the best synthetic turbulence implementation or to fine-tune the implementation used but to give an impression of how these routines can potentially improve the results.

4 Results and discussion

This section will describe the results of the test case described in Sect. 3. First, the performance of the open boundary implementation is evaluated using the coupled periodic–open boundary simulations. Second, the influence of input turbulence scales is described using the smoothed-input simulations. Third, the prospects of synthetic turbulence are explored.

4.1 Big Brother simulation

In this section the results of the Big Brother experiment are shown. In this setup the periodic boundary output is input into the simulation with open boundary conditions at the same spatial and temporal resolution. This setup allows us to investigate the definition and implementation of the boundary conditions. Any disturbances present in the simulation with open boundary conditions must be a direct result of the boundary implementation, as the periodic simulation supplies “perfect” boundary fields. It is a first necessary test that needs to be passed. The challenging areas are mainly the outflow (east) boundary and the north and south boundaries. At the outflow boundary, fields should leave the domain unperturbed and the area affected by reflections upstream of the outflow boundary should be minimal. The north and south boundaries are both in- and outflow boundaries and will therefore challenge the capability of the boundary conditions to switch from in- to outflow in time and space. The results from the simulation with open boundary conditions are compared to the reference case with periodic boundary conditions. We would like the mean field and the turbulence properties such as the length scales and energy distribution to be unaffected by the numerics of the boundary condition implementation. The two simulations do not have to match from a deterministic point of view, as the chaotic nature of the system will result in different placement of eddies between both simulations.

To investigate the sensitivity of the solution to the parameters of the open boundary implementation, the simulation is repeated for different sets of parameters. Each of the parameters is individually perturbed around the default values. The parameters and their values are shown in Table 3. Figure 4 shows the slab average profiles calculated over the last half hour of the simulation as a perturbation from the periodic profiles for potential temperature, eastward velocity, vertical potential temperature flux and eastward velocity variance. The profiles for the periodic simulation can be seen in Fig. 3. The black line represents the solution for the default values. Each colour represents a simulation where one of the parameters is perturbed, and the dashed or dotted line represents the perturbation value. Within the boundary layer, below 1000 m, the solution does not significantly depend on the values chosen for the parameters. All simulations are very close to the periodic simulation (within 1 %), indicating that the open boundary implementation does not have a significant impact on the solution.

At and above the inversion height, the simulation with a larger timescale for the Robin inflow conditions, τ₀=60 s, and the simulation without the buoyancy term in the top radiation boundary condition perform significantly worse than the other simulations. Without the buoyancy term in the top radiation boundary conditions, reflections from the top boundary result in distortions in the top layer of the simulation. Sometimes a sponge layer is implemented to dampen these types of reflections, but we do not need it here as, when used, the buoyancy term in the top radiation boundary condition solves the problem. The longer timescale for the Robin inflow conditions corresponds to a Robin boundary condition that is more weighted towards a Neumann boundary condition. Too long a timescale gives too much freedom at the inflow boundary and allows for waves to build up around the inversion layer. A shorter timescale such as used in the default settings therefore works better. The default timescale is not set to zero, which corresponds to Dirichlet conditions, because a slightly relaxed condition works better for simulations with lower mean background wind speeds.

For the integration length scale Δx^int and Δy^int, the simulation where they are set to the grid resolution shows the best results. This corresponds to Dirichlet boundary conditions for the boundary-normal velocity components. These settings work well for this setup because the boundary input is turbulent at the same resolution and from the same model. In other words, the simulation with open boundary conditions can find a solution that fits these boundary conditions. A larger integration length scale gives the LES more freedom and works better when the boundary input does not contain turbulence or is from a different model. The simulations have also been performed with a shorter time step of 2 s, and the results for all but the Robin boundary condition timescale remain the same. For the Robin boundary condition the optimum timescale is lower for a shorter time step, which requires further research. All the results shown from here on are obtained with the default settings.

Figure 4 Sensitivity analysis for the open boundary implementation parameters. Slab average profiles for simulations that have parameters perturbed around a default configuration (Table 3). The profiles are calculated over the last half hour of the simulation as a perturbation from the periodic profiles (Fig. 3). Left to right; potential temperature, eastward velocity, vertical potential temperature flux and eastward velocity variance.

Figures 5 and 6 show a top (xy) view at 110 m and a side (xz) view of the potential temperature respectively. The top view is shown as a perturbation with respect to the periodic slab average. The cross-sections are a snapshot after 6 h of simulation time. The top panel shows the results for the periodic simulation and the bottom panel for the simulation with open boundary conditions. The location of the xz cross-section within the xy cross-section (and vice versa) is shown by the dashed line. The slope of the solid line in the xz cross-section of the simulation with open boundary conditions corresponds to the ratio of the advective velocity scale ($U=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and convective velocity scale ($w^{*}=\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). Left (upstream) of this line, fields will be mainly dominated by information advected from the inflow boundary, whereas right of the line (downstream) the fields will be mainly influenced by convection originating from the surface boundary. The cross-sections are used to visually inspect the results to see if there are any discrepancies in the mean fields or turbulent structures. The simulations do not have to be similar from a deterministic point of view as the smallest differences at the boundaries would result in a different solution due to the chaotic nature of the system. The results of the open boundary simulation are very similar to the periodic simulation. The spatial scales and magnitude of the turbulent features resemble those of the periodic simulation. Up to 3 km from the inflow boundary (left), the turbulent features of the open boundary simulation are almost identical in shape and location to the periodic simulation, which shows that the turbulent boundary input fields at the inflow boundary are communicated well to the open boundary simulation. Further downwind they start to deviate as a result of the chaotic nature of the system. No clear disturbances at any of the boundaries are seen, and at the outflow boundary (right) the turbulent fields leave the domain without any significant reflections.

Figure 5 Horizontal cross-section of the potential temperature perturbation with respect to the periodic slab average at a height of 110 m for the periodic simulation (a) and open boundary simulation (b). The dotted line shows the location of the xz cross-section shown in Fig. 6.

Figure 6 Vertical cross-section of the potential temperature for the periodic simulation (a) and open boundary simulation (b). The dotted line shows the location of the xy cross-section shown in Fig. 5, and the slope of the solid line corresponds to the ratio of the advective velocity scale ($U=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and convective velocity scale ($w^{*}=\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$).

A more quantitative comparison of the influence of the open boundary conditions on the magnitude of the turbulent perturbations is obtained by calculating $\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{\sigma }}_y^{\mathrm{2}}\left(u\right)+{\mathit{\sigma }}_y^{\mathrm{2}}\left(v\right)+{\mathit{\sigma }}_y^{\mathrm{2}}\left(w\right)\right]$ for every time step and averaging it over the last half hour of the simulation. ${\mathit{\sigma }}_y^{\mathrm{2}}\left(\right)$ denotes the variance in the cross-wind (y) direction. This quantity is very close to the definition of turbulent kinetic energy (TKE) and will therefore be referred to as TKE from hereon. Figure 7 shows cross-sections of TKE for the periodic and open boundary condition experiments. The top panel shows the TKE for the periodic simulation and the bottom panel for the simulation with open boundary conditions. The dotted (dashed) grey contour lines mark the areas where the TKE values are smaller (larger) than the 2.5 % (97.5 %) percentile of the periodic simulation for that height. The slope of the solid black line corresponds to the ratio of the advective velocity scale ($U=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and the convective velocity scale ($w^{*}=\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and can be used as a measure of where the information from the surface boundary condition meets the information from the inflow boundary (left). The mean TKE values have similar magnitudes for both simulations. The simulation with open boundary conditions produces larger TKE values above the boundary layer and just before the outflow boundary (right). The increase in TKE above the boundary layer might be caused by the higher wind speeds present in the open boundary simulation (Fig. 4). The increased TKE values at the outflow boundary are the result of reflections and disappear 1 km upwind of the outflow boundary.

Figure 7 TKE profile derived from the cross-wind direction for the periodic (a) and open boundary (b) simulations. The slope of the solid line corresponds to the ratio of the advective velocity scale ($U=\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and convective velocity scale ($w^{*}=\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). The dotted (dashed) grey contour lines mark the areas where the TKE is smaller (larger) than the 2.5 % (97.5 %) percentile of the periodic simulation for that height.

To further quantify the differences between the simulations, we vertically integrate the TKE over the boundary layer (Fig. 8) along the cross-section shown in Fig. 7. We find that the magnitudes of TKE between the two simulations are very similar, indicating that the boundary conditions have virtually no influence on the Big Brother simulation, once again with the small exception of a slight accumulation of TKE at the outflow boundary.

Figure 8 TKE integrated over the boundary layer. The dashed lines show the mean and the mean ± 2 times the standard deviation.

A wavelet analysis of the potential temperature field is used to quantify the influence of the open boundary conditions on the power spectrum of the turbulence. Figure 9 shows a wavelet analysis for the periodic (top) and open boundary (bottom) simulations. A one-dimensional wavelet analysis is performed on an instantaneous xy slab after 6 h of simulation time. The wavelet analysis is done in the along-wind (x) direction. The results for each along-wind line are averaged over the cross-wind direction. A Morlet wavelet was used as the mother wavelet. The vertical axis shows the wavelength of the features on a logarithmic axis. The colours denote the wavelet power on a logarithmic scale. The shaded area indicates the cone of influence (COI); the COI describes the area that is potentially affected by boundary effects. These boundary effects result from the stretched wavelet extending beyond the edges of the domain, and results within the COI should therefore be ignored. The dotted (dashed) grey contour lines mark the areas where the wavelet energy is smaller (larger) than the 2.5 % (97.5 %) percentile of the periodic simulation for that wavelength. The wavelet analysis shows similar results for both simulations. As expected, the least energy is contained in the smallest wavelengths and the most energy is contained in features with a wavelength similar to the boundary layer height (≈10³ m). There are no clear differences visible between the periodic and open boundary wavelet analysis, which indicates that the open boundary implementation does not influence the turbulent power spectrum.

Figure 9 Wavelet analysis of the potential temperature at a height of 110 m for the periodic (a) and open boundary (b) simulations. The vertical axis shows the wavelengths of the features, the horizontal axis the distance from the inflow boundary and the colouring the energy present. The shaded area is the cone of influence and indicates the area that is potentially affected by boundary effects and results within should be ignored. The dotted (dashed) grey contour lines mark the areas where the wavelet energy is smaller (larger) than the 2.5 % (97.5 %) percentile of the periodic simulation for that wavelength.

From Figs. 5–9 it can be concluded that the influence of the open boundary implementation on the simulation is minimal. The slab-averaged fields, turbulent energy and spectral signature of the simulation are minimally perturbed by the implementation. Furthermore, the results of the sensitivity analysis show that solution is not sensitive to the values of the parameters as long as they are within a reasonable range and the buoyancy term in the top radiation boundary condition is used.

4.2 Smoothed-input simulations

This section will show and discuss the results of the smoothed-input simulations for different degrees of horizontal and temporal smoothing. This setup emulates the situation where the outer model provides boundary fields at a coarser spatial and/or temporal resolution than the LES. The panels in Figs. 10 and 11 show the same cross-sections as the bottom panels of Figs. 5 and 6 respectively for different degrees of smoothing. The horizontal axis of the panels shows the amount of smoothing in the temporal dimension and the vertical axis the amount of smoothing in the horizontal direction. The top-left cross-section is the result without smoothing and is the same as the bottom panels from Figs. 5 and 6. For low degrees of smoothing, σ_t≤30Δt and σ_x≤4Δx, the open boundary simulations resemble the periodic simulation and the solution is not significantly disturbed. For higher degrees of smoothing, wavelike structures emerge at the inflow boundary (left) that persist for up to 5 km into the domain. These structures become more prominent with increased smoothing. Horizontal smoothing (vertical axis) induces features that are aligned in the cross-wind direction. Temporal smoothing results in similar disturbances with the addition of some along-wind disturbances. The spatial–temporal smoothing does not affect the outflow boundary, where turbulent structures leave the domain unperturbed with no visual reflections.

Figure 10 Horizontal cross-section of the potential temperature perturbations with respect to the periodic simulation at a height of 110 m (similar to Fig. 5) for different degrees of smoothing. The horizontal axis of the panel shows the amount of smoothing in the temporal dimension and the vertical axis the amount of smoothing in the horizontal direction.

Figure 11 Vertical cross-section of the potential temperature (similar to Fig. 6) for different degrees of smoothing. The horizontal axis of the panel shows the amount of smoothing in the temporal dimension and the vertical axis the amount of smoothing in the horizontal direction.

Figure 12 TKE cross-section derived from the cross-wind direction (similar to Fig. 7) for different degrees of smoothing. The horizontal axis of the panel shows the amount of smoothing in the temporal dimension and the vertical axis the amount of smoothing in the horizontal direction.

Figure 12 shows the TKE cross-sections for the smoothed-input simulations. Smoothing the input reduces the turbulent scales present in the input data. This results in an area of reduced TKE downwind of the inflow boundary (left). The slope of the black line in Fig. 12 indicates the ratio between the advective and convective velocity scales ($U/w^{*}$). It is expected that upwind (left) of this line the solution will be predominately dominated by information advected from the inflow boundary, whereas downwind (right) of this line convection will take over. The area of reduced TKE values downwind of the inflow boundary increases with increased smoothing, and for large degrees of smoothing, the reduced TKE values extend much further than the line given by the $U/w^{*}$ ratio. For temporal smoothing of σ_t≥30Δt, a burst of TKE is present downwind of the reduced TKE area before settling to a TKE cross-section similar to that of the periodic simulation. This burst in TKE was also found by Muñoz-Esparza and Kosović (2018) and Kadasch et al. (2021). Our hypothesis is that the burst in TKE is a result of the clash between non-turbulent fields that are mainly governed by information supplied at the lateral inflow boundary and turbulent fields originating from surface convection. We believe that the sudden transition from non-turbulent flow to turbulent flow causes an overshoot in TKE. This phenomenon is also seen during the spinup time of (periodic) turbulent simulations. During the first hour, the turbulence in the boundary layer needs to build up. Only after this is developed is it capable of transporting the accumulated surface moisture and heat flux through the boundary layer, causing a peak not only in TKE but also in cloud fraction if clouds are formed on the top of the boundary layer (e.g. Siebesma and Cuijpers, 1995). In the worst cases (highest degrees of smoothing) it can take up to 6–7 km before the TKE settles to values similar to those of the periodic simulation. The TKE field near the outflow boundary is not affected by the smoothing. The wavelike structures seen in Fig. 10 are not visible in the TKE cross-sections, as they are aligned in the cross-wind direction, the same direction over which the TKE is calculated.

Once again, the results are quantified further by vertically integrating the TKE over the boundary layer along the cross-section for all simulations (Fig. 13). Each simulation is shown as a thin line, with the control (no smoothing) and representative simulations for strong temporal and/or spatial smoothing highlighted in colour. Compared to the Big Brother experiment (Fig. 8), deviations from the control (periodic) simulations are much larger, in particular at the inflow boundary but also elsewhere in the domain. Comparing Figs. 8 and 13 once again highlights that the limitations in the open boundary simulations are mostly introduced by the spatial and temporal smoothing of the boundary values and not by the implementation of the boundary conditions themselves.

Figure 13 TKE integrated over the boundary layer. Each line represents one of the simulations from Fig. 12, with the corners from the panels being highlighted in colour. The dashed blue lines represent the mean TKE from the periodic simulation and the mean ± 2 times the standard deviation.

The wavelet analysis for the smoothed-input simulations is shown in Fig. 14. For low horizontal and temporal smoothing, σ_x≤4Δx and σ_t≤30Δt, the influence on the results is small and the wavelet cross-section remains close to the periodic cross-section. For higher degrees of smoothing, an increase in energy for wavelengths of around 300 m is seen at the inflow boundary. The energy increase at these wavelengths represents the waves seen in Fig. 14. The energy increase is larger for increased smoothing. For high degrees of smoothing, a decrease in energy is seen for wavelengths around 1 km at the inflow boundary. This represents the lack of developed turbulence near the inflow boundary as a result of the missing turbulence in the input data. The energy distribution moves towards the periodic profile downstream of the inflow boundary, with the maximum energy moving towards turbulence of the scale of the boundary layer height. For the highest degrees of smoothing, this takes around 7 km. The smoothing does not influence the wavelet spectrum at the outflow boundary.

Figure 14 Wavelet analysis of the potential temperature at a height of 110 m (similar to Fig. 9) for different degrees of smoothing. The horizontal axis of the panel shows the amount of smoothing in the temporal dimension and the vertical axis the amount of smoothing in the horizontal direction.

Figure 15 Horizontal cross-section of the potential temperature perturbations with respect to the periodic simulation at a height of 110 m for different degrees of smoothing (similar to Fig. 10) with the addition of synthetic turbulence.

Figure 16 Vertical cross-section of the potential temperature for different degrees of smoothing (similar to Fig. 11) with the addition of synthetic turbulence.

Figure 17 TKE cross-section derived from the cross-wind direction for different degrees of smoothing (similar to Fig. 12) with the addition of synthetic turbulence.

The results analysed in Figs. 10–14 show that the input smoothing causes the solution to deteriorate. For high degrees of smoothing, turbulent structures are missing at the inflow boundary and cross-wind-oriented wavelike disturbances form. In the worst cases it can take up to 7 km before the turbulent intensity and spectral signal evolve towards values close to the results of the periodic simulation. These results are important to take into account when coupling LES models to regional weather models. The latter usually have a spatial resolution on the order of kilometres, and common output intervals are on the order of hours. This means that the ratios with the LES grid size and time step are at the bottom right of the shown panels. To avoid large ratios between the resolution of the input data and the LES model, repeated nesting can be used. With repeated nesting the LES can step down from the regional weather model resolution towards the desired resolution in steps with a determined refinement ratio. The results in this section suggest that a ratio of 4 between the spatial resolutions and a ratio of 30 between the temporal output and LES time step should not be exceeded. In practice this is often hard to achieve, especially the temporal constraint as weather model data are often saved on an hourly interval. Another approach is to artificially add finer turbulent scales to the input data. This can be done by turbulence recycling, dedicated turbulence simulations or synthetic turbulence (e.g. Tabor and Ahmadi, 2010).

4.3 Synthetic turbulence simulations

The previous section has highlighted significant issues at the inflow boundary when the boundary values are smoothed in space and/or time, resulting in a more laminar flow near that boundary. A potential approach to reduce these issues (Smirnov et al., 2001) is to add synthetic turbulence to the boundary values. The purpose of this section is to investigate how the results in our simulations are affected by doing so. The algorithm of Smirnov et al. (2001) is implemented and extended to give potential temperature perturbations as well (Sect. 4.3). Figures 15 and 16 show the cross-sections for potential temperature. Compared to Figs. 10 and 11 three things stand out. First, the addition of perturbations seems to remove the persistent wavelike structures at the inflow boundary (left). There are still disturbances present such as the perturbation at the inversion, but the persistent wavelike structure is gone. Second, the disturbances seem to disappear more quickly downstream. Third, the disturbances do not increase in magnitude with increased smoothing.

Figure 18 TKE integrated over the boundary layer. Each line represents one of the simulations from Fig. 17, with the corners from the panels being highlighted in colour. The dashed blue lines represent the mean TKE from the periodic simulation and the mean ± 2 times the standard deviation.

Figure 17 shows the TKE cross-sections for the simulations with synthetic turbulence. The addition of synthetic turbulence increases the TKE values directly downstream of the inflow boundary. The values are still below the developed-turbulence values of the periodic simulation. The synthetic turbulence does help to generate developed turbulence faster, which results in a smaller downstream area where the TKE is too low. The overshoot after the reduced TKE is also smaller in magnitude and area compared to Fig. 12. Furthermore, the overshoot no longer increases with increased smoothing in contrast to the results without added synthetic turbulence. The overshoot is similar in shape, magnitude and location for all smoothed simulations and is located near the line that represents the ratio of the convective and advective velocity scales, where the information from the inflow boundary meets with the information from the surface. This means that it can be predicted where the overshoot is and which part of the simulation should be ignored. All simulations settle to a profile close to the periodic simulation within 5 km of the inflow boundary. The simulations with high degrees of smoothing do have an area with too much TKE. The addition of synthetic turbulence does not seem to have an influence on the outflow boundary.

A quantitative comparison once again using the vertical integral of TKE over the boundary layer (Fig. 18) confirms the positive influence of adding synthetic turbulence, showing a much reduced discrepancy from the control simulation at the inflow boundary in comparison to the simulations without it (Fig. 13).

Figure 19 shows the wavelet analysis for the smoothed-input open boundary simulations. The increase in energy for wavelengths around 300 m is still visible at the inflow boundary, but the magnitude has been reduced by the turbulent perturbations. The decrease in energy for wavelengths around 1 km is no longer there. Furthermore, the wavelet profile converges much faster to the periodic profile and the results do not seem to worsen with increased smoothing. The addition of synthetic turbulence does not affect the outflow boundary.

Figure 19 Wavelet analysis of the potential temperature at a height of 110 m for different degrees of smoothing (similar to Fig. 14) with the addition of synthetic turbulence.

The results analysed in Figs. 15–19 show that the addition of synthetic turbulence on top of coarse input data can improve the simulation results. All of the inflow disturbances found in Sect. 4.2, as a result of coarse input data, were reduced in size and/or magnitude by the addition of synthetic turbulence. Furthermore, the location of the disturbances became predictable and their magnitude and size no longer increased with increased smoothing. The better performance when using synthetic turbulence may appear trivial. However, as we cannot add turbulence that is directly compatible with the LES solution, the synthetic turbulence could be dampened or generate artefacts near the inflow boundary. The fact that it does not shows the value of using it in our implementation.

5 Conclusions

This paper introduced an open boundary implementation for atmospheric large-eddy simulation models that was implemented in the Dutch Atmospheric Large-Eddy Simulation (DALES) model. The goal of this research was to give a detailed description of the implementation, investigate its performance, and show the influence of open boundary conditions and boundary input on the solution.

Radiation boundary conditions were implemented as an outflow condition for the boundary-normal velocity components at the lateral and top boundaries. At the top boundary, buoyancy was also taken into account, which negated the need to add a sponge layer in the upper parts of the domain. Neumann conditions were used for the other variables at outflow boundaries. For inflow boundaries a Robin boundary condition was derived for the cell-centred variables and tangential velocity components to allow for a smooth transition between in- and outflow boundaries, and a nudging condition was implemented for the boundary-normal velocity components.

Using a Big Brother-like setup, where a simulation with open boundary conditions was forced by an identical control simulation with periodic boundary conditions on the same spatial and temporal resolution, it was shown that the influence of the boundary implementation on the solution was minimal. Slab-averaged profiles showed that the mean profiles are conserved. Furthermore, cross-sections of the potential temperature field showed that the turbulent input data were communicated well through the inflow boundary and that the turbulent fields left the domain without reflections or perturbations at the outflow boundary. Cross-wind turbulent kinetic energy cross-sections showed that the energy in the turbulent perturbations were the same in the simulation with open boundary conditions and the control simulation with periodic boundary conditions. The energy spectrum of the perturbations was also unchanged, which was shown with a wavelet analysis.

To investigate the influence of the spatial and temporal resolution of the input data, the output of the periodic simulation was smoothed before feeding it to the simulation with open boundary conditions. Different degrees of spatial and temporal smoothing showed that a mismatch between input turbulent scales and model scales results in the generation of wavelike disturbances downstream of the inflow boundary. The disturbances grow in size and magnitude when the ratio between input and model scales grows. The lack of turbulence in the input data also results in an area of reduced turbulent kinetic energy downstream of the inflow boundary, where there is no developed turbulence. This area grew as the smoothing increased. For large degrees of smoothing it was found that the turbulent energy overshoots before settling to values similar to the periodic control simulation. For these reasons, it is advised to be careful when coupling a large-eddy simulation model with open boundary conditions to a coarser model. Repeated nesting can be used and is currently being explored to step down in multiple steps from coarse data to the desired resolution. The results of this research indicate that the refinement factor when nesting should not exceed 4 in the spatial dimension and 30 in the time dimension.

The potential of adding synthetic turbulence to the LBCs was explored, and the results show that it can help to reduce the disturbances found in size and magnitude and to speed up the process of obtaining developed turbulence by artificially reducing the gap between the input turbulent scales and model scales. The strong wavelike character of the disturbances were removed, and the length of the inflow area required for turbulence to develop was reduced. The disturbances and development area also became less dependent on the degree of smoothing, and the development area is given by the ratio of the advective and convective velocity scales. However, if possible, we would still advise keeping the spatial and temporal ratios between the input data and the LES below the earlier-mentioned values.

In summary, the implementation of open BCs described in this study provides a suitable framework for further investigating the use of the DALES model in “nested” mode. This provides a major advance in its utility as a science tool, as it increases its applicability to problems for which periodic BCs have strong limitations, such as over heterogeneous terrain. Spatial and temporal averaging of the boundary values, as is typical of embedding an LES into coarser-resolution mesoscale models, causes the results to deteriorate. The smoothing effects are much larger than those from the implementation of the open BCs themselves. Some of the deterioration can be overcome by adding synthetic turbulence at the inflow boundaries.