Spatial and temporal variability of turbulence dissipation rate in complex terrain

To improve parameterizations of the turbulence dissipation rate (ϵ) in numerical weather prediction models, the temporal and spatial variability of ϵ must be assessed. In this study, we explore influences on the variability of ϵ at various scales in the Columbia River Gorge during the WFIP2 field experiment between 2015 and 2017. We calculate ϵ from five sonic anemometers all deployed in a ∼4 km² area as well as from two scanning Doppler lidars and four profiling Doppler lidars, whose locations span a ∼300 km wide region. We retrieve ϵ from the sonic anemometers using the second-order structure function method, from the scanning lidars with the azimuth structure function approach, and from the profiling lidars with a novel technique using the variance of the line-of-sight velocity. The turbulence dissipation rate shows large spatial variability, even at the microscale, especially during nighttime stable conditions. Orographic features have a strong impact on the variability of ϵ, with the correlation between ϵ at different stations being highly influenced by terrain. ϵ shows larger values in sites located downwind of complex orographic structures or in wind farm wakes. A clear diurnal cycle in ϵ is found, with daytime convective conditions determining values over an order of magnitude higher than nighttime stable conditions. ϵ also shows a distinct seasonal cycle, with differences greater than an order of magnitude between average ϵ values in summer and winter.

1 Introduction

Numerical weather prediction models currently assume that the generation of turbulence within a grid cell is equal to the dissipation of turbulence ϵ within the same grid cell. While this assumption, which is appropriate for homogeneous and stationary flow (Albertson et al., 1997), can generally be considered valid when adopting a coarse grid (Lundquist and Chan, 2007; Mirocha et al., 2010), it breaks down when using models with finer horizontal resolution (Nakanishi and Niino, 2006; Krishnamurthy et al., 2011; Hong and Dudhia, 2012), as turbulence can be advected to a different grid cell before being dissipated. However, the scales at which the assumption of local equilibrium is no longer valid are currently not well understood, nor is how different atmospheric and topographic conditions can impact the development and decay of turbulent structures.

A more accurate representation of turbulence is crucially needed as it represents the fundamental process to transfer heat, momentum, and moisture in the atmospheric boundary layer (Garratt, 1994). Moreover, turbulence controls a wide range of processes with a direct effect on our socioeconomic activities: turbulence has impacts on forest fire development and propagation (Coen et al., 2013), it affects air traffic control with its influence on aviation meteorology and the dissipation of aircraft vortices (Gerz et al., 2005; Thobois et al., 2015), it determines the characteristics and impacts of pollutant dispersion (Huang et al., 2013), and it affects wind energy production and the lifetime of wind turbines themselves (Kelley et al., 2006). Moreover, the turbulence dissipation rate has been shown to have a primary role in the formation of frontal structures (Piper and Lundquist, 2004), the evolution of cyclones (Bister and Emanuel, 1998), and the development of flows in urban areas and other canopies (Baik and Kim, 1999; Lundquist and Chan, 2007). The precision of wind energy forecasting is also highly impacted by the accuracy of the representation of the turbulence dissipation rate. A recent sensitivity study (Yang et al., 2017; Berg et al., 2018) showed that up to 50 % of the variance in the turbine-height wind speed predicted by the Weather Research and Forecasting model (Skamarock et al., 2005) in complex terrain only depends on the accuracy of the parameterization of the turbulence dissipation rate.

Various techniques have been developed to calculate ϵ from different instruments. In general, all the proposed methods are based on the turbulence theory by Kolmogorov (1941), which represents the decay of turbulence eddies as an energy cascade in the inertial subrange, until the length scales are small enough for the turbulence kinetic energy to be dissipated by molecular diffusion in the viscous subrange. Turbulence dissipation can be calculated from sonic anemometers on meteorological towers (Champagne et al., 1977; Oncley et al., 1996) and super-high-frequency hot-wire anemometers suspended on tethered lifting systems (Frehlich et al., 2006; Lundquist and Bariteau, 2015) or flown on aircraft (Fairall et al., 1980) or UAVs (Lawrence and Balsley, 2013). Remote sensing instruments can provide additional insights into our understanding of turbulence dissipation by combining measurements at greater altitudes with their ease of deployment in complex terrain, despite their potential drawbacks of limited temporal frequency and inherent volume averaging (Frehlich and Cornman, 2002; Wang et al., 2016). Wind profiling radars (Shaw and LeMone, 2003; McCaffrey et al., 2017a), profiling lidars, and scanning lidars have all been successfully used to obtain turbulence measurements. For lidars, different approaches have been developed to retrieve ϵ: width of the Doppler spectrum (Smalikho, 1995; Banakh et al., 1995), line-of-sight velocity spectrum (Drobinski et al., 2000; O'Connor et al., 2010; Bodini et al., 2018), structure function (Frehlich, 1994; Banakh et al., 1996; Banakh and Smalikho, 1997; Smalikho et al., 2005; Frehlich et al., 2006; Wulfmeyer et al., 2016; Smalikho and Banakh, 2017), and range-gate filtering with a sub-grid-scale parameterization scheme (Krishnamurthy et al., 2010).

Here, we retrieve the turbulence dissipation rate from 11 instruments in a complex terrain region, thus building one of the widest observational assessments of ϵ to date. We explore how topography triggers the variability of ϵ at various temporal and spatial scales. We describe the WFIP2 field campaign in Sect. 2, and we define the characteristics of the sonic anemometers and wind profiling and scanning lidars that we use to estimate ϵ. We also describe the methods used to retrieve ϵ from the different instruments, and we further refine and extend a novel approach to derive ϵ from wind profiling lidars. In Sect. 3 we present the spatial variability of ϵ at both the microscale and mesoscale by comparing the estimates from multiple instruments in different locations, with particular attention to the impact that topography has on the spatial evolution of ϵ. In doing so, we also assess the climatology of turbulence dissipation in terms of both diurnal and seasonal cycles. Section 4 summarizes our results and suggests future work to further improve our understanding and representation of the turbulence dissipation rate in the boundary layer.

2 Data and methods

2.1 The WFIP2 field campaign

The Second Wind Forecast Improvement Project (WFIP2) (Shaw et al., 2019), which involved a field campaign (Wilczak et al., 2019) in the US Pacific Northwest between October 2015 and March 2017, was designed to improve numerical weather prediction model forecasts in complex terrain for wind energy applications. A large number of instruments was deployed in the Columbia River Gorge and basin in a region over 500 km wide. In this study, we focus on the evaluation of turbulence dissipation rate from instruments that span an approximately 300 km wide area. Two profiling lidars were located at the western and eastern edges of this region, at Troutdale (the only site on the western side of the Cascades) and Vansycle Ridge, respectively, with an additional scanning lidar located in Boardman (Fig. 1a). A region with a high density of instruments (HD Site in Fig. 1a), approximately ∼20 km wide, was located in the vicinity of the town of Wasco, from which we will analyze data from two wind profiling lidars, one scanning lidar (Fig. 1b), and the sonic anemometers on five meteorological towers (Fig. 1c). Extensive arrays of wind turbines are located on the northern side of the Columbia River and on the southwestern part of the studied region.

Figure 1 Map of the relevant instruments during the WFIP2 field campaign. The locations of the profiling lidars, scanning lidars, and meteorological towers used in this analysis are shown. (c) Elevation contour lines are shown every 10 m.

Multiple sonic anemometers were located on several meteorological towers at the Physics Site (Wilczak et al., 2019), which represented the finest array of instruments at WFIP2, aimed at having multiple measurements in an area similar in size to a grid cell of a high-resolution numerical weather prediction model. The site, covered by crop fields, is characterized by a moderately complex topography, with terrain elevation spanning from 405 m to 459 m a.s.l. (the elevations of the locations of the meteorological towers used in this study are reported in Table 1). Several wind turbines are located east of the Physics Site. The sonic anemometers used in this project provide 20 Hz measurements of the three components of the wind and virtual temperature at 10 m a.g.l., and they were operational from late March–early April 2016 to late April–mid-May 2017. To account for tower wake effects, data were excluded when the wind direction was within $\pm \mathrm{30}{}^{\circ }$ from the orientation of the tower boom (McCaffrey et al., 2017b). Less than 10 % of the data was excluded due to tower wake contamination.

Table 1 Elevation and period of data collection for the five 10 m sonic anemometers at the Physics Site, the four profiling lidars, and the two scanning lidars considered in this study, whose locations are shown in Fig. 1.

Data from the sonic anemometers were used to assess atmospheric stability calculated in terms of the Obukhov length L:

$$\begin{array}{}\text{(1)}& L=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v}}}\cdot u_{*}^{\mathrm{3}}}{k\cdot g\cdot \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}}},\end{array}$$

where θ_v is the virtual potential temperature (K) calculated from the virtual temperature measured by the sonic anemometers, $u_{*}=({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}$ is the friction velocity (m s⁻¹), k=0.4 is the von Kármán constant, g=9.81 m s⁻² is the acceleration due to gravity, and $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ is the kinematic sensible heat flux (m K s⁻¹). An averaging period of 30 min (De Franceschi and Zardi, 2003; Babić et al., 2012) has been used to apply the Reynolds decomposition and determine the fluxes. Based on the value of the Obukhov length, we classify neutral conditions for $L\le -\mathrm{500}$ m and L>500 m, unstable conditions for −500 m < L≤0 m, and stable conditions for 0 m $<L\le \mathrm{500}$ m (Muñoz-Esparza et al., 2012). Neutral conditions were infrequently recorded (less than 7 % of the time).

A WINDCUBE version 1 (v1) was located in Troutdale (12 m a.s.l.), about ∼20 km east of Portland, in a relatively flat region at the Portland–Troutdale Airport at the western edge of the Columbia River Gorge. The area is semi-urban, with some trees. This type of lidar (Aitken et al., 2012; Rhodes and Lundquist, 2013) measures line-of-sight velocity along four cardinal directions with a nominal zenith angle of 28^∘ and a temporal resolution of about 1 Hz along each beam direction. The measurements are taken every 20 m from 40 to 220 m a.g.l. The main technical specifications of the instrument are shown in Table 2.

A second WINDCUBE v1 and a WINDCUBE 200S scanning Doppler lidar were located at the Wasco Airport (456 m a.s.l.), in an area within the Columbia basin covered by short grass. The nearby region is characterized by moderately complex topography in the vicinity of the Columbia River. The WINDCUBE 200S performed a variety of planned position indicator (PPI), range-height indicator (RHI), and vertical stare scans within 15 min. Details on the scan patterns can be found in Choukulkar (2018). For this instrument we retrieve ϵ up to 300 m a.g.l..

A WINDCUBE version 2 (v2) was deployed on a low-grass surface on the top of Gordon Ridge (728 m a.s.l.), on the eastern side of the Cascades. A second v2 was deployed at Vansycle Ridge (542 m a.s.l.), in a site with grazed grass (Yang et al., 2013) about 20 km east of the Wallula Gap, where the Columbia River turns north, thus modifying the main topographic direction of the gorge. Compared to the WINDCUBE version 1, the v2 performs an additional line-of-sight velocity measurement along the vertical, and ∼4 s is required for the beam to complete the five-point scan strategy. The vertical resolution was of 20 m from 40 to 260 m a.g.l. (200 m a.g.l. for the v2 at Vansycle Ridge). Bodini et al. (2018) compared turbulence dissipation retrievals from colocated v1 WINDCUBEs and a WINDCUBE v2 and found a good agreement between the different instruments. Table 2 illustrates the major technical parameters of this lidar.

Finally, a Halo Streamline scanning Doppler lidar was deployed near a regional airport surrounded by farmland at Boardman (112 m a.s.l.). The long-range fiber-optic-based scanning Doppler lidar provides 3-D scanning capabilities and performed a wide range of scans covering the atmospheric boundary layer over a period of 15 min (Otarola, 2017). In this analysis, only the 5^∘ elevation angle scans with a scan rate of 1^∘ s⁻¹ were used to calculate the turbulence dissipation rate up to 120 m a.g.l. The other scans within the 15 min time period were not usable for turbulence calculations due to either fast scan rates (Frehlich et al., 2006) or low data availability.

For all the instruments, precipitation periods were excluded from the analysis based on measurements at two surface meteorological stations at the Wasco Airport and Troutdale (for the profiling lidar at that location).

Table 2 Main technical specifications of the lidars used in this study.

2.2 Turbulence dissipation rate from sonic anemometer

We estimate the turbulence dissipation rate from the sonic anemometers using the second-order structure function method, which has been demonstrated (Muñoz-Esparza et al., 2018) to provide ϵ retrievals with a lower error compared to the commonly used inertial-subrange energy spectrum method. The second-order structure function D_U of the horizontal velocity U at the position x is defined as a function of the spatial separation r as $D_U\left(r\right)\equiv <\left[U\right(x+r)-U(x){]}^{\mathrm{2}}>$, where $<\cdot >$ denotes an ensemble average. Within the inertial subrange, Kolmogorov's model (Kolmogorov, 1941; Frisch, 1995) relates the second-order structure function with the turbulence dissipation rate ϵ:

$$\begin{array}{}\text{(2)}& D_U\left(r\right)={\displaystyle \frac{\mathrm{1}}{a}}{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{r}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where a is the Kolmogorov constant, which we set equal to 0.52 (Paquin and Pond, 1971; Sreenivasan, 1995). By invoking Taylor's frozen turbulence hypothesis (Taylor, 1935), the spatial separation r can be written as temporal separation τ so that ϵ can be calculated as

$$\begin{array}{}\text{(3)}& \mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{U\mathit{\tau }}}{\left[a{D}_U\left(\mathit{\tau }\right)\right]}^{\mathrm{3}/\mathrm{2}}.\end{array}$$

We calculate ϵ every 30 s by fitting the Kolmogorov's theoretical model to the structure function calculated from the sonic anemometer data using a temporal separation between τ=0.1 and τ=2 s. From data inspection, measurements in the chosen time separation interval lie well within the inertial subrange, and therefore they fulfill the hypothesis of Kolmogorov's theory. Moreover, the high temporal resolution of the sonic anemometer suggests an adequate number of data points in this interval to obtain a robust estimate of the structure function.

2.3 Turbulence dissipation rate from wind profiling lidar

Measurements from wind Doppler lidars can extend our understanding of the variability of the turbulence dissipation rate thanks to their relatively easy deployment, even in prohibitive terrain conditions. Moreover, lidars can often provide measurements at higher altitudes compared to most meteorological towers, possibly out of the surface layer.

We follow the approach introduced by O'Connor et al. (2010) and refined by Bodini et al. (2018) to estimate ϵ from the variance of the line-of-sight velocity measured by the profiling lidars. Assuming locally homogeneous and isotropic turbulence, the one-dimensional spectrum S within the inertial subrange can be written as a function of the wave number k as

$$\begin{array}{}\text{(4)}& S\left(k\right)=a{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{k}^{-\mathrm{5}/\mathrm{3}},\end{array}$$

where a=0.52 is the one-dimensional Kolmogorov constant. By integrating Eq. (4) over the wavenumber space within the inertial subrange, the following expression can be found:

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}& {\displaystyle }=\underset{k}{\overset{k_{\mathrm{1}}}{\int }}S\left(k\right)\mathrm{d}k=-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}a{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}\left(k_{\mathrm{1}}^{-\mathrm{2}/\mathrm{3}}-k^{-\mathrm{2}/\mathrm{3}}\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}{\left({\displaystyle \frac{\mathit{ϵ}}{\mathrm{2}\mathit{\pi }}}\right)}^{\mathrm{2}/\mathrm{3}}\left(L_N^{\mathrm{2}/\mathrm{3}}-L_{\mathrm{1}}^{\mathrm{2}/\mathrm{3}}\right),\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}$ is the variance (averaged across the different beams) of the detrended line-of-sight velocity, and L₁ and L_N are the length scales that can be used instead of the wavenumbers by invoking Taylor's frozen turbulence hypothesis (Taylor, 1935). For a single sample, L₁ can be defined as

$$\begin{array}{}\text{(6)}& L_{\mathrm{1}}=Ut+\mathrm{2}z\sin \left({\displaystyle \frac{\mathit{\theta }}{\mathrm{2}}}\right),\end{array}$$

where U is the horizontal wind speed, t is the dwell time, θ is the half-angle divergence of the lidar beam, and z the height above ground level. The second term in Eq. (5) can typically be neglected as Doppler lidars generally have θ<0.1 mrad. For multiple samples, L_N=NL₁, where N is the number of samples used in the calculation.

The method relies on the fundamental assumption that the samples used in the calculation lie within the inertial subrange of turbulence. If longer samples are used, therefore including contributions from the outer scales, ϵ will be severely underestimated (Bodini et al., 2018). On the other hand, short samples will undermine the representativeness of the estimation of the turbulence contribution to variance (Lenschow et al., 1994), and a higher relative effect of instrumental noise (Lenschow et al., 2000) will also increase the error. Therefore, the choice of the sampling size N represents a crucial step to obtain accurate estimates of turbulent quantities, especially in stable conditions (Pichugina et al., 2008).

As shown in Bodini et al. (2018), the appropriate timescales for the lidar retrievals can be determined in different ways. When colocated sonic anemometers are available, the optimal values for N can be found by tuning the lidar method with the ϵ values derived from the sonic data. Another possibility is the use of spectral models (Kaimal et al., 1972; Panofsky, 1978; Olesen et al., 1984; Kristensen et al., 1989) to fit the experimental spectra from the lidar measurements and determine the extension of the inertial subrange from the fit (Tonttila et al., 2015).

In the WFIP2 case, no sonic anemometers colocated with the profiling lidars were available. Moreover, all the WINDCUBE lidars at WFIP2 operated in profiling mode using slant beams rather than in a purely vertical stare mode. Therefore, modeling the spectra of the line-of-sight velocity measured by these instruments is not trivial, as most of the spectral models are valid for either the purely horizontal or vertical components of wind speed, and projecting these models can lead to variance contamination (Newman et al., 2016). As a consequence, we further extend this method and we estimate the optimal sample length N to use in the retrieval of ϵ by determining the extension of the inertial subrange as the maximum in the curve representing a local regression of the spectrum of the line-of-sight velocity measured by the lidars. In doing so, we do not need to know the precise functional form for the spectrum of the measured radial velocity in an arbitrary slant direction. Using the dataset described in Bodini et al. (2018), with sonic anemometers colocated with lidars, we tested different local regression techniques, and we select the robust LOESS technique (Cleveland, 1979), with a span of 15 % of the total number of data points in each spectrum, which provided the best agreement (R²>0.95) with the ϵ values obtained from the fine-tuning with the estimates from the sonic anemometers. In the determination of the maximum of the local regression curve, we leave out frequencies greater than 0.05 Hz, which are most affected by instrumental noise (Frehlich, 2001). The distribution of sample size values we obtain is between 20 s (5th percentile) and 300 s (95th percentile). An example of the local regression of an experimental lidar spectrum at WFIP2 is shown in Fig. 2.

Figure 2 Example of the local regression of an experimental spectrum of the line-of-sight velocity measured by one of the four beams of the WINDCUBE v1 lidar at Wasco Airport at WFIP2. The red dashed line shows the maximum of the local regression curve. The orange dashed line shows the theoretical $-\mathrm{5}/\mathrm{3}$ slope of the spectrum in the inertial subrange.

Finally, the contribution due to instrumental noise needs to be considered. The observed variance ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}$ in Eq. (5) can be thought of as a combination of three different contributions, which can be considered as independent of one other (Doviak et al., 1993):

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{w}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{w}}^{\mathrm{2}}$ is the contribution from atmospheric turbulence at the scales the lidar can measure (Brugger et al., 2016), ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}$ is due to the instrumental noise, and ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ is related to the variation in the aerosol terminal fall velocity within the sampled volume, which can safely be ignored since the particle fall speed is typically very low (<1 cm s⁻¹). The contribution of instrumental noise ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}$ can be written as a function of the signal-to-noise ratio (SNR) (Pearson et al., 2009):

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{\Delta }{\mathit{\nu }}^{\mathrm{2}}\sqrt{\mathrm{8}}}{\mathit{\alpha }{N}_{\mathrm{p}}}}{\left(\mathrm{1}+{\displaystyle \frac{\mathit{\alpha }}{\sqrt{\mathrm{2}\mathit{\pi }}}}\right)}^{\mathrm{2}},\end{array}$$

where Δν is the signal spectral width, and α is the ratio of the lidar photon count to the speckle count (Rye, 1979), which can be calculated as a function of the bandwidth B as $\mathit{\alpha }=\frac{\mathrm{SNR}}{\sqrt{\mathrm{2}\mathit{\pi }}}\frac{B}{\mathrm{\Delta }\mathit{\nu }}$. The accumulated photon count N_p can be calculated as N_p=SNRnM, with n the number of lidar pulses that are averaged to get a profile and M the number of points sampled within a single range gate. Therefore, ϵ can be determined as

$$\begin{array}{}\text{(9)}& \mathit{ϵ}=\mathrm{2}\mathit{\pi }{\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}\mathrm{a}}}\right)}^{\mathrm{3}/\mathrm{2}}{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}}{L_N^{\mathrm{2}/\mathrm{3}}-L_{\mathrm{1}}^{\mathrm{2}/\mathrm{3}}}}\right)}^{\mathrm{3}/\mathrm{2}},\end{array}$$

with the accurate choice of the appropriate sample length N, as described.

2.4 Turbulence dissipation rate from scanning Doppler lidar

The turbulence dissipation rate from the scanning Doppler lidars is estimated using the azimuth structure function method (Frehlich et al., 2006; Krishnamurthy et al., 2011). The structure function from the radial velocity estimates can be used to retrieve the turbulence dissipation rate, the integral length scale, and the velocity variance, assuming a theoretical model for isotropic wind fields. In our approach, corrections for turbulence measurements have been considered to address the complications due to the inherent volumetric averaging of radial velocity over each range gate, the noise of the lidar data, and the assumptions required to estimate the effects of smaller scales of motion on turbulence quantities.

Both the scanning lidars have an azimuth scan rate of 1^∘ s⁻¹; the Halo Streamline has an accumulation time of 1 s, while the WINDCUBE 200S at Wasco has as a time of 0.5 s.

The structure function ${\widehat{D}}_{\mathrm{wgt}}$ of the mean Doppler lidar velocity perturbations, $\widehat{v^{\prime }}$, in the azimuth direction is given by

$$\begin{array}{ll}\text{(10)}& {\displaystyle }& {\displaystyle }{\widehat{D}}_{\mathrm{wgt}}(R,kR\mathrm{\Delta }\mathrm{\Phi },\mathit{\theta })={\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}-k}}\sum _{j=\mathrm{1}}^{N_{\mathrm{s}}-k}\left[\widehat{v^{\prime }}(R,(j-\mathrm{1})\mathrm{\Delta }\mathrm{\Phi },\mathit{\theta })\right.\\ {\displaystyle }& {\displaystyle }{\left.-\widehat{v^{\prime }}(R,(j+k-\mathrm{1})\mathrm{\Delta }\mathrm{\Phi },\mathit{\theta })\right]}^{\mathrm{2}}-\mathrm{2}{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}\left(R\right),\end{array}$$

where ΔΦ is the azimuth angular spacing between adjacent Doppler velocity estimates, and N_s is the number of velocity measurements for the sector scan. The estimation error is uncorrelated with the pulse-weighted velocity because each estimate is produced with different lidar pulses (assuming no multi-scattering effects); therefore, the velocity error variance ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{2}}\left(R\right)$ is only a function of the range gate (Krishnamurthy, 2008).

For homogeneous von Kármán turbulence over a two-dimensional plane, the following model (Hinze, 1959; Frehlich et al., 2006) for the structure function is valid:

$$\begin{array}{}\text{(11)}& D_{\mathrm{v}}(r,s)=\mathrm{2}{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}\left[\mathrm{\Lambda }\left({\displaystyle \frac{p}{L_{\mathrm{o}}}}\right)+{\mathrm{\Lambda }}_D\left({\displaystyle \frac{p}{L_{\mathrm{o}}}}\right)\left(\mathrm{1}-{\displaystyle \frac{r^{\mathrm{2}}}{p^{\mathrm{2}}}}\right)\right],\end{array}$$

where r denotes the distance along a fixed laser beam, $s=R({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}})$ is the transverse coordinate, $p=(r^{\mathrm{2}}+s^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ , L_o is the outer scale of turbulence, which is proportional to the integral length scale L_i, Λ(x) is the universal function, and

$$\begin{array}{}\text{(12)}& {\mathrm{\Lambda }}_D\left(x\right)={\displaystyle \frac{x^{\mathrm{4}/\mathrm{3}}}{{\mathrm{2}}^{\mathrm{1}/\mathrm{3}}\mathrm{\Gamma }(\mathrm{1}/\mathrm{3})}}{K}_{\mathrm{2}/\mathrm{3}}\left(x\right)=\mathrm{0.3}{x}^{\mathrm{2}/\mathrm{3}}{K}_{\mathrm{2}/\mathrm{3}}\left(x\right).\end{array}$$

Assuming that the averaged radial velocity can be written as a function of the instantaneous radial velocity and an effective spatial filter in terms of the pulse-weighting function and range-gate length of the lidar (Frehlich et al., 2006), the Doppler lidar azimuth structure function can be modeled as

$$\begin{array}{}\text{(13)}& D_{\mathrm{wgt}}(s,\mathit{\sigma },L_{\mathrm{o}})=\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}{G}_{\mathrm{a}}\left({\displaystyle \frac{s}{\mathrm{\Delta }\mathrm{p}}},\mathit{\mu },\mathit{\zeta }\right),\end{array}$$

where $s=R({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}})$, σ is the standard deviation of the transverse velocity fluctuations, and $G_{\mathrm{a}}(\mathit{\eta },\mathit{\mu },\mathit{\zeta })$ is the derived model based on weighted velocity estimates and the von Kármán model, as provided in Eq. (46) of Frehlich et al. (2006) and fully derived in Krishnamurthy (2008).

The parameters σ and L_o are estimated by minimizing the error between the lidar-derived structure function ${\widehat{D}}_{\mathrm{wgt}}(R,kR\mathrm{\Delta }\mathrm{\Phi },\mathit{\theta })$ and the model estimates ${\widehat{D}}_{\mathrm{wgt}}(s,\mathit{\sigma },L_{\mathrm{o}})$. The dissipation rate can then be estimated by (Hinze, 1959)

$$\begin{array}{}\text{(14)}& \mathit{ϵ}=\left(\mathrm{0.933}\right){\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{3}}}{L_{\mathrm{o}}}}.\end{array}$$

Although the assumption of homogeneous isotropic turbulence is not valid for every condition, the effect of anisotropy on the azimuth structure function is small (Krishnamurthy et al., 2011). Therefore, with an accurate choice of the scan angle and vertical resolution, the isotropic assumption can be relaxed in this algorithm for complex terrain applications. Using the selected scans described in the previous section, we retrieve ϵ from the WINDCUBE 200S and the Halo Streamline lidars every 15 min.

3 Results and discussion

The turbulence dissipation rate has been retrieved from the five sonic anemometers at the Physics Site, four profiling lidars, and two scanning lidars. This extensive network of measurements at WFIP2 allows for a unique assessment of the spatial and temporal variability, at various scales, of ϵ in complex terrain.

3.1 Microscale variability of turbulence dissipation rate in complex terrain

The analysis of the retrievals of turbulence dissipation rate from the five 10 m sonic anemometers, all located within a ∼4 km² area at the Physics Site (Fig. 1c), allows insight into the microscale variability of ϵ in the surface layer in complex terrain.

To gain first insights on the evolution of ϵ within the Physics Site, a portion of the time series of ϵ and correspondent wind speed from the five sonic anemometers can be analyzed (Fig. 3).

Figure 3 Time series from 24 June 2016 00:00 UTC to 27 June 2016 00:00 UTC comparing ϵ (a) and wind speed (b) from five sonic anemometers at 10 m a.g.l. at the Physics Site. Data have been smoothed with a 30 min running mean. Blue shaded areas show local nighttime conditions, while orange areas show local daytime periods.

The turbulence dissipation rate exhibits variability of at least 3 orders of magnitude over a diurnal cycle, with higher values generally observed during daytime conditions and lower values at night. However, the magnitudes observed in the diurnal cycle of ϵ show considerable variability among different days, with the minimum values during the night of calendar day 176 when high winds were recorded, similar to the maximum magnitudes observed during daytime convective conditions on day 178 when the wind was more quiescent. Moreover, although the five considered towers are all located within a ∼4 km² area, ϵ still shows considerable variability among the different sonic anemometers. This variability is particularly accentuated at night (especially for the night between calendar days 178 and 179), when ϵ varies more than an order of magnitude within the considered microscale region. This variability can be connected to the intermittent nature of the turbulence dissipation rate, for which a multifractal theory has been developed (Frisch, 1995).

Given this distinct variability of ϵ at different times of the day, the impact of atmospheric stability conditions can be additionally investigated throughout the ∼13 months of measurements at the Physics Site. To understand whether a systematic difference in the microscale variability of ϵ during different atmospheric stability conditions can be found, we calculated, at each time, the ratio between ϵ from each sonic anemometer and the average ϵ (at that time) from all five sonic anemometers. We then classified these ratios based on atmospheric stability, quantified as the median value of the Obukhov length from the five sonic anemometers. For each sonic anemometer, we estimate the variability of ϵ in different stability conditions in terms of the standard deviation of the distribution of these ϵ ratios, as reported in Table 3.

Table 3 Standard deviation of the distribution of the ratios between ϵ_i from each sonic anemometer and the average $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ from all five sonic anemometers for different atmospheric stability conditions.

For all five sonic anemometers, the standard deviation is higher during stable conditions compared to unstable conditions, with mean (across the five anemometers) values of 0.83 and 0.74, respectively. On average, in the surface layer, at the small spatial scales sampled within the Physics Site, ϵ shows a 12 % larger variability during nighttime stable conditions compared to daytime convective conditions.

Along with atmospheric stability, topographic features can have an impact on the variability of the turbulence dissipation rate, and the high-density array of meteorological towers at the Physics Site represents an ideal candidate to explore this relation at the microscale. Figure 4 shows the wind rose obtained from the 10 m sonic anemometer on the P03 meteorological tower at the Physics Site (the wind roses from the four other sonic anemometers are qualitatively similar to the one shown here and are reported in the Supplement). The prevailing wind directions at the Physics Site follow the dominant west–east direction of the Columbia River Gorge.

Figure 4 Wind rose computed from the data recorded by the 10 m sonic anemometer on the meteorological tower P03 at the Physics Site from 29 March 2016 to 15 May 2017.

As the wind at the site is almost always slightly southwesterly, it is interesting to study whether differences in turbulence dissipation rate can be found as the wind flows from the western to the eastern sides of the Physics Site. An analysis of the topography of the region reveals two distinct sets of terrain characteristics. The terrain on the west of the subgroup of towers on the western side of the Physics Site (towers P03 and P09) has slopes that reach 60 %, with average slopes larger than 6 %. In contrast, the remaining towers east of this cluster, which we will refer to as “eastern” (towers P04, P05, and P10), are surrounded by a terrain with more gentle slopes, which are on average less than 6 % and never exceed 25 %. We note that the far eastern side of the Physics Site includes an 80 m tower (Wilczak et al., 2019). We can first assess the topographic impact on the microscale variability of ϵ in terms of the distribution of the ratio between the mean ϵ from the groups of sonic anemometers on the two sides of the Physics Site (Fig. 5e).

Figure 5 Histogram of the ratio of the average ϵ retrieved from the three 10 m sonic anemometers on the eastern side of the Physics Site (towers P04, P05, and P10) to the average of ϵ retrieved from the two sonic anemometers on the western side of the Physics Site (towers P03, P09). The vertical dashed line shows the 1.0 ratio, which would indicate no difference, on average, in ϵ between the two sides of the site.

A systematic bias is observed in the values of ϵ on the two sides of the Physics Site, with the median value of turbulence dissipation on the eastern side being only 73 % of the median ϵ on the western side. These differences may be due to the drainage flows and channeling frequently observed at night at this site. The presence of steep topography increases the variability of the turbulence dissipation rate, even at small spatial scales (of the order of 2 km in this case).

To confirm this result, the correlation between ϵ retrievals from all the possible pairs of meteorological towers at the Physics Site can be studied (Fig. 6a).

Stations that are close by (separation <1 km) and on the same side of the Physics Site show high correlation coefficients (R>0.75). When considering pairs of stations on opposite sides of the Physics Site (with separations between 1 and 2 km), we find smaller correlations (R<0.7) for the turbulence dissipation rate, as is reasonable since the spatial separation between the towers increases. However, when looking at the correlation between the retrievals from the two sonic anemometers on the western side of the Physics Site, which have the highest separation (∼2.2 km), we still find a relatively high correlation coefficient (R>0.7). Larger separations do not represent the only dominant factor in determining a progressive reduction of the coefficient of correlation, as the specific interaction between the atmospheric flow and the topographic features in complex terrain seems to be capable of modifying the spatial evolution of correlation between turbulence dissipation at different locations.

The relationship between the correlation coefficient and separation can also provide confirmation of the larger variability of ϵ observed during stable conditions. When calculating the correlation coefficient between ϵ values classified in stable and unstable conditions, calculated in terms of the median value of the Obukhov length (Fig. 6b), we find systematically larger values of R during unstable conditions compared to stable conditions at every spatial separation. During quiescent stable conditions, the increased variability of ϵ even at the microscale determines a reduced correlation throughout the site. On the other hand, when considering the evolution of the correlation coefficient as a function of the elevation difference among the meteorological towers, no systematic trend can be found (plot shown in the Supplement).

Figure 6 Correlation coefficient R between log (ϵ) from different pairs of 10 m sonic anemometers at the Physics Site as a function of the separation between the single meteorological towers. (a) Different colors are used for pairs of towers both on the western side of the Physics Site, both on the eastern side, or on both sides. (b) Data points are classified as a function of atmospheric stability.

Finally, the temporal variability of the turbulence dissipation rate at the microscale can be assessed in terms of the annual cycle of ϵ. Increased daytime convection combined with stronger, on average, winds during the summer causes larger turbulent mixing, which in turn leads to higher values of dissipation rates compared to winter months. Figure 7a quantifies this process by showing how the median value of ϵ varies as a function of the month of the year for each of the five stations at the Physics Site. The annual cycle of wind speed is shown in the Supplement.

Figure 7 Median ϵ value for each calendar month and each considered sonic anemometer (a) and the correspondent interquartile range (b).

ϵ shows a clear annual cycle in the surface layer, with median ϵ values over an order of magnitude larger in summer than winter at all five locations considered within the Physics Site. As a consequence, the interquartile range of ϵ also reveals an annual cycle (Fig. 7b), with a larger range of variability in summer than winter, again with differences of orders of magnitude.

3.2 Mesoscale variability of turbulence dissipation rate in complex terrain

While the analysis of the heavily instrumented Physics Site provides a unique long-term dataset to explore the microscale variability of the turbulence dissipation rate in the surface layer, the four wind profiling lidars and the two scanning lidars allow for an evaluation of the variability of ϵ, at higher altitudes relevant for wind energy, in a region spanning ∼300 km.

The annual cycle in the turbulence dissipation rate found at the Physics Site can also be detected from the retrievals at higher altitude from the lidars at the mesoscale. Figure 8 shows the time series of ϵ from the different lidars at 100 m a.g.l., with a low-pass filter (15-day moving window) applied to filter out the high-frequency and diurnal fluctuations and focus on the seasonal trend. For the lidars at Gordon Ridge and at Vansycle Ridge, which were deployed for more than a year, two time series are plotted for the overlapping calendar days from different years. The time series of the seasonal cycle of wind speed for the different lidars is included in the Supplement.

Figure 8 Low-pass-filtered (with a 15-day moving average) time series of ϵ from the four considered profiling lidars and the two scanning lidars as a function of the calendar day at 100 m a.g.l.

The time series confirm that turbulence dissipation shows a distinct seasonal variability: ϵ is, on average, much higher during the summer, when strong convection increases turbulence production and consequently dissipation. Average ϵ values during winter are about 1 order of magnitude lower than what is observed in summer. Measurement records longer than a single year would be beneficial to filter out possible variations of ϵ linked with specific weather conditions, which, together with snow melting on the ground, possibly impacted the abrupt increase in average ϵ values at Gordon Ridge in the spring.

Moreover, the smoothed time series also reveals how the turbulence dissipation rate at Boardman and Gordon Ridge is, except for winter months, much larger than at the other locations, with the average time series at the other locations showing, on average, almost 1 order of magnitude lower values of ϵ. To explore why ϵ shows much larger values at these locations, Fig. 9 shows the wind roses and the correspondent roses of turbulence dissipation rate at 100 m a.g.l. for the WINDCUBE v2 and the Halo Streamline lidar.

Figure 9 Wind roses at 100 m a.g.l. from the WINDCUBE v2 at Gordon Ridge (a) and the Halo Streamline at Boardman (c) and the correspondent turbulence dissipation roses at the same altitude (b, d).

At Gordon Ridge, westerly winds are the prevailing pattern, with some northeasterly winds being the second most common situation. The highest values for ϵ are measured during westerly wind conditions, while cases with easterly winds rarely have $\mathit{ϵ}>{\mathrm{10}}^{-\mathrm{3}}$ m² s⁻³. When the wind flows from the west, the location of the WINDCUBE v2 lidar is at the easternmost edge of an area (∼6 km wide) with a particularly complex topography, where the Deschutes River (tributary of the Columbia) shapes a steep valley, with terrain slopes that locally exceed 70 % (see map in the Supplement).

Figure 10 Average diurnal climatology of ϵ for the summer (1 June–31 August) from the WINDCUBE v1 at Troutdale (a), the WINDCUBE v1 at the Wasco Airport (b), the WINDCUBE v2 at Gordon Ridge (c), and the WINDCUBE v2 at Vansycle Ridge (d). Average diurnal climatology for the winter (1 December to 28–29 February) from the WINDCUBE v2 at Gordon Ridge (e) and the WINDCUBE v2 at Vansycle Ridge (f). At this site, LST = UTC−8.

With the dominant southwesterly wind, the lidar at Boardman turns out to be located downwind (about 15 km) of a large wind farm. Wind farm wakes are associated with reduced wind speed and increased turbulence (Tennekes and Lumley, 1972), which can have important impacts on wind energy production downwind (Lissaman, 1979; Nygaard, 2014). Wind speed deficits from wind farm wakes have been observed using SAR (Christiansen and Hasager, 2005; Hasager et al., 2006), radars (Hirth et al., 2015), and aircraft measurements (Platis et al., 2018) up to 25 km downwind of the plants. Systematic turbulence measurements that far downwind of wind farms have not yet been made. However, turbulence dissipation measurements 2–3 rotor diameters in the wake of a single turbine (Lundquist and Bariteau, 2015) showed an elevated level of ϵ. Therefore, the increased dissipation aloft observed at Boardman is likely due to the increased turbulence aloft in the wind farm wake. Wind roses and turbulence dissipation roses for the other lidars are included in the Supplement.

The seasonal variability of turbulence dissipation can be additionally investigated by considering the differences in the average daily conditions of ϵ throughout the year. Figures 10 and 11 show the average diurnal climatology of the turbulence dissipation rate at the various locations of the four profiling lidars and the two scanning lidars, respectively. The left column shows the average climatology for the summer, calculated as average conditions from 1 June to 31 August. For the profiling lidars at Gordon Ridge and Vansycle Ridge and the scanning lidar at Wasco, which were also deployed during winter months, the panels on the right show the average daily cycle for the winter using ϵ retrievals from 1 December to the end of February.

Figure 11 Average diurnal climatology of ϵ for the summer (1 June–31 August) from the WINDCUBE 200S at Wasco (a) and the Halo Streamline at Boardman (b). Average diurnal climatology for the winter (1 December to 28–29 February) from the WINDCUBE 200S at Wasco (c). At this site, LST = UTC−8.

For all the lidars, we neglect the heights at which less than 15 % of data within the considered season are available (the complete data availability is shown in the Supplement). In all the locations, turbulence dissipation rate shows a clear diurnal cycle, with higher values during daytime convective conditions and lower values at night, with differences greater than 1 order of magnitude, especially in summer. The intercomparison between the plots from the different lidars also confirms the impact of topography in determining much higher average values of ϵ at Gordon Ridge compared to what is recorded at the other locations. In particular, daytime summer values are about 1 order of magnitude higher than what is found from the other lidars. At Boardman, large average values of dissipation are found aloft at night. In fact, the increased turbulence in the wind farm wake can be further advected during nighttime stable periods, when stronger stratification is found in the boundary layer. The comparison between the summer climatologies (left panels) with the winter ones (right panels) reveals how larger values of ϵ are found during the summer compared to what is found in the wintertime diurnal climatology, when daytime ϵ values are about 2 to 3 orders of magnitude lower, and with a much weaker difference between daytime and nighttime average conditions. It is reasonable to expect that the increased diurnal convection during the summer months determines much stronger turbulent mixing, which in turn causes higher values of turbulence dissipation.

4 Conclusions

Although turbulence is a fundamental transport mechanism in the atmospheric boundary layer, current numerical weather prediction models are limited in their representation of turbulence, for which a local equilibrium between the production and dissipation (ϵ) of turbulence is assumed. The error introduced by the parameterization of ϵ has been shown to be responsible for up to 50 % of the variance of hub-height wind speed predicted by models. Detailed study of observations in the surface layer has great potential for reducing the uncertainty in our understanding of the turbulence dissipation rate. Although methods to retrieve ϵ, at least from in situ measurements, have been known for decades, comprehensive analysis of the spatial and temporal variability of ϵ using data from instruments covering wide regions had not been fully explored to date. In this study we have presented an extensive assessment of the variability, both in space and time, of the turbulence dissipation rate in complex terrain at both the microscale and mesoscale using measurements from both in situ and remote sensing instruments. The impact of topography and other forcings, like large wind farms, on the variability of ϵ has been captured at the different sampled scales.

The turbulence dissipation rate has been calculated from five 10 m sonic anemometers, four wind profiling lidars, and two scanning lidars deployed at the WFIP2 field campaign in the Columbia River Gorge and basin from fall 2015 to spring 2017. The sonic anemometers were all located in an area with an extension of approximately 2 km × 2 km, and they therefore allow for an assessment of the variability of ϵ in the surface layer at the microscale. More homogeneous turbulence across the investigated region is caused by convective mixing during the day. On the other hand, considerable differences (up to 1 order of magnitude) in ϵ are found at night when comparing retrievals of ϵ from the different meteorological towers. On average, ϵ is 12 % more variable during nighttime stable conditions than during unstable convective conditions. Systematic differences emerged from ϵ measured on the western and eastern sides of the Physics Site, the former being located downwind of terrain with larger slopes compared to the latter, thus suggesting the possible impact of terrain slope in triggering the variability of ϵ. The change in correlation between ϵ in different locations is not fully determined purely by spatial separation, as topographic features maintain an importance in influencing it. Therefore, the representation of the turbulence dissipation rate in complex terrain, especially during nighttime stable conditions, needs to be extremely localized to fully capture the turbulence variability in the surface layer.

The variability of ϵ at the mesoscale can be analyzed from the 100 m altitude retrievals from the four wind profiling lidars and the two scanning lidars, which were deployed over a region ∼300 km wide. For the profiling lidars, the retrieval approach proposed in Bodini et al. (2018) has been further refined here and tested to derive ϵ without the need for in situ measurements colocated with the lidars. The profiling lidar located at the topographically complex Gordon Ridge site systematically detected ϵ values which, on average, were over 1 order of magnitude higher than what was measured by the profiling lidar deployed in the gentler Troutdale, Wasco Airport, and Vansycle Ridge sites. The dominant westerly winds at the site resulted in the location of this lidar being on the downwind edge of an orographic complex, therefore experiencing a strong increase in turbulence production and consequently dissipation. Similarly, the scanning lidar located at Boardman showed higher values of ϵ due to increased turbulence in the wake of a wind farm.

The extensive duration of the WFIP2 field campaign has allowed for the evaluation of the annual cycle of ϵ: the increased convective mixing in summer determines higher values of ϵ compared to the typically more quiescent winter conditions, with an average difference that can reach 1 order of magnitude, both at the microscale and at the mesoscale, in the surface layer and above. We have determined the impact of this seasonal cycle on the average diurnal climatology of ϵ. Overall, ϵ is, on average, up to 3 orders of magnitude higher in summer compared to winter. The diurnal cycle, with higher values of ϵ during daytime convective conditions and lower values at night, is much stronger during the summer when diurnal differences in ϵ values are about 2 orders of magnitude, while the reduced daytime convection during wintertime leads to a more uniform average daily climatology, with less than 1 order of magnitude of difference between daytime and nighttime values of ϵ.

Future work can explore and compare the variability of ϵ from other datasets in different topographic conditions, as well as in the offshore environment (Peña et al., 2009; Canadillas et al., 2010; Türk and Emeis, 2010). Assessing the spatial and temporal variability of ϵ within a typical grid cell of a mesoscale model will provide further insights into the validity of sub-grid-scale ϵ parameterization schemes during various atmospheric stability conditions. As this variability appears to be dependent on several different atmospheric and topographic factors, complex techniques are likely needed to provide accurate spatial representations of ϵ over a mesoscale grid. Sophisticated tools such as physics-driven machine-learning techniques (Sharma et al., 2011; Xingjian et al., 2015; Alemany et al., 2018; Gentine et al., 2018) are paving the path to accurately capture the microscale variability of ϵ in mesoscale models.