Response to Reviewer 3 Comments : Blowing snow detection from ground-based ceilometers : application to East Antarctica

Blowing snow detection from ground-based ceilometers: application to East Antarctica Alexandra Gossart1, Niels Souverijns1, Irina V. Gorodetskaya2,1, Stef Lhermitte3,1, Jan T.M. Lenaerts4,1,5, Jan H. Schween6, Alexander Mangold7, Quentin Laffineur7, and Nicole P.M. van Lipzig1 1Department of Earth and Environmental Sciences, KU Leuven, Leuven, Belgium 2Centre for Environmental and Marine Sciences, Department of Physics, University of Aveiro, Aveiro, Portugal 3Department of Geosciences and Remote Sensing, Delft University of Technology, Delft, the Netherlands 4Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, The Netherlands 5Departement of Atmospheric and Oceanic Sciences, University of Colorado, Boulder CO, USA 6Institute of Geophysics and Meteorology, Koeln University, Koeln, Germany 7Royal Meteorological Institute of Belgium, Brussels, Belgium

Abstract. Blowing snow impacts Antarctic ice sheet surface mass balance by snow redistribution and sublimation. However, numerical models poorly represent blowing snow processes, while direct observations are limited in space and time. Satellite retrieval of blowing snow is hindered by clouds and only the strongest events are considered. Here, we develop a blowing snow detection (BSD) algorithm for ground-based remote-sensing ceilometers in polar regions and apply it to ceilometers at Neumayer III and Princess Elisabeth (PE) stations, East Antarctica. The algorithm is able to detect (heavy) blowing snow layers reaching 30 m height. Results show that 78 % of the detected events are in agreement with visual observations at Neumayer III station. The BSD algorithm detects heavy blowing snow 36 % of the time at Neumayer (2011Neumayer ( -2015 and 13 % at PE station (2010)(2011)(2012)(2013)(2014)(2015)(2016). Blowing snow occurrence peaks during the austral winter and shows around 5 % interannual variability. The BSD algorithm is capable of detecting blowing snow both lifted from the ground and occurring during precipitation, which is an added value since results indicate that 92 % of the blowing snow is during synoptic events, often combined with precipitation. Analysis of atmospheric meteorological variables shows that blowing snow occurrence strongly depends on fresh snow availability in addition to wind speed. This finding challenges the commonly used parametrizations, where the threshold for snow particles to be lifted is a function of wind speed only. Blowing snow occurs predominantly during storms and overcast conditions, shortly after precipitation events, and can reach up to 1300 m a.g.l. in the case of heavy mixed events (precipitation and blowing snow together). These results suggest that synoptic conditions play an important role in generating blowing snow events and that fresh snow availability should be considered in determining the blowing snow onset.

Introduction
Understanding the Antarctic ice sheet (AIS) response to atmospheric and oceanic forcing is crucial given its large potential impact on sea level rise (Rignot and Thomas, 2002;Rignot and Jacobs, 2002;Rignot et al., 2011;Shepherd et al., 2012). AIS mass balance is governed by the difference between surface mass balance (SMB) and solid ice discharging into the ocean. Solid precipitation is the only source term for the SMB. Maltwater runoff and surface sublimation are terms removing mass at the surface of the AIS. In addition, the sublimation of the suspended snow particles adds to this removal. The fourth process is the wind-induced erosion or redeposition of transported snow particles from one location to another (Takahashi et al., 1988). Snow particles can be dislodged from the snow surface, picked up by the wind and lifted from the ground into the near-surface atmospheric layer. Drifting snow events are shallower than blowing snow events. Drifting snow typically stays below 2 m height whereas blowing snow can reach heights of several hundreds of meters. The transport involves a mix of suspension and saltation transport modes (Leonard et al., 2011), with a dominance of saltating particles (Bagnold, 1974) in the case of drifting snow, and suspended particles in blowing snow layers (Mellor, 1965). Despite its importance, the role of blowing snow on AIS SMB is currently poorly quantified. If we consider the ice sheet in its whole, the contribution of blowing snow is rather small: around 0-6 % (Loewe, 1970;Déry and Yau, 2002;Lenaerts and van den Broeke, 2012). However, blowing snow is crucial for the regional SMB (Gallée et al. , 2001;Déry and Yau, 2002;Lenaerts and van den Broeke, 2012;Groot Zwaaftink et al., 2013) through the displacement and relocation of the snow particles (Déry and Tremblay, 2004). This phenomenon occurs approximatively on 70 % of the Antarctic continent during winter (Palm et al., 2011). In addition, drifting and blowing snow sublimation contributes substantially to SMB (Kodama et al., 1985;Takahashi et al., 1992;Thiery et al., 2012;Dai and Huang, 2014). This process can even be more effective to remove mass than surface sublimation (van den Broeke et al., 2004). The combination of blowing snow sublimation and transport is estimated to remove from 50 to 80 % (van den Broeke, 1997;Frezzotti et al., 2004;van den Broeke et al., 2008;Scarchili et al., 2010) of the accumulated snow on coastal areas. Moreover, removal of the snow by the wind can locally lead to the formation of blue ice areas (Takahashi et al., 1988;Bintanja et al., 1995), which have a lower albedo and therefore enhance surface melt, and could affect ice shelf stability and collapse (Lenaerts et al., 2017). Blowing snow also plays a role in determining snow surface characteristics (Déry and Yau, 2002), affecting surface energy balance (Yamanouchi and Kawaguchi, 1985;Mahesh et al., 2003;Lesins et al., 2009).
Many studies have focused on a minimum wind speed as a threshold to dislodge snow particles, depending on the snow surface properties (Budd et al., 1966). Schmidt (1980Schmidt ( , 1982 explained that cohesion between snow particles requires higher wind speeds or a higher impacting force of particles on the snowpack. In addition, the presence of liquid water in the snow and enhanced snow metamorphism with the higher atmospheric temperatures in the summer induce varying wind thresholds throughout the year (Bromwich, 1988;Li and Pomeroy, 1997).
Currently, simulations of the AIS SMB are highly uncertain since both precipitation and blowing snow processes are poorly constrained and probably lead to inconsistencies between the atmospheric modeled precipitations and the measured snow accumulation value (Frezzotti et al., 2004;van de Berg et al., 2005;Scarchili et al., 2010;Groot Zwaaftink et al., 2013;Gorodetskaya et al., 2015). In addition, strong blowing snow also hampers ground detection from satellites, and biases can be induced in efforts to study the Antarctic surface elevation due to the presence and radiative properties of blowing snow (Mahesh et al., 2002(Mahesh et al., , 2003. Efforts have been made to retrieve blowing snow from satellite data, but, while it offers a large area coverage, the detection is limited to clear-sky conditions and blowing snow layers thicker than 30 m (Palm et al., 2011) and make use of a wind threshold criterion. Moreover, ground validation remains essential to evaluate satellites measurements. A number of measurement campaigns have been organized in various regions of the AIS, using different types of devices: nets, mechanical traps and rocket traps, photoelectric and acoustic sensors, or piezoelectric devices (Leonard et al., 2011;Barral et al., 2014;Amory et al., 2015). However, custom-engineered sensors are rather expensive and scarce (Leonard et al., 2011), and both the remoteness of the continent and the harshness of the climate are limitations to widespread use of these devices.
In this study we propose a new method to detect blowing snow by the use of ground-based remote-sensing ceilometers. Ceilometers are robust cloud-base height detection devices. Frequently used in airports and designed to report visibility for pilots, the backscatter signal of these ground-based low-power lidars contain further information, widely used for scientific purposes regarding boundary-layer investigation (Marcowicz et al., 1997;Eresmaa et al., 2006;Heese et al., 2010;Thomas, 2012). They have been used to detect the vertical extent of aerosol layers below 5 km, mixing height layers  and the detection of the early stage of radiation fog . Several algorithms have been developed to detect cloud-base height in specific areas: at the polar regions using the polar threshold algorithm (Van Tricht et al., 2014) or at temperate latitudes with the temporal height tracking algorithm (Martucci et al., 2010) and the standard Vaisala algorithm (Flynn, 2004). Ceilometer networks are also developed as a potential to cover larger regions (Illingworth et al., 2015). Over the Antarctic continent, the environmental conditions require that research stations are usually equipped with robust instruments that are able to withstand low temperatures and high winds. Ceilometers can be operated autonomously and continuously in environmental conditions between −40 and +60 • C, with up to 100 % relative humidity and 50 m · s −1 wind speeds (Vaisala User's guide, 2006). Compared to lidars, ceilometers have numerous advantages, such as eyesafe operation, low first range gate and relatively low price, making it one the most abundant cloud detection device on the ice sheets (Van Tricht et al., 2014;Wiegner et al., 2014).
The goal of this paper is to present a new methodology for blowing snow detection (BSD) using the ceilometer attenuated backscatter profile and to estimate the frequency of blowing snow at Neumayer III and Princess Elisabeth (PE) stations. Subsequently, we apply the BSD algorithm to infer various blowing snow statistics (occurrence, depth) and investigate meteorological conditions during blowing snow  events. We conclude by examining the applicability of the BSD algorithm to other Antarctic sites.

Ceilometers
Ceilometers consist of a single-wavelength, eye-safe active laser transmitter that emits pulses in the vertical direction and an avalanche diode receiver that collects the pulse signal. The laser pulse backscattered by molecules, aerosols, precipitation and cloud particles present in the atmosphere at height z is detected by the ceilometer receiver. Typically, the backscatter intensity depends on the concentration or size of particles in the air, but the ceilometer receiver also detects noise induced by the device's electronics and the background light. The lidar equation enables us to get the return signal strength from the emitted laser pulse (Münkel et al., 2006): (1) The attenuated backscatter profile at range a, β att (sr −1 ·m −1 ), is a product of the true backscatter coefficient β at height z, taking into account the two-way attenuation of the lidar due to the transmittance of the atmosphere (τ 2 ). A height normalization is applied to the retrieved signal to remove the excessive decrease in backscatter intensity. The detected signal is reported at the center of the 10 m range gate (i.e., for a signal measured between 50 and 60 m, the value of the range gate will be attributed to a height of 55 m (range bin 6)). The ceilometer measures continuously and the standard output, β att , is displayed in a time-height cross section, with a 10 m vertical resolution and 15 s temporal resolution (Fig. 1). The cloud-base height is the standard ceilometer output determined from the backscattered signal. In addition, the quantitative information that can be derived from the ceilometer measurements is the instantaneous magnitude of the signal received by the diode that provides information on the backscattering properties of the atmosphere at determined heights (Wiegner et al., 2014;Madonna et al., 2015). Other properties such as optical depth, size and density of the scatterers would require to know the lidar ratio, and a reliable estimate of lidar ratio is complicated (Wiegner et al., 2014). In addition, this is only possible if the ceilometer is calibrated, which is very challenging since the signal-to-noise ratio (SNR) has to be large enough in the troposphere (Wiegner et al., 2014), and this is not done in the present study. This implies that quantification of blowing snow displacement, and the determination of blowing snow properties such as particles density, shape or number cannot be derived from the ceilometer attenuated backscatter signal at Neumayer III and PE stations.

The cloud-precipitation observatory at PE station
The PE station is located on the Utsteinen ridge in Dronning Maud Land (DML), East Antarctica ( Fig. 2 and Table 1). A cloud and precipitation observatory was set up on the roof of the station (approx. 10 m above the ridge) during the summer season of 2009-2010 and is still operational under the Hydrant and Aerocloud projects (www.aerocloud.be). The observatory contains an automatic weather station (AWS) and a set of ground-based remote-sensing instruments. The observatory was designed to be operated year-round, including the winter period when PE is unmanned. The station and the set of instruments are controlled remotely via a satellite connection.
The Vaisala CL-31 ceilometer was installed on the roof of the station in December 2009 and is operational at present (see Table 2 and Fig. 3). Several outages of the energy provision system limit the data mainly to Antarctic summer season (December to March is best represented). Only during the year 2015 continuous measurements were obtained.
The Metek vertical-profiling precipitation radar (microrain radar, MRR), set up in 2010, enables us to retrieve snowfall rates using the return from the vertical-profiling Doppler radar operating at a frequency of 24 GHz. The raw Doppler spectra is post-processed following Maahn and Kollias (2012) to calculate radar reflectivity profiles which are then linked to snowfall rates using the newly developed Ze-Sr relation for PE by Souverijns et al. (2017). A full description of MRRs can be found in Klugmann et al. (1996) and the radar set up at PE is described in Gorodetskaya et al. (2015). The monitoring of the instruments set up on the roof of the station is done via a webcam. For a specifications of the instruments, see also Gorodetskaya et al. (2013) and Gorodetskaya et al. (2015).
The AWS was set up 300 m from the station for recording meteorological parameters, broadband radiative fluxes and snow height changes (Gorodetskaya et al., 2013). It is designed to work continuously in remote locations, enabling studies of mass balance and radiative fluxes. The AWS was established in February 2009 and replaced by a new station in December 2015, both designed by the Institute for Marine and Atmospheric Research, University of Utrecht (Utrecht, the Netherlands). The station provides hourly mean data of near-ground and air temperature, air pressure, wind speed and direction, relative humidity and radiative fluxes (downwards and upwards short-and longwave radiation) (for details on sensors, see Table S1 in the Supplement). Postprocessing of the data includes a treatment for relative humidities as described by Anderson (1994) for humidity with respect to ice and a correction for the relative humidities above 100 % following van den Broeke et al. (2004). The temperature gradient is computed as the difference between the 2 m and surface temperatures over the distance between the sensors (Gorodetskaya et al., 2013).

Neumayer III research station
Neumayer III research station is located on the Ekström ice shelf, in the northeastern Weddell Sea ( Fig. 2 and Table 1). Researchers are present year-round at the station and it is equipped with various instruments. Measurements include upper air soundings, ozone soundings, radiation measurements and weather observations. Weather measurements are carried out since 1981 at Neumayer, and the station is the weather forecasting center for DML. The synoptic observations at Neumayer III include 2 (König-Langlo, 2011, 2013 and 10 m air temperature, air pressure, wind vector at 2 and 10 m height, 2 m dew point temperature, cloud presence, information on the type and height of clouds, horizontal visibility, and past and present weather including snowdrift and whiteout (for a description of the sensors, see Table S2). The measurements are carried out daily every 3 h except for 03:00 and 06:00 UTC. In this paper we use the visual observations of blowing snow, classified into nine categories (S8 code) according to the Word Meteorological Organization (WMO) coding system (see Table S3). The visual observations regarding blowing snow are performed as follows (detection procedure from Gert König-Langlo, personal communication, 2016): "if the wind exceeds 5 m · s −1 , the observer goes out about 100 m windward from the research station and observes the snow surface. No target is used to detect blowing snow against, and during winter (in complete darkness) a small flashlight is used. The distinction between blowing and drifting snow is made according to the Table 1. Meteorological conditions at Princess Elisabeth and Neumayer III stations. For extended climatology, see Gorodetskaya et al. (2013) for PE station and König-Langlo and Loose (2007)    height of the blowing snow layer in relation to the eye level: drifting snow below the eye level and blowing snow above. Further, if the blowing snow layer is not too dense, one can distinguish blowing snow with or without precipitation by an additional observation from the roof of the station." The set of instruments present at Neumayer III station includes a Vaisala ceilometer CL-51, set up on the roof of the station and operating continuously since 15 January 2011 (see Table 2 and Fig. 3). The blowing snow record (König-Langlo, 2014, 2016 at Neumayer station is analyzed together with the atmospheric measurements available from the synoptic observations. An overview of the climatic conditions is given in Table 1. The data are freely available interactively from https://www.pangaea.de/ (Konig-Langlo, 2015, 2016.
3 Data treatment and BSD algorithm

Preprocessing
We average every 15 s β att profile over 1 h using a running mean to create mean attenuated backscatter profiles at every time step and avoid the variability due to turbulence and hardware noise. Figure 4 shows the resulting β att for 24 April 2016 at 09:30 UTC, based on the average of 240 profiles (120 preceding and 120 following 09:30 UTC). An additional reason for the integration of the signal over longer time periods is that it improves the SNR. No additional SNR correction is performed on the raw data, as we found that a temporal SNR higher than 0.3 would remove parts of the blowing snow signal (Gorodetskaya et al., 2015).
There are two sources of noise and artifacts affecting the ceilometer backscatter signal: the hardware of the Vaisala ceilometers and the internal processing of the data (Kotthaus et al., 2016). Firstly, a heater is incorporated in the device to stabilize the laser temperature in cold environments. This heater is placed close to the laser transmitter and the periodic turning on (when a minimum temperature is reached by the instrument) and off (when the laser temperature is high enough) of the heater introduces a small periodic variation in the stability of the emitted signal (and therefore of the detected signal). This effect is stronger in the first range bins, closest to the device, and the hourly running mean enables us to smooth out most of this signal variation. Secondly, the internal processing of the signal includes a built-in correction for the partial overlap of the laser in the first range bins. This overlap is due to the coaxial configuration of the laser: the same lens is used for the emitted and the received signals, made possible by the use of mirrors (Spirnhirne, 1993;Vande Hey, 2015). The total overlap is only reached at the seventh range bin (65 m) for the CL-31 (Vande Hey, 2015; Kotthaus et al., 2016). However, the partial overlap in the near-ground range bins does not imply that the minimum detection range is at 65 m only; in the case the signal returned by the close-range scatterers is large enough (which is mainly the case during blowing snow), it will be recorded even before the overlap onset (Vande Hey, 2015). Lastly, the CL-31 backscatter profile is constrained in the lowest bins by a builtin function to correct for unrealistically high values resulting from window obstruction. However, this correction likely introduces artifacts in the signal in the first range bins. As a result of the periodic switching on and off of the heater and the low overlap in the first range gate, the reported value of the combined β att signal in the lowermost range bin is systematically and unrealistically higher than the signal in the next bins (Vaisala, personal communication, 2016). We therefore exclude the signal reported in the lowermost range bin in our analysis and start investigating the profile from the second range bin, 15 m above the CL-31 and CL-51 ceilometers onwards.
Moreover, artifacts have been observed in the ceilometer profiles at both stations (also visible in Fig. 6). There is a discontinuity in all profiles between the fourth and fifth range bins (35 and 45 m). This discontinuity is also visible in profiles where the instrument is hooded, which are supposedly mimicking full atmospheric attenuation and thus recording the background noise produced by the hardware and electronics. Many authors have reported artifacts in the lowest range bins (below 70 m height) that are usually excluded during processing for boundary-layer investigation (Wiegner et al., 2014). This local minimum is also reported by Sokol (2014) at the fifth range bin (45 m) during the whole duration of his campaign, as well as by Martucci et al. (2010) and Tsaknakis et al. (2011). Kotthaus et al. (2016) stated that these are likely due to the correction applied by Vaisala to prevent unrealistic values in the lower bins, related to the obstruction of the window and the internal noise. In the case of Vaisala instruments, the output is already corrected with a correction function, unknown to the user, which cannot be modified (Wiegner et al., 2014). This has to be kept in mind when using the profile information to detect blowing snow.

The BSD algorithm
Studies investigating the boundary-layer properties based on ceilometer β att make use of both properties of the signal (shape and intensity) to evaluate the presence and extent of a particular layer, e.g., in order to determine the height of the mixing layer (Wiegner et al., 2014). For such an analysis, five methods have been developed (Emeis et al., 2008), including a threshold method and a gradient method . In the first case, the mixing height is attained when the intensity of the signal drops below a fixed threshold value (Münkel and Rasanen, 2004). The second method  considers the minimum of the first or second derivative of the backscattering profile as top of the mixing layer (Sicard et al., 2004). To detect the occurrence of blowing snow, Palm et al. (2011) used a combination of both types of methods on the CALIOP (satellite-borne) backscatter. First, the intensity of the backscatter in the bin closest to the detected ground return must exceed a certain threshold. Second, the decrease of the profile of the signal with height indicates the presence of blowing snow: the concentration of particles close to the ground is much higher than in the overlying layers (Takeuchi, 1980;Schmidt, 1982;Palm et al., 2011). This is associated with a sharp vertical gradient, where the β att profile decreases strongly in the very first range bins. In addition, a wind speed threshold is applied (3 m · s −1 at 10 m).
The approach used for the BSD algorithm is similar, but there is no wind speed criterion in our analysis. In addition, the ceilometer is ground based, allowing the detection of blowing snow during overcast conditions. The algorithm method is displayed in Fig. 5. To detect blowing snow, the intensity of the backscatter signal at the lowest usable bin must exceed a certain threshold (Sect. 3.3), and the intensity of the signal must decrease in the next range bins indicating a particles density greater in the lower levels than at layers directly above. As previously highlighted, clean air molecules cannot be distinguished because the signal associated with it is smaller than the noise generated by the hardware (Wiegner et al., 2014;Kotthaus et al., 2016) and by the background light (Vande Hey, 2015), polluting the signal in the lowest bins. To distinguish the presence of scatterers (aerosols, blowing snow particles, cloud particles, etc.) present in the atmosphere from these artifacts, we need to investigate the signal intensity representative for cloudless conditions, i.e., the average β att of the lowest usable range bin received by the ceilometer during scatter-free conditions. Clear-sky days are manually selected for the whole period using the daily quicklooks (Fig. 1) and are days when the quicklook background is uniform and without precipitation or clouds and when the time series of the signal in the lowest usable range bin is stable around a low value (corresponding to hardware and background noises) to avoid low-level disrupting signal. Next, we compute the 99th percentile of all clear-sky β att signal in this range bin as threshold value (for calculation, see Sect. 3.3). As such, it is representative of the presence of scatterers ex-ceeding the value for clear sky. Since the noise is instrumentdependent, individual preprocessing and thresholds have to be defined for each instrument the BSD algorithm is applied to.
After comparing the backscatter signal to the clear-sky threshold, the BSD algorithm investigates the shape of the β att profile. A regular clear-sky ceilometer profile (signal intensity vs. height) does not show intense vertical variations; in the infrared, the transmission term is close to 1 and decreases only slightly with height. This implies that any important variation in the β att signal can be attributed to the particles backscatter. The profiles for blowing snow in Fig. 6 show a typical sharp decrease until bin 8-10 (75-95 m height), above which the signal keeps decreasing steadily (blue line): this is the signature of clear-sky blowing snow. The red profile, in contrast, shows a reincrease in intensity around the 15th bin (145 m height), overlying the blowing snow signal: this indicates the presence of scatterers interpreted as precipitation or clouds (denoted by the arrow on the graph). If there is no blowing snow while precipitation and/or clouds are present, the profile does not decrease prior to the increase at higher levels (black line in Fig. 6). The algorithm therefore investigates the shape of the profile in order to detect blowing snow. A condition is set that a blowing snow profile implies that the mean of the overlying bins 3 to 7 (25 to 65 m) must be lower than the signal in the lowest usable range bin (15 m). In this way, the discontinuity, as described in Sect. 3.1. (visible in Figs. 1 and 6 between 35 and 45 m) does not affect our retrievals. In order to detect blowing snow occurring during clouds or precipitation, the profile shape is analyzed to identify a second increases in the signal intensity above the seventh bin (65 m height). A clear differentiation between clouds or precipitation cannot be made on the basis of the ceilometer alone, but the presence of clouds and/or precipitation can be identified. This analysis is carried out for both blowing snow and the absence of blowing snow measurements. The information retrieved from the MRR (hourly precipitation rates) is collocated to ceilometer BSD to determine the time (in hours) since the last precipitation event at PE station.
Inherent to this profile-based method, the detection of blowing snow during precipitating events is limited to cases when the blowing snow signal is preserved close to the ground. In the case of precipitation associated with storms, there is always blowing snow due to the high wind displacing the snow, and no distinction between precipitation and blowing snow is possible, as the ceilometer signal is entirely attenuated near the surface (Gorodetskaya et al., 2015). It is thus not possible to get signal in the overlying bins, and the profile of the backscatter intensity might not decrease upwards. Such intense precipitating events mixed with blowing snow are identified as having a second bin signal higher than 1000 ·10 −5 · km −1 · sr −1 (threshold adapted from Gorodetskaya et al., 2015). In those cases, the events are classified as a intense mixed event, and the profile analysis is omitted by the algorithm.
In addition to the detection of blowing snow, the BSD algorithm quantifies the height of the layer. This is done as follows: if the profile decreases steadily (indication of absence of precipitation), the range gate at which the intensity of β att drops under the clear-sky threshold value is the top of the layer. Anything above this height is considered clear sky. If there is precipitation or a cloud during the blowing snow event, the shape of the backscatter profile does not decrease monotonously but shows an increase in higher levels. In that case, the range gate at which the profile increases again is the top of the blowing snow layer and the base of the cloud and/or precipitation (the green line around the seventh bin in Fig. 6, for the black and the red profiles). Layer height definition is illustrated in Fig. S3 (Supplement).

Application of the BSD algorithm to different stations
The BSD algorithm is designed to detect blowing snow events reaching heights of minimum 15 m and is based on the Vaisala CL-31 located at PE station, for the period 2010-2016. It is applicable to other ceilometers: we applied the BSD algorithm to backscatter data from the Vaisala CL-51 ceilometer at Neumayer station for the years 2011-2015.
The time (15 s) and height resolution (10 m) are the same for both instruments. We can therefore apply the BSD algorithm in the same fashion to both datasets with the only difference being the β att near-surface threshold (first step in the BSD algorithm used to identify the presence of blowing snow). We obtain a threshold of 21 ·10 −5 · km −1 · sr −1 for the CL-31 ceilometer at PE, based on 127 clear-sky days out of a total of 1064 days. The threshold at Neumayer is of 32.5 ·10 −5 · km −1 · sr −1 , based on 125 clear-sky days out of 1444 days.

Frequency of blowing snow
In order to investigate the type of blowing snow detected by the BSD algorithm, we compare it to visual observations at Neumayer, carried out routinely at 09:00, 12:00, 15:00, 18:00, 21:00, and 24:00 UTC. All ceilometer measurements are considered over 1 h, corresponding to the time at which visual observations are carried out. We identify a blowing snow event when blowing snow is present in at least 80 profiles (20 min). The WMO visual observations are categorized in six classes of blowing and/or drifting snow events, ranging in intensity and whether there is precipitation or not (Table S3). Before we start the comparison, it should be noted that visual observations are difficult to perform, and the error associated with it is not quantified. Therefore, in this part we refer to the number of measurements that match or mis-match between the BSD algorithm and visual observations rather than using the visual observations as "ground truth". The total number of measurements, N, is the total number of visual observations performed during which the ceilometer is also measuring, independently of whether there is blowing snow or not (N = 10 854; six measurements per day for the 2011-2015 period). The match ratio (Eq. 2) is the total agreement between visual and BSD algorithm detections over N; with N BS both when both the ceilometer and the observer detect blowing snow, and N BS none when neither the ceilometer nor the observer detects blowing snow. Mismatches occur when only one of the methods detects blowing snow, while the other does not: N BS ceilo if blowing snow is only reported by the BSD algorithm, but not reported by the visual observations (commission error), and N BS vis when only the visual observer records blowing snow, referred to as omission error.
The results (Tables 3 and S4) show a very good match for the detection of blowing snow events, and the optimum, 78 %, is reached for events classified as all blowing snow with or without precipitation: (2) The lowest match (70 %) is found when all blowing and drifting snow is taken into account: the number of visually detected events strongly increases since more categories are included, whereas the number of detections by the BSD algorithm is fixed. In 21 % of the time, the visual observer reports something (blowing or drifting snow) that is not detected by the BSD algorithm. This is related to the fact that the ceilometer points upwards and is elevated at a height of 17 m above the surface, which prevents it from detecting shallow layers of drifting snow. A fraction of 84 % visually observed heavy blowing snow events is detected by the BSD algorithm: In this case, we consider visual observations reported as "heavy blowing snow" only. For 95 % of the N BS ceilo events not reported as heavy blowing snow by the observer, intensities of the backscatter signal are below 1000 ·10 −5 · km −1 · sr −1 ; it is therefore likely that those events are classified as "slight" or "moderate" by the visual observer instead of being considered heavy. For the N BS vis , 54 % do not attain the threshold indicating the presence of scatterers and in 46 % of the cases the ceilometer attenuated backscatter profile does not decrease with height.
We also compare the skills of the BSD algorithm to different evaluation metrics (Allouche et al., 2006) (the equations for each of the metrics are presented in the Supplement). The accuracy, highest for the category collecting all blowing snow events, is the proportion of correctly detected events. To take into account omission errors, sensitivity is used and the best score is attained by both heavy blowing snow categories (with and without precipitation). Specificity reflects commission errors, and the categories encompassing most events (blowing and drifting snow) perform best. Finally, since the N BS none is larger than the other categories, the matches are likely biased, and we therefore use the true skill statistics (TSS; (Allouche et al., 2006)), which is a method enabling us to measure the overall accuracy and correct for the accuracy expecting to occur by chance, which also accounts for commission and omission errors, while being independent of prevalence in the data. TSS range from −1 to +1, where values under zero indicate no better performance than random, and the closest the result to 1 the better the agreement. This metric clearly indicates that heavy blowing snow is the best detected category of events.
To investigate to which extent the BSD algorithm is limited by precipitation, we compare matches and mismatches for each of the categories with and without precipitation. The value for N BS both doubles and even nearly triples when including events occurring during precipitation while N BS ceilo decreases by nearly 50 % and N BS vis increases by the same amount. This indicates that the BSD algorithm is not impeded by the presence of precipitation: the commission errors in the non-precipitation category were mostly blowing snow events that are encompassed when taking into account events occurring together with precipitation.
Moreover, we gather that the ceilometer algorithm is not limited to heavy blowing snow but also detects a number of events referenced under "moderate", or "slight blowing snow events", and even occasionally "drifting snow". This is revealed by the fact that N BS ceilo reduces as we consider less intense and shallower types of events (Table S4).
The frequency is calculated here by reporting the sum of all hours during which blowing snow occurs at Neumayer Table 3. Detection numbers and scores of the different categories of observations. The first four columns give the numbers for all four categories: N BS both stands for blowing snow detected by both the algorithm and the visual observations; N BS none indicates both methods agree that there is no blowing snow; N BS ceilo and N BS vis represent detections by the algorithm and the observer only, respectively (the corresponding percentages are presented in Table S4). The four last columns give the scores; TSS is true skill statistics. B stands for blowing and D for drifting snow; prec is precipitation. The total number of measurements is 10 584.  (2007), who report 20 % of drifting and 40 % drifting and blowing snow for the 1981-2006 period. However, there is an interannual variability that reaches ±5 %, also observed by Lenaerts et al. (2010). The pattern visible in Fig. 7 is common for blowing snow over Antarctica: a seasonal cycle peaking during the Antarctic winter (March-November) and displaying lower values for the rest of the year (Mahesh et al., 2003;Lenaerts et al., 2010;Scarchili et al., 2010;Palm et al., 2011;Amory et al., 2017). The overall blowing snow frequency is computed at PE for the 2010-2016 period and reaches 13 %. This lower blowing snow frequency at PE can be explained by the location of the station: the station is shielded from the katabatic winds by the Utsteinen mountain range, making it a quieter zone between the flows diverged to the sides of the station (Parish and Bromwich, 2007), while Neumayer III station is located on the ice shelf and experiences higher wind speeds (see Table 1) and is more exposed to storms. In addition, the limited availability of Antarctic winter data (due to power failures at the station) leads to an underestimation of the blowing snow frequency as mostly extended summer period was used, and only one winter is taken into account. The frequencies measured by the BSD algorithm are larger than those retrieved by satellite method: Palm et al. (2011) give a range of 0-10 % blowing snow for both locations. This can be related to the number of blowing snow events occurring together with clouds and precipitation that were missed by the satellite and to the different spatial and temporal dimensions of the different methods. Of all blowing snow detected events, 67 % is mixed with intense events at Neumayer III and 43 % at PE station. Cloudless blowing snow is very rare at Neumayer III station (8 % of the events), while it reaches 30 % at PE station.

Time since last precipitation and blowing snow occurrence
We investigate the time lag between the last precipitation event and the onset of blowing snow events at PE station. The majority of blowing snow occurs during or within a day after a precipitation event (nearly 60 and over 80 % of the blowing snow occurrences, respectively). There is a clear drop for larger time lags (Fig. 8a). This is, however, not so obvious if we normalize the distribution of blowing snow events, taking into account the total number of ceilometer measurements within each time lag after precipitation (Fig. 8b). A possible explanation is that there are less measurements as we go in time, as very long periods without any precipitation (300 h, more than 10 days) are rather rare at PE station, and that blowing snow occurred during those measurements. This can also be linked to the fact that the blowing snow particles detected by the BSD algorithm might originate from another location where there is precipitation, while snowfall is not reported by the precipitation radar at the station itself.

Blowing snow and meteorological regimes
The near-surface atmosphere changes associated with blowing snow events are investigated for both stations, and detailed means and standard deviation are displayed in Tables S5, S6 and S7. We investigate how blowing snow relates to weather regimes, derived from the hierarchical cluster analysis using PE AWS data following Gorodetskaya et al. (2013) which defines the weather regimes at PE station: "cold katabatic", "warm synoptic" and "transitional synoptic". The cold katabatic regime is characterized by slower wind speeds and lower relative humidity, reduced incoming longwave radiation, a slight surface pressure increase and a substantial temperature inversion. Warm synoptic conditions involve higher wind speeds and specific humidity, strongly positive anomalies of incoming longwave radiation. The surface pressure is slightly lower than, and the temperature inversion is strongly reduced compared to average conditions. Finally, average wind speeds, humidity and incoming longwave radiation, as well as slightly lower surface pres-  sure, are observed during the transitional regime, when the situation evolves from synoptic to katabatic or the other way around (Gorodetskaya et al., 2013). We therefore investigate the specific meteorological conditions (near-surface temperature inversion, relative humidity, surface temperature, wind speed and direction, in-and outgoing longwave fluxes and the time since the last precipitation event) during blowing snow events.
For all three categories of blowing snow events, the 2 m wind direction shows a preferential easterly or northeasterly orientation at both Neumayer and PE, while the absence of blowing snow is characterized by a wider spectrum of wind directions (Figs. 9 and 10). Positive anomalies in wind speed and relative humidity occur during blowing snow events.
Cyclonic events are a common feature at Neumayer (König-Langlo and Loose, 2007), bringing easterly winds during which most of the drifting and blowing snow occur. Also at PE, most of the blowing snow events (N = 1643, 92 %) are associated with the warm synoptic and transitional regimes, when moist air is brought from the ocean, that precipitate inland (Gorodetskaya et al., 2013). Thiery et al. (2012) also showed that at PE drifting snow sublimation occurs mostly during transitional regimes. These regimes occur 41-48 % of the time (Gorodetskaya et al., 2013. Very few blowing snow events occur in cloudless cold conditions (cold katabatic regime), when the northerly winds blows from the interior towards the coast (N = 139; 8 %).
Intense mixed events (Fig. 5) occur together with strong northeasterly winds -87 • to the north at 10 m · s −1 (PE) and 65 • to the north at 13 m · s −1 (Neumayer III) -warmer surface temperatures and higher relative humidity. These are the signature of storms associated with synoptic events, during which the turbulent mixing reduces the vertical temperature gradient (Gorodetskaya et al., 2013). The majority (60 %) of the blowing snow events occur during storms or overcast conditions (with cloud and/or precipitation). These mixed events have generally a short time lag since the last precipitation event and reach high atmospheric levels. Dry blowing snow has a mean wind direction of 120 • to the north at PE and 77 • at Neumayer III, lower wind speeds (6-7 m · s −1 ) and a greater temperature inversion. The mean time lag since the last precipitation event at PE (23 h) indicates that these events most likely occur shortly after a storm and that cloudless blowing snow (8 %) is mostly associated with katabatic winds. Apart from these factors, sastrugi might also have an impact on blowing snow (Amory et al., 2017) but are not measured here.

Depth of the blowing snow layer
The height of the blowing snow layer (algorithm explained in Sect. 3.2) varies according to different parameters: the wind speed and the size and density of the snow particles. In addition, the presence of clouds and precipitation also influences the vertical extent of the layer. Blowing snow layer depths at Neumayer III and PE show a predominance of shallow layers (55 and 75 % thicknesses below 100 m, respectively; Fig. 11). However, there is little inter-event variability in blowing snow layer height at both stations. At both stations, blowing snow layers with the highest vertical extend occur during blowing snow mixed with precipitation. Mean blowing snow layer height during precipitating event reaches 331 m, while clear-sky mean blowing snow layer depth is limited to 78 m at PE. The values found for both stations are consistent with the mean blowing snow layer height detected by ground-based lidar at South Pole (Mahesh et al., 2003), although somewhat lower. The thickness of the blowing snow layer detected by the BSD algorithm is probably underesti- Figure 11. Distribution of the height of the blowing snow layer at (a) PE station (b) at Neumayer station, blowing snow accompanied with precipitation in blue, blowing snow without precipitation in red. mated in the case of heavy blowing snow events due to total attenuation of the signal before reaching the top of the layer.
We further tested the hypothesis that the height of the blowing snow layer is related to wind speed at PE station. While there is no correlation (also found by Mahesh et al., 2003), the height of the blowing snow layer is related to the time since last precipitation (Fig. 12). The height of the blowing snow layer can reach up to 1000 m within 24 to 48 h after precipitation, and 95 % of the blowing snow layers thicker than 500 m occur shortly after the last precipitation event. Blowing snow events taking place much later after the precipitation event are limited to a vertical extend lower than 100 m thick.

Applicability of the algorithm
The BSD algorithm developed for the Vaisala CL-31 ceilometer at PE was applied to the Vaisala CL-51 ceilometer at Neumayer III station. Comparing the BSD algorithm detections to visual observations at Neumayer showed a good agreement and the ability of the BSD algorithm to detect (heavy) blowing snow events, under both dry and precipitating conditions. Satellite detections of blowing snow, although covering the whole continent, are limited to clear-sky conditions. The BSD algorithm, however, is able to detect blowing snow during most of the storms, which is an improvement as the majority of blowing snow occurs together with cloud and/or precipitation. When we limit the analysis to (heavy) blowing snow, the algorithm detects 78 % of the events indicated by the observer. However, there are cases where the ceilometer does not detect events classified as heavy events by the observer. In such cases, it has to be kept in mind that blowing and drifting snow observations are extremely challenging, with a potentially large but unknown error on the observations. Furthermore, the hourly time filtering applied leads to commission and omission errors (short-lived events are likely removed from the running mean). However, such events induce much smaller blowing snow transport rates than strong events, and we suspect that omitting them will only reduce blowing snow transport rates by a small percentage. A limitation of the BSD algorithm is that both ceilometers are set up on the roof of the station, 17 m at Neumayer III and 12 m above the ground in the main wind direction at PE. In addition, 15 m have to be added to account for the discard of the first range bin. There, ceilometers will report the most significant blowing snow events (higher than 30 m) and most drifting snow and shallow blowing snow events are not detected. If setting up a ceilometer with the aim of measuring blowing snow, the device should be placed as close to the ground as possible to also retrieve shallower blowing snow events. Ceilometers can retrieve the presence of blowing snow, but other properties such as optical depth, size and density measurement are only possible if the ceilometer is calibrated, which is very challenging and not done in this paper. The BSD algorithm can be applied to any ceilometer located in Antarctica, but we recommend using a bin width of 10 m for operating ceilometers to detect blowing snow, which is the case at PE and Neumayer III. Since the Vaisala CT25K at Halley station uses a 30 m vertical resolution, it was not used in this study. Gallée et al. (2001) stated that snowpack properties mainly determine snow erosion: dendricity, density, sphericity and particles size regulate the availability of snow for transportation. These parameters change with metamorphism and impact the threshold friction velocity and thus the minimum wind speed required for particles uplift from the ground. Here, we do not apply any wind speed threshold to the detection of blowing snow, whereas some modeling studies assume a drifting snow dependency on temperature and wind speed (Giovinetto et al., 1992;Yau, 1999, 2002;Yang et al., 2010). Palm et al. (2011), for instance, used a minimum wind speed criterion to detect blowing snow from satellite backscatter, potentially leaving out some events.

Wind speed vs. snow availability
We find that the presence of freshly fallen snow has a great impact on blowing snow occurrence and blowing snow layer height. As postulated by Mahesh et al. (2003), the end of a large snow storm with high wind speeds could still hold snow particles suspended in the air, even if the wind speed has already dropped to lower speeds than those required to dislodge the particles from the ground at the onset of the blowing snow event. Conversely, if no particles are available for the wind to pick up, blowing snow might not occur even though the wind speeds are high. The large majority of blowing snow events occur under synoptic disturbances (92 % at PE and Neumayer III) rather than katabatic conditions. These disturbances are also associated with higher wind speeds and are often accompanied with precipitation. In those cases, snow is available for transport. At PE, the explanation for the limited occurrence of blowing snow under katabatic conditions might lie in the fact that the station is shielded by the Sør Rondane Mountains, as well as the limited availability of fresh snow and the reduced turbulence during those events compared to synoptic conditions, maintaining particles aloft. This, together with the reduced number of blowing snow events occurring under katabatic winds, might indicate that the effect of katabatic winds on blowing snow occurrence has been overestimated at both PE and Neumayer stations and that synoptic events bringing fresh snow is most likely a determining factor for blowing snow at those locations.

Conclusions
Various observations, models and satellite studies have been performed to quantify and investigate blowing snow on the Antarctic continent. We present here our novel BSD algorithm, designed to retrieve blowing snow events, but not drifting snow, from ground-based remote-sensing ceilometers.
The algorithm has proven to be reliable in detecting blowing snow at Neumayer station in up to 78 % of the cases when compared to visual observations. The presence of precipitation does not substantially limit the retrieval by the ceilometer. This is an improvement to satellite detection, which is limited to clear-sky conditions and therefore missing a great part of the blowing snow as more than half of the blowing snow happens during a storm at PE and Neumayer III stations. Blowing snow detected by the BSD algorithm occurs 36 % of the time at Neumayer station and 13 % at PE station, with an interannual variation of ±5 % and seasonal cycle that peaks during the Antarctic winter. We further conclude that most of the blowing snow events happen during or shortly after precipitation, brought to the continent by the easterly winds associated with synoptic systems. The availability of fresh snow mainly determines the onset of blowing snow, and the available fresh snow can be lifted to higher heights than during katabatic conditions at PE and Neumayer stations. This highlights again the limitation of wind speed thresholds, when applied to blowing snow retrieval methods. The properties of the snow particles, as well as the availability of fresh snow, need to be taken into account in order to accurately initiate blowing snow in models. Since ceilometers are low-cost robust instruments, and often deployed at stations for the purpose of aircraft operations, our newly developed algorithm opens opportunities for long-term monitoring networks of consistent blowing snow observations. These can further be used to evaluate satellite retrieval and combined to produce blowing snow products over the ice sheets.
Code availability. The algorithm is freely available upon request to alexandra.gossart@kuleuven.be Data availability. Data from Neumayer station are freely available on the Pangaea portal (König-Langlo, 2011, 2013, 2014, 2016 and data from the instruments at Princess Elisabeth station are available upon request (www.aerocloud.be).