Low-level mixing height detection in coastal locations with a scanning Doppler lidar

. Mixing layer height (MLH) is one of the key parameters in describing lower tropospheric dynamics and capturing its diurnal variability is crucial, especially for interpreting surface observations. In this paper we introduce a method for identifying MLH below the minimum range of a scanning Doppler lidar when operated at vertical. The method we propose is based on velocity variance in low-elevation-angle conical scanning and is applied to measurements in two very different coastal environments: Limassol, Cyprus, during summer and Loviisa, Finland, during winter. At both locations, the new method agrees well with MLH derived from turbulent kinetic energy dissipation rate proﬁles obtained from vertically pointing measurements. The low-level scanning routine frequently indicated non-zero MLH less than 100 m above the surface. Such low MLHs were more common in wintertime Loviisa on the Baltic Sea coast than during summertime in Mediterranean Limassol.


Introduction
Turbulent atmospheric mixing is a key process in the lower troposphere for climate, weather and air quality.Turbulent mixing is regarded as a significant player in aerosol microphysical processes (e.g.Nilsson et al., 2001;Wehner et al., 2010;Hirsikko et al., 2013) and in cloud microphysics (e.g.Pinsky et al., 2008).Representing turbulent mixing in numerical models requires an understanding of the variability in space and time, the length scales involved and the processes that are responsible: friction, surface heating, shear (Baklanov et al., 2011).Stably stratified layers and diurnal cycles require particular attention (Holtslag et al., 2013).From an air quality perspective, the height of the layer that is in constant contact with the surface, i.e. mixing layer height (MLH) is a critical parameter governing the dispersion of air pollutants (e.g.White et al., 2009).The continuous monitoring of the atmospheric mixing profile covering the lowest few kilometres is not a straightforward task and there are only a few long-term data sets (e.g.Harvey et al., 2013;Schween et al., 2014).In many studies the mixed layer is characterised indirectly, often in terms of MLH inferred from aerosol backscatter profiles (e.g.Baars et al., 2008;Korhonen et al., 2014) or temperature profiles (e.g.Beyrich and Leps, 2012).However, indirect methods in MLH estimation may potentially suffer from erroneous interpretation, especially for stably stratified layers and during the initial early morning phase or afternoon collapse of a convective boundary layer (e.g.Pearson et al., 2010;Schween et al., 2014).
In situ measurements of turbulent mixing in the lower troposphere have been conducted successfully on a variety of platforms.For instance, turbulent mixing measurements have been carried out with radiosondes (Harrison and Hogan, 2006), tethered balloons (Siebert et al., 2003), various aircraft (Muschinski and Wode, 1998;Khelif et al., 1999) and recently small unmanned aerial vehicles (Martin et al., 2011).Deployed in-aircraft in situ sensors can yield mixing information with an unsurpassed resolution both spatially and temporally (Muschinski and Wode, 1998;Martin et al., 2014).The drawback with these methods is that they are restricted to short-term campaigns.Only radiosonde measurements are possible routinely, but turbulent sensors are not yet part of the standard operational package.Furthermore, the temporal resolution for routine operational launches is rather coarse; typically a maximum of four per day.Mastbased sonic anemometers provide excellent turbulent measurements with high temporal resolution, but masts taller than 100 m are rare; therefore remote sensing techniques are currently the only viable option for long-term continuous monitoring through the entire lower troposphere.
Several remote sensing techniques, such as Doppler sodar (e.g.Beyrich, 1997;Seibert et al., 2000;Emeis et al., 2008), Doppler lidar (e.g.Harvey et al., 2013;Schween et al., 2014) and radar wind profiler (e.g.Bianco et al., 2008;Emeis et al., 2008), enable continuous measurements of the vertical wind velocity (w) profile with high time resolution.Subsequently, these measurements can be processed to provide vertical profiles of vertical velocity variance, σ 2 w (e.g.Pearson et al., 2010), or turbulent kinetic energy (TKE) dissipation rate (e.g.O'Connor et al., 2010) with resolutions better than a few minutes.However, no single remote sensing instrument has been able to cover the full range of MLHs from instrument level up to the top of the atmospheric boundary layer.Doppler sodars can cover the low range from 10 m up to a few hundred metres, possibly reaching 1 km in good conditions (e.g.Emeis et al., 2008), but, in many environments, the daytime convective MLH exceeds the range of sodars.On the other hand, vertically pointing Doppler lidars and radar wind profilers are usually sensitive enough to reach the top of the atmospheric boundary layer and beyond but cannot see closer than an instrument-specific range, typically 100-200 m (Bianco et al., 2008;Pearson et al., 2009;Srinivasulu et al., 2012), which hampers the detection of low-level only mixing.
Scanning Doppler lidars can partially overcome a minimum height limitation by retrieving radial wind velocity measurements at a low elevation angle (Banta et al., 2006).In this paper we present a method whereby scanning Doppler lidars can identify the presence of turbulent mixing from the instrument level up, making it possible to cover the full range of potential MLHs with an appropriate selection of scan types from one instrument.From the various possible low-elevation angle scanning patterns, we selected vertical azimuth display (VAD) scans as the basis for low-level MLH detection.The main reason for choosing VAD scans is that it simultaneously provides the horizontal wind profile (e.g.Browning and Wexler, 1968) and can be utilised in studying the surface effects on the wind field in the vicinity of the Doppler lidar.The method for identifying turbulent mixing and obtaining MLH is described in Sect. 2. By applying the VAD-based MLH detection to two data sets from very different environments, i.e. summertime Mediterranean coast and wintertime Baltic Sea coast, in Sect.3, we show that this method compares well with MLH inferred from vertically pointing measurements.With this method we frequently identified MLHs below the lowest usable range gate of the Doppler lidar operating in vertical mode at both locations.

Measurements
In this study we utilise Doppler lidar measurements from two locations: Limassol, Cyprus, and Loviisa, Finland (Fig. 1).
Measurements were carried out with a Halo Photonics Streamline scanning Doppler lidar (Pearson et al., 2009).The Halo Photonics Streamline is a 1.5 µm pulsed Doppler lidar with a heterodyne detector that can switch between co-and cross-polar channels (Pearson et al., 2009).Standard operating specifications are given in Table 1, and the minimum range of the instrument is 90 m.The accumulation time per ray varies according to scan type and location.At Limassol, measurements were performed at the Cyprus University of Technology campus (34.6756 • N, 33.0403 • E; 15 m a.s.l.) from 22 August to 15 October 2013.Solar noon at Limassol is 09:40 UTC.The measurement site was on a rooftop 600 m NE from the Mediterranean Sea shoreline (Fig. 1).The Limassol campaign took place in typical Mediterranean summer conditions with surface temperatures ranging from +15 to +35 • C, very low cloud cover and three rain showers.
At Loviisa, the measurement campaign took place on the Fortum Power and Heat Oy nuclear power plant site on Hästholmen island (60.3660 • N, 26.3500 • E; 24 m a.s.l.) from 10 December 2013 to 17 March 2014.Solar noon at Loviisa is 10:20 UTC.Hästholmen island is approximately 2000 m long in the SE-NW direction and 500 m wide in the SW-NE direction (Fig. 1).From the measurement platform at the centre of the island, the distance to shoreline in the SW direction was 300 m and in the NE direction was 200 m.The Loviisa measurements were representative of winter conditions in the Baltic Sea region, characterised by surface temperatures ranging from −25 to +10 • C, frequent bursts of rain and snow and few cloud-free conditions.
At both locations, a fixed scanning routine was operated continuously throughout the campaign.At Limassol the scanning schedule consisted of two VAD scans (e.g.Brown-  ing and Wexler, 1968) at 30 • and 10 • elevation angle every 20 min (Table 2), a three-beam Doppler beam swing every 20 min as an alternative method for retrieving the vertical profile of horizontal wind (e.g.Lane et al., 2013) and a range height indicator scan every 10 min.Besides scanning, 15 out of every 20 min were available for vertically pointing measurements.Vertically pointing data were recorded at 3 s integration time and every second beam was measured with the cross-polar receiver (Pearson et al., 2009).At Loviisa, the terrain allowed VAD scans to be operated at a lower elevation angle than at Limassol and thus the scan schedule at Loviisa consisted of a 24-azimuthal-direction VAD scan at 15 • elevation angle and a 72-azimuthaldirection VAD scan at 4 • elevation angle (Table 2).At Loviisa, vertically pointing measurements were initially operated at 8 s integration time, but on 11 February 2014 the integration time was increased to 16 s.Again, every second vertical beam was measured with the cross-polar receiver as at Limassol.Due to the longer integration time for each beam when scanning, there were 13 min free per 30 min for vertically pointing measurements.At both locations the focus of the Doppler lidar telescope was set to 2000 m.In this study we utilise only the VAD scans and co-polar vertically pointing measurements.

Atmos
Velocity measurement uncertainty is directly related to the instrument sensitivity, so there is a potential trade-off between achieving the high temporal resolution suitable for investigating turbulent conditions at the expense of measurement sensitivity and uncertainty.The integration times for each particular scan at each location were selected to achieve the highest temporal resolution possible while retaining sufficient sensitivity for each individual measurement; at Loviisa, a much lower atmospheric aerosol loading required a longer integration time per individual ray.Additional optimisation was required when implementing the scan schedule so that the relevant scans were acquired while still providing enough vertically pointing coverage.

Estimating turbulent mixing from a VAD
In a smooth homogeneous wind flow field, radial velocities in a VAD follow a sinusoidal curve.A vertical profile of the horizontal wind vector is obtained from the sinusoidal fit (Fig. 2a) to the VAD at each elevation level (Browning and Wexler, 1968).Deviations from this ideal shape (cf.Fig. 2a) can originate from several processes such as wind field divergence or deformation (Browning and Wexler, 1968)  surement in a VAD, the observed radial velocity (V R ) at a single range gate can be expressed analogous to the Reynolds formulation in statistical turbulence theory: where V wind is the radial component of homogeneous horizontal wind and R is the deviation from the homogenous horizontal wind.The horizontal wind component V wind is then derived by fitting a sinusoidal curve to V R as a function of the azimuthal angle (Fig. 2a) for each elevation level in a VAD, with R being the residuals of these fits.Contrary to the Reynolds formulation, however, the deviation term R in a VAD contains more components than the purely turbulent fluctuations.These additional components arise from fact that, over the relatively large volume of atmosphere covered by the VAD scan, there are often nonturbulent changes in the horizontal wind speed and direction on scales of hundreds of metres to kilometres that are also contained in R in Eq. ( 1) when V wind is estimated from a fit over the full VAD.In quiescent conditions, the non-turbulent contributions may even form most of the residual.This is evident when correlating the residuals of two consecutive range gates as in Fig. 2b; a high correlation indicates that the residuals are dominated by flow patterns with length scales that are large compared to the 30 m radial resolution of the instrument.
The contribution of some non-turbulent processes (such as translation, divergence, relative vorticity, stretching and shearing deformation) to R can be estimated through a Taylor expansion (Browning and Wexler, 1968).However, within the first 100 m above the surface, the focus of this study, R can also contain significant contributions from local distortions.Typical surface-induced changes would be wind flows around a building or wind channelling around an island.In most environments quantifying this kind of local effects with reasonable accuracy would require such detailed knowledge of the surface at the measurement location that determining turbulent contribution directly from Eq. ( 1) is impractical.Therefore, we consider the change in V R from one range gate at distance r to the next range gate at (r +30 m) for one radial measurement: (2) In this way, the problem of determining the non-turbulent changes in the wind field within the VAD volume has been reduced to the question of whether these changes are significant on the scale of the range resolution of the Doppler lidar.
The changes in the wind field that can be represented with the Taylor expansion occur over length scales that are significantly larger than the Doppler lidar range resolution of 30 m (e.g.Browning and Wexler, 1968) and thus their contribution to R can be considered negligible.Similarly, if the localised effects in the wind field can be considered smooth at the Doppler lidar range resolution scale (here 30 m), i.e. the surface is reasonably homogeneous or the measurement is not close to the surface, then we can assume also their contribution to R to be very small.
With these assumptions we can now estimate where V turb is the contribution from turbulent mixing that we are interested in and i.e. we have assumed that the instrumental measurement uncertainty for the two range gates is uncorrelated.Thus the proxy variable for identifying turbulent mixing is the variance of the difference of residuals from two consecutive elevation levels in a VAD subtracted by the variance of the corresponding measurement uncertainty (cf.Fig. 2c):

Atmos
where the subscript i refers to individual radials and n is the number of radials within a VAD.The measurement uncertainty variance σ 2 υ is estimated from the δ 2 i calculated for every ith radial: where median is taken over all radials at each range gate.The instrument uncertainty in velocity δ is primarily a function of the signal-to-noise ratio (SNR) (Pearson et al., 2009) and is calculated for every individual radial velocity measurement in a VAD.Based on pointing accuracy tests performed during both campaigns towards hard targets at known direction and distance, we consider the pointing accuracy error negligible especially when compared to the uncertainty arising from the measurement itself.The choice of median in Eq. ( 6) is to avoid outliers skewing the distribution.In general, most of the points in a particular VAD display similar sensitivity, but the presence of cloud or precipitation can increase SNR, whereas certain radials may be completely or partially obscured by buildings or trees.Measurements in the presence of precipitation have to be carefully evaluated as the signal can be dominated by the terminal fall velocity of the drop, and large drops do not necessarily track the turbulent motion of the air.Here, we only calculate σ 2 VAD when SNR > 0.0025 for at least 20 out of 24 points at any given elevation level.In terms of σ 2 υ the SNR limit of 0.0025 is equivalent to a threshold of 1.58 m 2 s −2 using the instrument specifications given in Table 1.This is derived by calculating δ from SNR according to O'Connor et al. ( 2010) and then applying Eq. ( 6).Note that this assumes that the turbulence is isotropic in nature.
Abrupt changes in surface roughness also have an impact on the turbulent properties of the wind field (Garrat, 1990).For instance, the 4 • elevation angle VAD at Loviisa is so close to the surface that the radial wind field is clearly affected by the change of surface roughness moving from sea to land.To minimise the effect of surface roughness changes at Loviisa we considered only a 55 • wide sector (i.e. 12 azimuthal directions) upwind of the island to derive MLH from the 4 • elevation angle VAD.For 10 • and higher elevation angle VADs we used the full 360 • to determine σ 2 VAD and MLH.

Turbulent kinetic energy dissipation rate
The dissipation rate of TKE was determined from the vertical velocity measurements according to the method described by O' Connor et al. (2010).This method utilises the velocity variance over a specific number of samples, from which the dissipation rate is derived using appropriate advective length scales obtained from the vertical profiles of horizontal wind.The implicit constraint in this method is that the advective length scales for calculating the variance should remain within the inertial subrange.Typically, this means that the total time available for collecting samples for one variance profile should not exceed about 3 min and that the number of samples available per dissipation rate profile therefore depends on the integration time for an individual ray.
For vertical profiles with 3 min resolution at Limassol, optimal operating conditions provided 60 vertical velocity samples per dissipation rate profile.At Loviisa, a much lower atmospheric aerosol loading required a longer integration time per individual ray to obtain sufficient sensitivity; therefore only 11 vertical velocity samples per 3 min resolution dissipation rate profile were available.Vertical resolution was 30 m, and dissipation rate values were determined only when the relative uncertainty in the variance was less than 1, i.e. observed variance at least twice the theoretical contribution from noise.

Comparison to vertical wind speed variance
We first investigate how σ 2 VAD correlates with the vertical velocity variance σ 2 w calculated directly from the vertically pointing time series.Note that we do not expect σ 2 VAD to be equivalent to σ 2 w since the effective measurement volumes can encompass very different length scales.Figure 3 shows that, at Limassol, σ 2 VAD from both 30 and 10 • elevation VADs correlates reasonably well with σ 2 w .At Loviisa the correlation is not quite as good for the 15 • VAD (Fig. 3), which we partly attribute to the narrow range of observed values available; note that for high-latitude Finland in winter there is minimal diurnal influence on the turbulent mixing.In addition, the reduction in instrument sensitivity due to the low aerosol loading results in fewer data points at Loviisa and very few low σ 2 w values during the campaign.At Loviisa, the correlation between the 15 and 4 • elevation VADs is rather poor (Fig. 3d), but the data are close to the 1 : 1 line and show a similar scatter to data from Limassol.Altogether, the relationship between σ 2 VAD and σ 2 w is reasonably linear especially for σ 2 w < 0.1 m 2 s −2 (Fig. 3), which is the most important range for determining the MLH.The reasonably close agreement with σ 2 VAD values at different elevation angles (Fig. 3b  and d) suggests that the relationship between σ 2 VAD and σ 2 w is independent of the VAD elevation angle (Fig. 3) and thus the VAD elevation angle is not a critical parameter for considering σ 2 VAD as a proxy for turbulent mixing.

Comparing MLH from a VAD and vertically pointing measurements
The diurnal variation of TKE dissipation rate calculated from vertically pointing measurements is plotted together with σ 2 VAD for Limassol in Figs. 4 and 5 and for Loviisa in Figs. 6  and 7.In a qualitative sense the diurnal dissipation rate and σ 2 VAD profiles agree.The greatest difference is the instrument sensitivity for VADs below 15 • in elevation: the signal clearly peters out at a lower altitude than for the vertically pointing data.This is due to two reasons: firstly, at low elevation angles the range, or path length, to a given altitude is obviously much further than when vertically pointing and the relative rate of attenuation in the vertical plane is correspondingly higher; secondly, the integration time for vertically pointing rays can be increased without a major impact on the scan schedule.
We used the simplest possible MLH detection scheme -a constant threshold value -to assess the usefulness of σ 2 VAD as a measure of the MLH.For σ 2 VAD , the MLH was diagnosed    the changes in the VAD-based MLH due to ±0.01 m 2 s −2 changes in the σ 2 VAD threshold are small.For applications where accurate MLH detection is of critical importance, a more sophisticated MLH detection scheme than the flat threshold used here would be appropriate.The flat threshold is clearly a reasonable and robust initial estimate, and the selected threshold of 0.05 m 2 s −2 gives the best agreement with the MLH inferred from the TKE profile with a threshold of 10 −4 m 2 s −3 .

Frequency of low mixing level heights at Limassol and Loviisa
The vertically pointing TKE dissipation rate data close to surface indicate no significant mixing, implying MLH below the lowest vertical range gate at Limassol 42 % of the time and at Loviisa 62 % of the time.The VAD-based MLH estimates show that, at Limassol, in 58 % of the cases when MLH must be below the vertically pointing altitude limit, there is a shallow mixed layer at the surface (Fig. 9).At Loviisa the VADs indicate a shallow mixed layer at the surface on 87 % of the cases when the MLH must be below the vertically pointing altitude limit (Fig. 9).At Limassol the MLH exhibits a clear diurnal cycle with low-altitude mixing levels occurring almost exclusively during night time (Fig. 9b).This agrees with radiosonde observations of mixing level heights at coastal Mediterranean locations during summer (e.g.Seidel et al., 2012).At Loviisa, however, very low MLHs are also common during daytime (Fig. 9d), typical for cold conditions with a stably stratified atmosphere and minimal surface heating (e.g.Liu and Liang, 2010).

Conclusions
We have shown for two very different environments that a low-elevation-angle Doppler lidar VAD scan can be used to identify the presence of turbulent mixing in the atmosphere.Although the method was developed with measurements at coastal locations, the method should be applicable in any environment where VADs can be performed.Furthermore, the VAD-based proxy for turbulence can be used to identify the MLH.If scanning at a very low elevation angle is feasible at the measurement location, the VAD-based MLH can detect the presence or absence of mixing from the instrument level up.However, at elevation angles lower than 10 • , the impact of surface roughness changes across the VAD volume must be taken into consideration.
Comparison of MLHs from vertically pointing data and VADs shows reasonably good agreement, especially considering the simplicity of the MLH detection scheme used here.However, the rapid increase in radius for VADs at low elevation angles limits the altitude range of MLH retrievals from VADs.Therefore, to cover the full range of MLHs from ground level up, a combination of vertically pointing and VAD measurements are most suitable.In this manner, turbulent mixing can be identified from the surface up to heights of 1 km or more continuously with good time resolution.At the same time, VADs can be used to retrieve the wind profile, which in turn can be used to e.g.identify wind shear generated mixing.
Finally we have demonstrated that very shallow MLHs can be present during the majority of the time when vertically pointing measurements indicate no mixing; i.e.MLH is below the lowest usable range gate in vertically pointing measurement mode.At Limassol, representing Mediterranean summer time conditions, such low MLHs occurred only during the night; at Loviisa, in Baltic Sea wintertime conditions, very low MLHs were also common during the day.

Figure 1 .
Figure 1.Location of Loviisa and Limassol measurement sites.In the 20 km × 20 km topographic map inserts for Loviisa (National Land Survey of Finland, 2014) and Limassol (United States Geological Survey, 2014) the location of the lidar is indicated by a red dot.
Figure 2. (a) Radial wind speed and sinusoidal fit on 24 August 2013 at 03:57 (UTC) for two consecutive altitudes from a 10 • elevation angle VAD scan.(b) Scatterplot of the residuals of the fit presented in panel (a).(c) Histogram of the difference of residuals at 148 and 153 m a.s.l.altitudes from a 10 • elevation angle VAD scan for a calm period (03:57 UTC) and a turbulent period (11:59 UTC) with corresponding variances indicated in the legend.The σ 2 υ (Eq.6) is 0.03 at 03:57 UTC and 0.04 at 11:59 UTC.

Table 2 .
Scan settings at Limassol and Loviisa.At Loviisa the 15 • elevation angle VAD integration time was decreased from 10 to 7 s on 11 February 2014.

Table 3 .
Comparison of VAD and vertically pointing MLH estimate.For Limassol data the 25th, 50th and 75th percentile of the VAD-based MLH are presented for the cases when the nearest-neighbour TKE-based MLH is 120 m a.s.l.For Loviisa data, the 25th, 50th and 75th percentile of the VAD-based MLH from the 15 • elevation angle VAD are presented for the cases when the nearest-neighbour TKE-based MLH is 159 m a.s.l.The 25th, 50th and 75th percentile of the VAD-based MLH from the 4 • elevation angle VAD at Loviisa are presented for the cases when the nearest-neighbour 15 • elevation angle VAD-based MLH is 55 m a.s.l.