Self-nested large-eddy simulations in PALM model system v21.10 for offshore wind prediction under different atmospheric stability conditions

Large-eddy simulation (LES) resolves large-scale turbulence directly and parametrizes small-scale turbulence. Resolving micro-scale turbulence, e.g., in wind turbine wakes, requires both a sufficiently small grid spacing and a domain large enough to develop turbulent flow. Refining a grid locally via a nesting interface effectively decreases the required computational time compared to the global grid refinement. However, interpolating the flow between nested grid boundaries introduces another source of uncertainty. Previous studies reviewed nesting effects for a buoyancy-driven flow and observed a secondary circulation in the two-way nested area. Using a nesting interface with a shear-driven flow in LES, therefore, requires additional verification. We use PALM model system 21.10 to simulate a boundary layer in a cascading self-nested domain under neutral, convective, and stable conditions and verify the results based on the wind speed measurements taken at the FINO1 platform in the North Sea.

We show that the feedback between parent and child domains in a two-way nested simulation of a non-neutral boundary layer alters the circulation in the nested area, despite spectral characteristics following the reference measurements. Unlike the pure buoyancy-driven flow, a non-neutral shear-driven flow slows down in a two-way nested area and accelerates after exiting the child domain. We also briefly review the nesting effect on the velocity profiles and turbulence anisotropy.

1 Introduction

Large-eddy simulation (LES) allows performing a detailed process study for areas and situations where we lack appropriate field measurements. For this reason, LES is widely used for high-fidelity simulations of wind flows in wind energy applications. When considering the turbulent flow, the grid resolution should be sufficiently high to resolve the relevant turbulence scales (Wurps et al., 2020). Increased grid resolution comes at the cost of gradually increased computational time. The overall computational time can be reduced by refining a grid locally through the nesting interface. While improving the grid resolution, a nesting interface introduces new uncertainties in the simulation. Such nesting effects are documented for buoyancy-driven flows, with the strongest influence observed for the two-way nesting mode (Moeng et al., 2007; Hellsten et al., 2021). A buoyancy-driven flow develops a secondary circulation and decreased velocity inside the nested area – the effect becomes prominent for the data averaged over several hours. However, buoyancy-driven flows are characterized by near-zero wind speed, while the wind energy research primarily deals with wind speeds of 5–25 m s⁻¹. Therefore, shear-driven LES with the nesting interface requires additional verification.

We use the Fortran-based LES code PALM 21.10 (Maronga et al., 2020) to simulate wind flow with a speed of 12.5 m s⁻¹ at the reference height of 119 m for three stability conditions: true neutral (NBL), convective (CBL), and stable (SBL) boundary layers. The initial velocity and turbulence intensity profiles are defined to match 1 h averages of the sonic anemometer time series as processed by Nybø et al. (2019). The domain is simulated for a non-nested grid and nested grids with one-way or two-way nesting modes. The resulting turbulence statistics are then compared with the measurements to evaluate the model's performance.

2 Data

The reference measurements contain wind speed directional components u, v, and w recorded with sonic anemometers during the Offshore Boundary-Layer Experiment at FINO1 (OBLEX-F1) campaign in 2015–2016 in the North Sea. The meteorological mast is installed on the FINO1 platform located in the North Sea at 54^∘00^′53.5^′′ N, 6^∘35^′15.5^′′ E, 45 km to the north of the German island of Borkum.

Table 1 Aggregated statistics of 1 h sonic anemometer time series.

The sonic anemometers were installed at the meteorological mast at 40, 60, and 80 m. The measurements were processed by Nybø et al. (2019) and organized into 1 h time series of 1 Hz frequency. Each processed series corresponds to different pairs of a stability condition and mean wind speed at the reference height of 119 m. This height was chosen as an outlook into future wind turbine development and corresponds to a hub height of the DTU reference 10 MW turbine (Bak et al., 2013). The reference height unifies different stability conditions under the assumption of a similar flow speed. Due to the computational time restrictions, we simulate only those conditions where the horizontal wind speed reaches approximately ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{119}}=\mathrm{12.5}$ m s⁻¹ at the reference height (Table 1).

The wind speed ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{119}}$ at the reference height was estimated from the measurement data. Since the measurements are originally available only for three levels, the mean wind speed profile was approximated by Nybø et al. (2020) by fitting the logarithmic law

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{u}\left(z\right)=u_{\mathrm{F}\mathrm{1}}\left[{\displaystyle \frac{\ln \left(\frac{z}{z_{\mathrm{0}}}-\mathit{\psi }\right)}{\ln \left(\frac{z_{\mathrm{F}\mathrm{1}}}{z_{\mathrm{0}}}-\mathit{\psi }\right)}}\right],\end{array}$$

where the wind speed u_F1 measured at FINO1 is taken for the highest available level z_F1=80 m, and the stability correction function ψ is defined as in Stull (1988):

$$\begin{array}{}\text{(2)}& \mathit{\psi }=\left\{\begin{array}{ll}\mathrm{0}& \text{ – NBL},\\ -\mathrm{2}\ln \frac{\mathrm{1}+x}{\mathrm{2}}-\ln \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}+\mathrm{2}\arctan x-\frac{\mathit{\pi }}{\mathrm{2}}& \text{ – CBL},\\ \mathrm{4.7}\mathit{\zeta }& \text{ – SBL},\end{array}\right.\end{array}$$

where $x=(\mathrm{1}-\mathrm{15}\mathit{\zeta }{)}^{\mathrm{1}/\mathrm{4}}$. The stability parameter ζ is derived from the height above the surface z and Obukhov length L as

$$\begin{array}{}\text{(3)}& \mathit{\zeta }={\displaystyle \frac{z}{L}}.\end{array}$$

The roughness length z₀ in Eq. (1) is, therefore, a fitting parameter to be found. The estimation is based on the assumption that the boundary layer extends beyond 119 m so that the logarithmic law can be applied to the mean wind profile. During the simulation, we attempt to match the mean wind profile, including the estimated wind speed at 119 m and turbulence intensity calculated for levels 40, 60, and 80 m.

3 Methodology

3.1 PALM LES model

We perform a free-flow large-eddy simulation (LES) using the Fortran code PALM developed at Leibniz Universität Hannover (Maronga et al., 2020). PALM utilizes a staggered Arakawa C grid: the velocity components are defined at the grid cell edges and are shifted by a half-grid spacing; the scalar variables are defined at the center of a grid cell. The subgrid-scale fluxes are resolved via the Deardorff 1.5-order closure model.

By default, PALM solves prognostic equations for the velocity components u, v, and w and potential temperature θ. If the stability condition is set to true neutral, the temperature is considered constant, and the corresponding equation is not solved. Buoyancy terms are also not considered in a true neutral simulation

A nested simulation in PALM consists of at least one child domain inside a parent domain. Each child domain can simultaneously be a parent domain for another child domain, thus forming a cascading self-nested structure. The top-level parent domain is further referred to as the root domain to make a distinction from inner parent domains. Overall, PALM supports simulation of one root domain and up to 63 child domains.

The nesting algorithm is constructed in a way to optimize computational time for multiple child domains (Hellsten et al., 2021). The nested domains communicate via interpolation which is performed just before the pressure-correction step, so that the time-consuming pressure solver is run only once per the time step. The solution at the nested boundaries of a parent domain – velocity components and scalar quantities, e.g., temperature and humidity – is linearly interpolated to all nested boundaries, except the bottom surface, as boundary conditions. The bottom surface is always located at a zero level as in the root domain and utilizes Dirichlet or Neumann boundary conditions as prescribed in the corresponding child domain input files.

After the interpolation, the prognostic equations are solved for a child domain. In the case of cascading nesting, the procedure is repeated until the solution is found for all nested domains at the current step. In a one-way nesting case, the simulation proceeds to the pressure-correction step, so the solution in parent domains remains unaffected by the solution in child domains. In a two-way nesting case, PALM uses an anterpolation scheme – a term suggested by Sullivan et al. (1988) and first described by Clark and Farley (1984). The technical details behind the implementation in PALM are explained in Hellsten et al. (2021). Each child domain anterpolates its solution via first-order integration to the respective parent domain before the pressure-correction step. Therefore, the two-way nested solution remains similar in the nested area, while the one-way nested solution may eventually diverge for parent and child domains.

3.2 Precursor and main LES run parameters

One of the ways PALM can simulate a turbulent flow is a precursor scheme, which does not require complex dynamic input data and effectively reduces the domain size required for turbulence development (Witha et al., 2014). First, a small precursor domain is simulated with cyclic boundaries until the flow reaches a steady state. The resulting mean wind speed and temperature profiles are then copied over the larger main domain to set up an initial non-cyclic flow with a developed turbulence. Provided that the main run is simulated with the same forcing as the precursor, the mean profiles in the main run remain stationary.

The size of the precursor domain is usually smaller than for the main run, and the y-shift procedure is performed at left/right cyclic boundaries to avoid non-physical regularity of the flow (Munters et al., 2016). The y-shift procedure is also applied in the main run for an additional disruption of regularity. Using the precursor scheme also ensures that an idealized input flow remains the same within a stability case regarded.

The grid characteristics of the root and innermost child domain in the PALM simulation were selected to closely match the SOWFA simulation in Nybø et al. (2020). The ratio between the parent and child domains' grid spacing, thus, would reach 8 (from 10 to 1.25 m for NBL and CBL cases) or 4 (from 5 to 1.25 m for SBL case). As shown by Hellsten et al. (2021), the discrepancy with a fine-grid simulation in PALM increases if the grid spacing ratio is 4 or higher. Therefore, we add intermediate child domains and reduce the grid spacing by a factor of 2 until the desired refinement is reached. Hence, NBL and CBL simulations contain three child domains, while the SBL simulation has two (Tables 2 and 3, Fig. 1).

Table 2 Grid parameters for NBL and CBL nested domains (Fig. 1a).

Table 3 Grid parameters for SBL nested domains (Fig. 1b).

Figure 1 Nested domains schematic. (a) NBL and CBL domains and (b) SBL domains.

We perform one-way and two-way nested simulations. To evaluate the nesting effect, we also simulate domains without nested grids using the same precursor flow. Due to high computational time and memory requirements, we only simulate non-nested domains for the grid spacing of Δ_x = 10 and 5 m.

Table 4 Input parameters of the precursor runs.

The precursor profiles undergo development during a simulation and thus may deviate from the initial profiles. The precursor's input parameters are then selected so that the resulting steady-state profiles of mean wind speed and turbulence intensity follow the values estimated from the measurements, particularly the wind speed at the reference height. The Coriolis force is switched off; hence the required wind speed and turbulence intensity profiles in the precursor run are enforced by a combination of the parameters: the initial mean wind ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{0}}$, the pressure gradient forcing $\mathrm{d}p/\mathrm{d}x$, and the roughness length z₀. The NBL case is run as the true neutral flow with no heat flux. The CBL case is defined via the positive heat flux $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}$ in addition to the parameters mentioned above. The SBL case uses surface cooling over time $\mathrm{d}{T}_{\mathrm{s}}/\mathrm{d}t$ instead of the heat flux (Wurps et al., 2020). NBL and SBL cases start with zero temperature gradient; the CBL case has an initial temperature gradient of 1 K (100 m)⁻¹. The surface temperature T_s is varied to match the conditions observed during the reference meteorological measurements at FINO1. The precursor domain characteristics and input parameters are listed in Tables 2–4.

During the precursor simulation, the initial profiles are altered due to the influence of pressure forcing and heat fluxes. The resulting precursor profiles are provided in Table 5; the same profiles are used to initialize the main run.

Table 5 Steady state of the precursor runs – turbulent inflow for the main run.

We run main simulations for 3 h with a dynamic time step selected by the model. The simulation is then continued for another hour with the fixed time step of Δt = 0.05 s to obtain a high-frequency output. Then, we probe time series of each wind speed component at the center of the innermost child domain and the corresponding points of the parent domain (Fig. 1). The high-frequency time series are further used to compare turbulence statistics with the measurements. Spatial averages (cross-sectional flows, profiles) are calculated for 10 min periods.

3.3 Turbulence characteristics

We evaluate the model performance based on turbulence characteristics: power spectrum, coherence, co-coherence, and phase. The coherence represents a correlation between time series a(t) and b(t) at two points separated by a certain distance δ and is calculated as follows

$$\begin{array}{}\text{(4)}& {\text{Coh}}_{ab}={\displaystyle \frac{S_{ab}}{\sqrt{S_{aa}{S}_{bb}}}},\end{array}$$

where S_aa and S_bb are the spectral densities of a(t) and b(t), while S_ab is the cross-spectrum of the same series.

The co-coherence represents the real part of the coherence

$$\begin{array}{}\text{(5)}& {\text{Co}}_{ab}=\mathit{Re}{\text{Coh}}_{ab}=\mathit{Re}{\displaystyle \frac{S_{ab}}{\sqrt{S_{aa}{S}_{bb}}}}.\end{array}$$

The phase ϕ_ab shows the level of synchronicity between time series a(t) and b(t)

$$\begin{array}{}\text{(6)}& {\mathit{\varphi }}_{ab}=\text{arctan}{\displaystyle \frac{\mathit{Re}{\text{Coh}}_{ab}}{\mathit{Im}{\text{Coh}}_{ab}}}.\end{array}$$

Since the measurement time series are available only for three levels, 40, 60, and 80 m, the spectra are calculated and compared at h=80 m for the total horizontal $U=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}$ and vertical w wind speed. The co-coherence is calculated for two vertical separations of δ = 20 m (between levels 60 and 80 m) and δ = 40 m (between levels 40 and 80 m). The sampling frequency for the LES time series matches the output frequency $f_{\mathrm{s}}^{\mathrm{LES}}=\mathrm{1}/\mathrm{0.05}$ s = 20 Hz, and the segment length is chosen as 60 s. The sampling frequency for the measurement time series is lower, $f_{\mathrm{s}}^{\mathrm{mast}}=\mathrm{1}/\mathrm{0.1}$ s = 10 Hz, although the segment length is left the same.

3.4 Flow characteristics for load analysis

We also review flow characteristics relevant to the turbine performance analysis: power law coefficient and turbulence anisotropy.

The power law is commonly applied to assess wind resources at the hub height from near-surface wind speed measurements.

$$\begin{array}{}\text{(7)}& U\left(z\right)={\stackrel{\mathrm{‾}}{U}}_{\mathrm{10}}{\left({\displaystyle \frac{z}{\mathrm{10}}}\right)}^{\mathit{\alpha }},\end{array}$$

where ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{10}}$ is the wind speed at z=10 m and α is the power law exponent. The power law exponent is sensitive to atmospheric conditions and is usually approximated with a constant; e.g., $\mathit{\alpha }=\mathrm{1}/\mathrm{7}$ is applicable to neutral onshore sites but not other stabilities (Touma, 1977). Often, the approximations do not reflect seasonal and diurnal variations in mean wind profiles (Bratton and Womeldorf, 2011; Jung and Schindler, 2021). Hence, simulating a long time series with the LES gives a possibility to study wind profiles in detail.

The anisotropic turbulence naturally develops in a simulation with an anisotropic grid resolution (Haering et al., 2019) but may also occur in isotropic grids, such as those used in this study. The anisotropic turbulence affects wind turbine loads, particularly fatigue loads. Therefore, it is important to evaluate its strength in the simulation (Dimitrov et al., 2017). We estimate turbulence anisotropy by comparing spectra of velocity components for the normalized frequency $f_n=fz/U_z$, where z=80 m and U_z is the horizontal velocity at this level. We compute ratios $S_{vv}/S_{uu}$ and $S_{ww}/S_{uu}$ for all regarded cases at f_n≈1. The closer both ratios are to the theoretical value of $\mathrm{4}/\mathrm{3}=\mathrm{1.333}$, the more isotropic the simulated turbulence is (Weiler and Burling, 1967; Smedman et al., 2003).

4 Results

4.1 Nesting effects

All LESs are run at 1024 cores for each case with a time step of Δt = 0.05 s; the required simulation times for each scenario are summarized in Table 6. Since the domains vary in size and number of grid points, we compare not the total CPU time but the CPU time per second of the simulated time. The non-nested coarse domain (Δ_x = 10 m) is not computationally demanding, regardless of the stability case. However, the required CPU time gradually increases if the grid spacing is reduced globally for the whole domain. As could be seen for the NBL case, the CPU time per second of the simulated time increases from 5.1 s for Δ_x=10 m to 31.7 s for Δ_x=5 m, respectively. Refining the grid locally by adding child domains increases the CPU time compared to the coarse reference non-nested grid (Δ_x = 10 m). Still, the nested simulation finishes faster than the globally refined non-nested simulation (Δ_x = 5 m), while allowing better a local grid refinement up to Δ_x = 1.25 m.

Table 6 CPU time in seconds used per second of simulated time. All simulations run at 1024 cores with a time step of Δt = 0.05 s.

Both NBL and CBL simulations have the same domain structure and grid spacing (Table 2). However, CBL simulations require more CPU time compared to the respective NBL (true neutral) simulations due to solving the temperature equation. SBL simulations use CPU time comparable to NBL simulations due to having one child domain less and a smaller root domain size – and thus a lower overall number of the grid points (Table 3).

Two-way nested simulations require additional ∼ 2–3 s of the CPU time per simulated time step to anterpolate the child domain solution back to the parent domain. This results in about 10 % increase in the CPU time compared to one-way nesting.

It should be noted that, unless obtaining high-frequency time series is the main goal of a simulation, the time step can be gradually increased for non-nested runs in order to speed up the computation. The computational time will, nevertheless, increase in a similar proportion with the global grid refinement. The time step in nested runs is still limited by the lowest grid spacing in child domains. For example, the dynamic step in the regarded configuration does not exceed 0.075 s to satisfy Courant–Friedrichs–Lewy condition.

Depending on the simulation conditions, LES produces different results in the nested area. If the true neutral case is defined in PALM explicitly via setting a corresponding flag, the one-way and the two-way nested simulations behave similarly with respect to grid spacing and feedback between domains (Fig. 2). Switching on the true neutral flag means that the temperature equation and buoyancy terms are not considered in the calculations. As long as those terms are introduced for non-neutral simulations, the two-way nested simulation results in a decreased flow speed in child domains.

Figure 2 NBL, flow at the reference height of 119 m for different wind speed components: (a) one-way nesting and (b) two-way nesting.

Since the child domains anterpolate their solution back to the parent domain, the area of reduced flow speed spreads to the root domain. While the effect is less prominent for the instantaneous fields, it becomes apparent in the 10 min averaged flow (Fig. 3). The induction of downward vertical wind in two-way nested simulations was already described by Hellsten et al. (2021) for the 5 h averaged buoyancy-driven flow in PALM. Hellsten et al. (2021) argued that the effect of the secondary circulation described by Moeng et al. (2007) was caused solely by the insufficient domain size and explained it with the different grid spacing and subsequent divergence of the vertical heat flux in the parent and child domains. The researchers hypothesized that the secondary circulation was an inevitable side effect of the two-way nesting solution due to the better resolution of the turbulence mixing in child domains. In the case of the shear-driven flow, we observe that the slowing effect is more prominent and develops faster. The effect emerges in the beginning of the simulation within 20 min – an approximate time required for the precursor flow to pass the main run domain. In addition, some of the quantities of a shear-driven flow, mainly the vertical velocity w, are not uniformly distributed inside the child domains (Fig. 4).

Figure 3 SBL, flow at the reference height of 119 m for different wind speed components: (a) one-way nesting and (b) two-way nesting.

Figure 4 The 10 min average profiles, SBL two-way nested case. (a) Sampling points; (b) the mean flow is slowed down in the nested area; (c) the vertical flow near the entrance of the nested area remains weak but becomes stronger as the flow passes through the nested area.

4.2 Subgrid scales

LES resolves scales larger than the grid spacing directly but approximates smaller scales. In a well-resolved flow, the unresolved (subgrid) scales should not exceed the resolved ones. This relation holds for all simulations performed, implying that the grid spacing of Δ=10 m is already small enough for the given flow (Fig. 5). The grid refinement does not strongly affect momentum fluxes, except for the CBL case (Fig. 5b), where turbulent eddies are generally larger than in the NBL and SBL cases. The effect from the nesting mode is also the most pronounced in CBL simulations (Fig. 5b). The resolved $\stackrel{\mathrm{‾}}{wu}$ and $\stackrel{\mathrm{‾}}{wv}$ fluxes remain stationary in the one-way nesting mode but decrease over time in the two-way nesting mode and eventually merge.

Figure 5 Comparison of resolved and subgrid-scale momentum fluxes for different stability simulations and nesting modes.

Figure 6 Comparison of near-surface resolved and subgrid-scale momentum fluxes for different stability simulations and nesting modes.

The subgrid-scale fluxes consistently remain near zero for all levels except near-surface cells, where the turbulence intensity is expected to be high due to the surface influence (Fig. 6). Consequently, the near-surface subgrid-scale fluxes are comparable to resolved-scale fluxes. However, the subgrid-scale fluxes at lower levels tend to zero faster as the grid spacing is refined. Unlike the one-way nesting mode, the resolved fluxes in the two-way nesting mode show a non-monotonic behavior near the surface in the intermediate child domains. The effect is observed in all two-way simulations, including true neutral conditions. Therefore, it cannot be solely caused by the flow difference in the nested and non-nested areas, despite the flux profiles being time and spatial averages. The occurring non-monotonic behavior can be rather attributed to the way PALM performs anterpolation from a child to the parent domain.

4.3 Turbulence characteristics

Since the flow is driven by the pressure gradient instead of the Coriolis force, the flow is aligned with the x axis, and the wind direction remains nearly constant. The fluctuations of the lateral component v are stronger for the measurement time series. Therefore, we compare turbulence statistics of the horizontal wind speed u from the LES results to the total horizontal flow in the measurements $U=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}$ and omit the lateral component v for the LES data.

Figure 7 Spectra for the horizontal velocity u at the height z=80 m. (a) NBL case, (b) CBL case, and (c) SBL case.

Figure 8 Spectra for the vertical velocity w at the height z=80 m. (a) NBL case, (b) CBL case, and (c) SBL case.

In one-way nested simulations, the parent domain does not receive feedback from the child domain. Consequently, the spectral characteristics of non-nested domains with the grid spacing of Δ_x=10 m (NBL and CBL) and 5 m (SBL) match the characteristics of the corresponding domain in a one-way nesting simulation (Figs. 7 and 8). The individual spectra of the nested domains lie apart from each other but show improvement as the grid spacing is reduced. The inertial subrange resolved by LES widens as the grid becomes more refined; however, it is not fully resolved despite the grid spacing being reduced to Δ_x = 1.25 m.

The two-way nesting mode ensures feedback between the nested domains. Therefore, the root and child domain spectra lie closer to each other and to the one-way spectra of the most refined child domain (Δ_x = 1.25 m). Despite the exchange between domains in the two-way nested case, the spectral characteristics do not coincide perfectly. The inertial subrange being shorter for Δ_x=10 m than for the refined domains implies that the grid resolution is the limiting factor, and the solution for the root domain cannot be improved further even in the two-way nesting case.

Despite the NBL case being simulated as a true neutral condition, it showed good agreement with the measurements on par with the CBL case. The result suggests that it is possible to omit a weak heat flux in neutral cases to save computational time and avoid secondary circulation in the two-way nesting mode.

The SBL simulations largely overestimate the energy contained in low-frequency eddies. The inertial subrange of the corresponding measurement time series also starts at higher frequencies, unlike in the NBL and CBL cases. The LES does not fully resolve high frequencies despite gradually reduced grid spacing. Hence the overall agreement for the SBL case is worse than for NBL and CBL. When comparing available measurement profiles for the specific period of SBL time series, we did not observe anomalies or irregularities, such as reported by Kettle (2014), which could be studied as a possible cause of a discrepancy. The existing studies on SBL simulations with PALM (Beare et al., 2006; Wurps et al., 2020) do not compare simulated spectra against measurements but evaluate other aspects, such as fluxes and grid resolution influence. Hence, simulating SBL in PALM may require additional studies focusing on turbulence characteristics.

In order to match the SBL spectra shape, we performed a short SBL simulation with lower forcing, which led to a decreased turbulence intensity but stronger mean profile shear. The results are provided in Appendix.

The coherence, co-coherence, and phase are plotted against the reduced frequency:

$$\begin{array}{}\text{(8)}& f_{\mathrm{r}}={\displaystyle \frac{f\mathit{\delta }}{\stackrel{\mathrm{‾}}{u}}},\end{array}$$

where f is the original frequency, δ is the vertical separation distance, and $\stackrel{\mathrm{‾}}{u}$ is the mean wind speed of the two regarded levels, 60 and 80 m for δ=20 m or 40 and 80 m for δ=40 m.

The coherence and co-coherence calculated for NBL and CBL coarse domains (Δ_x = 10 m) and δ=20 m show strong deviation from the measurements for the one-way and non-nested simulations at f_r>1 (Figs. 9a and 10a). The tendency to the coherence/co-coherence value of 0.5 suggests that the time series at points separated by δ=20 m remain partially correlated in the coarse grid, which is not the case for the corresponding measurements. While the most refined child domain (Δ_x = 1.25 m) shows a good match between the LES and measurement series (Figs. 9b and 10b), the agreement already improves for Δ_x = 5 m, and the correlation falls to zero for f_r>0.5. The SBL case shows better agreement for the root domain because of the lower initial grid spacing Δ_x = 5 m. Nevertheless, the coherence is noticeably overestimated for low f_r compared to the measurements (Fig. 9a, b). The time series are generally uncorrelated for the vertical separation of δ=40 m both for the LESs and measurements (Figs. 9c, d and 10c, d). However, the NBL case does not capture the high coherence value at f_r=0 observed in the measurements.

Figure 9 Coherence for the horizontal velocity u and different stability cases. (a) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 20 m. (b) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 20 m. (c) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 40 m. (d) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 40 m.

Figure 10 Co-coherence for the horizontal velocity u and different stability cases. (a) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 20 m. (b) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 20 m. (c) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 40 m. (d) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 40 m.

The phase plots are in line with the coherence. The time series are in phase for f_r<0.1, where the coherence is above zero. The effect is strong for the low vertical separation of δ=20 m (Fig. 11a, b) and is in good agreement with the measurements. The phase becomes more chaotic as the vertical separation distance increases to δ=40 m (Fig. 11c, d), while the time series become less correlated (Figs. 9c, d and 10c, d).

Figure 11 Phase plot for the horizontal velocity u and different stability cases and domains. (a) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 20 m. (b) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 20 m. (c) Root domain (Δ_x = 10 m for NBL and CBL, Δ_x = 5 m for SBL), vertical separation δ = 40 m. (d) Innermost child domain (Δ_x = 1.25 m, all cases), vertical separation δ = 40 m.

4.4 Other flow characteristics

4.4.1 Power law

In general, the power law coefficient follows the known trend, also observed in the measurement profile fits (Table 7): high value in the stable layer and low value in the convective layer (Touma, 1977). The discrepancy between exact values of α in measurement and simulated fits is primarily caused by the different way of obtaining U₁₀. For sonic data, U₁₀ is calculated from the previously estimated profile Eq. (1). The LES returns the full mean profile on the pre-defined grid, so U₁₀ can be interpolated to the level of z=10 m. U₁₀ derived from LES data consistently deviates from measurements U₁₀ by 10 %–20 %, thus affecting the estimation of the power law exponent.

Table 7 Estimated power law coefficient.

The estimated power law coefficient α shows little variation for the NBL and CBL domains of the same refinement but implies high sensitivity of the SBL profiles. Considering higher shear in the SBL profiles, the grid refinement may affect the estimation of U₁₀ more strongly than lower shear NBL and CBL profiles.

4.4.2 Turbulence anisotropy

The anisotropy estimation captures only general trends seen in the measurements with the nesting modes being radically different between each other (Fig. 12). Since the inertial subrange resolved in a one-way nested root domain is slightly shorter than of a two-way root domain (Figs. 7–8), f_n≈1 may fall outside of the resolved subrange and provide a less precise estimation. The two-way nested cases approach closer to the anisotropy seen in the measurement, although the anisotropy strength may not match the value seen in the measurement data. The divergence is particularly strong for the SBL simulation, primarily caused by the differences in power density spectra discussed in Sect. 4.3.

Figure 12 Comparison of anisotropy across the regarded stability and nesting cases. The color map is centered at the value $\mathrm{4}/\mathrm{3}=\mathrm{1.333}$.

5 Conclusions

We performed nested LES of three stability cases for the horizontal mean wind speed of 12–13 m s⁻¹ at the reference height of 119 m. The simulations were verified by comparing turbulence characteristics to the corresponding measurement time series. The comparison showed that the grid spacing of Δ_x = 10 m was insufficient for NBL and CBL simulations; the spectral and coherence characteristics had improved their agreement with the measurements after the spacing was reduced to Δ_x = 5 m via nesting or a refined non-nested domain simulation. The inertial subrange was not fully resolved despite further refinement and remained narrower than for the measurement time series even at Δ_x = 1.25 m.

We confirmed that the nesting mode does not affect the true neutral simulation, unlike when the temperature equation is solved along with other prognostic equations for CBL and SBL conditions. In the case of CBL or SBL, the flow inside the child domain differed for the one-way and two-way nesting. The two-way nested simulation produced a secondary circulation resulting in a decreased velocity and increased turbulence intensity in the child domains. Due to a strong horizontal shear, the irregularities in lateral and vertical velocity profiles were spread non-uniformly; e.g., the downward flow was stronger at the exit of the nested domain. The horizontal flow accelerated after leaving the nested area so that the mass conservation law was not violated eventually. Unlike the existing research on buoyancy-driven flows, the two-way nesting effects in a shear-driven flow emerged in the first hour of the LES and did not dissipate as the simulation proceeded for 3 more hours.

In theory, the two-way nesting is a good option to refine the grid in the area of interest of a non-homogeneous flow, e.g., wind turbine wakes, as the feedback between parent and child domain allows accounting for the irregularities after the flow exits the nested area. However, the fast development of a secondary circulation in the shear-driven flow limits the two-way nesting application strictly to the true neutral condition. The one-way nested simulation did not add anomalies to the flow; each child domain only refined the grid spacing and resolved small turbulence scales. We, therefore, recommend using the one-way nesting mode for the wind turbine wake simulation. In the case when the two-way nesting mode is preferable, only a true neutral setup does not produce secondary circulation.

Appendix A: SBL simulation with reduced forcing

We performed a test simulation of an SBL precursor for the same wind speed but weaker pressure gradient (−0.0001 Pa m⁻¹ instead of −0.0005 Pa m⁻¹) and slightly stronger surface cooling (−0.3 K s⁻¹ instead of −0.2 K s⁻¹). As a result of the decreased forcing, the developed profiles deviated from the reference measurements and showed stronger shear but lower turbulence intensity (Fig. A1). Due to the computational time constraints we simulate only a non-nested main run for a comparison of spectral characteristics. We observe a better agreement with the measurements spectra (Fig. A2), especially in the w component, whose spectrum does not follow a $-\mathrm{5}/\mathrm{3}$ theoretical slope. Therefore, we are able to match only one of two – either SBL profiles or SBL spectra – and observe a strong discrepancy in another.

Figure A1 Precursor run profiles with original and reduced pressure forcing. (a) Horizontal flow mean profile and (b) turbulence intensity profile.

Figure A2 Main run spectra with original and reduced pressure forcing. (a) Horizontal velocity spectrum and (b) vertical velocity spectrum.