Methodology for deriving the telescope focus function and its uncertainty for a heterodyne pulsed Doppler lidar

Doppler lidars provide two measured parameters, radial velocity and signal-to-noise ratio, from which winds and turbulent properties are routinely derived. Attenuated backscatter, which gives quantitative information on aerosols, clouds, and precipitation in the atmosphere, can be used in conjunction with the winds and turbulent properties to create a sophisticated classification of the state of the atmospheric boundary layer. Calculating attenuated backscatter from the signal-to-noise ratio requires accurate knowledge of the telescope focus function, which is usually unavailable. Inaccurate assumptions of the telescope focus function can significantly deform attenuated backscatter profiles, even if the instrument is focused at infinity. Here, we present a methodology for deriving the telescope focus function using a co-located ceilometer for pulsed heterodyne Doppler lidars. The method was tested with Halo Photonics StreamLine and StreamLine XR Doppler lidars but should also be applicable to other pulsed heterodyne Doppler lidar systems. The method derives two parameters of the telescope focus function, the effective beam diameter and the effective focal length of the telescope. Additionally, the method provides uncertainty estimates for the retrieved attenuated backscatter profile arising from uncertainties in deriving the telescope function, together with standard measurement uncertainties from the signal-to-noise ratio. The method is best suited for locations where the absolute difference in aerosol extinction at the ceilometer and Doppler lidar wavelengths is small.

1 Introduction

Coherent Doppler lidar systems are capable of providing accurate radial Doppler velocities at high temporal and spatial resolution and have been employed across a wide range of scientific and operational fields. Meteorological applications include the retrieval of turbulent properties to determine the strength, location, and source of mixing in the atmospheric boundary layer and, with many systems having scanning capability, the retrieval of winds. Information on the targets responsible for the radial Doppler velocities measured by the Doppler lidar (e.g. aerosol, cloud, precipitation) would greatly aid the interpretation of both the velocities and the products derived from them, but this requires quantitative use of the signal power received by the instrument.

The performance of a Doppler lidar depends on the signal-to-noise ratio, SNR, of the system, as SNR determines the radial velocity uncertainty (Rye and Hardesty, 1993; Pearson et al., 2009). The outgoing laser beam can be focused to improve the SNR at ranges close to the focal length (Pearson et al., 2002), and this is often used to improve the Doppler lidar velocity data quality and data availability, particularly in the atmospheric boundary layer. The optimal choice of focus will depend on the atmospheric conditions at the deployment location (Hirsikko et al., 2014).

Knowledge of how the choice of instrument parameters, such as the effective focal length of the telescope, impacts the SNR profile is necessary in order to obtain profiles of the attenuated backscatter coefficient (Zhao et al., 1990). A comprehensive overview of the theoretical considerations in determining the performance of coherent Doppler lidar systems was given by Frehlich and Kavaya (1991), who provided analytical expressions for deriving the expected signal measured by the coherent detector for a given target for a range of instrument configurations, including analytical expressions for the telescope focus function (also termed coherent responsivity). Most analytical expressions assume ideal Gaussian beams, which may not always be appropriate (Hill, 2018); hence experimental approaches have also been used to determine the impact of beam aberrations (Hu et al., 2013).

The profile of the attenuated backscatter coefficient has the potential to be used in real time by weather forecasters (Illingworth et al., 2019), as it can be used in the same manner as for ceilometers. This includes the detection of liquid, supercooled-liquid, mixed-phase, and ice clouds (Hogan et al., 2003; Van Tricht et al., 2014; Tonttila et al., 2015), aerosol layer and mixing-height determination (Flentje et al., 2010; Kotthaus and Grimmond, 2018), and the retrieval of precipitation parameters (Lolli et al., 2018).

In addition to providing velocity estimates for wind and turbulence, the inclusion of the profile of the attenuated backscatter coefficient is advantageous for Doppler lidar boundary layer classification schemes (Tucker et al., 2009; Harvey et al., 2013; Manninen et al., 2018) by enhancing the discrimination between aerosol, cloud, and precipitation, and it can be used for tracking elevated aerosol plumes (Hannon et al., 1999). The combination of attenuated backscatter profiles from coherent Doppler lidars with other profiling instruments permits additional retrievals; for example, together with a ceilometer (Westbrook et al., 2010b), or with a cloud radar (Träumner et al., 2010), it can yield drizzle drop size and precipitation rate. There is also an advantage to obtaining attenuated backscatter and Doppler velocity measurements in the same measurement volume, since this will simplify the calculation of cloud or aerosol mass fluxes (Engelmann et al., 2008).

Therefore, an accurate profile of the attenuated backscatter coefficient requires confidence in the parameters used to generate the telescope focus function. The parameters may not be known a priori or may differ from what is assumed, and incorrect values can result in artefacts and very large biases in the attenuated backscatter coefficient. We present a methodology for deriving the parameters of the telescope focus function experimentally from co-located Doppler lidar and ceilometer observations, together with the uncertainties in the function parameters. The ceilometer, for which the overlap function is known, provides our reference attenuated backscatter profiles. This methodology is relevant for coherent Doppler lidars designed for meteorological applications with maximum ranges suitable for observing the full extent of the boundary layer and beyond. Note that a calibration constant may still need to be determined and applied after implementing the calculated telescope focus function to retrieve the profile of the attenuated backscatter coefficient (Westbrook et al., 2010a; Chouza et al., 2015).

The theoretical description of the telescope focus function is outlined in Sect. 2. In Sect. 3, we introduce the instruments and the methodology for deriving the parameters of the telescope focus function experimentally. An iterative least-squares regression using weighted mean square error (MSE) is used to find the best solution for the telescope focus function, where the weights represent the measurement uncertainties in both instruments. The use of long time periods (1 year or more) also provides an estimate of the uncertainties in the parameters for the telescope focus function, which can then be propagated through to uncertainties in the retrieved attenuated backscatter coefficients. The methodology is applied to different instruments in multiple locations in Sect. 4, and the validation of the method is presented in Sect. 5.

2 Theory

2.1 Telescope focus function

Following Frehlich and Kavaya (1991), the coherent Doppler lidar equation can be expressed as

$$\begin{array}{}\text{(1)}& \mathrm{SNR}\left(R\right)={\displaystyle \frac{\mathit{\eta }cE}{\mathrm{2}h\mathit{\nu }B}}{\displaystyle \frac{A_{\mathrm{e}}\left(R\right)}{R^{\mathrm{2}}}}{\mathit{\beta }}^{\prime }\left(R\right),\end{array}$$

where SNR is the signal-to-noise ratio, varying as a function of range, R, from the instrument; β^′ is the attenuated backscatter coefficient; η is the detector quantum efficiency; c is the speed of light; E is the beam energy; h is Planck's constant; ν is the optical frequency; B is the receiver bandwidth; and A_e is the effective receiver area.

For a monostatic system emitting a circular Gaussian beam, using a circular aperture, and having matched filters, the effective receiver area is given by Frehlich and Kavaya (1991) and Henderson et al. (2005):

$$\begin{array}{}\text{(2)}& A_{\mathrm{e}}\left(R\right)={\displaystyle \frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}\left(\mathrm{1}+{\left(\frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}\mathit{\lambda }R}\right)}^{\mathrm{2}}{\left(\mathrm{1}-\frac{R}{f}\right)}^{\mathrm{2}}+{\left(\frac{D}{\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}}\right)}^{\mathrm{2}}\right)}},\end{array}$$

where D is the 1∕e² effective diameter of a Gaussian beam, λ is the laser wavelength, f is the effective focal length of the telescope for the transmitter and receiver, and ρ₀ is a turbulent parameter, also termed transverse field coherence length.

Collecting the range-dependent terms, we obtain a unitless telescope focus function:

$$\begin{array}{}\text{(3)}& T_{\mathrm{f}}\left(R\right)={\displaystyle \frac{A_{\mathrm{e}}\left(R\right)}{R^{\mathrm{2}}}},\end{array}$$

which is also termed the coherent responsivity (Frehlich and Kavaya, 1991).

The profile of the attenuated backscatter coefficient is then obtained by rearranging Eq. (1):

$$\begin{array}{}\text{(4)}& {\mathit{\beta }}^{\prime }\left(R\right)={\displaystyle \frac{\mathrm{2}h\mathit{\nu }B}{\mathit{\eta }cE}}{\displaystyle \frac{\mathrm{SNR}\left(R\right)}{T_{\mathrm{f}}\left(R\right)}}.\end{array}$$

Figure 1 Telescope focus functions for (a) varying f with D=70 mm, (b) varying D with f=1000 m, and (c) varying D with f being infinity.

Figure 1a shows how T_f(R) depends on the telescope focal length, f, and Fig. 1b how T_f(R) depends on D. Both figures show that the apparent focus – i.e range to the T_f(R) maximum – is always closer than f and that decreasing D shortens the apparent focus. This makes estimation of the parameters by eye in T_f(R) prone to errors, since the apparent focus cannot be translated into f without knowledge of D.

Figure 1c shows that even if the telescope is focused at infinity, knowledge of the D is essential to derive attenuated backscatter coefficient profiles. While the gradient of T_f(R) may be independent of D at the near and far ranges, the relative magnitude is not, and the potential variation is high in the range of the profile that is commonly of most interest.

2.2 Uncertainty in the attenuated backscatter coefficient

Assuming that the parameters T_f(R) and SNR are independent, and have uncertainties that can be described as Gaussian, the relative random uncertainty in the attenuated backscatter coefficient is

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{{\mathit{\beta }}^{\prime }}=\sqrt{{\mathit{\sigma }}_{\mathrm{S}}^{\mathrm{2}}+{\mathit{\sigma }}_{T_{\mathrm{f}}}^{\mathrm{2}}},\end{array}$$

where σ_S is the relative uncertainty in the Doppler lidar SNR, and ${\mathit{\sigma }}_{T_{\mathrm{f}}}$ is the relative uncertainty in T_f(R). An expression for deriving σ_S is given by Manninen et al. (2018), and we describe our method for obtaining ${\mathit{\sigma }}_{T_{\mathrm{f}}}$ in Sect. 4.2.

3 Application to data

There are three range-dependent unknowns in Eq. (2): f, D, and ρ₀. We first assume that we can neglect ρ₀ and then describe a method for estimating f and D, together with their uncertainties, which can then be propagated to obtain the uncertainty in the attenuated backscatter coefficient. The impact of the ρ₀ parameter is discussed in Sect. 4.3.

3.1 Instruments

We used measurements taken from the U.S. Department of Energy Atmospheric Radiation Measurement (ARM, Mather et al., 2016) observatories. We selected five sites with co-located ceilometer and Doppler lidar instruments: Southern Great Plains, US (SGP); tropical west Pacific, Darwin, Australia (Darwin); Barrow, Alaska, US (North Slope of Alaska, NSA); Graciosa, Azores (Graciosa); and Ascension Island, Atlantic, UK (Ascension).

The Doppler lidars operated by ARM comprise both Halo Photonics StreamLine and StreamLine XR versions. These are commercially available heterodyne pulsed systems capable of full hemispheric scanning and operated at a temporal resolution of 1–2 s (see Table 1). The focus for the StreamLine version can be set by the operator, whereas the StreamLine XR has the focus set by the manufacturer; however ARM has had some instruments upgraded from their original specification.

The ceilometer at all sites was a Vaisala CL31 ceilometer, which has a coaxial design and full overlap before 100 m and a temporal resolution of 30 s (more specifications given in Table 2).

Table 1 Halo Photonics StreamLine and StreamLine XR heterodyne Doppler lidar specifications. Values in parentheses refer to the specification of the Doppler lidar during the first period in Darwin.

Table 2 Vaisala CL31 ceilometer specifications.

3.2 Methodology

3.2.1 Telescope focus function parameter estimation

The methodology for deriving the parameters of the telescope focus function compares profiles from a co-located Doppler lidar and ceilometer using an iterative least-squares regression to find the best solution. The method follows the process diagram given in Fig. 2.

Figure 2 Process diagram of the telescope focus function parameter estimation.

Before input, the Doppler lidar SNR data had a background correction applied to reduce bias (Manninen et al., 2016). Both ceilometer and Doppler lidar data were averaged to a common 30 min, 30 m vertical resolution grid. If the Doppler lidar vertical resolution was larger than 30 m (as was the case for one period from Darwin), linear interpolation was used to match resolutions. After averaging, data below a minimum threshold of −22.2 dB (Manninen et al., 2018) were discarded. The threshold is based on the expected noise floor for the instruments considered here (Halo StreamLine and StreamLine XR) and should probably be modified for different instruments.

The data were then filtered to select only those portions of the profiles that are considered reliable for comparison. Ceilometer data below 195 m were discarded to ensure that only data with full overlap were used.

Due to the wavelength difference between the Doppler lidar and the ceilometer, it cannot be assumed that the atmospheric backscattering properties are the same at both wavelengths. However, we are only interested in the profile shape, not the absolute values, so profiles from the Doppler lidar and the ceilometer can be compared as long as they contain only one type of scatterer, and one which can be assumed to be distributed homogeneously throughout the portion of the profile used for comparison. Hence, the portion of a profile selected for comparison should contain only one aerosol layer, no clouds, and no precipitating hydrometeors.

We removed clouds by identifying the range gate 150 m below the cloud base detected by the ceilometer and excluding all data beyond this. Elevated aerosol layers and precipitating hydrometeors were filtered out by identifying layers using a convolution of the ceilometer profile with a haar-wavelet to detect changes in the gradient. The base of the second layer was identified where the gradient was increasing over two range gates, and all data above this were discarded. This process also eliminates noisy profiles with low SNR.

The T_f(R) parameter estimation is performed on a profile-by-profile basis for each profile where the filtering process leaves eight successive range gates present. From Eq. (4), dividing the Doppler lidar SNR profile with the appropriate T_f(R) will generate a Doppler lidar attenuated backscatter profile whose shape should match that of the ceilometer attenuated backscatter profile.

We use a brute-force approach to iterate through a range of reasonable f and D values, generating a corresponding T_f(R) and Doppler lidar attenuated backscatter profile for each combination of values. The ceilometer profile and resulting Doppler lidar attenuated backscatter profiles are normalised so that the integral value of the unfiltered portion is unity. We then use a least-squares regression using weighted MSE to find the best solution (smallest MSE), where the weights represent the measurement uncertainties in both instruments.

Collecting results over many profiles results in a bivariate distribution; the peak of this distribution is chosen as the best estimate of f and D, and hence the best estimate of T_f(R), using Eq. (3).

3.2.2 Outlier removal

Occasionally, data of poor quality pass the filtering step in Fig. 2. The most common issues are noisy ceilometer data and a bias in the Doppler lidar SNR profiles. If not screened, these occasional profiles result in significantly altered T_f(R) estimation. Any noise in the ceilometer data is magnified by the profile length often being relatively short, and hence large uncertainty in even a single range gate can skew the regression. Doppler lidar SNR bias will impact the normalisation process, changing the T_f(R) selected by the method due to the now incorrect profile shape. Due to the non-linearity of the T_f(R) parameter estimation process, these issues result in regression solutions wildly inconsistent with the estimates based on good data. These outliers, which do not fall within the normal uncertainty observed in good data, are then removed from the bivariate distribution of solutions before calculating the uncertainty estimates.

We used the median absolute deviation, MAD (Huber and Ronchetti, 2009; Leys et al., 2013), to distinguish outliers in the bivariate distribution of estimated f and D. MAD can be calculated using

$$\begin{array}{}\text{(6)}& \mathrm{MAD}=b\mathrm{med}\mathit{\left\{}\right|x_i-\mathrm{med}\mathit{\left\{}{x}_i\mathit{\right\}}\left|\mathit{\right\}},\end{array}$$

where b=1.4826 when the distribution excluding the outliers is normal. However, the distribution of f and D may not meet this criterion due to the non-linearity of T_f(R) and the computational T_f(R) estimation process. We expect the distributions of D and f⁻² to be close to normal and will use f⁻² rather than f to determine outliers. Additionally, the peak of the bivariate distribution may not always coincide with the medians of the univariate D and f⁻² distributions, and, hence, we use a modified form of Eq. (6),

$$\begin{array}{}\text{(7)}& \mathrm{MAD}=b\mathrm{med}\mathit{\left\{}\right|x_i-\mathrm{peak}\mathit{\{}{x}_i,y_i\mathit{\}}\left|\mathit{\right\}}.\end{array}$$

We selected three MADs as the threshold for flagging outliers:

$$\begin{array}{}\text{(8)}& \sqrt{{\left({\displaystyle \frac{f_i^{-\mathrm{2}}-\mathrm{med}\mathit{\left\{}{f}_i^{-\mathrm{2}}\mathit{\right\}}}{{\mathrm{MAD}}_{f^{-\mathrm{2}}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{D_i-\mathrm{med}\mathit{\left\{}{D}_i\mathit{\right\}}}{{\mathrm{MAD}}_D}}\right)}^{\mathrm{2}}}\ge \mathrm{3}.\end{array}$$

Assuming the distribution excluding the outliers to be normal, three MADs correspond to 3 standard deviations of the underlying distribution. In cases where f is at infinity, all estimates with a finite f will be flagged as outliers.

Table 3 Best estimates of f and D together with their uncertainties for Doppler lidars at five ARM sites. Total estimates refers to the number of ceilometer and Doppler lidar profiles suitable for comparison. Good estimates refers to the number of estimates remaining after outlier filtering with MAD.

Figure 3 Distributions of T_f(R) parameter estimates from (a) Darwin from 21 June 2011 to 22 July 2012, (b) SGP from 1 January 2015 to 2 May 2016, and (c) SGP from 3 May 2016 to 5 June 2017. Outliers filtered using MAD ≥3 are marked in black.

4 Results

4.1 Parameter estimation

We applied the T_f(R) estimation method to Doppler lidars at five ARM atmospheric observatories. Figure 3a shows the distribution of f and D calculated for the Doppler lidar operating at Darwin in northern Australia between 21 June 2011 and 22 July 2012. This Doppler lidar is a StreamLine and the distribution of f is positively skewed, as explained in Sect. 3.2.2. The distribution displays a slightly wider peak than expected for a normal distribution.

Figure 3b shows the distribution of f and D for the StreamLine Doppler lidar operating at SGP from 1 January 2015 to 2 May 2016. The distribution close to the peak is really tight, while the outliers have substantial spread. Many of the poor estimates responsible for the outlier spread occur during January and February in both years, while for the rest of the period the estimates are remarkably consistent. On 3 May 2016 the Doppler lidar at SGP was changed to a StreamLine XR, and Fig. 3c shows the distribution of f and D from 3 May 2016 to 5 June 2017. The change in instrument version, from StreamLine to StreamLine XR, is clearly seen in the change in D, whereas the best estimate for f did not change. However, inspecting the data by eye would suggest a significantly shorter apparent focus, and the T_f(R) calculated using the best estimates for f and D also exhibits a significantly shorter apparent focus. Consequently, the StreamLine XR in SGP has been noted to suffer from poor SNR at the boundary layer top.

The bivariate distributions of f and D show notable variations in how tight they are around the peak, likely a result of differences in data quality between the instruments. The best estimates of f and D and their uncertainties for all sites are presented in Table 3. The Doppler lidar measurements at Darwin were split into two periods, as there was a 2-month break in the measurements between these two periods. We performed the T_f(R) parameter estimation separately for both periods. The best estimates from these periods differ from each other, which is expected as some adjustments were made to the instrument. The telescope focal length for this instrument is directly adjustable by any operator, while the beam diameter is set by the manufacturer and is not modifiable by the operator. We note that the D estimates are the same for these two periods within the margin of error calculated.

For the sites and instruments selected here, only the Doppler lidar at NSA had f set to infinity. In fact, all StreamLine lidars have D in the vicinity of 25 mm, whereas D for the StreamLine XR versions is about half this. Nevertheless, the variation between instruments of the same version is not negligible and should be taken into account when calculating T_f(R) and then attenuated backscatter.

The final step to obtain attenuated backscatter profiles is to apply a calibration constant, which can be achieved using the liquid cloud calibration method (Westbrook et al., 2010a; O'Connor et al., 2004).

Figure 4 (a) Doppler lidar attenuated backscatter coefficient assuming a generic T_f(R), (b) corrected Doppler lidar attenuated backscatter coefficient, and (c) ceilometer attenuated backscatter coefficient, from Darwin on 28 May 2012.

The parameters f and D calculated for period 1 in Darwin have been used to derive T_f(R) and the results applied in Fig. 4. This shows the utility of the method, able to provide reliable Doppler lidar attenuated backscatter profiles in Fig. 4b that show no overcorrection below 1 km and display similar in-cloud values to the ceilometer in Fig. 4c. It is expected that the aerosol attenuated backscatter coefficients will differ due to the different scattering properties of aerosol at the different wavelengths; the scattering properties of cloud droplets remain similar at the two wavelengths (O'Connor et al., 2004; Westbrook et al., 2010a).

Figure 5 Distributions of the Monte Carlo simulation (MCS) input values used for calculating ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$. Values are obtained from (a) resampling, (b) assuming $N(f^{-\mathrm{2}},D)$, and (c) assuming N(f,D). All distributions contain 5337 samples.

Table 4 ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ uncertainty envelopes generated using MCS with three different sets of input values.

4.2 Uncertainty

A computational method was used to calculate the uncertainty in the estimated T_f(R) as it is a non-linear function of f and D. We used Monte Carlo simulation (MCS) (Morgan and Henrion, 1990), where a distribution of input values is fed into a model, here the effective receiver area Eq. (2), and the uncertainty is obtained from the distribution of the output. The input values can be created either from observed statistics or by bootstrapping, i.e. resampling the data. We created three different sets of input values for our MCS:

1. Resampling the individual estimates of f and D provided directly by the T_f(R) estimation method (i.e. those displayed in Fig. 3) after excluding outliers.

2. Generating the values from the statistics presented in Table 3, assuming that D and f⁻² are normally distributed and independent, $N(f^{-\mathrm{2}},D)$.

3. Generating the values from statistics presented in Table 3, assuming D and f to be normally distributed and independent, N(f,D).

For each set of input values, the relative uncertainty in T_f(R) is calculated as

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)={\displaystyle \frac{{\mathit{\sigma }}_{T_{\mathit{\varphi }}}\left(R\right)}{T_{\mathrm{f}}\left(R\right)}},\end{array}$$

where ${\mathit{\sigma }}_{T_{\mathit{\varphi }}}\left(R\right)$ is expressed in terms of the mean-squared deviation of the MCS-simulated telescope focus function, T_ϕ(R), from the best estimate of T_f(R),

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{T_{\mathit{\varphi }}}\left(R\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}{\mathrm{\Sigma }}_{i=\mathrm{1}}^N\left(T_{{\mathit{\varphi }}_i}\right(R)-T_{\mathrm{f}}(R\left)\right)},\end{array}$$

to avoid underestimating the uncertainty resulting from the asymmetry in T_ϕ(R). This also allows us to estimate the impact of refractive turbulence on the uncertainty estimate.

Examples of the three input parameter distributions are presented in Fig. 5. Resampling (Fig. 5a) is the most accurate method as it does not require assumptions about the parameter distributions and their independence. We recommend resampling as the primary method for uncertainty calculation. Using the $N(f^{-\mathrm{2}},D)$ distribution (Fig. 5b) produces a set of input values that appear to be a reasonable approximation, except that the distribution is not as tight around the peak. Using the N(f,D) distribution (Fig. 5c) produces a set of input values that tend to overemphasise shorter values of f and underemphasise higher values. We also note that the central bin of the resampled distribution contains 50 % more samples than the central bin of the statistically generated distributions do. We presume that this is a consequence of the variation in SNR not being necessarily normally distributed.

Figure 6 Relative telescope focus function uncertainties, ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$, generated using MCS with three different sets of input values for (a) Darwin from 21 June 2011 to 22 July 2012 and (b) NSA from 30 July 2014 to 31 December 2017.

Figure 6a displays ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for Darwin showing the range dependence of the uncertainty, with much larger uncertainties for ranges close to either side of the focus (f=590 m). The profile of uncertainties obtained with each set of MCS input values exhibit a similar shape, with $N(f^{-\mathrm{2}},D)$ being closer to resampling than N(f,D) in the near field.

Figure 6b displays ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for the Doppler lidar at NSA, which has f set to infinity; therefore ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ is only dependent on the uncertainty in D. Note the reduced uncertainties around 200–400 m, which are expected when examining Fig. 1c.

The largest value of ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ provides the uncertainty envelope value for each site, which is presented in Table 4. Resampling provides values ranging from 0.12 for the updated instrument at SGP to 0.32 at NSA. MCS values created using $N(f^{-\mathrm{2}},D)$ provided similar values, whereas MCS using N(f,D) often provided much larger uncertainties.

Figure 7 Impact of turbulent parameter, ρ₀, on the telescope focus function, T_f(R), and relative uncertainties, ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$, for different ${C_n}^{\mathrm{2}}$ profiles. (a) Selected profiles of ${C_n}^{\mathrm{2}}$ with range; (b) T_f(R) and (c) ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for Darwin from 21 June 2011 to 22 July 2012; (d) T_f(R) and (e) ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for NSA from 30 July 2014 to 31 December 2017.

4.3 Impact of refractive turbulence

So far we have neglected the potential impact of turbulence on T_f(R) arising from the refractive turbulent parameter, ρ₀, in Eq. (2). An expression for ρ₀ is given in Frehlich and Kavaya (1991),

$$\begin{array}{}\text{(11)}& {\mathit{\rho }}_{\mathrm{0}}\left(R\right)={\left[H{k}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{R}{\int }}{C_n}^{\mathrm{2}}\left(z\right)(\mathrm{1}-z/R{)}^{\mathrm{5}/\mathrm{3}}\mathrm{d}z\right]}^{-\mathrm{3}/\mathrm{5}},\end{array}$$

where H=2.914383, $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$, and ${C_n}^{\mathrm{2}}\left(z\right)$ is the refractive turbulence at range z. We chose three profiles with constant ${C_n}^{\mathrm{2}}\left(z\right)$ and a realistic vertical profile based on the most turbulent case presented by Roadcap and Tracy (2009). Figure 7 shows the impact that these different profiles have on T_f(R) and the resulting resampling calculation of ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for two Doppler lidar instruments with different T_f(R). Values of ${C_n}^{\mathrm{2}}$ up to ${\mathrm{10}}^{-\mathrm{14}}{\mathrm{m}}^{-\mathrm{2}/\mathrm{3}}$ have negligible impact on T_f(R), and even the realistic profile only showed a slight increase in ${\mathit{\sigma }}_{T_{\mathrm{f}}}\left(R\right)$ for the instrument with a focus set closer than infinity. This suggests that the impact of turbulence can be safely neglected for low values of ${C_n}^{\mathrm{2}}$, and for most applications, it can also be neglected when operating in the vertical. Hence, turbulence has no significant impact on the methodology described here for deriving the parameters f and D and their uncertainties from vertical profiles but can be included for completeness.

However, it is clear that the turbulent impact should not be ignored when measuring at low elevation angles close to the horizon, where a profile similar to ${C_n}^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{13}}{\mathrm{m}}^{-\mathrm{2}/\mathrm{3}}$ may be possible. Figure 7 shows that such a profile has a major impact on T_f(R), especially in the far range. In these cases, the parameters f and D obtained from vertical measurements are still applicable, but ρ₀(R) must also be calculated or estimated in order to derive the profile of attenuated backscatter, β^′(R).

Figure 8 Doppler lidar attenuated backscatter coefficient and apparent lidar ratio, ηS, from (a) Darwin on 8 May 2012 and (b) NSA on 20 August 2014.

5 Validation

The liquid cloud calibration method (O'Connor et al., 2004; Westbrook et al., 2010a) is used to determine a calibration constant by integrating attenuated backscatter profiles containing fully attenuating liquid clouds, which have a well-constrained apparent lidar ratio, ηS, where η is a multiple scattering factor and S is the lidar ratio. In the absence of multiple scattering, ηS can be assumed to be independent of the height of the cloud.

This calibration method can be used to evaluate the estimated T_f(R) for Doppler lidar by checking whether the attenuated backscatter profiles obtained for the Doppler lidar after applying T_f(R) indeed provide similar ηS values for liquid clouds at different heights.

Figure 8 shows examples of Doppler lidar attenuated backscatter profiles after calibration and the derived apparent lidar ratio at two sites, Darwin and NSA. These sites have different values of f: Darwin has f=590 m and NSA has f set to infinity. For both cases, liquid clouds are present throughout the day, with altitudes varying from 2 to 6 km. When fully attenuating liquid clouds are present, the apparent lidar ratio is close to the expected value of 20 sr, regardless of the height of the cloud, thus confirming that the method of estimating T_f(R) is valid.

Limitations

Table 3 shows that the proportion of data that can be used for the T_f(R) parameter estimation varies significantly from site to site. Over a third of the available profiles from SGP are used, whereas only 0.3 % pass the filtering for Ascension. The lack of suitable profiles at Ascension is explained by the almost constant low cloud cover at this site, with very few profiles having a sufficient number of successive range gates.

Data quality is also a limiting factor, so at sites with very low aerosol optical depth, AOD, such as NSA, the Doppler lidar SNR decreases so rapidly that again there are few profiles having a sufficient number of successive range gates. Low AOD also impacts the performance of the ceilometer, with 48 % of the estimates at NSA discarded as outliers even after the initial filtering was performed. While the outlier removal can separate the good and the poor estimates, the largest uncertainty in D was at NSA. We attempted to perform the T_f(R) parameter estimation on a Doppler lidar from an ARM campaign in Cape Cod but could not obtain reliable estimates due to the low SNR of the ceilometer data.

The T_f(R) parameter estimation method is suitable only in situations where there is minimal difference in atmospheric extinction within the aerosol layer between the two instrument wavelengths of 910 and 1500 nm. Using AERONET (Aerosol Robotic Network) AOD measurements co-located at the ARM atmospheric observatories, the median difference in AOD at 870 and 1640 nm varied between 0.016 and 0.027, which should correspond closely to what might be expected for the difference between the ceilometer and Doppler lidar. Very occasional periods of notable AOD differences were observed at some sites, but including these profiles in time series extending beyond a year will have negligible impact on the T_f(R) parameter best estimate. However, there were breaks in the AOD measurements, and some periods experiencing a significant differential extinction may have gone unnoticed. An additional filter using AERONET AOD measurements to remove profiles experiencing significant differential extinction could be included in Fig. 2 for those sites where this may be an issue.

6 Conclusions

We have developed a method for deriving the telescope focus function and its uncertainty for pulsed heterodyne Doppler lidars, and we applied the method to Halo Photonics StreamLine and XR Doppler lidars. The method compares profiles of the Doppler lidar SNR to profiles of the attenuated backscatter coefficient from a co-located ceilometer, producing estimates for two parameters of the T_f(R): the effective focal length for the telescope, f; and the 1∕e² effective diameter of a Gaussian beam, D. This method was developed because it also provides uncertainties in f, D, and T_f(R), which are necessary for quantitative use of the subsequently derived attenuated backscatter profiles. The method can be used to check the manufacturer specifications for these parameters, calculate them if not known, and also check their stability over time. The method does not necessarily require a permanently co-located ceilometer as the estimates of f and D can be made from a short time series with a co-located ceilometer; the Doppler lidar can then be moved to another location without a co-located ceilometer using the T_f(R) estimated previously, assuming that there has been no modification to the instrument parameters.

The method was applied to data from Doppler lidars with different configurations deployed at five ARM sites. Relative uncertainties in f for these instruments ranged from 6 % to 17 %, with the median uncertainty being 10 %; the relative uncertainty in D ranged from 2 % to 12 %, with a median of 3 %. The uncertainty in T_f(R) was calculated using Monte Carlo simulation, using three methods to prepare the input values. We recommend the direct resampling method, but reasonable results were obtained using statistically derived input values assuming a normal distribution. The envelope of relative uncertainties in T_f(R) ranged from 13 % to 32 % and depended on both the instrument configuration and the instrument location. We also show that, even for a Doppler lidar with the focus set at infinity, uncertainty remains in estimating T_f(R) arising from the uncertainty in D. The method was validated by calculating the apparent lidar ratio of fully attenuating liquid clouds from T_f(R) corrected profiles of Doppler lidar attenuated backscatter.

The impact of turbulence on T_f(R) was also investigated and was found to have no significant impact on the methodology described here for deriving the parameters f and D and their uncertainties from vertical profiles. However, the turbulent impact should not be ignored when measuring at low elevation angles close to the horizon, as it can modify T_f(R) considerably, especially in the far range. In these cases, the parameters f and D obtained from vertical measurements are still applicable, but the turbulent contribution to T_f(R) should be included when deriving the attenuated backscatter coefficient.

The T_f(R) estimation method is suitable only for conditions where the differential extinction at the two wavelengths of the Doppler lidar and the ceilometer is small, which can be identified, for example, using AOD from co-located AERONET observations.