Boundary layer height above the Great Barrier Reef studied using drone and Mini-Micropulse LiDAR measurements

2.1. Drone sampling

Drone-based meteorological profiles were obtained between late February and early April 2023. The flights, conducted using a DJI Matrice 300 (M300) quadcopter drone, reached altitudes of up to 1.5 km. Therefore, all flights were beyond visual line of sight and approved by the Australian Civil Aviation Safety Authority. Flights were conducted daily at 09:00, 12:00 and 15:00 hours Australian eastern standard time. Temperature, pressure and relative humidity were measured on the drones using InterMet Systems iMet-XQ2 UAV (iMet) sensors. For a more in-depth explanation of the drone sampling methodology, refer to Eckert et al. (2023). The 3-D position of the mounted iMet sensors and flight patterns (ascent, mission peak and descent) are determined using the drone’s flight logs rather than relying on barometric measurements or the built-in global navigation satellite system recordings from the atmospheric sensors.

Potential temperature Θ (K) is the temperature that an unsaturated parcel of dry air would have if it were brought adiabatically and reversibly from its initial state to a standard pressure level (typically 1000 hPa), as shown in Eqn 1. Because the difference between the potential temperature of unsaturated moist air and that of dry air remains small, no distinction needs to be made for most practical purposes (World Meteorological Organization 1966):

where R ÷ c_P = 2 ÷ 7, R represents the ideal gas constant for dry air (8.314 J K⁻¹ mol⁻¹) and c_p represents the specific heat capacity at constant pressure (29.14 kJ kmol⁻¹ K⁻¹). To determine the PBLH from the drone Θ profile, using the 1.5-theta method, we first calculate the mean Θ in the lowest 300 m of the profile (Θ_300 m), as Θ varied little in range 0–300 m and it was expected that the PBLH would be significantly above 300-m altitude during the daytime in this location. The PBLH is then identified as the first altitude at which Θ = Θ_300 m + 1.5 K (Nielsen-Gammon et al. 2008; Hu et al. 2010).

2.2. MP LiDAR sampling

The MP LiDAR used at One Tree Island (Droplet Measurement Technologies; https://www.dropletmeasurement.com/product/mini-mpl/) was located on the upper part of the beach in front of the main research station. The laser is green (532-nm wavelength), eye-safe and operating at a frequency of 2500 Hz. The integration time was 30 s and vertical resolution was 30 m as recommended by the manufacturer. We utilise data from both co-polarised and cross-polarised channels. The MP LiDAR scans were taken at 90, 45 and 0° relative to the horizon. Horizontal scans were used employed for calibration, and we use only the vertical scans for PBLH analysis in this work.

Several corrections must be applied to the MP LiDAR signal before proceeding, as outlined in Campbell et al. (2002). Firstly, the background signal and after-pulse effect must be subtracted. The background signal accounts for any ambient light at 532 nm, whereas the after-pulse signal arises from detector saturation due to the laser flash, at the very beginning of a sampling period. Raw MP LiDAR signals are also divided by an overlap function that accounts for the drop off in near-field signal due to the combination of a narrow receiver field of view and a wide transmitting aperture. The overlap range of MP LiDAR signals is typically ~5 km. Overlap and after-pulse calibrations were carried out routinely throughout the campaign period following the procedures outlined in Welton et al. (2000, 2002). The dead time correction, resulting from saturation caused by photon avalanche events, was unaltered from that supplied by the manufacturer.

The raw MP LiDAR signal P as a function of range r comprises contributions from aerosol scattering, Rayleigh scattering, background (mostly solar) photons at 532 nm, and instrumental noise:

where β denotes backscatter coefficients, σ denotes extinction coefficients, subscripts R and A indicate Rayleigh and aerosol scattering, N_B denotes the 532-nm background atmospheric contribution and E relates to the laser energy. The raw MP LiDAR signal is a combination of co-polarised and cross-polarised return, calculated according to Córdoba-Jabonero et al. (2021) and Flynn et al. (2007):

The normalised relative backscatter (NRB) is calculated from the raw signal by subtracting N_B and the after-pulse signal (A(r)) then dividing by the overlap function O(r) and finally multiplying by the square of the range (r²). Key sources of uncertainty associated with NRB calculation include the choice of after-pulse and overlap calibration functions and laser energy fluctuations (Welton and Campbell 2002). Overall uncertainty in the NRB retrieval is typically less than 5% (Welton et al. 2000).

2.3. MP LiDAR cloud filtering

To ensure we sampled the PBLH in clear-sky conditions, we implemented a cloud-filtering algorithm based on the raw LiDAR signals and a co-located all-sky cloud camera (Solmirus All Sky Imaging System M1v, https://solmirus.com/asis-m1v). To validate cloud detection from the MP LiDAR signals, we looked for cloud-masked pixels in a circular region with 5° radius from the centre of each all-sky camera image (Fig. 2). Given that the receiving field of view of the MP LiDAR is 220 mrad (0.01°), this results in a very conservatively strict cloud-flagging. The cloud-masking in the all-sky images were determined by the all-sky camera manufacturer’s (Solmirus) internal algorithm, with cloud-flagged pixels shown in white in the top row of Fig. 2.

Fig. 2.

Cloud detection using the all-sky camera and Mini-Micropulse (MP) LiDAR. Top row shows all-sky camera images from 27 February 2023 at OTI for scenes that are (a) predominantly cloud obscured, (b) slightly cloud obscured and (c) cloud-free. Colours indicate the ratio of green to blue pixels with cloud-flagged pixels in white. Red circles indicate the 5° radius region in the centre of each image, relevant for validating the MP LiDAR cloud-detection algorithm. Bottom row shows MP LiDAR raw count vertical profiles, and the gradient of this signal, at the same time as each camera image. Red dashed lines indicate flagging of a cloud-base in the MP LiDAR signal. Note that cloud-flagging was carried out through the depth of the troposphere (MP LiDAR signal 0–10-km altitude) but only the lowest 3 km of the profile is shown for clarity here.

Here, we develop and test a threshold-based MP LiDAR detection algorithm that, as in Wang and Sassen (2001), is based on the raw LiDAR signal rather than the range-corrected signal or calibrated NRB. The algorithm begins by calculating the derivative of the entire raw LiDAR signal. The raw MP LiDAR signal typically decreases with altitude, and a cloud or aerosol layer base may be present when the signal slope increases with altitude (gradient becomes positive). Example vertical profiles of the raw MP LiDAR signal and gradient on 27 February 2023 are shown in the bottom panels of Fig. 2. To exclude single-layer noise spikes, we only flag a potential layer base if the gradient exceeds zero in three successive layers (over 90-m altitude). We then examine the characteristics of the signal and the gradient in the 20 altitude bins above each flagged layer base. Guided by threshold characteristics in Wang and Sassen (2001), we tuned the following threshold values by comparing to the cloud camera images. A cloud layer is flagged if (a) the layer’s maximum raw signal exceeds 0.4 counts μs⁻¹ and the minimum gradient is less than −0.5 counts μs⁻¹ m⁻¹, (b) the layer peak’s raw signal is more than twice the raw signal at the base, and the altitude is below 4 km or (c) the layer’s peak’s raw signal is more than 1.2 times the raw signal at the base, and the altitude is above 4 km. Once a cloud is identified in a profile, the algorithm goes on to search the remainder of the profile and can retrieve multiple cloud layers in one profile. The cloud-flagging is sensitive to predominantly cloudy scenes such as Fig. 2a where >50% of the all-sky camera’s centre pixels are marked as cloud, and to scenes were as little as 2% of the centre pixels are cloud marked (Fig. 2b). Sky scenes with no cloud-marked pixels, such as Fig. 2c, do not result in a cloud detection in the MP LiDAR profile.

2.4. MP LiDAR PBLH calculation

We used 30-min averaged MP LiDAR NRB profiles to calculate the PBLH, encompassing the time taken to complete a drone ascent and descent. The simplest method employed for calculating the PBLH from cloud-free MP LiDAR NRB profiles is the gradient method. The PBLH_Grad was determined to be the altitude of the most negative NRB signal gradient. An example drone potential temperature profile, including the 1.5-theta PBLH, is shown in Fig. 3a. Alongside this, an example vertical NRB profile, taken on 6 March 2023 (Fig. 3b), shows the gradient of that signal, indicating the PBLH detection at the minimum value of the gradient. Signal-to-noise in MP LiDAR NRB profiles decreases with altitude, so to isolate the minimum gradient applicable to the PBLH rather than noise, we took the minimum gradient between 0- and 3-km altitude. As there are no sources of uncertainty other than those inherent to the NRB calculation, we take the gradient method PBLH uncertainty to be 5%, following the upper estimate of typical NRB uncertainty in Welton et al. (2000).

Fig. 3.

(a) An example potential temperature profile from the drone at midday on 6 March 2023, showing the planetary boundary layer height (PBLH) determined by the 1.5-theta method. (b) An example normalised relative backscatter (NRB) profile taken at the same time by the Mini-Micropulse (MP) LiDAR, showing the fitted ideal profile shape and associated ideal profile method-derived PBLH. (c) Gradient of the MP LiDAR NRB signal and associated gradient method PBLH. (d) Covariance transform of the Haar wavelet function and the NRB signal, showing the wavelet method PBLH. (e) Greyscale image showing the mean cloud-free MP LiDAR NRB in each daylight hour (local time) throughout the campaign. Overlayed are the campaign mean PBLH estimates from each MP LiDAR method. The means of all estimated drone PBLH values at 09:00, 12:00 and 15:00 hours are represented by yellow stars. Error bars are the standard deviations of the means.

The IP method aims to fit an idealised four-parameter equation for the shape of the NRB signal, following Steyn et al. (1999) and as outlined in Li et al. (2017) and Dang et al. (2019). The IP is given as follows:

The error function denoted by ε is defined in Dang et al. (2019). The fit parameter Z_m is the PBLH and s is the thickness of the transition region between the boundary layer and the free troposphere (the entrainment zone). The fit parameters B_m and B_u refer to the backscatter below and above the entrainment zone, respectively. The standard error on the fit for parameter Z_m is calculated and added in quadrature to the NRB uncertainty to produce the overall IP method PBLH uncertainty. The shape of B(z) is shown in magenta in Fig. 3b, where the IP function has been fitted to the measured NRB signal. The entrainment zone defined by s in Eqn 4 is indicated in Fig. 3b by the magenta-shaded region.

The wavelet method determines the PBLH by covariance transform between a wavelet profile and the LiDAR NRB profile. Here we use the Haar wavelet shape h (Gamage and Hagelberg 1993) and covariance transform w(s, b) (Brooks 2003):

where z is the profile altitude, s is the dilation of the wavelet shape and b is the translation of the wavelet. The lower and upper bounds of the LiDAR profile NRB(z) are given by z₀ and z₀ respectively. Although other wavelet shapes can be used, the Haar wavelet has been shown to be robust for PBLH detection in Baars et al. (2008), Dang et al. (2019) and Wang et al. (2021). The value of b when the minimum of w(a, b) is found, is taken as the PBLH. This relies on a suitable value of s, which physically represents the thickness of the entrainment zone, as for the IP method, and as indicated for an example NRB profile in Fig. 3b. On inspecting the NRB profiles, the entrainment zone was determined to typically be ~200 m thick, hence a value of s = 0.2 km is used throughout for wavelet method PBLH detection. We make a conservative estimate of the uncertainty of this approach by calculating the standard deviation of Haar wavelet PBLH values retrieved with s = 0.1, 0.2, 0.3 and 0.4 km. The standard error calculated from this is added in quadrature to the 5% NRB uncertainty to produce an overall wavelet method PBLH uncertainty. The shape of the covariance transform function applied to the Haar wavelet shape and the NRB profile at midday on 6 March 2023, is shown in Fig. 3c.

3.1. Results overview

Fig. 3 illustrates the key aim of this work; the detection of the PBLH in drone and MP LiDAR vertical profiles. The profiles and PBLH retrievals are from midday on 6 March 2023. The Θ and NRB profiles are typical of those throughout the campaign, with the PBLH evident by inspection as both an increase in Θ, and decrease in NRB, ~1-km altitude. In this example, the three MP LiDAR methods agree within 100 m of each other and the drone 1.5-theta PBLH. An increase in Θ with altitude above the PBLH (Fig. 3a) indicates a stable boundary layer. To assess the atmospheric stability across the campaign, we calculate the mean Θ value in the layers 100 and 400 m thick, above the calculated 1.5 theta PBLH. Only 2% of the profiles show a Θ decrease in the 100 m above the PBLH. None show a Θ decrease on average over the 400 m above the PBLH. These results highlight the consistent stability of the boundary layer during the campaign.

The campaign average diurnal NRB profile is shown in Fig. 3b, between 0 and 1.5 km altitude. This represents the mean NRB vertical profile in each daylight hour of the day. We show campaign average results rather than daily PBLH evolution examples because the frequency of broken clouds over One Tree Island meant there were no completely clear days. Apart from some noise in the early afternoon, above 1 km, a consistent diurnal cycle emerges with backscatter increasing in magnitude and altitude throughout the morning, before decreasing in magnitude and altitude through the afternoon. Despite this, none of the MP LiDAR PBLH methods show strong variations in their diurnal cycle. No consistent increase in PBLH through the middle of the day is observed, when solar heating is strongest, indicating a lack of convective boundary layer development. This indicates that the boundary layer is well mixed and stable throughout the day. Previous boundary layer measurements in the region, at Heron Island (~30 km north-west of One Tree Island) in 2012 using radiosondes, found evidence of greater diurnal boundary layer variability (MacKellar et al. 2013) – variations in boundary layer structure therein were attributed to convection associated with the reef. The lack of PBLH diurnal variability over the campaign implies that the diurnal MP LiDAR backscatter variability in Fig. 3b is unlikely to be explained by vertical transport from the surface. Instead, the daytime aerosol enhancement may be caused by horizontal transport of aerosols into the sampling domain, possibly including aerosols regenerated from cloud droplet evaporation (Bott 1999).

The mean PBLH using the IP method is consistently higher than other MP LiDAR methods. The gradient method consistently produces the lowest mean PBLH values, whereas the wavelet method is in between. The overall campaign mean (±s.d.) PBLH values are 910 ± 200 m for the IP method, 860 ± 290 m for the wavelet method and 810 ± 260 m for the gradient method. No evidence of convective boundary layer development is observed in the mean drone PBLH values in Fig. 3b either, with the maximum PBLH at 09:00 hours (830 ± 190 m), followed by slight decreases at 12:00 hours (790 ± 170 m) and 15:00 hours (780 ± 240 m). The PBLH values here are within the range of summertime marine boundary layer heights reported at Heron Island in MacKellar et al. (2013) and McGowan et al. (2022).

3.2. Comparison of drone and MP LiDAR PBLH

Fig. 4a shows a timeseries of PBLH retrievals from the drone (data from both iMet sensors and both ascent and descent flights are shown at each sampling time) the MP LiDAR (gradient method), overlayed on the overall MP LiDAR backscatter timeseries. Data gaps, due to cloud filtering of the MP LiDAR results, are evident as white segments in the NRB timeseries and missing MP LiDAR gradient PBLH retrievals. The backscatter-derived PBLH from the gradient method is plotted as a continuous timeseries (as long as data are available), which indicates the high-frequency variability of the PBLH retrieval, and can be up to 200 m over a 2-h period. The timeseries shows frequent overlap between the MP LiDAR gradient method and the drone within the respective uncertainties. Broad temporal features are captured by each method, such as lower average PBLH values for 10–17 March than in the rest of the comparison period. Data from a weather station mounted at 10-m altitude on One Tree Island (Lufft GmbH WS800–UMB Smart Weather Sensor, data not shown) show that periods of lower PBLH correspond to low surface wind speed (typically <20 km h⁻¹) from the northerly sector and high relative humidity (>75%). Periods of higher PBLH tend to occur with trade winds from the southerly and easterly sectors, often >20 km h⁻¹ and relative humidity of 60–75%. These broad temporal characteristics of the timeseries agree between MP LiDAR PBLH methods.

Fig. 4.

(a) Timeseries of planetary boundary layer height (PBLH) estimates from the drone and Mini-Micropulse (MP) LiDAR (gradient method) throughout the campaign, overlayed on the MP LiDAR normalised relative backscatter (NRB) timeseries. (b), (c) and (d) Scatterplot comparisons of MP LiDAR v. drone PBLH estimates for gradient, IP and wavelet methods, respectively. Dashed lines indicate 1:1, solid lines are the linear fit to all data; all other colours in this figure correspond to measurement times 09:00, 12:00 and 15:00 hours, consistent with Fig. 3. Error bars for drone measurements are the standard deviation of the four PBLH retrievals in each timestep; error bars for MP LiDAR results are calculated as described in the Methods section.

Here, we compare two distinctly different quantities as an indicator of PBLH: a thermodynamically derived measure and an aerosol backscatter-derived measure. These quantities may not necessarily be the same. Under decoupled boundary layer conditions, or in the presence of elevated aerosol pollution layers, for example, the interpretation of backscatter in relation to the thermodynamic PBLH is more complicated (Bohlmann et al. 2018; Alexander and Protat 2019). Nevertheless, Alexander and Protat (2019) and Bohlmann et al. (2018) also show that under clean marine conditions and when the boundary layer is well-mixed, the thermodynamic- and backscatter-derived PBLH will likely be similar. We use a statistical comparison of MP LiDAR PBLH methods compared to the drone 1.5-theta results to demonstrate the consistency of the similarity between the thermodynamic and backscatter-derived PBLH in the One Tree Island dataset. We used orthogonal distance regression (ODR) analysis, considering the relative uncertainties associated with the individual MP LiDAR methods and with averaging the four drone results at each sampling time. Fig. 4b–d show scatterplots with linear ODR fits for all data (black lines) and for each sampling hour (coloured lines), for the gradient, IP and wavelet methods respectively. Comparison data existed for 37 days, encompassing 111 discrete sampling times. Following the cloud filtering, 57 comparison points remain. The cloud filtering does not disproportionately remove data from any one sampling time, with each having a similar number of total comparison points: 18 for 09:00 hour, 21 for 12:00 hours and 18 for 15:00 hours. The uncertainties vary between MP LiDAR methods. The 5% uncertainty ascribed to the gradient method, ~40 m on average, is the lowest uncertainty of the three methods and just larger than the 30-m vertical resolution of the MP LiDAR. The uncertainty associated with the IP method is also low at ~6% (on average ~55 m). The uncertainty associated with the Haar wavelet method is much larger at ~18% (on average ~ 150 m).

The scatterplots show moderate scatter in each comparison. Considering the entire dataset, Pearson’s R coefficients are 0.68 for the IP method, 0.62 for the gradient method and 0.59 for the wavelet method. Considering all sampling times, R values range from 0.36 (wavelet method, 12:00 hours) to 0.92 (gradient method, 09:00 hours). The R values are consistently highest for the 09:00 hours comparison. No one MP LiDAR method has consistently higher correlation; however, the gradient and IP methods arguably show a tighter data density around the 1:1 line than the wavelet method. The gradient method’s overall R value is likely detrimentally influenced by five examples of the MP LiDAR result being much higher.

The ODR regression slopes (Fig. 4) are all greater than 1. This indicates that all MP LiDAR methods tend to locate the boundary layer lower than the drone at low altitude and higher than the drone at high altitude. For the whole data comparison, the slope is 1.9 for the gradient method, 1.5 for the IP method and 1.4 for the wavelet method. For the gradient and IP methods, the 09:00 hours comparison slope is closest to 1, with 12:00 and 15:00 hours comparison slopes being variable. The wavelet method shows the most consistent slope between sampling times, with range 1.1–1.5. Owing to the regression slopes all being >1, y-intercept values (not shown) for the ODR linear fits range are all negative. Here, we reiterate that the IP and gradient method data, apart from a few discrete poor performing examples, appear more closely clustered to the 1:1 line (Fig. 4b, c) than for the wavelet method (Fig. 4d), despite the wavelet method slopes typically being closer to 1. The higher slope values for IP and gradient methods are likely influenced these methods’ lower uncertainty compared to the wavelet method, boosting the relative uncertainty of the x-ordinate and making the slope steeper.

For the whole dataset, the NMB is −7.4% for the wavelet method, −6.5% for the IP method and only −1.1% for the gradient method. The NMB is negative for all methods and nearly all sampling times, indicating that the drone 1.5-theta PBLH value is on average lower than the MP LiDAR estimate. The generally good agreement between the drone and all MP LiDAR methods is indicated by the fact that the absolute value of the NMB for all comparisons except one (12:00 hours, wavelet method) is <10%, which means differences are typically <100 m.