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 km2 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.


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 15 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 not valid anymore are currently not well understood, as well as how different atmospheric and topographic conditions can impact the development and decay of turbulent structures. 20 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 socio-economic 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, turbulence dissipation rate has been shown to have a primary role in the formation of frontal structures (Piper and Lundquist,5 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 turbulence dissipation rate. A recent sensitivity study (Yang et al., 2017;Berg et al., 2019) 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 parametrization of turbulence 10 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 15 (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), flown on aircrafts (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 a limited temporal frequency and their inherent volume averaging (Frehlich and Cornman, 2002;Wang et al., 20 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 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 25 a subgrid-scale parametrization scheme (Krishnamurthy et al., 2010).
Here, we retrieve turbulence dissipation rate from eleven instruments in a complex terrain region, thus building one of the widest observational assessment of to date. We explore how topography triggers the variability of at various temporal and spatial scales. We describe the WFIP2 field campaign in Section 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 30 different instruments, and we further refine and extend a novel approach to derive from wind profiling lidars. In Section 3 we present the spatial variability of at both the microscale and mesoscale by comparing the estimates from multiple instruments in different locations, with a 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 turbulence dissipation rate in the boundary layer.

The WFIP2 field campaign
The Second Wind Forecast Improvement Project (WFIP2) , which involved a field campaign (Wilczak 5 et al., 2019) in the U.S. 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 which span an approximately 300km wide area. Two profiling lidars were located at the western and eastern edges of this region, at Troutdale and Vansycle Ridge, respectively, with an additional scanning lidar located 10 in Boardman (Figure 1, panel a). A region with a high-density of instruments (HD Site in Figure 1, panel a), 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 and one scanning lidar ( Figure 1, panel b) and the sonic anemometers on five meteorological towers (Figure 1, panel c).
Multiple sonic anemometers were located on several meteorological towers at the Physics Site , which represented the finest array of instruments at WFIP2, aimed at having multiple measurements in an area similar in size to 15 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 ASL (the elevation of the locations of the meteorological towers used in this study are reported in Table 1). Extensive arrays of wind turbines are located on the northern side of the Columbia River and on the south-western part of the studied region, as well as to the east of the Phyiscis Site. The sonic anemometers used in this project provide 20 Hz measurements of the three components of the wind and virtual 20 temperature at 10m AGL, 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 ±30 • from the orientation of the tower boom (McCaffrey et al., 2017b). Less than 10% of the data were excluded due to tower wake contamination.
Data from the sonic anemometers were used to assess atmospheric stability, calculated in terms of the Obukhov length L: 25 where θ v is the virtual potential temperature (K), calculated from the virtual temperature measured by the sonic anemometers; u * = (u w 2 +v w 2 ) 1/4 is the friction velocity (m s −1 ); k = 0.4 is the von Kármán constant; g = 9.81m s −2 is the acceleration due to gravity; and w θ v is the kinematic sensible heat flux (m K s −1 ). An averaging period of 30 minutes (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 ≤ −500 m and L > 500 m; unstable conditions 30 for −500 m < L ≤ 0 m; and stable conditions for 0 m < L ≤ 500 m (Muñoz-Esparza et al., 2012). Neutral conditions were infrequently recorded (less than 10% of the times).  Table 2. (RHI) and vertical stare scans within 15 minutes. Details of the scan patterns can be found in Choukulkar (2018). For this instrument we retrieve up to 300 m AGL.

10
A WINDCUBE version 2 (v2) was deployed on a low-grass surface on the top of Gordon's Ridge (728 m ASL). A second v2 was deployed at Vansycle Ridge (542 m ASL), in a site with grazed grass (Yang et al., 2013) about 20 km east of the Table 1. Elevation and period of data collection of 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 Figure 1. Finally, a Halo Streamline scanning Doppler lidar was deployed near a regional airport surrounded by farmland at Boardman (112m ASL). The long range fiber-optic based scanning Doppler lidar provides 3D scanning capabilities, and performed a wide range of scans covering the atmospheric boundary layer over a period of 15 minutes (Otarola, 2017). In this analysis, only the 5 • elevation angle scans with a scan rate of 1 • s −1 were used to calculate turbulence dissipation rate up to 120 m AGL. The 10 other scans within the 15-minute time period were not usable for turbulence calculations due to either fast scan rates 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).
2.2 Turbulence dissipation rate from sonic anemometer 15 We estimate 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 average. Within the inertial subrange, Kolmogorov's model (Kolmogorov, 1941) relates the second-order structure function with the turbulence dissipation rate : where a is the Kolmogorov constant, which we set equal to 0.52 (Paquin and Pond, 1971;Sreenivasan, 1995). By invoking 5 Taylor's frozen turbulence hypothesis (Taylor, 1935), the spatial separation r can be written as temporal separation τ , so that can be calculated as: We calculate every 30s by fitting the Kolmogorov's theoretical model to the structure function calculated from the sonic anemometer data using a temporal separation between τ = 0.1s and τ = 2s. From data inspection, measurements in the chosen 10 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.

Turbulence dissipation rate from wind profiling lidar
Measurements from wind Doppler lidars can extend our understanding of the variability of turbulence dissipation rate thanks 15 to the 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 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 where a = 0.52 is the one-dimensional Kolmogorov constant. By integrating (4) over the wavenumber space within the inertial 5 subrange, the following expression can be found: where σ 2 v is the variance of the detrended line-of-sight velocity, and L 1 and L N are the length scales which can be used instead of the wavenumbers by invoking Taylor's frozen turbulence hypothesis (Taylor, 1935). For a single sample, L 1 can be defined 10 as 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 AGL. The second term in (6) can typically be neglected as Doppler lidars generally have θ < 0.1 mrad. For multiple samples, where N is the number of samples used in the calculation. 15 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 20 turbulent quantities, especially in stable conditions (Pichugina et al., 2008).
As shown in Bodini et al. (2018), the appropriate time scales for the lidar retrievals can be determined in different ways.
When co-located 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 25 extension of the inertial subrange from the fit (Tonttila et al., 2015).
In the WFIP2 case, no sonic anemometers co-located 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 contamina- 30 tion (Newman et al., 2016). As a consequence, we further extend this method and we estimate the optimal sample length N to  Figure 2. Example of 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 dashed line shows the theoretical −5/3 slope of the spectrum in the inertial subrange.
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 co-located 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, 5 which provided the best agreement (R 2 > 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.05Hz, which are most affected by instrumental noise (Frehlich, 2001). An example of local regression of an experimental lidar spectrum at WFIP2 is shown in Figure 2.
Finally, the contribution due to instrumental noise needs to be considered. The observed variance σ 2 v in (5) can be thought 10 as a combination of three different contributions, which can be considered as independent of one other (Doviak et al., 1993): where σ 2 w is the contribution from atmospheric turbulence at the scales the lidar can measure (Brugger et al., 2016), σ 2 e is due to the instrumental noise, and σ 2 d 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 −1 ). The contribution of instrumental 15 noise σ 2 e can be written as a function of the signal-to-noise ratio (SNR) Pearson et al. (2009): where ∆ν is the signal spectral width; α 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 α = SN R within a single range gate. Therefore, can be determined as with the accurate choice of the appropriate sample length N , as described.

Turbulence dissipation rate from scanning Doppler lidar
Turbulence dissipation rate from the scanning Doppler lidars is estimated using the azimuth structure function method (Frehlich 5 et al., 2006;Krishnamurthy et al., 2011). The structure function from the radial velocity estimates can be used to estimate 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 10 of 1 • s −1 , and an accumulation time of 1 s, which determine an azimuth spacing ∆Φ = 1 • .
The structure functionD wgt of the mean Doppler lidar velocity perturbations,v , in the azimuth direction is given bŷ 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 15 is produced with different lidar pulses (assuming no multi-scattering effects); therefore, the velocity error variance σ 2 e (R) 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: 20 where r denotes the distance along a fixed laser beam, s = R(φ 1 − φ 2 ) is the transverse coordinate, p = (r 2 + s 2 ) 1/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 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 25 azimuth structure function can be modeled as where s = R(φ 1 − φ 2 ), σ is the standard deviation of the transverse velocity fluctuations and G a (η, µ, ζ) is the derived model based on weighted velocity estimates and the von Kármán model, as provided in Equation (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 functionD wgt (R, kR∆Φ, θ) and the model estimatesD wgt (s, σ, L o ). The dissipation rate can then be estimated by (Hinze, 1959) 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 10 15 minutes.

Results and Discussion
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.   each time, the ratio between from each sonic anemometer and the average (at that time) from all the 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. For all five sonic anemometers, the standard deviation is higher during stable conditions compared to unstable conditions, with mean (across the five anemometers) values 5 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 turbulence dissipation rate, and the high-density array of meteorological towers at the Physics Site represents and 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 10 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. As the wind at the site is almost always slightly south-westerly, 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. The five meteorological towers at the Physics Site can be divided into the sub-group on the western side of the Physics

15
Site (towers P03 and P09) and the remaining ones east of this cluster, which we will refer to as "eastern" (towers P04, P05, and P10). We note that the far eastern side of the Physics Site includes an 80-m tower . 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 ( Figure 5). To confirm this result, the correlation between retrievals from all the possible pairs of meteorological towers at the Physics Site can be studied ( Figure 6, panel a). Stations which 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 turbulence dissipation rate, as reasonable since the spatial separation between the towers increases. However, when looking at the correlation between the retrievals 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 seem to be capable of modifying the spatial evolution of correlation between turbulence dissipation at different locations.
The relationship between correlation coefficient and separation can also provide a confirmation of the larger variability of 5 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 ( Figure 6, panel b), 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 10 meteorological towers, no systematic trend can be found (plot shown in the Supplement).
14 Atmos. Chem. Phys. Discuss., https://doi.org /10.5194/acp-2018-1131 Manuscript under review for journal Atmos. Chem. Phys.  Finally, the temporal variability of turbulence dissipation rate at the microscale can be assessed in terms of the annual cycle of . The increased daytime convection combined with stronger, on average, winds during the summer cause larger turbulent mixing, which in turns leads to higher values of dissipation rates compared to winter months. Figure 7  with median values over an order of magnitude larger in summer than winter, at all the five locations considered within the Physics Site. As a consequence, the inter-quartile range of also reveals an annual cycle (Figure 7, panel b), with a larger range of variability in summer than winter, again with differences of orders of magnitude.

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 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 turbulence dissipation rate found at the Physics Site can also be detected from the retrievals, at higher 5 altitude, from the lidars at the mesoscale. Figure 8 shows the time series of from the different lidars, 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's 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. The time series confirm that turbulence dissipation shows a distinct seasonal   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, which can have important impacts on wind energy production downwind (S. 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 10 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 turbulence dissipation rate at the various locations of the four profiling lidars and the two scanning lidars, respectively. The left column values are about two to three 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 20 stronger turbulent mixing, which in turns causes higher values of turbulence dissipation.

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, where a local equilibrium between production and dissipation ( ) of turbulence is assumed. The error introduced by the parametrization of has been shown to be responsible for up to 50% 25 of the variance of hub-height wind speed predicted by models. The detailed study of observations in the surface layer has a great potential for reducing the uncertainty in our understanding of 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 turbulence dissipation rate in complex terrain at both the 30 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. 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 x 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 the convective mixing during the day. On the other hand, considerable differences (up to one 5 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, in correspondence with the dominant westerly wind pattern, thus suggesting the possible impact of topography in triggering the variability of . The change of correlation between in different locations is not fully determined purely by spatial separation, as topographic features maintain an 10 importance in influencing it. Therefore, the representation of 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 here been further refined and tested to derive without the need of in situ measurements colocated with the lidars. The profiling lidar located at the topographically complex Gordon's Ridge site systematically detected values which, on average, were over one order of magnitude higher than what was measured by the profiling lidar deployed in 5 the gentler Troutdale, Wasco Airport and Vansycle Ridge sites. The dominant westerly winds at the site resulted in the location of this lidar to be 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 the 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 10 convective mixing in summer determines higher values of compared to the typically more quiescent winter conditions, with an average difference which can reach one 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 three 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, where diurnal differences 15 in values are of about two orders of magnitude, while the reduced daytime convection during wintertime leads to a more uniform average daily climatology, with less than one 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). Once this systematic assessment of 20 the variability of turbulence dissipation has been completed from different regions, all the insights on the spatial and temporal variability of should be incorporated into numerical weather prediction model representations of turbulence. As this variability appears to be dependent on several different atmospheric and topographic factors, machine learning techniques, which have already been successfully applied to a broad range of complex atmospheric problems (Sharma et al., 2011;Xingjian et al., 2015;Alemany et al., 2018;Gentine et al., 2018), could provide a reliable representation of .

25
Data availability. The data of the sonic anemometers and wind Doppler lidars at the WFIP2 field campaign are publicly available at https: //a2e.energy.gov/data.
Author contributions. JKL, LKB and MP helped designing and carrying out the field measurements. NB analyzed the data from the sonic anemometers and the profiling lidars and made the figures, in close consultation with JKL. RK analyzed the data from the scanning lidars.
NB wrote the paper, with significant contributions from JKL and RK. All the coauthors contributed to the refining the paper text.