A Halo Photonics Stream Line XR Doppler lidar has been deployed for the Indianapolis Flux Experiment (INFLUX) to measure profiles of the mean horizontal wind and the mixing layer height for quantification of greenhouse gas emissions from the urban area. To measure the mixing layer height continuously and autonomously, a novel composite fuzzy logic approach has been developed that combines information from various scan types, including conical and vertical-slice scans and zenith stares, to determine a unified measurement of the mixing height and its uncertainty. The composite approach uses the strengths of each measurement strategy to overcome the limitations of others so that a complete representation of turbulent mixing is made in the lowest km, depending on clouds and aerosol distribution. Additionally, submeso nonturbulent motions are identified from zenith stares and removed from the analysis, as these motions can lead to an overestimate of the mixing height. The mixing height is compared with in situ profile measurements from a research aircraft for validation. To demonstrate the utility of the measurements, statistics of the mixing height and its diurnal and annual variability for 2016 are also presented. The annual cycle is clearly captured, with the largest and smallest afternoon mixing heights observed at the summer and winter solstices, respectively. The diurnal cycle of the mixing layer is affected by the mean wind, growing slower in the morning and decaying more rapidly in the evening with lighter winds.
The mixing height (MH) is a crucial value for numerical weather prediction (NWP) modeling and air quality monitoring. Interpretation of in situ measurements of trace gases requires knowledge of the depth through which the trace gases are well mixed (Lin et al. 2003; Locatelli et al. 2015). Quantification of emissions using mesoscale inversions is particularly sensitive to accurate representation of the MH (Lauvaux and Davis 2014). The MH also affects the profiles of temperature, moisture, and momentum in the lower atmosphere, which are imperative for accurate forecasts from NWP output. The parameterization of planetary boundary layer (PBL) processes and the vertical extent of mixing in NWP models affects predictions of severe storms (Cohen et al. 2015); cloud cover, which is important for solar power production (Cintineo et al. 2014); and the mean wind through entrainment at the top of the mixed layer (Tennekes and Driedonks 1981). Additionally, a small daytime MH may lead to strong heating over a shallow depth and have been observed during heat waves (Kunkel et al. 1996). Thus, accurate measurements of the MH are necessary for air quality studies and validation of NWP output.
However, measuring the MH both continuously and accurately is difficult. Radiosondes provide an instantaneous snapshot of the temperature, moisture, and wind profile, which can be used to estimate the MH using various parcel methods (Seibert et al. 2000); however, the vertical profile may not be representative of the mean MH in general, especially if the radiosonde ascends through a pronounced updraft or downdraft. Radar wind profilers have been shown to accurately measure the MH during convective conditions (Angevine et al. 1994), but they are generally unable to measure the MH at night. Similarly, ceilometers and other backscatter lidars have shown promise in continuously measuring the MH during the day (Münkel et al. 2007; Haeffelin et al. 2012), but they have difficulty distinguishing the nocturnal mixing layer (ML) from the residual layer at night and during morning and evening transition periods (Schween et al. 2014).
Doppler lidar observations have been used to estimate the MH, most often using either backscatter or turbulence information from vertical stares (e.g., Hogan et al. 2009; Barlow et al. 2011; Huang et al. 2017). Tucker et al. (2009) evaluates the accuracy of various techniques and finds that vertical velocity variance generally provides the most accurate measure of the MH. However, because of minimum range issues and the presence of nonturbulent motions (e.g., gravity waves), vertical velocity variance profiles alone are insufficient to identify the MH. Vakkari et al. (2015) shows that using the variability along a low-angle conical scan reduces the minimum range, so the MH may be identified when it is below the minimum range of a vertical stare; these MH estimates from conical scans generally agreed with MH estimates from vertical measurements when there was overlapping coverage. Pichugina et al. (2008) uses vertical-slice scans to measure the horizontal velocity variance, which may be used to quantify mixing all the way down to the surface to detect a shallow MH (Pichugina and Banta 2010). However, the shallow conical and vertical-slice scans alone may not be able to measure turbulence through a deep enough layer to penetrate the top of the mixed layer. To overcome the challenges and limitations of each individual scan strategy, a composite approach leveraging the advantages of each scanning strategy is necessary to provide continuous measurements of the MH.
To take advantage of all measured quantities from different scan types, a fuzzy logic approach is introduced here. This technique blends all the data together from multiple scans to determine a unified measurement of the MH and the uncertainty of the estimate. Bianco and Wilczak (2002) used a similar approach to combine several variables measured by wind profiling radar to find the top of the convective PBL. The algorithm presented here has been initially developed to determine the MH for the Indianapolis Flux Experiment (INFLUX), which is a multi-institution collaborative project to measure the greenhouse gas emissions from the metropolitan area of Indianapolis, Indiana (Davis et al. 2017). For INFLUX, a Halo Photonics Stream Line Doppler lidar has been deployed to suburban Indianapolis to measure the mean horizontal winds, turbulence, and the MH.
The paper is organized in the following manner. The distinction between the PBL and MH, as well as measurement of each, is discussed in section 2. The experiment, instrument, scanning strategies, and measurements are described in section 3. The operation of the algorithm itself is thoroughly described in section 4. An intercomparison between the Doppler lidar MH and in situ measurements is shown in section 5. A brief climatology of the MH over Indianapolis in 2016, demonstrating the utility of these measurements, is provided in section 6. The strengths and limitations of this algorithm, along with areas for further research and improvement, are discussed in section 7, and a summary of the algorithm and results is proved in section 8.
2. Definition and measurement of the mixing height
While the PBL and the mixing layer are closely related, the two are not always identical or interchangeable with each other. According to the American Meteorological Society (2017), the PBL is defined as “the bottom layer of the troposphere that is in contact with the surface of the earth.” Thus, the PBL is the layer of the air directly influenced and responsive to surface forcings (Stull 1988). The mixing layer is the depth of air near the ground where pollutants or other passive tracers are vertically dispersed by convection or mechanical turbulence within about an hour (Seibert et al. 2000), where the MH is the top of this layer. During vigorously convective conditions in the unstable PBL, the mixing layer can be characterized as “well mixed,” in which passive tracers are quickly mixed and relatively constant with height throughout the entire mixing layer. During well-mixed time periods, the mixing layer and PBL are often considered identical. Conversely, the PBL height is often ill-defined during the morning or evening transition periods and stable conditions, as the depth to which surface forcings influence the atmosphere is ambiguous and has been identified in multiple ways that are not always in agreement (see Vickers and Mahrt 2004).
Various definitions for the PBL height are largely the result of different ways in which the PBL is measured depending on the type of instrument used. The PBL height can be determined in different ways from radiosonde vertical profiles of temperature, humidity, and wind (Seidel et al. 2010), such as the use of a parcel method during strong convection (Holzworth 1964) or analysis of surface inversions and the wind profile during stable conditions (Joffre et al. 2001). Sodar measurements of the PBL are largely based on the backscatter, which is proportional to the temperature structure function parameter ; thus, sodar PBL heights are sensitive to local heat and momentum fluxes and gradients in temperature and wind speed (Beyrich 1997). Similarly, PBL heights from radar wind profilers are largely based on peaks in the refractive index structure function parameter , which is typically large at the top of the convective PBL, where a humidity gradient is often present (Cohn and Angevine 2000). Lidar backscatter has been used extensively to measure PBL heights (e.g., Menut et al. 1999; Davis et al. 2000; Banks et al. 2015; Huang et al. 2017). A more thorough examination of the methods for measuring the PBL height and MH specifically, including a discussion of the advantages and limitations of different instrumentations and techniques, is provided by Seibert et al. (2000).
Since the PBL is often poorly defined, especially under stable conditions, it is a difficult quantity to measure continuously. However, the MH is a clearly defined quantity that can be quantified continuously with the appropriate measurements. Herein, the presented algorithm is designed to measure the MH as defined by Seibert et al. (2000, p. 1002): “The mixing height is the height of the layer adjacent to the ground over which pollutants or any constituents emitted within this layer or entrained into it become vertically dispersed by convection or mechanical turbulence within a time scale of about an hour.” Since these MH measurements are primarily used for air quality applications, this definition is refined during well-mixed convective time periods to be the mean height that pollutants are dispersed. Using this definition, Doppler lidar measurements of turbulent quantities and other contextual information can be used to directly determine the MH.
3. Doppler lidar deployment and measurements
The INFLUX began in 2010 with the primary goal of measuring citywide greenhouse gas emissions at a high spatial (1 km) and temporal (weekly or finer) resolution. To make these measurements, a suite of instrumentation has been installed in and around the Indianapolis metropolitan area to measure trace gas quantities, sensible and latent heat fluxes, and meteorological quantities, in addition to other episodic measurements, such as research aircraft observations. To complement these measurements, a Halo Photonics Stream Line Doppler lidar was deployed in Indianapolis in 2013 primarily to provide the mean wind profile in the lower part of the atmosphere and turbulence information, including the MH. A more thorough overview of the INFLUX project in general is provided by Davis et al. (2017).
For INFLUX, the lidar has been installed on the roof of a building four stories ( m) above ground level in Lawrence, Indiana, 17 km to the northeast of downtown Indianapolis (see map in Fig. 1). With a southwesterly mean wind direction, the lidar is situated climatologically downwind of the urban core. The site is near the edge of the metropolitan area, with a state park to the north and west. The lidar has been running nearly continuously since August 2013 with the exception of a 6-month gap at the end of 2015, when the system was being upgraded to a Stream Line XR. Since the system upgrade enhanced its sensitivity, the scanning strategy was modified to take advantage of this improvement. Data presented here were taken after the upgrade; thus, only a description of the operation of the Stream Line XR is provided. Details about the lidar design are given in Pearson et al. (2009), and principal system parameters are provided in Table 1. Since 2016, the Doppler lidar has operated with 48-m range gates to a maximum range of 10 km. The system operates at 10 kHz, and 5000 pulses are accumulated so that range-resolved line-of-sight and backscatter intensity estimates are recorded at 2 Hz.
The lidar has been operating in a scan sequence that repeats every 20 min continuously. The scan pattern consists of conical [plan position indicator (PPI)] scans at elevation angles ϕ (above the horizon) of 3°, 10°, 35.3°, and 60°; vertical-slice [range–height indicator (RHI)] scans to the south and east; a zenith stare (lasting 10 min during the day, 4 min at night); and quasi-horizontal stares at to the south and east at night. The scans, their key measured variables, and their height coverage are summarized in Table 2. Based on the scan geometry, the vertical resolution of the measurements is generally dense near the ground ( m) and reduces with increasing altitude. With sufficient aerosol loading, the combined scans provide coverage from a few meters above the ground to the top of and sometimes above the ML. Thus, minimum range issues typically associated with vertically pointed lidars are not pertinent.
The measurements taken by the various scans are used to produce vertical profiles of the mean wind, turbulence, and signal intensity. Profiles of the mean zonal u, meridional υ, and vertical w winds are produced by applying the velocity–azimuthal display (VAD) technique (Browning and Wexler 1968) using radial velocity measurements at different azimuthal angles θ. The measured wide-band signal-to-noise (SNR) can be used to estimate profiles of the relative aerosol backscatter. This is done by normalizing SNR profiles by a range function, which is how the SNR varies in range from the lidar in constant aerosol loading. The range function is calculated as the mean SNR for each separate range gate, using measurements from low ϕ (<3°) angles during convective periods when aerosols can be assumed to be well mixed. While the range function itself varies based on atmospheric conditions, including refractive turbulence (Frehlich and Kavaya 1991; Belmonte and Rye 2000) and depolarization of the backscatter (Sakai et al. 2000), the normalization of SNR produces range-corrected intensity (RCI) profiles that are adequate to identify gradients at the MH. While profiles of RCI may be made from any scan type (RHIs, PPIs, etc.), only RCI profiles from vertical and shallow stares are used here. Additionally, the variance of the SNR () can be useful to identify the MH, since variance is often large at the MH as clean free troposphere air is entrained into the convective ML (Hooper and Eloranta 1986). Thus, profiles of have also been calculated from the stares.
Each of the scans used in INFLUX provides some information on the amount of turbulent mixing in the atmosphere. From the PPI scans, the variance of the residuals of the VAD fit () at each range gate is directly proportional to the variances of u, υ, and w (, , and , respectively) (Eberhard et al. 1989). Values of are more sensitive to for high ϕ scans and more affected by and at low ϕ. At , the is the turbulence kinetic energy (TKE; Eberhard et al. 1989). Profiles of and can be made at user-defined bins (5 m here) from RHI scans to the east and south, respectively (Banta et al. 2006; Pichugina et al. 2008). The vertical stares provide continuous measurements of w at 2 Hz, allowing for a straightforward time series analysis to calculate . Similarly, the shallow quasi-horizontal stares allow for calculation of and for the two look angles by assuming that the contributions of is negligible at . An SNR threshold of −23 dB was used to remove highly noisy radial velocity data, after which uncorrelated noise in the calculated velocity variances was removed for the PPI scans and stares. This is done by applying a structure function fit to the autocovariance of w in the time series and residuals of the VAD fit along a range bin around the PPI, similarly to the technique described by Lenschow et al. (2000). In addition to removing instrument noise, this structure function fit corrects for averaging effects over the pulse volume and accumulation time (Bonin et al. 2016). With this noise removal technique, it is assumed that the lower frequency portion of the inertial subrange is resolved and that turbulence is frozen for the time and space for which the structure function fitting is applied (typically <10 s, corresponding to the 30° sector of PPI). When a portion of the inertial subrange is not captured, such as during stable conditions or far ranges in PPI scans where the distance between adjacent beams becomes large, the value of may be overestimated and the atmospheric velocity variance undermeasured. Further details on these turbulence measurement techniques and their operation, noise removal, and relative accuracies can be found in Bonin et al. (2017).
4. Algorithm description
a. Fuzzy logic overview
To take advantage of all the measured quantities from the different scanning strategies, all of the observations need to be combined together in some manner for a unified determination of the MH. Fuzzy logic is naturally suitable for combining different measured variables together for determining the physical properties of a sample. Additionally, fuzzy logic avoids depending on a single threshold that is valid for all conditions, which results in large sensitivity to the MH estimate (Schween et al. 2014). Within the field of atmospheric science and remote sensing, fuzzy logic has been used extensively, especially within the radar community. For polarimetric weather radar, fuzzy logic has been used for hydrometeor classification (Vivekanandan et al. 1999; Liu and Chandrasekar 2000), identification of nonprecipitation echoes and radar artifacts (Gourley et al. 2007; Mahale et al. 2014), and discrimination between stratiform and convective precipitation (Yang et al. 2013), to name a few uses. Additionally, fuzzy logic has been used to measure the convective PBL depth using quantities measured by a radar wind profiler (Bianco and Wilczak 2002). While fuzzy logic has been used extensively in the atmospheric sciences with radar measurements, this is the first time to the authors’ knowledge that fuzzy logic has been applied to Doppler lidar.
Fuzzy logic is essentially the mapping of multiple input variables to determine the quality or characteristic of a measurement (Mendel 1995). The input variables are related to an output characteristic through membership functions, which vary from zero to one. A membership value of one indicates that a measurement is a member of a certain classification. Membership values from different inputs are aggregated together through a weighted mean, after which the aggregate is defuzzified to infer a characteristic of the measurement. While fuzzy logic can be used to determine many different possible classifications (e.g., as in hydrometeor identification; Vivekanandan et al. 1999), herein we use it to determine whether turbulent mixing is present in a measurement. The vertical extent of turbulent mixing is used to identify the MH. Since the scan pattern described in section 3 takes 20 min to complete and all of the data from these scans are used together, the MH measurement is intrinsically a 20-min average.
b. Nonturbulent motion detection
Before the MH can be determined, it is imperative to identify and separate nonturbulent fluctuations in the mean wind (, , and ) from those that are turbulent. Mahrt (2014) provides a thorough discussion of nonturbulent submeso motions, some of which are resolved by Doppler lidar, such as drainage flows (e.g., Post and Neff 1986) and internal wave motions (e.g., Banakh and Smalikho 2016). These submeso motions and other flow features can contaminate turbulence measurements (Bonin et al. 2017). Drainage and heterogeneous flows in complex terrain particularly affect turbulence measurements taken with PPIs and RHIs, as the spatial variability is used to calculate turbulence quantities. Since Indianapolis is practically flat (terrain varies by m over a 6-km-radius circle from the lidar), local flows are not apparent in the PPI and RHI scans used here. However, waves with periodic alternating rising and sinking motions have been frequently observed during INFLUX, particularly within the stable PBL and at the top of a developing convective ML as is shown in Fig. 2. These wavelike motions need to be identified so that they are not misconstrued as turbulent mixing, which would result in the MH being overestimated.
For purposes here it is not necessary to parse out the exact velocity variance from atmospheric turbulence from wave motions, which would require wavelet analysis or multispectral resolution techniques (Vickers and Mahrt 2003; Viana et al. 2010). Instead, it is only necessary to determine whether turbulent mixing is occurring within a layer. Within the PBL turbulent energy is apparent across a wide range of frequencies as it cascades from larger scales to smaller scales, as is reflected in the velocity spectrum (Kaimal et al. 1976). Energy from pure wave motions is confined to the frequencies of the waves in the packet. Since observed lower-atmospheric waves have periods on the order of minutes to tens of minutes (e.g., Finnigan et al. 1984; Viana et al. 2009; Toms et al. 2017), simply evaluating the w variance at high frequencies or applying a high-pass filter removes effects from waves alone, while some variance from turbulent motions existing within the inertial subrange remains.
For each range gate in the zenith stare, values of in total (across all frequencies) and at high frequencies (denoted by the subscript HF), defined here as the frequency Hz (period less than 1 min), are calculated from the w spectrum above the noise floor. Since the noise can be assumed to be white (Frehlich and Yadlowsky 1994), the noise floor is , as determined using the structure function fit from the autocovariance of w, is equally distributed across all frequencies in the power spectrum. Example profiles of and total during a period when waves were observed above the MH are shown in Fig. 2 to demonstrate how profiles of can be used to differentiate between turbulent and nonturbulent motions. The power spectral density (PSD), used to calculate , within the turbulent and wave layers is also provided in Fig. 2 to demonstrate the differing spectral characteristics. Since the distributions of the velocity power spectrum are a function of the wind speed U (Kaimal et al. 1976), values of are scaled by U in a given layer. A layer of air is flagged as a “wave” if two criteria are met: 1) is small (<0.0025 m s−1) and 2) m2 s−2 ( for that variable, as defined in section 4c). These criteria are chosen to ensure 1) any turbulent mixing is small and 2) the submeso motions are large enough to affect the retrieval of the MH. The scaling by U also ensures that large convective eddies when winds are weak are not misconstrued as wave motions. For layers that are flagged as waves, velocity variances from other scans (RHIs, PPIs) are disregarded in all further analysis.
c. Turbulent layer identification
After submeso motions have been detected using vertical stares, data from all the scans can be combined. To make all of the turbulence quantities comparable, the measurements are fuzzified, or transformed into values that vary between zero and one, according to their membership function. A membership value of one indicates that the measurement is part of a turbulently mixed classification, while zero indicates that a measurement is not. An example of the half-trapezoidal-shaped membership functions that are used here is shown in Fig. 3. The membership functions have two parameters— and —that determine for what values the function increases above zero and reaches a maximum of one. The values of and , which are given in Table 3, are centered and based on threshold values used for identifying the MH in previous studies (i.e., Tucker et al. 2009; Barlow et al. 2011; Vakkari et al. 2015). These values vary depending on the input quantity and what type of scan was used for the measurement. For example, even though both RHIs and shallow stares measure the same quantities ( and ), contamination from can be easily quantified and removed using established techniques from the stare but not from the RHIs. Thus, the range between and is larger for the RHIs. Similarly, and for from the VAD residuals is dependent on ϕ, since higher scans are less sensitive to and but more dependent on , which is typically smaller than the horizontal variances.
Since measurement heights from different scan types vary depending on the geometry, the fuzzified values are linearly interpolated to a common 5-m-height grid. These values are aggregated by taking the mean of the fuzzified values. While some versions of fuzzy logic use weighted means, all inputs are weighted equally here. This aggregate is then used to identify a first guess for the MH . This value for is the lowest height where the aggregate is . The threshold of 0.5 is chosen because it naturally indicates that half of the inputs indicate that turbulent mixing is occurring. Likewise, by using the lowest height where the threshold is met ensures that the turbulent layer is connected to the surface.
Examples of the different measures of turbulence and the aggregate produced by combining these measurements are shown in Fig. 4 for data collected on 17 October 2016. This day was chosen to demonstrate the utility of combining and using measurements from the different scans for a complete representation of the entire PBL, with zenith stares measuring to the top of the PBL and low elevation scans filling in the measurement gap below the minimum range of the zenith stares. Additionally, nonturbulent motions above the MH, such as shown in Fig. 2, are apparent on this day and flagged as such so that they are not misconstrued as turbulent motions when incorporated into the fuzzy aggregate. These waves are apparent in Fig. 4a from 1400 to 2300 UTC between 1500 and 2000 m where is >0.1 m2 s−2, but the is simultaneously small, indicating no turbulence in the layer. The aggregated mean membership function from all the measures of turbulence is shown in Fig. 4f, from which is determined. Generally, the aggregated mean is smooth in height and time at the lower heights (<500 m), since it is well covered by all the different scan types; thus, more inputs are used in the average. At higher heights the aggregate can be discontinuous in time, as often only zenith stare measurements are used as inputs. These observations during the day are subject to sampling limitations, as the vertical extent of the measured is closely related to the maximum height of an updraft, causing fluctuations in w, during the 10-min stare.
d. Use of other indicators for MH
During well-mixed periods, vertical profiles of other measured quantities that are well mixed themselves (i.e., roughly constant with height above the surface layer and below the MH) can be used to improve the MH measurement. Specifically herein, vertical profiles of the mean horizontal wind and RCI are used to refine the measurement of the MH. Large gradients of these quantities are often apparent at the interface of the well-mixed layer and the free troposphere. Additionally, cleaner and drier air being entrained into the mixing layer frequently results in a large at the MH (Menut et al. 1999; Hennemuth and Lammert 2006). Herein, these aforementioned profiles of mean wind, RCI, and are referred to as “indicators” of turbulent mixing; while these quantities are not a direct measure of turbulence itself, their profiles can show where the air is well mixed. These indicators are often a more accurate measure of the mean MH (defined as the height where quantities are well mixed) during strongly convective conditions than profiles of turbulence quantities. To demonstrate this, profiles of , , and calculated from 10-min time series from zenith stares are shown in Fig. 5. Perturbation in the instantaneous MH caused by updrafts and downdrafts are visualized in Figs. 5a,b. Strong updrafts extend above the mean MH, resulting in instantaneous enhancements in RCI and fluctuations in w. Statistics calculated from this period, provided in Fig. 5c, show that can be large () above the mean MH, as a result of these updrafts that extend above the mean height where aerosols are well mixed.
A Haar wavelet (Haar 1910) covariance transform on a vertical profile can be used to detect sharp changes in a given quantity. Haar wavelet covariance transforms have been extensively applied to backscatter profiles from lidars to determine the MH during convective conditions (e.g., Cohn and Angevine 2000; Davis et al. 2000; Brooks 2003). Similarly, a Haar wavelet transform is applied to the profiles of RCI here to detect gradients in the aerosol content at the MH. The dilation of the Haar wavelet used here is 200 m. A wavelet transform can also be employed on profiles of u and υ to detect changes at the top of a well-mixed layer. An example of u, υ, and the vector-summed Haar wavelet (dilation = 200 m, as well) transform operated on profiles of each component separately is shown in Figs. 6a,b. In the wavelet transform of the wind profile, pronounced peaks are apparent both at the MH (1250 m) and near the surface, where friction causes the wind speed to sharply decrease near the ground. These peaks in the wavelet transform of both u, υ, and RCI are used to refine the measurement of the MH.
The location of up to the five largest local peaks in the profiles of and wavelet transforms of both u, υ, and RCI are selected for further analysis. While these peaks are often associated with the MH, their location can be far away ( km) from the turbulent layer. This frequently happens at night, when the mixing layer is shallow but the top of the residual layer, where aerosols are well mixed from the previous day, extends well above the MH. Thus, there is a large gradient in RCI at the top of the residual layer that is not associated with the mixing layer. This can be seen in Fig. 7c, in which the wavelet transform of RCI shows a pronounced peak at 1000–1500 m between 0000 and 1200 UTC when the MH remains below 300 m. To prevent these peaks from being misconstrued as the MH, peaks are considered further only if they are located near (specifically ).
The peaks are incorporated into this fuzzy logic technique by forming a second-generation aggregate. In contrast to the membership functions used in section 4c, which are based on the local value of the quantity, membership functions here are dynamically created for each time step considering the location and size of the peaks. This process is done for each type of peak (i.e., and the wavelet transform of u, υ, and RCI) separately. An example membership function for a corresponding wind profile is shown in Fig. 6c. A membership value of one, indicating solidly within the turbulent mixing layer, is assigned to all heights in the profile below the lowest picked peak. A membership value of zero, indicating above the MH, is assigned to all heights above the highest peak. In between picked peaks, the membership function decreases in steps. The size of the step is the magnitude of the peak divided by the sum of all the peak magnitudes. If no peak is located near in a given profile, then no membership function is created and the input is not used in the second-generation aggregate.
The second-generation aggregate is produced by taking a weighted mean of the membership values used in the first-generation aggregate and the membership values described above for u, υ, , and RCI profiles. Since these indicators are a more direct measure of the depth where aerosols and momentum are well mixed, their membership values are weighted twice as much as those of the direct measurements of turbulence in section 4c, which capture the maximum updraft height and not necessarily the mean mixing depth. An example of the second-generation aggregate and the indicators of mixing are shown in Fig. 7. During well-mixed convective conditions between 1400 and 2200 UTC, the MH can generally be tracked using the three shown indicators. By using these three indicators combined, the limitations of any one are mitigated. For example, at 1500 UTC the clearly shows the continuous growth of the mixing layer, while the wavelet transform on RCI detects only the larger gradient associated with intermittent clouds at 1 km.
Similar to how is found, the final determination of the MH is the lowest height where the second-generation aggregate falls below 0.5, as shown in Fig. 7d. Again, this ensures that the turbulent layer is connected to the surface. In Fig. 8, the difference between the first- and second-generation aggregates is shown, along with and . This shows the impact of incorporating the indicators into the fuzzy logic. At night (0000–1200 UTC), the indicators hardly impact the measurement, since the layer is not well mixed. Conversely, the impact of these indicators is clearly visible during the day, when the measurement of is more continuous than and is generally lower than . Both of these features are related to the fact that turbulence profiles are not sufficient alone to properly measure the convectively well-mixed layer. Thermal plumes extend beyond the mean height where quantities are well mixed, resulting in turbulent motions above the mean MH. This leads to a small high bias in compared to the true MH, which varies in time depending on the altitude to which the convective plume reaches. Additionally, the second-generation aggregate is typically larger above and smaller below , indicating that the second-generation aggregate decreases more slowly with height. This difference is primarily due to the fact that more inputs, each giving a slightly different result, are combined into the second-generation aggregate.
e. Reliability metrics
Since the algorithm produces continuous measurements of regardless of the quality of the input data, it is necessary to create metrics useful in assessing the quality and accuracy of the output. Uncertainty estimates for are created for each time step to reflect the degree of agreement between the various inputs. The lower and upper uncertainty bounds for are set to the lowest height where the second-generation aggregate falls below 0.75 and 0.25, respectively; these bounds are rather arbitrary and could be changed depending on the level of uncertainty between the different inputs desired. When the inputs are in general agreement with each other, such as in the example in Fig. 7d, the uncertainty in is small. Generally, the uncertainty is largest during the afternoon, when the PBL is deep, and during the evening transition period, when the MH can quickly decrease between different scans within the 20-min scan cycle.
The uncertainty bounds themselves are subject to the quality of the input measurements. If the observations themselves are affected by rain or by low aerosol loading, then both the and its uncertainty estimates may be unreliable. Thus, other flags are necessary to assess the validity of the measurement. Rain generally causes a large as a result of differing fall speeds measured based on a distribution of hydrometeor sizes, often leading to an anomalously large . A flag for precipitation is created if the mean column-averaged m s−1 during a vertical stare. This flag is not perfect, as drizzle and snow have smaller fall velocities and may not be detected; additionally, virga may affect the measurement but may not be detected, since it is not affecting the entire column.
Flags are also produced to indicate whether the measurement of is an upper or lower limit if there are no lidar observations below or above , respectively. If the MH is very shallow (<20 m) or nonexistent, no turbulence may be detected by the lidar even at its lowest measurement height. In this case a flag is created, indicating that the reported is an upper limit; that is, there is no measurable mixing at or above , but no data are available below to determine the exact MH or even the existence of a mixing layer. Conversely, it is possible that the lidar inputs indicate continuous mixing up to the reported , but no lidar data are available above to confirm that the true MH is not higher. This happens when the mixing layer is cloud topped and when aerosol loading is low. A cloud-topped mixing layer is indicated if is collocated with a cloud, which are detected in an RCI profile when the wavelet transform is less than −0.5, indicating a large increase in RCI with height. This can be seen in Fig. 7c between 1500 and 1900 UTC at m, where a cloud is collocated with . When aerosol loading is low and the MH is sufficiently deep, the lidar may not be sensitive enough to detect the MH. These periods are identified and flagged as a lower limit for if the layer is turbulently mixed from the surface to the highest usable data (range gate below height where becomes m2 s−2), but there are no coincident gradients in RCI or the other indicators. If is not flagged as a lower limit, upper limit, or cloud topped, then the measurement is “good.”
5. Comparison with in situ measurements
Other field measurements for INFLUX provide occasional opportunities to validate measurements from the Doppler lidar. Occasionally, research aircraft flights are conducted around the Indianapolis area to quantify greenhouse gas emission rates from the urban area (e.g., Heimburger et al. 2017). During some of these flights, the aircraft performs helical ascents and descents over the lidar site, allowing aircraft and lidar profiles to be compared. The aircraft is outfitted with a suite of sensors, including a cavity ring-down spectrometer for in situ measurements of carbon dioxide, methane, and water vapor (Crosson 2008), and a microbead thermistor for temperature observations. Profiles of these measurements can be used to determine the MH, based on where trace gases are well mixed and on the potential temperature . Unfortunately, there are an insufficient number of aircraft profiles to perform a statistical comparison of the MH, since the lidar has been upgraded in the beginning of 2016, after which the measurements are more trustworthy because of the increased sensitivity.
Given the limited dataset, aircraft profile measurements are intercompared with the corresponding Doppler lidar in Fig. 9. Data collected during a helical ascent and a descent, with a circle radius of km, over the lidar site provide two independent profiles in quick succession. During this period the lidar remained nearly constant at 1550 m. Depending on the quantity used to identify the MH, the MH could be interpreted to be located anywhere from 1350 m (from the methane and water vapor profiles during the ascent) to 1770 m (from the carbon dioxide profile during the descent). A strong capping inversion often located near the MH is evident near 1500 m during both flights. The lidar and its uncertainty are located in the middle of the indicated MH from the different profiles, validating the measurement. The large change in the profiles between the ascent and descent, taken minutes apart, demonstrates the uncertainty in using quasi-instantaneous soundings to determine the MH. Any individual profile may not be representative of the mean conditions, especially under strongly convective conditions when the apparent MH at a given location and time may be affected by a local upward or downward perturbation from an updraft or downdraft. This representativeness issue of individual profiles from soundings, specifically of radiosondes, has been discussed by Hennemuth and Lammert (2006). Conversely, the lidar uses profile observations collected over the 20-min scan cycle to produce a single measurement of the MH, yielding a more robust estimate of the mean .
While a case is presented here to demonstrate the validity of measured by Doppler lidar, a statistical comparison of against another in situ or remote sensor would provide more insight into any systematic biases, limitations, and the average error of . At the Indianapolis site, there are no routine measurements available nearby (within km) for a robust statistical analysis. In the future, a more comprehensive analysis is planned through comparison with PBL heights derived from a radar wind profiler and radiosonde measurements.
6. Mixing height statistics in 2016
Since has been measured nearly continuously for an extended time (>1 year), its diurnal and seasonal variability can be uniquely investigated from a statistical and climatological standpoint. This brief analysis highlights the primary advantages of automated determination of continuously using the described method for large datasets. Before the following statistics were calculated, 1.5% of the data were removed when precipitation was detected and flagged. For 2016, the daily afternoon and nocturnal mean have been calculated and are shown in Fig. 10. The annual cycle in the afternoon mean is clearly visible, as it reaches a peak around the summer solstice and reaches a minimum around the winter solstice. The nocturnal shows an opposite cycle, as it is generally deeper in winter than in summer. This can be attributed to the fact that the near-surface winds are typically stronger in winter (not shown), which generally maintains a mechanically mixed PBL (Sun et al. 2012; Bonin et al. 2015; Mahrt et al. 2015), and the sky is cloudier, limiting surface cooling. In addition to the seasonal cycle, there is large variability in the afternoon and nocturnal on a day-to-day basis as they respond to synoptic scale forcings.
With being measured every 20 min, the diurnal cycle of can be analyzed as well. Since the magnitude of the afternoon can vary greatly (–3000 m) depending on the season and synoptic pattern, it is necessary to normalize the values by the daily maximum , , to make the diurnal evolution of comparable among different days. Composites of the daily cycle are centered on sunrise and sunset to understand how the evolves around the morning and evening transitions. Using observations from all of 2016, separate composites are made for low ( m s−1), moderate (5 m s−1 m s−1), and high ( m s−1) wind speeds in the lowest 1 km, as shown in Fig. 11. In interpreting these plots, the composites are the most meaningful in the hours around sunrise and sunset; since the length of daylight varies between 9 and 15 h per day during the year, the composites are a combination of different times of day further from the normalization time. Thus, data 9 h after sunrise are not identical to 3 h before sunset.
The diurnal cycle of is clearly visible in Fig. 11. Generally, begins to increase shortly after sunrise as heating begins and quickly decreases just before sunset as turbulence decays. The nocturnal is generally deeper with increasing wind speeds, as turbulence is mechanically generated and better sustained. Subtle but notable differences are also evident in the evolution in the morning and evening depending on U. When the wind speed is low, grows slowly after sunrise for a few hours compared to windier conditions, after which its depth rapidly increases. This is likely attributed to the fact that a pronounced surface-based temperature inversion often forms over land when winds are weak (Bonin et al. 2015), and most of the initial heating in the morning would be used to erode the inversion over a shallow layer. This inversion does not often form when the mean wind speed is larger (Bonin et al. 2015); thus, heating in the early morning can be dispersed through a deeper layer, causing to increase shortly after sunrise. The timing of the afternoon and evening transition is also dependent on U. On average, reaches a maximum 4.5 h before sunset when m s−1, while reaches a maximum later at 2.5 h before sunset when m s−1. Thus, the evening transitional period is generally delayed when winds are stronger, as mechanically generated turbulence impedes stabilization (Blay-Carreras et al. 2014).
The presented algorithm has many advantages over other individual approaches to determine the MH. This composite technique combines information from all the scan angles, mitigating detection issues when the MH is below the minimum range of a zenith pointed lidar while also maintaining sensitivity to detect a deep MH from the zenith measurements. The fuzzy logic method combines multiple techniques used separately in previous studies to obtain a unified determination of the MH, even though each measurement alone (Haar wavelet on RCI, , etc.) may give a different estimate of the MH. The uncertainty in the MH can be related to how the different measures converge onto the same value of the MH. Since the Doppler lidar is directly able to measure mixing processes, it is straightforward in distinguishing gradients in RCI as the MH or top of the residual layer. While this distinction is possible with backscatter lidars, it is not trivial (de Bruine et al. 2017). Additionally, the weighting of the inputs can be modified based on definition and application of the MH measurements. While gradients in RCI and the horizontal wind are weighted heavily here given the air quality application, as they indicate the mean height to which quantities, including pollutants, are mixed, the direct measures of turbulence could be assigned a larger weighting for a different objective.
Despite these advantages over other traditional approaches to determine the MH, there are still challenges and limitations of this technique that will need to be addressed in the future. The method in which turbulence is measured from the PPI and RHI scans inherently relies on homogeneity in the mean flow (Bonin et al. 2017). While this assumption can be often safely made in Indianapolis, where the terrain is flat, it is often not valid in areas with complex terrain. This heterogeneity leads to inaccurate measures of turbulence and ultimately to poor measurements of . Thus, heterogeneity in the mean flow needs to be distinguished from turbulence features in PPI and RHI scans. Vakkari et al. (2015) propose a method to measure mixing from PPIs in complex terrain, which may be useful when applying this composite method in areas with heterogeneous flow.
While the results here have been shown only from a Halo Photonics Stream Line XR Doppler lidar, this technique is currently being similarly applied to measurements from a Leosphere Windcube 200S, the high-resolution Doppler lidar (Grund et al. 2001), and a custom-built lidar currently under development. With a similar scanning pattern as used here, this algorithm should be applicable to measurements from Doppler lidar systems as well. As the fuzzy logic technique is applied to the different systems operating with varying parameters (pulse width, accumulation time, pulse repetition frequency, etc.) and scan cycles, the membership functions in Table 3 are adjusted. For example, the upper and lower limits for the membership function for are generally lowered for a longer pulse width or accumulation time, as smaller turbulent motions are averaged out by the Doppler lidar. This fuzzy logic technique can also be used with limited inputs, such as only data from zenith stares. When this strategy is implemented, the fuzzy logic combines information from only and RCI to determine a unified measure of but with the same minimum range limitations as any technique from a zenith stare alone.
A new fuzzy logic–based composite technique to determine the MH from Doppler lidar data has been presented. The method is able to take advantage of the strengths of each measurement and scan strategy to overcome the limitations of others. For instance, shallow conical scans can be used to take measurements at high vertical resolution near the surface and compliment zenith stares that can take measurements several kilometers into the atmosphere. Thus, the MH can be detected autonomously under a wide variety of conditions, such as when mixing is only a few tens of meters to when it is several kilometers deep. Within the algorithm nonturbulent submeso motions are spectrally identified so that they are not misconstrued as turbulence, leading to an overestimate of the MH.
Measurements collected for INFLUX are used to demonstrate the algorithm. While sample detailed data from only one day are shown here for brevity, time–height cross sections of wind, SNR, and with the for all of 2016 through the present are available online (https://www.esrl.noaa.gov/csd/groups/csd3/measurements/influx16/) to show the algorithm operation for a variety of conditions. Measurements of the MH from Indianapolis in 2016 are used to investigate and quantify both the annual and diurnal variability of the MH. The afternoon MH was largest around the summer solstice; however, the nocturnal MH tended to be deeper in winter. The timing of the growth and decay of the MH is closely related to the magnitude of the mean wind. The deepening of the MH in the morning is delayed but ultimately more rapid, and the decay is quickened when winds are lighter.
This algorithm can be and is being adapted to different types of Doppler lidars operating different scanning strategies. Efforts are underway to further validate from the Doppler lidar against other measurements, such as from radiosondes and radar wind profilers. These measurements of the MH are being used within the larger context of INFLUX to measure greenhouse gas emissions and to validate numerical weather prediction model output.
We thank Scott Sandberg and Ann Weickmann for their time and effort in deploying, maintaining, and developing supplementary equipment to assist in the automated operation of the lidar in Indianapolis. We thank our hosts at Ivy Tech in Lawrence, Indiana, whose staff generously provided assistance when needed. We also thank James Whetstone for his continued support of INFLUX. Funding for this work was provided by the NIST Greenhouse Gas Measurements Research Program via Interagency Agreement (IAA) M121271 and NOAA’s Earth System Research Laboratory.