A New Methodology to Characterise the Radar Bright Band Using Doppler Spectral Moments from Vertically Pointing Radar Observations

The detection and characterisation of the radar Bright Band (BB) are essential for many applications of weather radar quantitative precipitation estimates, such as heavy rainfall surveillance, hydrological modelling or numerical weather prediction data assimilation. This study presents a new technique to detect the radar BB levels (top, peak and bottom) for Doppler radar spectral moments from the vertically pointing radars applied here to a K-band radar, the MRR-Pro (Micro Rain Radar). The methodology includes signal and noise detection and dealiasing schemes to provide realistic vertical Doppler velocities of precipitating hydrometeors, subsequent calculation of Doppler moments and associated parameters and BB detection and characterisation. Retrieved BB properties are compared with the melting level provided by the MRR-Pro manufacturer software and also with the 0 °C levels for both dry-bulb temperature (freezing level) and wet-bulb temperature from co-located radio soundings in 39 days. In addition, a co-located Parsivel disdrometer is used to analyse the equivalent reflectivity of the lowest radar height bins confirming consistent results of the new signal and noise detection scheme. The processing methodology is coded in a Python program called RaProM-Pro which is freely available in the GitHub repository.


Introduction
Precipitating hydrometeors undergo various processes as they fall, including water vapour condensation, coalescence, break-up or evaporation for liquid water and ice nucleation, riming, aggregation or accretion for the solid phase [1]. One of the most important processes occurs as falling particles cross the 0 • C isotherm level, also called melting level, where solid water particles begin to melt and eventually transform completely into liquid particles [2,3]. The atmospheric layer where this process takes place is known as the melting layer and may produce a characteristic radar signature, the so-called radar Bright Band (hereafter BB), a term originated from the local maxima caused by high reflectivity values visible in the equivalent reflectivity vertical profile [4]. The BB is caused by differences in the dielectric constants, shape and terminal fall speeds of liquid and solid hydrometeor precipitating particles, which lead to abrupt changes of the radar backscattered power within the BB. The most evident BB signatures are produced under stratiform cold rain conditions [5,6] as updrafts, and vertical mixing present in convective precipitation do not provide the proper conditions for BB formation.
The presence of a BB in volumetric operational weather radar observations may produce local overestimations of rainfall amounts, which has led to the development of different procedures to detect and correct BB effects [2,[7][8][9]. This is particularly important for events with rapidly changing characteristics, for example, with quick transitions from An OTT Parsivel-2 disdrometer [29], hereafter Parsivel, co-located with the MRR-Pro, provided precipitation particle size and fall speed spectra at the radar level. These parameters allow comparisons between the MRR-Pro and Parsivel for different variables, such as rainfall rate or radar reflectivity. Finally, on the same roof is located the Barcelona Remote Sens. 2021, 13, 4323 3 of 20 radio sounding station (WMO code 08190), which performs two soundings a day (at 00 and 12 UTC) that were also used in this study.

Data Acquisition
The data files generated by the MRR-Pro manufacturer software are used as input files for the processing with RaProM-Pro. These files are in netcdf format and contain basic configuration settings of the data acquisition, the raw data (so-called spectral reflectivity, or spectrum raw according to the manufacturer) and derived parameters, such as radial velocity spectra or an estimate of existence of the melting layer. Figure 1 schematically shows a selection of the content of the files and also indicates which variables are used in the proposed methodology (shown in green), which are, essentially, configuration parameters and raw data, from which derived parameters can be calculated.

Interval of velocity
Δv m·s −1 0.19 An OTT Parsivel-2 disdrometer [29], hereafter Parsivel, co-located with the MRR-Pro, provided precipitation particle size and fall speed spectra at the radar level. These parameters allow comparisons between the MRR-Pro and Parsivel for different variables, such as rainfall rate or radar reflectivity. Finally, on the same roof is located the Barcelona radio sounding station (WMO code 08190), which performs two soundings a day (at 00 and 12 UTC) that were also used in this study.

Data Acquisition
The data files generated by the MRR-Pro manufacturer software are used as input files for the processing with RaProM-Pro. These files are in netcdf format and contain basic configuration settings of the data acquisition, the raw data (so-called spectral reflectivity, or spectrum raw according to the manufacturer) and derived parameters, such as radial velocity spectra or an estimate of existence of the melting layer. Figure 1 schematically shows a selection of the content of the files and also indicates which variables are used in the proposed methodology (shown in green), which are, essentially, configuration parameters and raw data, from which derived parameters can be calculated.  The parameters used by RaProM-Pro include, for each time step, a matrix s(n, i) with the spectral raw values (labelled by the manufacturer as "spectrum raw") for all vertical levels n and Doppler frequencies i. The matrix s(n, i) contains the ratio between emitted and received power after the Fourier Transform is computed by the radar and represent the intensity of the echo backscattered by the precipitation particles. The spectral values are processed with the information provided in different arrays, such as "range" (list of Remote Sens. 2021, 13, 4323 4 of 20 heights above sea level) or "transfer function", which allows for correcting them according to their distance from the radar. Finally, a calibration constant unique to each radar unit is also provided.
To check the consistency of the new methodology proposed, additional parameters provided by the manufacturer-see [30]-are also used for comparison. These parameters include the SNR (signal to noise ratio), or the Doppler radial velocity and spectrum width and four different versions of radar reflectivity. These versions of radar reflectivity and equivalent radar reflectivity are computed in two stages. Firstly, without considering rainfall attenuation effects (in the so-called attenuated version of these variables, Za and Zea). Then, after the calculation of the drop size distribution (N), an estimate of the rainfall attenuation is calculated considering the Path Integrated rain Attenuation (PIA), which then allows the non-attenuated version of the reflectivity and equivalent reflectivity to be computed (Z, Ze). Finally, three additional parameters are considered for each height level: an estimate of the Melting Layer (ML)-expressed as a probability, a value between 0 and 1-the Liquid Water Content (LWC) and the Rainfall Rate (RR).

Processing Method
The processing software provided by the MRR-Pro manufacturer performs reasonably well in most meteorological conditions. However, in some cases, the original de-aliasing method provides limited results, as illustrated in Section 3.2. In order to develop a new de-aliasing scheme, spectral reflectivity has to be computed, so a new approach is also considered for the signal and noise processing described in this section.
The proposed processing method starts from the transformation of spectral raw data values read from the netcdf matrix S, for each level n and Doppler bin i, to their physical value given by spectral reflectivity (η), as described in two-steps in the following Equations (1) and (2): where CC is the calibration constant, TF(n) is the transfer function, δr is the height resolution, n is the number of height gates and i is the number of Doppler bins. From the spectral reflectivity, it is possible to calculate several physical parameters, such as hydrometeor velocity, equivalent radar reflectivity and the precipitation type classification, as described by [31]. The processing method consists of the following four main stages ( Figure 2): (1). removal of noise and peaks detection from the raw signal, (2). dealiasing of the spectrum to improve the detection of the vertical velocity, (3). computation of attenuation path integrated (PIA) factors and (4). calculation of radar parameters using the corrected spectrum and the BB characterisation. The results are saved in a netcdf output file. Stages (2) and (4) are particularly novel. Note that Stage (1) must be performed before Stage (2) as spectral reflectivity (η), required for dealiasing, is not available in the manufacturer's netcdf file.
The parameters used by RaProM-Pro include, for each time step, a matrix s(n,i) with the spectral raw values (labelled by the manufacturer as "spectrum raw") for all vertical levels n and Doppler frequencies i. The matrix s(n,i) contains the ratio between emitted and received power after the Fourier Transform is computed by the radar and represent the intensity of the echo backscattered by the precipitation particles. The spectral values are processed with the information provided in different arrays, such as "range" (list of heights above sea level) or "transfer function", which allows for correcting them according to their distance from the radar. Finally, a calibration constant unique to each radar unit is also provided.
To check the consistency of the new methodology proposed, additional parameters provided by the manufacturer-see [30]-are also used for comparison. These parameters include the SNR (signal to noise ratio), or the Doppler radial velocity and spectrum width and four different versions of radar reflectivity. These versions of radar reflectivity and equivalent radar reflectivity are computed in two stages. Firstly, without considering rainfall attenuation effects (in the so-called attenuated version of these variables, Za and Zea). Then, after the calculation of the drop size distribution (N), an estimate of the rainfall attenuation is calculated considering the Path Integrated rain Attenuation (PIA), which then allows the non-attenuated version of the reflectivity and equivalent reflectivity to be computed (Z, Ze). Finally, three additional parameters are considered for each height level: an estimate of the Melting Layer (ML)-expressed as a probability, a value between 0 and 1-the Liquid Water Content (LWC) and the Rainfall Rate (RR).

Processing Method
The processing software provided by the MRR-Pro manufacturer performs reasonably well in most meteorological conditions. However, in some cases, the original de-aliasing method provides limited results, as illustrated in Section 3.2. In order to develop a new de-aliasing scheme, spectral reflectivity has to be computed, so a new approach is also considered for the signal and noise processing described in this section.
The proposed processing method starts from the transformation of spectral raw data values read from the netcdf matrix S, for each level n and Doppler bin i, to their physical value given by spectral reflectivity (ƞ), as described in two-steps in the following Equations (1) and (2): where CC is the calibration constant, TF(n) is the transfer function, δr is the height resolution, n is the number of height gates and i is the number of Doppler bins. From the spectral reflectivity, it is possible to calculate several physical parameters, such as hydrometeor velocity, equivalent radar reflectivity and the precipitation type classification, as described by [31]. The processing method consists of the following four main stages ( Figure 2): (1). removal of noise and peaks detection from the raw signal, (2). dealiasing of the spectrum to improve the detection of the vertical velocity, (3). computation of attenuation path integrated (PIA) factors and (4). calculation of radar parameters using the corrected spectrum and the BB characterisation. The results are saved in a netcdf output file. Stages (2) and (4) are particularly novel. Note that Stage (1) must be performed before Stage (2) as spectral reflectivity (ƞ), required for dealiasing, is not available in the manufacturer's netcdf file.

Signal and Noise Detection
The signal and noise separation is performed considering the algorithm proposed by [32], similarly to what was described by [30], but considering two steps. The first step consists of comparing the ratio of the squared mean spectral reflectivity and its variance with a specific threshold or limit fixed for a given integration time, given by Equation (3): where Limit equals the time resolution chosen (T i ). This step is applied iteratively while the condition is verified. Each iteration implies evaluating a peak candidate, and if (3) is fulfilled, then the peak is discarded. The signal remains until the condition is false and will be considered background noise. More details of the implementation of this first step are detailed in [31]. The second step adds a new condition where the spectral reflectivity peak divided by the mean of the spectrum must be equal to or greater than a threshold value equal to 1.3 as shown in Equation (4): in order to be considered a real signal. Note that (4) is applied after verifying (3) so that both conditions must be satisfied. Then, the next step is the noise determination. The SNR is calculated using its definition expressed in dB, according to the manufacturer's documentation, given by Equation (5): It is noted that SNR values provided by the manufacturer are substantially lower than those obtained with RaProM-Pro; in particular, they contain negative SNR values, i.e., signal below the noise level. This is a consequence of the different schemes applied for signal determination by the manufacturer and the methodology proposed, despite other derived variables presenting very similar values. An example is shown in Figure 3 comparing values obtained by RaProM-Pro and the manufacturer for equivalent reflectivity, SNR and Doppler velocity (Figure 3a-c respectively). SNR values present systematic differences around 20 dB but very similar values for the other variables, except for a few vertical velocity outliers due to the different dealiasing methods discussed below.

Signal and Noise Detection
The signal and noise separation is performed considering the algorithm proposed by [32], similarly to what was described by [30], but considering two steps. The first step consists of comparing the ratio of the squared mean spectral reflectivity and its variance with a specific threshold or limit fixed for a given integration time, given by Equation (3): where Limit equals the time resolution chosen (Ti). This step is applied iteratively while the condition is verified. Each iteration implies evaluating a peak candidate, and if (3) is fulfilled, then the peak is discarded. The signal remains until the condition is false and will be considered background noise. More details of the implementation of this first step are detailed in [31]. The second step adds a new condition where the spectral reflectivity peak divided by the mean of the spectrum must be equal to or greater than a threshold value equal to 1.3 as shown in Equation (4): in order to be considered a real signal. Note that (4) is applied after verifying (3) so that both conditions must be satisfied. Then, the next step is the noise determination. The SNR is calculated using its definition expressed in dB, according to the manufacturer's documentation, given by Equation (5): It is noted that SNR values provided by the manufacturer are substantially lower than those obtained with RaProM-Pro; in particular, they contain negative SNR values, i.e., signal below the noise level. This is a consequence of the different schemes applied for signal determination by the manufacturer and the methodology proposed, despite other derived variables presenting very similar values. An example is shown in Figure 3 comparing values obtained by RaProM-Pro and the manufacturer for equivalent reflectivity, SNR and Doppler velocity (Figures 3a, 3b and 3c respectively). SNR values present systematic differences around 20 dB but very similar values for the other variables, except for a few vertical velocity outliers due to the different dealiasing methods discussed below.  An additional analysis is performed for radar reflectivity comparing the lowest valid radar height bin (from 150 to 200 m above radar level) and the co-located Parsivel disdrometer ( Figure 4) considering 1 min sampling periods. Both radar processing schemes compare very well with Parsivel values, with slight discrepancies that may be explained by instrumental differences-see [33]. More details about the signal and noise detection scheme can be found in Appendix A.
An additional analysis is performed for radar reflectivity comparing the lowest valid radar height bin (from 150 to 200 m above radar level) and the co-located Parsivel disdrometer ( Figure 4) considering 1 min sampling periods. Both radar processing schemes compare very well with Parsivel values, with slight discrepancies that may be explained by instrumental differences-see [33]. More details about the signal and noise detection scheme can be found in Appendix A.

Dealiasing
Spectral reflectivity aliasing occurs when the target returns a signal outside the unambiguous range interval. A systematic method to correct aliasing in MRR-2 was proposed by [34] and was implemented with some modifications by [31]. The two methods are based on the estimated velocity parameters calculated from equivalent radar reflectivity in [35]. Here we propose a different approach, where only signal continuity between vertical levels is used, instead of the parameters estimated in [35]. According to the manufacturer's documentation, the radar manufacturer processing is able to detect upward movements of precipitation particles, but in some cases, this detection is not possible, and velocities are aliased. Figure 5 shows an example where the manufacturer velocity spectrum ( Figure 5a) shows a suspicious pattern between 4000 and 5000 m, potentially caused by aliasing. By extending or unfolding the spectra to both sides (Figure 5b), the vertical continuity of the spectra allows a consistent dealiased spectra profile to be selected ( Figure 5c).

Dealiasing
Spectral reflectivity aliasing occurs when the target returns a signal outside the unambiguous range interval. A systematic method to correct aliasing in MRR-2 was proposed by [34] and was implemented with some modifications by [31]. The two methods are based on the estimated velocity parameters calculated from equivalent radar reflectivity in [35].
Here we propose a different approach, where only signal continuity between vertical levels is used, instead of the parameters estimated in [35]. According to the manufacturer's documentation, the radar manufacturer processing is able to detect upward movements of precipitation particles, but in some cases, this detection is not possible, and velocities are aliased. Figure 5 shows an example where the manufacturer velocity spectrum ( Figure 5a) shows a suspicious pattern between 4000 and 5000 m, potentially caused by aliasing. By extending or unfolding the spectra to both sides (Figure 5b), the vertical continuity of the spectra allows a consistent dealiased spectra profile to be selected ( Figure 5c).
An additional analysis is performed for radar reflectivity comparing the lowest valid radar height bin (from 150 to 200 m above radar level) and the co-located Parsivel disdrometer ( Figure 4) considering 1 min sampling periods. Both radar processing schemes compare very well with Parsivel values, with slight discrepancies that may be explained by instrumental differences-see [33]. More details about the signal and noise detection scheme can be found in Appendix A.

Dealiasing
Spectral reflectivity aliasing occurs when the target returns a signal outside the unambiguous range interval. A systematic method to correct aliasing in MRR-2 was proposed by [34] and was implemented with some modifications by [31]. The two methods are based on the estimated velocity parameters calculated from equivalent radar reflectivity in [35]. Here we propose a different approach, where only signal continuity between vertical levels is used, instead of the parameters estimated in [35]. According to the manufacturer's documentation, the radar manufacturer processing is able to detect upward movements of precipitation particles, but in some cases, this detection is not possible, and velocities are aliased. Figure 5 shows an example where the manufacturer velocity spectrum ( Figure 5a) shows a suspicious pattern between 4000 and 5000 m, potentially caused by aliasing. By extending or unfolding the spectra to both sides (Figure 5b), the vertical continuity of the spectra allows a consistent dealiased spectra profile to be selected ( Figure 5c).  The dealiased spectral reflectivity allows, in this case, to detect upward movements of precipitation particles between 4000 and 5000 m. Figure 6 shows the corresponding time-height display of this case where RaProM-Pro detects upward movements, unlike the manufacturer original output, which indicates high downward values. More challenging cases, for example, with convective precipitation and strong windshear might not be de-tected by the new proposed scheme, which was designed to deal with typical BB conditions (see Appendix B for more details).
The dealiased spectral reflectivity allows, in this case, to detect upward movements of precipitation particles between 4000 and 5000 m. Figure 6 shows the corresponding time-height display of this case where RaProM-Pro detects upward movements, unlike the manufacturer original output, which indicates high downward values. More challenging cases, for example, with convective precipitation and strong windshear might not be detected by the new proposed scheme, which was designed to deal with typical BB conditions (see Appendix B for more details). After the dealiasing is applied, different Doppler moments from the spectral reflectivity are computed, including the equivalent radar reflectivity (dBZ), the Doppler velocity (m/s), the spectral width (m/s), the skewness and the kurtosis (Equations (6)-(10)): where λ is the radar wavelength, |K| 2 is the dielectric factor, in this case, liquid water and Δv is the Nyquist velocity. Note that the radar reflectivity does not yet consider the possible effects of rainfall attenuation, which is computed in the next subsection.

Attenuation Calculation
Weather radars operating in attenuated frequencies, such as the K-band, may be affected by rainfall attenuation, impacting specific parameters such as radar reflectivity (Z), liquid water content (LWC) and rain rate (RR). Attenuation is calculated to determine the amount of signal loss integrated along a path (in height) by absorption and scattering by precipitating particles. PIA values are computed following an iterative process described in [36], shown schematically in Figure 7. After the dealiasing is applied, different Doppler moments from the spectral reflectivity are computed, including the equivalent radar reflectivity (dBZ), the Doppler velocity (m/s), the spectral width (m/s), the skewness and the kurtosis (Equations (6)-(10)): where λ is the radar wavelength, |K| 2 is the dielectric factor, in this case, liquid water and ∆v is the Nyquist velocity. Note that the radar reflectivity does not yet consider the possible effects of rainfall attenuation, which is computed in the next subsection.

Attenuation Calculation
Weather radars operating in attenuated frequencies, such as the K-band, may be affected by rainfall attenuation, impacting specific parameters such as radar reflectivity (Z), liquid water content (LWC) and rain rate (RR). Attenuation is calculated to determine the amount of signal loss integrated along a path (in height) by absorption and scattering by precipitating particles. PIA values are computed following an iterative process described in [36], shown schematically in Figure 7.
Essentially, drop size distributions N'(D, n), at each level n, are calculated considering an attenuation factor (the PIA) multiplied by the previous (attenuated) drop size distribution Na(D, n), computed from the Doppler spectra assuming Mie scattering conditions [37,38]. As these calculations are only valid for liquid precipitation particles falling at terminal fall speeds, an additional procedure that provides a hydrometeor classification type for each bin height [31] is applied so that attenuation can be used consistently only for liquid precipitation. Essentially, drop size distributions N'(D, n), at each level n, are calculated considering an attenuation factor (the PIA) multiplied by the previous (attenuated) drop size distribution Na(D, n), computed from the Doppler spectra assuming Mie scattering conditions [37,38]. As these calculations are only valid for liquid precipitation particles falling at terminal fall speeds, an additional procedure that provides a hydrometeor classification type for each bin height [31] is applied so that attenuation can be used consistently only for liquid precipitation.  [36]. N'(D, n) and Na(D, n) are, respectively, the drop size distribution and the attenuated drop size distribution; ke is the specific rain attenuation and σe is the single-particle extinction coefficient, calculated with the Mie theory.
As described in Figure 7, the maximum PIA value is 10, because for higher values, the scheme may not work properly. If PIA reaches the value of 10, the manufacturer processing stops calculating it for higher range bins. However, RaProM-Pro assigns a constant value of 10 for internal processing reasons. The final output parameter of PIA in RaProM-Pro is simply the PIA value expressed in dB (11), called DBPIA: Figure 7. PIA calculation flow chart adapted from [36]. N'(D, n) and Na(D, n) are, respectively, the drop size distribution and the attenuated drop size distribution; k e is the specific rain attenuation and σ e is the single-particle extinction coefficient, calculated with the Mie theory.
As described in Figure 7, the maximum PIA value is 10, because for higher values, the scheme may not work properly. If PIA reaches the value of 10, the manufacturer processing stops calculating it for higher range bins. However, RaProM-Pro assigns a constant value of 10 for internal processing reasons. The final output parameter of PIA in RaProM-Pro is simply the PIA value expressed in dB (11), called DBPIA: The DBPIA calculation is included in RaProM-Pro processing despite it not being applied in the BB determination procedure

Bright Band Calculation
The new methodology proposed to detect the BB is based on [39], plus a novel approach considering the vertical variation of the skewness computed from the spectrum Doppler velocity at each height. Skewness provides information about the asymmetry of the fall velocity distribution and indicates that in the BB, snowflakes or ice particles have started to melt [3,4,22]. The change of shape and aerodynamics of the solid particles as they melt modifies the averaged Doppler velocity and the spectrum shape. This change can be observed in the velocity distribution provided by the spectrum reflectivity at each height, where the maximum value changes from being tilted to the right to being tilted to the left, which implies a change of sign of the skewness. Figure 8 shows an example observed by the MRR-Pro, highlighting the different spectra shape above, within and below the BB calculated by RaProM-Pro.
The DBPIA calculation is included in RaProM-Pro processing despite it not being applied in the BB determination procedure

Bright Band Calculation
The new methodology proposed to detect the BB is based on [39], plus a novel approach considering the vertical variation of the skewness computed from the spectrum Doppler velocity at each height. Skewness provides information about the asymmetry of the fall velocity distribution and indicates that in the BB, snowflakes or ice particles have started to melt [3,4,22]. The change of shape and aerodynamics of the solid particles as they melt modifies the averaged Doppler velocity and the spectrum shape. This change can be observed in the velocity distribution provided by the spectrum reflectivity at each height, where the maximum value changes from being tilted to the right to being tilted to the left, which implies a change of sign of the skewness. Figure 8 shows an example observed by the MRR-Pro, highlighting the different spectra shape above, within and below the BB calculated by RaProM-Pro.  The change in the shape of the spectra is clearly visible in Figure 8, where the progressive appearance of raindrops at the expense of melted solid particles (Figure 8b) modifies the Doppler velocity spectra, widening it to the right due to higher fall speeds. This leads to a symmetric or slightly right-skewed spectrum distribution, which implies a change of the skewness from negative to positive values. Determining the height where the skewness sign changes is thus a key feature to obtain the height of the BB.
The method proposed is detailed in the flowchart shown in Figure 9, which describes, first, the BB detection approach, based on [39], and then the BB characterisation, which computes the BB top and the BB bottom. The remaining BB feature, the BBpeak, is the level located between the BB bottom and the BB top, where the skewness is maximum and should be close to the melting level. An additional checking is performed to remove BB detections of virga cases, simply verifying that precipitation reaches the ground.
The procedure is applied to each MRR-Pro vertical profile (in our case, available every 10 s). Then the results (BB top, BB peak and BB bottom) are smoothed temporally, considering a generalised exponential moving average [40], allowing a more continuous signal of BB characteristics, but keeping the original temporal resolution. Note that the current implementation of the scheme detects the lowest BB present in a vertical profile. The method proposed is detailed in the flowchart shown in Figure 9, which describes, first, the BB detection approach, based on [39], and then the BB characterisation, which computes the BB top and the BB bottom. The remaining BB feature, the BBpeak, is the level located between the BB bottom and the BB top, where the skewness is maximum and should be close to the melting level. An additional checking is performed to remove BB detections of virga cases, simply verifying that precipitation reaches the ground.
The procedure is applied to each MRR-Pro vertical profile (in our case, available every 10 s). Then the results (BB top, BB peak and BB bottom) are smoothed temporally, considering a generalised exponential moving average [40], allowing a more continuous signal of BB characteristics, but keeping the original temporal resolution. Note that the current implementation of the scheme detects the lowest BB present in a vertical profile.

Results
The methodology presented in the previous section is illustrated in Figure 10, displaying both radiosonde data ( Figure 10a) and MRR-Pro data (Figure 10b) for a clear BB case. The figure shows the sounding profiles of dry-bulb, wet-bulb and dew point temperature (Figure 10a) recorded the 31 March 2020 at 12 UTC and radar equivalent

Results
The methodology presented in the previous section is illustrated in Figure 10, displaying both radiosonde data ( Figure 10a) and MRR-Pro data (Figure 10b) for a clear BB case.
The figure shows the sounding profiles of dry-bulb, wet-bulb and dew point temperature (Figure 10a) recorded the 31 March 2020 at 12 UTC and radar equivalent reflectivity Ze from MRR-Pro and Parsivel observations; the latter was plotted at the lowest height level and resampled at a 10 s resolution to match MRR-Pro observations, from 12 to 15 UTC (Figure 10b). The sounding plot panel explicitly shows that the 0 • C wet-bulb temperature is relatively lower than the freezing level (0 • C dry-bulb temperature) due to low saturation. The wind profile is also plotted, showing west wind components below the BB, which is consistent with the precipitation fall streaks visible from the BB. Figure 10b shows consistent reflectivity values of ground measures from Parsivel and the first lowest valid height bin from the radar, for example, the alternating maxima and minima reflectivity columns. The 12 UTC melting levels and 0 • C wet-bulb temperature levels obtained with the sounding match the detected BB top and BB bottom well, respectively. This example illustrates well the ability of the new methodology to provide a temporal evolution of the BB details with a 10 s resolution.

Results
The methodology presented in the previous section is illustrated in Figure 10, displaying both radiosonde data ( Figure 10a) and MRR-Pro data (Figure 10b) for a clear BB case. The figure shows the sounding profiles of dry-bulb, wet-bulb and dew point temperature (Figure 10a) recorded the 31 March 2020 at 12 UTC and radar equivalent reflectivity Ze from MRR-Pro and Parsivel observations; the latter was plotted at the lowest height level and resampled at a 10 s resolution to match MRR-Pro observations, from 12 to 15 UTC (Figure 10b). The sounding plot panel explicitly shows that the 0 °C wet-bulb temperature is relatively lower than the freezing level (0 °C dry-bulb temperature) due to low saturation. The wind profile is also plotted, showing west wind components below the BB, which is consistent with the precipitation fall streaks visible from the BB. Figure  10b shows consistent reflectivity values of ground measures from Parsivel and the first lowest valid height bin from the radar, for example, the alternating maxima and minima reflectivity columns. The 12 UTC melting levels and 0 °C wet-bulb temperature levels obtained with the sounding match the detected BB top and BB bottom well, respectively. This example illustrates well the ability of the new methodology to provide a temporal evolution of the BB details with a 10 s resolution. The results are presented in the following subsections, considering three different statistical comparisons. The first one is performed comparing the manufacturer's BB product and the proposed methodology. The second and third subsections compare, respectively, the new methodology and the original manufacturer product, with radio sounding observations.

Manufacturer ML Height vs. RaProM-Pro BBpeak Height
The previous example is the starting point to assess differences between the manufacturer's BB product and the new proposed methodology. Figure 11 shows radar reflectivity and Doppler velocity profiles processed with each methodology, also including the ML manufacturer's product and the proposed BB product. In this case, the melting layer heights computed by the manufacturer are always plotted between the estimated BB top and BB bottom, which provide a more complete description of the BB. Additionally, RaProM-Pro offers a more detailed description of the precipitation field (better depiction of contours and inclusion of additional weather echoes) around 15 UTC, particularly above the BB (from 4000 to 5000 m) thanks to the new noise and peak detection methodology. The results are presented in the following subsections, considering three different statistical comparisons. The first one is performed comparing the manufacturer's BB product and the proposed methodology. The second and third subsections compare, respectively, the new methodology and the original manufacturer product, with radio sounding observations.

Manufacturer ML Height vs. RaProM-Pro BB peak Height
The previous example is the starting point to assess differences between the manufacturer's BB product and the new proposed methodology. Figure 11 shows radar reflectivity and Doppler velocity profiles processed with each methodology, also including the ML manufacturer's product and the proposed BB product. In this case, the melting layer heights computed by the manufacturer are always plotted between the estimated BB top and BB bottom, which provide a more complete description of the BB. Additionally, RaProM-Pro offers a more detailed description of the precipitation field (better depiction of contours and inclusion of additional weather echoes) around 15 UTC, particularly above the BB (from 4000 to 5000 m) thanks to the new noise and peak detection methodology. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 20 A quantitative analysis is provided by comparing the melting layer (ML) height provided by the manufacturer and the BBpeak height calculated with RaProM-Pro, given by: The ML product is calculated using an Artificial Intelligence approach [41], providing probability values; a probability greater than 0.75 is considered a reliable ML height determination, being the maximum probability of the selected height for the melting level. The BB peak has been described in Section 3. Figure 12 shows a histogram of Diff_ML for 39 days, selecting a time window of ±1 h from the sounding launch, resulting in around 2600 cases. It displays a nearly symmetrical distribution pattern with a single-mode centred in the second negative class (from −50 to −100 m), indicating that the manufacturer's ML height is slightly lower than the BBpeak (the averaged value of Diff_ML is −89 m and the standard deviation is 180 m). More than 80% of cases do not exceed 200 m, and the tails of the distribution fall quickly, which despite some differences, reach values higher than 400 m. A quantitative analysis is provided by comparing the melting layer (ML) height provided by the manufacturer and the BB peak height calculated with RaProM-Pro, given by: The ML product is calculated using an Artificial Intelligence approach [41], providing probability values; a probability greater than 0.75 is considered a reliable ML height determination, being the maximum probability of the selected height for the melting level. The BB peak has been described in Section 3. Figure 12 shows a histogram of Diff_ML for 39 days, selecting a time window of ±1 h from the sounding launch, resulting in around 2600 cases. It displays a nearly symmetrical distribution pattern with a single-mode centred in the second negative class (from −50 to −100 m), indicating that the manufacturer's ML height is slightly lower than the BB peak (the averaged value of Diff_ML is −89 m and the standard deviation is 180 m). More than 80% of cases do not exceed 200 m, and the tails of the distribution fall quickly, which despite some differences, reach values higher than 400 m. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 20 Figure 12. Histogram of differences between the melting level height provided by the manufacturer and the Bright Band peak height from RaProM-Pro. The height bin resolution is equal to the histogram class bin, 50 m.

Sounding Observations vs. RaProM-Pro BB Levels
The 0 °C dry-bulb temperature level (freezing level, hereafter h0) and the 0 °C wetbulb temperature level (hereafter hw,0) are compared here with the BBtop, BBpeak and BBbot heights. We consider the same 39 days studied in the previous subsection, with the same time intervals of ±1 h from the sounding launch time to minimise spurious differences caused by rapidly changing BB heights, leading to around 5327 cases.
The evaluation is performed using the parameters Diff_top, Diff_bot, Diff_peak and Diff_Tw defined as: where BBtop, BBbot and BBpeak are the heights from the BB top, BB bottom and BB peak, respectively. Note that a priori, Diff_top should be close to 0 as solid particles begin to melt when they reach the melting level, Diff_bot should be greater than 0, as it takes some time to completely melt all solid particles and, by definition, Diff_peak, should be between Diff_top and Diff_bot. Regarding the expected value of Diff_Tw, the recent study of [42] indicated that BBpeak heights were very similar to hw,0, so a value close to 0 would be consistent with that result. Figure 13a shows histograms of Diff_top and Diff_bot, indicating similar patterns but centred, respectively, below and above h0: BBtop mode is between 150 and 200 m, and BBbot is between −250 and −2000 m. The fact that BBtop occurs mostly above h0 (so it is not close to 0 as initially expected) can be explained by the tendency of solid particles to increase aggregation just above the melting level, producing larger snowflakes and reducing the

Sounding Observations vs. RaProM-Pro BB Levels
The 0 • C dry-bulb temperature level (freezing level, hereafter h 0 ) and the 0 • C wetbulb temperature level (hereafter h w,0 ) are compared here with the BB top , BB peak and BB bot heights. We consider the same 39 days studied in the previous subsection, with the same time intervals of ±1 h from the sounding launch time to minimise spurious differences caused by rapidly changing BB heights, leading to around 5327 cases.
The evaluation is performed using the parameters Diff_top, Diff_bot, Diff_peak and Diff_Tw defined as: where BB top , BB bot and BB peak are the heights from the BB top, BB bottom and BB peak, respectively. Note that a priori, Diff_top should be close to 0 as solid particles begin to melt when they reach the melting level, Diff_bot should be greater than 0, as it takes some time to completely melt all solid particles and, by definition, Diff_peak, should be between Diff_top and Diff_bot. Regarding the expected value of Diff_Tw, the recent study of [42] indicated that BB peak heights were very similar to h w,0 , so a value close to 0 would be consistent with that result. Figure 13a shows histograms of Diff_top and Diff_bot, indicating similar patterns but centred, respectively, below and above h 0 : BB top mode is between 150 and 200 m, and BB bot is between −250 and −2000 m. The fact that BB top occurs mostly above h 0 (so it is not close to 0 as initially expected) can be explained by the tendency of solid particles to increase aggregation just above the melting level, producing larger snowflakes and reducing the number of smaller particles [6]. This would lead to a change in the skewness spectrum, which would be detected by the proposed methodology as the BB top . On the other hand, the mode value of Diff_bot found is reasonable, compared with previous existing studies, such as [4]. Figure 13b shows that the Diff_peak histogram presents a similar pattern to Diff_top and Diff_bot, but, as expected, the mode value, more pronounced (corresponding to a more leptokurtic distribution), is centred between the previous two modes, close to 0. Finally, Diff_T w presents a slightly thicker mode, with a maximum between −50 and 0 m, which is just one class below the Diff_peak mode.
Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 20 number of smaller particles [6]. This would lead to a change in the skewness spectrum, which would be detected by the proposed methodology as the BBtop. On the other hand, the mode value of Diff_bot found is reasonable, compared with previous existing studies, such as [4]. Figure 13b shows that the Diff_peak histogram presents a similar pattern to Diff_top and Diff_bot, but, as expected, the mode value, more pronounced (corresponding to a more leptokurtic distribution), is centred between the previous two modes, close to 0. Finally, Diff_Tw presents a slightly thicker mode, with a maximum between −50 and 0 m, which is just one class below the Diff_peak mode. Figure 13. (a). Histograms of the differences between Bright Band top and Bright Band bottom heights and the sounding-derived freezing levels Dif_top (blue) and Dif_bot (green), respectively, and (b). histograms of differences between Bright Band peak height and the sounding-derived height of zero wet-bulb temperature, Diff_Tw (red) and differences between Bright Band peak height and the freezing level, Diff_peak (yellow). Values are analysed around ±1.5 h from the sounding launch time.

Sounding Observations vs. Manufacturer ML Levels
In this subsection, the same analysis performed in Section 4.2 is applied to the manufacturer's ML product for 1749 cases detected in the same time period considered above. However, in this case, no BB top nor bottom are considered; only a BB peak is given here by the maximum probability of the ML height product (exceeding 75% as mentioned earlier), denoted as MLmax. We computed the differences between the MLmax and soundingderived zero dry and wet-bulb temperature heights (Diff_peak_Man and Dif_Tw_Man respectively), given by: The distributions of these variables are displayed in Figure 14, similarly to Figure  13b. Both variables show similar patterns and a common main mode corresponding to the same class of height differences (−150 to −100 m) and are relatively wide (three to four different height classes). Secondary modes exceeding 5% relative frequencies are found around −300 and +400 m for Diff_peak_Man and around −500 m for Dif_Tw_Man. Figure 13. (a). Histograms of the differences between Bright Band top and Bright Band bottom heights and the sounding-derived freezing levels Diff_top (blue) and Diff_bot (green), respectively, and (b). histograms of differences between Bright Band peak height and the sounding-derived height of zero wet-bulb temperature, Diff_T w (red) and differences between Bright Band peak height and the freezing level, Diff_peak (yellow). Values are analysed around ±1.5 h from the sounding launch time.

Sounding Observations vs. Manufacturer ML Levels
In this subsection, the same analysis performed in Section 4.2 is applied to the manufacturer's ML product for 1749 cases detected in the same time period considered above. However, in this case, no BB top nor bottom are considered; only a BB peak is given here by the maximum probability of the ML height product (exceeding 75% as mentioned earlier), denoted as ML max . We computed the differences between the ML max and sounding-derived zero dry and wet-bulb temperature heights (Diff_peak_Man and Diff_Tw_Man respectively), given by: The distributions of these variables are displayed in Figure 14, similarly to Figure 13b. Both variables show similar patterns and a common main mode corresponding to the same class of height differences (−150 to −100 m) and are relatively wide (three to four different height classes). Secondary modes exceeding 5% relative frequencies are found around −300 and +400 m for Diff_peak_Man and around −500 m for Diff_Tw_Man.

Discussion
According to the results shown in the previous section, the proposed BB detection and characterisation method provides some advantages.
Firstly, the new dealiasing scheme presents an improvement in some cases where the manufacturer standard processing fails. Despite the fact that the new dealiasing cannot handle very complex cases, such as those found on convective precipitation with intense turbulence, environments favourable to BB conditions, with moderate updrafts, are reasonably well identified, improving the original manufacturer's software capabilities.
Secondly, the BB detection proposed, when compared to radiosounding derived zero dry and wet-bulb temperatures, provide narrower difference height distributions compared to those obtained with the manufacturer's ML product. This suggests that, despite both schemes performing similarly, the proposed methodology gives lower differences compared to the observed freezing level. Moreover, the new method gives explicit information about the BB top and bottom, information not available from the ML manufacturer's product, and the number of detections is considerably higher (5327 vs. 1749 in the period examined).
Despite improvements in the dealiasing approach, limitations of the proposed method include the inability to correct fully folded velocity profiles and also convection with strong windshear. However, these conditions do not typically produce BBs. On the other hand, strong windshear with stratiform precipitation, leading to tilted precipitation streaks, might be a problem for the BB scheme because it examines single radar vertical profiles. Moreover, the proposed scheme cannot handle either multiple BB cases as only the lowest BB can be detected. In any case, the new method based on the vertical variability of the Doppler speed skewness provides a good basis for the further development of more sophisticated BB detection methods.

Summary and Conclusions
The work presents a new methodology to process spectral raw reflectivity data from a K-band vertically pointing Doppler radar, implemented for the Metek MRR-Pro system, and called RaProM-Pro, which is freely available. RaProM-Pro can be used complementary with the manufacturer's software and provides additional features, such as an improved signal and noise detection scheme, an advanced dealiasing method and a new Bright Band product (including top, peak and bottom levels).

Discussion
According to the results shown in the previous section, the proposed BB detection and characterisation method provides some advantages.
Firstly, the new dealiasing scheme presents an improvement in some cases where the manufacturer standard processing fails. Despite the fact that the new dealiasing cannot handle very complex cases, such as those found on convective precipitation with intense turbulence, environments favourable to BB conditions, with moderate updrafts, are reasonably well identified, improving the original manufacturer's software capabilities.
Secondly, the BB detection proposed, when compared to radiosounding derived zero dry and wet-bulb temperatures, provide narrower difference height distributions compared to those obtained with the manufacturer's ML product. This suggests that, despite both schemes performing similarly, the proposed methodology gives lower differences compared to the observed freezing level. Moreover, the new method gives explicit information about the BB top and bottom, information not available from the ML manufacturer's product, and the number of detections is considerably higher (5327 vs. 1749 in the period examined).
Despite improvements in the dealiasing approach, limitations of the proposed method include the inability to correct fully folded velocity profiles and also convection with strong windshear. However, these conditions do not typically produce BBs. On the other hand, strong windshear with stratiform precipitation, leading to tilted precipitation streaks, might be a problem for the BB scheme because it examines single radar vertical profiles. Moreover, the proposed scheme cannot handle either multiple BB cases as only the lowest BB can be detected. In any case, the new method based on the vertical variability of the Doppler speed skewness provides a good basis for the further development of more sophisticated BB detection methods.

Summary and Conclusions
The work presents a new methodology to process spectral raw reflectivity data from a K-band vertically pointing Doppler radar, implemented for the Metek MRR-Pro system, and called RaProM-Pro, which is freely available. RaProM-Pro can be used complementary with the manufacturer's software and provides additional features, such as an improved signal and noise detection scheme, an advanced dealiasing method and a new Bright Band product (including top, peak and bottom levels).
The study illustrates the advantages of RaProM-Pro, such as the capability to detect weaker signals or the robust detection of updrafts thanks to a novel dealiasing scheme. The derived parameters, such as radar reflectivity, mostly match the values provided by the manufacturer (R 2 of 0.992) and is also consistent with independent observations from colocated disdrometer data. The new Bright Band product, based on changes of the skewness of spectral Doppler velocities, compares favourably with the manufacturer's Melting Level product and also with collocated radio sounding observations, both qualitatively in selected examples and quantitatively, as revealed by a study considering 39 days.
Based on the current results, future work is planned to perform a long-term study of BB features in more detail, including BB occurrence, height, thickness and atmospheric conditions (dry and moist BBs).
The methodology can be used for both research and operational applications and could be adapted to other vertically pointing radars. RaProM-Pro is written in Python and is freely available at the GitHub repository.  Data Availability Statement: The proposed methodology coded as a Python program is available at the GitHub repository. Radiosonde data are available from the Meteorological Service of Catalonia (meteo.cat) and MRR-Pro and Parsivel data are available from the authors upon request.
Acknowledgments: Radiosonde data of the Meteorological Service of Catalonia were used in this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This Appendix provides additional information about the new processing methodology regarding noise and signal detection introduced in Section 3.1. Based on 3 h of precipitation data recorded from 12 to 15 UTC 31 March 2020 (Figure 11), 138,240 points (height gates) were examined. Figure A1 shows a mask of three possible cases regarding the signal and noise detection of each method: gates (pixels) with the signal detected by both methods, gates with noise detected by both methods and pixels detected only by one of the methods. The study illustrates the advantages of RaProM-Pro, such as the capability to detect weaker signals or the robust detection of updrafts thanks to a novel dealiasing scheme. The derived parameters, such as radar reflectivity, mostly match the values provided by the manufacturer (R 2 of 0.992) and is also consistent with independent observations from co-located disdrometer data. The new Bright Band product, based on changes of the skewness of spectral Doppler velocities, compares favourably with the manufacturer's Melting Level product and also with collocated radio sounding observations, both qualitatively in selected examples and quantitatively, as revealed by a study considering 39 days.
Based on the current results, future work is planned to perform a long-term study of BB features in more detail, including BB occurrence, height, thickness and atmospheric conditions (dry and moist BBs).
The methodology can be used for both research and operational applications and could be adapted to other vertically pointing radars. RaProM-Pro is written in Python and is freely available at the GitHub repository.  Data Availability Statement: The proposed methodology coded as a Python program is available at the GitHub repository. Radiosonde data are available from the Meteorological Service of Catalonia (meteo.cat) and MRR-Pro and Parsivel data are available from the authors upon request.
Acknowledgments: Radiosonde data of the Meteorological Service of Catalonia were used in this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This Appendix provides additional information about the new processing methodology regarding noise and signal detection introduced in Section 3.1. Based on 3 h of precipitation data recorded from 12 to 15 UTC 31 March 2020 (Figure 11), 138,240 points (height gates) were examined. Figure A1 shows a mask of three possible cases regarding the signal and noise detection of each method: gates (pixels) with the signal detected by both methods, gates with noise detected by both methods and pixels detected only by one of the methods.  It is clear that both methods perform similarly, as most signal and noise detections are identical; only a small fraction, close or at the precipitation contours, is detected differently, being RaProM-Pro a bit more sensitive (about 2.8% more signal detection than the manufacturer's method, as listed in Table A1).

Appendix B
This Appendix provides more details about the dealiasing scheme proposed. The dealiasing scheme is based on the original work by [43,44], and it has been tested for 39 2 h events with precipitation and radiosounding data with the aim to apply it for BB detection. Other more challenging situations, such as convective precipitation with windshear or strong turbulence, where typically BB is not present, may not produce good results. Figure A2 shows one of these cases, recorded on 27 July 2019, displaying the Doppler vertical velocity provided by the manufacturer and the new dealiasing scheme. Figure A3 (analogous to the simpler case shown in Figure 5) illustrates the steps of the dealiasing applied. Two basic concepts are considered in the scheme: Doppler bin clustering (of precipitation and non-precipitation blocks) and vertical continuity of precipitation blocks.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 o It is clear that both methods perform similarly, as most signal and noise detecti are identical; only a small fraction, close or at the precipitation contours, is detected ferently, being RaProM-Pro a bit more sensitive (about 2.8% more signal detection t the manufacturer's method, as listed in Table A1).

Appendix B
This Appendix provides more details about the dealiasing scheme proposed. The dealiasing scheme is based on the original work by [43,44], and it has been tes for 39 2 h events with precipitation and radiosounding data with the aim to apply it BB detection. Other more challenging situations, such as convective precipitation w windshear or strong turbulence, where typically BB is not present, may not produce go results. Figure A2 shows one of these cases, recorded on 27 July 2019, displaying the D pler vertical velocity provided by the manufacturer and the new dealiasing scheme. Figure A3 (analogous to the simpler case shown in Figure 5) illustrates the step the dealiasing applied. Two basic concepts are considered in the scheme: Doppler clustering (of precipitation and non-precipitation blocks) and vertical continuity of p cipitation blocks.   For each height, two Doppler bin groups are identified: blocks of continuous precipitation bins and blocks of non-precipitation bins (hereafter gaps) ( Figure A3a). Then, the spectra are extended to both sides of the Nyquist velocity interval, i.e., considering stronger fall speeds (adding to the right of the original spectra the spectra immediately above the central one) or upward speeds (adding to the left of the original spectra the spectra immediately below the central one). A new grouping of Doppler bins is applied to the extended spectra, providing new precipitation blocks and gaps ( Figure A3b). Now three options are possible to select the dealiased velocity profile, and we assume that only one, for each height, is valid. Starting from the second lowest valid level to higher ones, the selection criteria is that the average velocity of the level considered is closest to the average velocity of the level below ( Figure A3c).
In this convective case, the results seem reasonable for aliasing found from 1500 to ca. 2500 m ASL, but it is not so clear for heights above 2500 m ASL. For each height, two Doppler bin groups are identified: blocks of continuous precipitation bins and blocks of non-precipitation bins (hereafter gaps) ( Figure A3a). Then, the spectra are extended to both sides of the Nyquist velocity interval, i.e., considering stronger fall speeds (adding to the right of the original spectra the spectra immediately above the central one) or upward speeds (adding to the left of the original spectra the spectra immediately below the central one). A new grouping of Doppler bins is applied to the extended spectra, providing new precipitation blocks and gaps ( Figure A3b). Now three options are possible to select the dealiased velocity profile, and we assume that only one, for each height, is valid. Starting from the second lowest valid level to higher ones, the selection criteria is that the average velocity of the level considered is closest to the average velocity of the level below ( Figure A3c).
In this convective case, the results seem reasonable for aliasing found from 1500 to ca. 2500 m ASL, but it is not so clear for heights above 2500 m ASL.