the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Open boundary conditions for atmospheric largeeddy simulations and their implementation in DALES4.4
Franciscus Liqui Lung
Christian Jakob
A. Pier Siebesma
Fredrik Jansson
Open boundary conditions were developed for atmospheric largeeddy simulation (LES) models and implemented into the Dutch Atmospheric LargeEddy 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.
 Article
(23807 KB)  Fulltext XML
 BibTeX
 EndNote
Largeeddy 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 limitedarea 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 largescale 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 largescale 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 boundarynormal 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 convectiondominated flows. Incompressible LES models such as PALM (Maronga et al., 2015, 2020) and MesoNH (Lac et al., 2018) and the fully compressible WRFLES (Skamarock et al., 2021) have the option to use radiation boundary conditions for the boundarynormal velocity components based on one of the previously mentioned implementations. For the other variables, (homogeneous) Neumann BCs, which specify the boundarynormal 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)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 WRFLES to conduct twoway 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 oneway and a twoway 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 nonturbulence 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 oneway and a twoway nested setup with WRFLES 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 WRFLES domain and the NCEP data set. However, they found that the cloud fields can be strongly modulated by mesoscale organization, especially in highwind conditions where the clouds align with the mean wind direction. They found little benefit of twoway nesting over oneway nesting. Talbot et al. (2012) coupled WRFLES within WRF in a realistic oneway 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 nearsurface 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 differentresolution mesoscale simulations. They test their results both with and without the addition of the cell perturbation method of MuñozEsparza 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 highestresolution mesoscale model. Heinze et al. (2017) used a oneway 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)^{2} campaign and conclude that when it comes to smallscale 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 LargeEddy 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, nonturbulenceresolving LAM. The LBCs are developed to minimize reflections and the area needed for turbulence to develop and to allow for potential future oneway 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 oneway nesting with coarser LAMs, the open LBCs will be developed such that they allow timevarying 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 nonturbulenceresolving 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.
This section will describe the implementation of the open boundary conditions in the Dutch Atmospheric LargeEddy 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 ArakawaC 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}), subfilterscale 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 boundarynormal 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 boundarynormal 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 Boundarynormal velocity components
The boundary condition for the boundarynormal velocity components depends on whether the cell is an in or outflow cell. An inflow cell for the boundarynormal 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 outwardpointing boundarynormal 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.
In Eq. (1) u_{n} is the boundarynormal velocity component; U the phase speed of the disturbances; $\partial /\partial n$ the boundarynormal 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' thirdorder Runge–Kutta scheme (Heus et al., 2010). The spatial derivative is discretized using a firstorder upwind scheme.
For nondispersive 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 MesoNH (Lafore et al., 1998; Lac et al., 2018) and WRFLES (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}}$.
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 firstorder upwind scheme (Eq. 2).
2.1.2 Inflow
For inflow cells, the boundarynormal 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 thirdorder Runge–Kutta scheme used by DALES (Heus et al., 2010).
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 heightintegrated mass flux through each boundary does not change in time. This limits, however, the functionality for timevarying 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.

The input boundarynormal 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}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}S=\mathrm{0}\end{array}$$ 
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}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}S=\underset{{S}^{\text{int}}}{\iint}\mathit{\rho}{\mathit{u}}^{\mathrm{B}}\cdot \widehat{\mathit{n}}\phantom{\rule{0.125em}{0ex}}\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{\u0303}}{u}}_{\mathrm{n}}}{\partial t}=\frac{\partial u}{\partial t}\mathit{\u03f5}$ as the tendency from either Eq. (1) or Eq. (4) minus the correction term. Within DALES, the tendencies for the boundarynormal 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.
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).
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 nonhomogeneous Neumann boundary condition for the modified pressure (defined in Heus et al., 2010) $\frac{\partial \mathit{\pi}}{\partial n}=\mathit{\u03f5}$ 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 boundarynormal 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 Boundarytangential velocity components and cellcentred variables
This section will discuss the boundary conditions for the cellcentred variables and the tangential velocity components. These variables are not computed at the boundary. Instead, ghost cells are used together with a secondorder central discretization to determine the behaviour of the variable at the boundary. The implementation is different for in and outflow boundaries. For the cellcentred 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 boundarynormal 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.
In Eq. (8) ψ is any of the cellcentred 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 slabaveraged vertical derivative into account,
in which 〈〉 denotes a slab average. The decision to use homogeneous Neumann boundary conditions for all but the boundarynormal 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 boundarynormal velocity components results in an illposed 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. MesoNH 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.
We also impose the condition that the boundary value is nudged towards a given input value ψ^{B} over a timescale τ.
Combining these two constraints gives
which can be rewritten in the form of a Robin boundary condition:
The behaviour of Eq. (13) is determined by the value of u_{n}τ. Dirichlet and homogeneous Neumann conditions correspond to different limits.
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 τ:
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 τ_{0}. The third condition allows us to set the Robin inflow condition to Dirichlet inflow conditions when the nudging timescale is set to τ_{0}=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
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 subfilterscale turbulence kinetic energy (SFSTKE) 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 τ=τ_{0}. 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).
At the top of the domain the slabaveraged vertical gradient is taken into account. The Robin boundary condition at the top of the domain is given by Eq. (18).
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 nonhomogeneous anisotropic (near)divergencefree velocity perturbations for a given covariance matrix $\stackrel{\mathrm{\u203e}}{{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,
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{\u203e}}{{{\mathit{\theta}}^{\prime}}^{\mathrm{2}}}$ and correlated to w^{R} for a given $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\theta}}^{\prime}}$.
The RFG algorithm is easy to implement and is computationally inexpensive. A downside of the RFG algorithm is that it produces a Gaussianmodellike 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.
The test case setup used in this research is summarized in Fig. 2 and consists of a Big Brotherlike 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.
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{\u203e}}{{w}^{\prime}{\mathit{\theta}}^{\prime}}}_{\mathrm{s}}=\mathrm{0.115}\phantom{\rule{0.125em}{0ex}}\mathrm{K}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, a zero surface momentum flux of ${u}^{*}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ and a geostrophic forcing in the east–west direction corresponding to ${u}_{g}=\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$. The simulation is initialized with an east–west velocity of $U=\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$; a north–south velocity of $V=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}\mathrm{K}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}$ above. This corresponds to a convective velocity scale of ${w}^{*}=\mathrm{1.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}\mathrm{km}$ with a horizontal resolution of $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{60}\phantom{\rule{0.125em}{0ex}}\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 SFSTKE scheme described in Heus et al. (2010). For the advection of all variables, DALES' secondorder 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 secondorder 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.
The Big Brotherlike 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 nostress 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 indepth 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.
In practice, the open boundary conditions will often be used to couple the LES to a coarserresolution model, such as a mesoscale weather model. To study the impact of coarseresolution (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 oneway nesting setup and moves from the LES being nested in a turbulenceresolving model to a nonturbulenceresolving 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 smoothedinput open boundary simulations are repeated with the addition of synthetic turbulence. The perturbations are created using heightdependent covariance matrices for u and θ obtained from the differences between the smoothed and nonsmoothedinput 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 realcase setup, but it allows us to see how the algorithm would perform in a bestcase scenario. The purpose of these simulations is not to find the best synthetic turbulence implementation or to finetune the implementation used but to give an impression of how these routines can potentially improve the results.
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 smoothedinput 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, τ_{0}=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 boundarynormal 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.
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 crosssections 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 crosssection within the xy crosssection (and vice versa) is shown by the dashed line. The slope of the solid line in the xz crosssection of the simulation with open boundary conditions corresponds to the ratio of the advective velocity scale ($U=\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$) and convective velocity scale (${w}^{*}=\mathrm{1.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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 crosssections 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.
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 crosswind (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 crosssections 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}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$) and the convective velocity scale (${w}^{*}=\mathrm{1.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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.
To further quantify the differences between the simulations, we vertically integrate the TKE over the boundary layer (Fig. 8) along the crosssection 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.
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 onedimensional wavelet analysis is performed on an instantaneous xy slab after 6 h of simulation time. The wavelet analysis is done in the alongwind (x) direction. The results for each alongwind line are averaged over the crosswind 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^{3} 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.
From Figs. 5–9 it can be concluded that the influence of the open boundary implementation on the simulation is minimal. The slabaveraged 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 Smoothedinput simulations
This section will show and discuss the results of the smoothedinput 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 crosssections 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 topleft crosssection 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 crosswind direction. Temporal smoothing results in similar disturbances with the addition of some alongwind disturbances. The spatial–temporal smoothing does not affect the outflow boundary, where turbulent structures leave the domain unperturbed with no visual reflections.
Figure 12 shows the TKE crosssections for the smoothedinput 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 crosssection similar to that of the periodic simulation. This burst in TKE was also found by MuñozEsparza and Kosović (2018) and Kadasch et al. (2021). Our hypothesis is that the burst in TKE is a result of the clash between nonturbulent 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 nonturbulent 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 crosssections, as they are aligned in the crosswind 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 crosssection 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.
The wavelet analysis for the smoothedinput 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 crosssection remains close to the periodic crosssection. 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.
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 crosswindoriented 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 crosssections 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 17 shows the TKE crosssections 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 developedturbulence 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 smoothedinput 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.
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.
This paper introduced an open boundary implementation for atmospheric largeeddy simulation models that was implemented in the Dutch Atmospheric LargeEddy 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 boundarynormal 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 cellcentred variables and tangential velocity components to allow for a smooth transition between in and outflow boundaries, and a nudging condition was implemented for the boundarynormal velocity components.
Using a Big Brotherlike 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. Slabaveraged profiles showed that the mean profiles are conserved. Furthermore, crosssections 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. Crosswind turbulent kinetic energy crosssections 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 largeeddy 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 earliermentioned 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 coarserresolution 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.
The current version of DALES is available from the project website – https://github.com/dalesteam/dales (last access: 25 March 2024) – under the GNU General Public License. The exact version of the model used to produce the results used in this paper is archived on Zenodo (https://doi.org/10.5281/zenodo.10046420), as are input data and scripts to run the model and produce the plots for all the simulations presented in this paper (Liqui Lung et al., 2023).
FLL conceptualized the paper, did the formal analysis and visualization, implemented the methodology and software, and wrote and edited the draft. CJ and APS supervised during the project, reviewed the draft and supported in the conceptualization. FJ supported the software implementation and reviewed the draft.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We acknowledge the use of ECMWF's computing and archive facilities; the funding provided by the Australian Research Council Centre of Excellence for Climate Extremes (CE170100023); and the support of the Ruisdael Observatory, scientific research infrastructure which is (partly) financed by the Dutch Research Council (NWO, grant no. 184.034.015)
This research is funded by the Australian Research Council Centre of Excellence for Climate Extremes (grant no. CE170100023) and by the Ruisdael Observatory scientific research infrastructure, which is (partly) financed by the Dutch Research Council (NWO, grant no. 184.034.015).
This paper was edited by Sylwester Arabas and reviewed by two anonymous referees.
Carpenter, K. M.: Note on the paper “Radiation conditions for the lateral boundaries of limitedarea numerical models”, Q. J. Roy. Meteor. Soc., 108, 717–719, https://doi.org/10.1002/qj.49710845714, 1982. a
Craske, J. and Van Reeuwijk, M.: Robust and accurate open boundary conditions for incompressible turbulent jets and plumes, Comput. Fluids, 86, 284–297, https://doi.org/10.1016/j.compfluid.2013.06.026, 2013. a, b, c
Deardorff, J.: Numerical Investigation of Neutral and Unstable Planetary Boundary Layers, J. Atmos. Sci., 29, 91–115, https://doi.org/10.1175/15200469(1972)029<0091:NIONAU>2.0.CO;2, 1972. a
Denis, B., Laprise, R., Caya, D., and Côté, J.: Downscaling ability of oneway nested regional climate models: The BigBrother Experiment, Clim. Dynam., 18, 627–646, https://doi.org/10.1007/s0038200102010, 2002. a, b, c
Giometto, M., Christen, A., Meneveau, C., Fang, J., Krafczyk, M., and Parlange, M.: Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface, Bound.Lay. Meteorol., 160, 425–452, https://doi.org/10.1007/s1054601601576, 2016. a
Hedley, M. and Yau, M.: Radiation Boundary Conditions in Numerical Modeling, Mon. Weather Rev., 116, 1721–1736, https://doi.org/10.1175/15200493(1988)116<1721:RBCINM>2.0.CO;2, 1988. a, b, c, d, e, f
Heinze, R., Dipankar, A., Henken, C. C., Moseley, C., Sourdeval, O., Trömel, S., Xie, X., Adamidis, P., Ament, F., Baars, H., Barthlott, C., Behrendt, A., Blahak, U., Bley, S., Brdar, S., Brueck, M., Crewell, S., Deneke, H., Di Girolamo, P., Evaristo, R., Fischer, J., Frank, C., Friederichs, P., Göcke, T., Gorges, K., Hande, L., Hanke, M., Hansen, A., Hege, H.C., Hoose, C., Jahns, T., Kalthoff, N., Klocke, D., Kneifel, S., Knippertz, P., Kuhn, A., van Laar, T., Macke, A., Maurer, V., Mayer, B., Meyer, C. I., Muppa, S. K., Neggers, R. A. J., Orlandi, E., Pantillon, F., Pospichal, B., Röber, N., Scheck, L., Seifert, A., Seifert, P., Senf, F., Siligam, P., Simmer, C., Steinke, S., Stevens, B., Wapler, K., Weniger, M., Wulfmeyer, V., Zängl, G., Zhang, D., and Quaas, J.: Largeeddy simulations over Germany using ICON: a comprehensive evaluation, Q. J. Roy. Meteor. Soc., 143, 69–100, https://doi.org/10.1002/qj.2947, 2017. a, b, c, d, e
Heus, T., van Heerwaarden, C. C., Jonker, H. J. J., Pier Siebesma, A., Axelsen, S., van den Dries, K., Geoffroy, O., Moene, A. F., Pino, D., de Roode, S. R., and VilàGuerau de Arellano, J.: Formulation of the Dutch Atmospheric LargeEddy Simulation (DALES) and overview of its applications, Geosci. Model Dev., 3, 415–444, https://doi.org/10.5194/gmd34152010, 2010. a, b, c, d, e, f, g, h, i, j, k
Huang, S., Li, Q., and Wu, J.: A general inflow turbulence generator for large eddy simulation, J. Wind Eng. Ind. Aerod., 98, 600–617, https://doi.org/10.1016/j.jweia.2010.06.002, 2010. a
Kadasch, E., Sühring, M., Gronemeier, T., and Raasch, S.: Mesoscale nesting interface of the PALM model system 6.0, Geosci. Model Dev., 14, 5435–5465, https://doi.org/10.5194/gmd1454352021, 2021. a, b, c, d
Klemp, J. and Wilhelmson, R.: The Simulation of ThreeDimensional Convective Storm Dynamics, J. Atmos. Sci., 35, 1070–1096, https://doi.org/10.1175/15200469(1978)035<1070:TSOTDC>2.0.CO;2, 1978. a, b, c, d, e, f
Kurppa, M., Hellsten, A., Auvinen, M., Raasch, S., Vesala, T., and Järvi, L.: Ventilation and Air Quality in City Blocks Using LargeEddy Simulation–Urban Planning Perspective, AtmosphereBasel, 9, 65, https://doi.org/10.3390/atmos9020065, 2018. a
Lac, C., Chaboureau, J.P., Masson, V., Pinty, J.P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.P., Lebeaupin Brossier, C., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V., and Wautelet, P.: Overview of the MesoNH model version 5.4 and its applications, Geosci. Model Dev., 11, 1929–1969, https://doi.org/10.5194/gmd1119292018, 2018. a, b, c, d, e
Lafore, J.P., Stein, J., Asencio, N., Bougeault, P., Ducrocq, V., Duron, J., Fischer, C., Héreil, P., Mascart, P., Masson, V., Pinty, J.P., Redelsperger, J.L., Richard, E., and Arellano, J.: The MesoNH Atmospheric Simulation System. Part I: Adiabatic formulation and control simulations, Ann. Geophys., 16, 90–109, https://doi.org/10.1007/s005850050582, 1998. a, b
Lilly, D.: The presentation of smallscale turbulence in numerical simulation experiments, 281, National Center for Atmospheric Research: Boulder, CO, USA, https://doi.org/10.5065/D62R3PMM, also featured in Proc. IBM Sci. Comput. Symp. on Environmental Science 1967, 195–210, 1966. a
Liqui Lung, F. P. A., Siebesma, A. P., Jansson, F. R., Arabas, S., Axelsen, S. L., Attema, J., Azizi, V., Beets, C., Boeing, S. J., de Bruine, M., Chylik, J., Cuijpers, H., van Dorp, P., van der Dussen, J., Duynkerke, P., van Heerwaarden, C., Heus, T., Janssens, M., Jonker, H., Moene, A., Nieuwstadt, F., Ouwersloot, H., van den Oord, G., Pourquie, M., de Roode, S., Neggers, R., Pedruzo, X., Sikma, M., van Stratum, B., Vila, J., and van Zanten, M.: franslql/dales: DALES 4.4 with open boundary conditions, (v4.4_openBC), Zenodo [code and data set], https://doi.org/10.5281/zenodo.10046420, 2023. a
Maronga, B., Gryschka, M., Heinze, R., Hoffmann, F., KananiSühring, F., Keck, M., Ketelsen, K., Letzel, M. O., Sühring, M., and Raasch, S.: The Parallelized LargeEddy Simulation Model (PALM) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments, and future perspectives, Geosci. Model Dev., 8, 2515–2551, https://doi.org/10.5194/gmd825152015, 2015. a, b, c, d, e
Maronga, B., Banzhaf, S., Burmeister, C., Esch, T., Forkel, R., Fröhlich, D., Fuka, V., Gehrke, K. F., Geletič, J., Giersch, S., Gronemeier, T., Groß, G., Heldens, W., Hellsten, A., Hoffmann, F., Inagaki, A., Kadasch, E., KananiSühring, F., Ketelsen, K., Khan, B. A., Knigge, C., Knoop, H., Krč, P., Kurppa, M., Maamari, H., Matzarakis, A., Mauder, M., Pallasch, M., Pavlik, D., Pfafferott, J., Resler, J., Rissmann, S., Russo, E., Salim, M., Schrempf, M., Schwenkel, J., Seckmeyer, G., Schubert, S., Sühring, M., von Tils, R., Vollmer, L., Ward, S., Witha, B., Wurps, H., Zeidler, J., and Raasch, S.: Overview of the PALM model system 6.0, Geosci. Model Dev., 13, 1335–1372, https://doi.org/10.5194/gmd1313352020, 2020. a, b, c
Mazzaro, L., MuñozEsparza, D., Lundquist, J., and Linn, R.: Nested mesoscaletoLES modeling of the atmospheric boundary layer in the presence of underresolved convective structures, J. Adv. Model. Earth Sy., 9, 1795–1810, https://doi.org/10.1002/2017MS000912, 2017. a, b
Mehta, D., Zuijlen, A., Koren, B., Holierhoek, J., and Bijl, H.: Large Eddy Simulation of wind farm aerodynamics: A review, J. Wind Eng. Ind. Aerod., 133, 1–17, https://doi.org/10.1016/j.jweia.2014.07.002, 2014. a
Mirocha, J., Kosović, B., and Kirkil, G.: Resolved Turbulence Characteristics in LargeEddy Simulations Nested within Mesoscale Simulations Using the Weather Research and Forecasting Model, Mon. Weather Rev., 142, 806–831, https://doi.org/10.1175/MWRD1300064.1, 2014. a
Moeng, C.H., Dudhia, J., Klemp, J., and Sullivan, P.: Examining TwoWay Grid Nesting for Large Eddy Simulation of the PBL Using the WRF Model, Mon. Weather Rev., 135, 2295–2311, https://doi.org/10.1175/MWR3406.1, 2007. a, b, c, d
MuñozEsparza, D. and Kosović, B.: Generation of Inflow Turbulence in LargeEddy Simulations of Nonneutral Atmospheric Boundary Layers with the Cell Perturbation Method, Mon. Weather Rev., 146, 1889–1909, https://doi.org/10.1175/MWRD180077.1, 2018. a
MuñozEsparza, D., Kosovic, B., Mirocha, J., and Beeck, J.: Bridging the Transition from Mesoscale to Microscale Turbulence in Numerical Weather Prediction Models, Bound.Lay. Meteorol., 153, 409–440, https://doi.org/10.1007/s1054601499569, 2014. a
MuñozEsparza, D., Kosović, B., van Beeck, J., and Mirocha, J.: A stochastic perturbation method to generate inflow turbulence in largeeddy simulation models: Application to neutrally stratified atmospheric boundary layers, Phys. Fluids, 27, 035102, https://doi.org/10.1063/1.4913572, 2015. a
Ol'shanskii, M. A. and Staroverov, V. M.: On simulation of outflow boundary conditions in finite difference calculations for incompressible fluid, Int. J. Numer. Meth. Fl., 33, 499–534, https://doi.org/10.1002/10970363(20000630)33:4<499::AIDFLD19>3.0.CO;27, 2000. a
Orlanski, I.: Simple boundary condition for unbounded flows, J. Comput. Phys., 21, 251–269, https://doi.org/10.1016/00219991(76)900231, 1976. a, b, c, d, e, f, g, h, i
Sani, R. and Gresho, P.: Résumé and remarks on the Open Boundary Condition Minisymposium, Int. J. Numer. Meth. Fl., 18, 983–1008, https://doi.org/10.1002/fld.1650181006, 1994. a, b, c, d
Schalkwijk, J., Jonker, H., Siebesma, A., and Meijgaard, E.: Weather Forecasting Using GPUBased LargeEddy Simulations, B. Am. Meteorol. Soc., 96, 715–723, https://doi.org/10.1175/BAMSD1400114.1, 2015. a
Siebesma, A. P. and Cuijpers, J. W. M.: Evaluation of Parametric Assumptions for Shallow Cumulus Convection, J. Atmos. Sci., 52, 650–666, https://doi.org/10.1175/15200469(1995)052<0650:EOPAFS>2.0.CO;2, 1995. a
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D., and Huang, X.: A Description of the Advanced Research WRF Model Version 4.3, Tech. Rep. NCAR/TN556+STR, National Center for Atmospheric Research, Boulder, CO, US, https://doi.org/10.5065/1dfh6p97, 2021. a, b, c
Smirnov, A., Shi, S., and Celik, I.: Random Flow Generation Technique for Large Eddy Simulations and ParticleDynamics Modeling, J. Fluid. Eng.T. ASME, 123, 359–371, https://doi.org/10.1115/1.1369598, 2001. a, b, c, d
Sommerfeld, A.: Partial differential equations in physics, Academic press, https://doi.org/10.1016/B9780126546583.X50010, 1949. a, b, c
Sommeria, G.: ThreeDimensional Simulation of Turbulent Processes in an Undisturbed Trade Wind Boundary Layer, J. Atmos. Sci., 33, 216–241, https://doi.org/10.1175/15200469(1976)033<0216:TDSOTP>2.0.CO;2, 1976. a
Sullivan, P. P., Moeng, C.H., Stevens, B., Lenschow, D. H., and Mayor, S. D.: Structure of the Entrainment Zone Capping the Convective Atmospheric Boundary Layer, J. Atmos. Sci., 55, 3042–3064, https://doi.org/10.1175/15200469(1998)055<3042:SOTEZC>2.0.CO;2, 1998. a
Tabor, G. and Ahmadi, M.: Inlet conditions for large eddy simulation: A review, Comput. Fluids, 39, 553–567, https://doi.org/10.1016/j.compfluid.2009.10.007, 2010. a
Talbot, C., BouZeid, E., and Smith, J.: Nested Mesoscale LargeEddy Simulations with WRF: Performance in Real Test Cases, J. Hydrometeorol., 13, 1421–1441, https://doi.org/10.1175/JHMD11048.1, 2012. a, b
Wesseling, P.: Principles of computational fluid dynamics, vol. 29, Springer Science & Business Media, https://doi.org/10.1007/9783642051463, 2009. a
Zhu, P., Albrecht, B., Ghate, V., and Zhu, Z.: Multiplescale simulations of stratocumulus clouds, J. Geophys. Res., 115, D23201, https://doi.org/10.1029/2010JD014400, 2010. a, b, c