Melting Layer Detection and Observation with the NCAR Airborne W-Band Radar

: A melting layer detection algorithm is developed for the NCAR 94 GHz airborne cloud radar (HIAPER CloudRadar, HCR). The detection method is based on maxima in the linear depolarization ratio and a discontinuity in the radial velocity ﬁeld. A melting layer ﬁeld is added to the radar data, which provides detected, interpolated, and estimated altitudes of the melting layer and the altitude of the 0 ◦ C isotherm detected in model temperature data. The icing level is deﬁned as the lowest melting layer, and the cloud data are ﬂagged as either above (cold) or below (warm) the icing level. Analysis of the detected melting layer shows that the offset between the 0 ◦ C isotherm and the actual melting layer varies with cloud type: in heavy convection sampled in the tropics, the melting layer is found up to 500 m below the 0 ◦ C isotherm, while in shallow clouds, the offset is much smaller or sometimes vanishes completely. A relationship between the offset and the particle fall speed both above and below the melting layer is established. Special phenomena, such as a lowering of the melting layer towards the center of storms or split melting layers, were observed.


Introduction
The mission of the Earth Observing Laboratory (EOL) at the National Center for Atmospheric Research (NCAR) is to develop and deploy observing facilities and provide expertise and data services to aid the scientific community in observational field campaigns. One of the instruments deployed by EOL on the NSF/NCAR High-performance Instrumented Airborne Platform for Environmental Research (HIAPER) Gulfstream V aircraft is the HIAPER Cloud Radar (HCR, [1]), a 94 GHz W-band radar. Since 2015, HCR standard radar data products from three major field campaigns have been released to the scientific community [2]. The HCR team is committed to developing new data products to assist scientists in accomplishing their research goals and to increase the benefit of HCR data to the scientific community. One of these data products is a recently developed estimate of the altitude of the melting layer.
Different physical mechanisms govern the cloud processes above and below the melting layer, and knowledge of its altitude is, therefore, essential for the correct interpretation of cloud radar signals. If the actual melting layer altitude is unknown, temperature data are often used as a proxy. However, because it takes some time for the melting process to take effect once a frozen particle reaches positive temperatures on its downward trajectory, the melting layer is generally found below the 0 • C isotherm. How far below not only depends on the properties of the particles, such as their type or size, and external factors, such as temperature gradient and air motion, but is also determined by cloud physical and dynamical processes, as will be shown in our analysis (Section 4.1). Using 0 • C temperature data instead of the actual melting layer altitude will therefore lead to inaccuracies in data interpretation.
To the careful observer of radar data, the location of the melting layer is often obvious with characteristic signatures in measured parameters, such as a bright band in the reflectivity (DBZ) and linear depolarization ratio (LDR) fields, or an abrupt change in radial velocity (VEL). However, for scientists using an automated as opposed to a purely visual approach to radar data analysis, a quantitative estimate of the altitude of the melting layer from radar data is highly desirable. Therefore, we developed software that detects and flags the melting layer altitude in HCR data.
Using radar data for melting layer detection, especially the polarimetric variables, has proven very successful. Ground-based radars in the S-band [3,4], C-band [5], and K-band [6] were successfully used for bright band and melting layer detection. Bright band detection is also part of the convective/stratiform partitioning algorithms for the space-borne Ka-and Ku-band radars of the Global Precipitation Measurement (GPM) mission [7], and the W-band radar of the CloudSat mission [8].
The primary goal of the proposed melting layer data product is to allow users to distinguish easily between warm and cold clouds below and above the melting layer. It is not aimed at studies of the melting layer itself and is not intended to distinguish between the top and the bottom of the melting layer as is the case with other more specialized algorithms, e.g., [9].
In this study, we give a detailed description of the melting layer detection algorithm and the resulting data product. The algorithm, as developed, was applied to HCR data from three different field campaigns ranging in location from the tropics to the Southern Ocean (see Section 2.1 for details). We analyze the results of the algorithm in a qualitative and quantitative way and compare them to supplementary data sets. Noteworthy observations of special melting layer features are highlighted and discussed.

HCR Data
HCR has been deployed in three major field campaigns. The Cloud Systems Evolution in the Trades (CSET) project took place from 1 July to 12 August 2015, and consisted of 16 research flights (RFs) between California and Hawaii [10]. It was designed to describe and explain the evolution of the boundary layer aerosol, cloud, and thermodynamic structures within the north Pacific trade winds. The CSET HCR data are available at https://doi. org/10.5065/D6CJ8BV7 [11]. The Southern Ocean Clouds, Radiation, Aerosol Transport Experimental Study (SOCRATES) was conducted to further our understanding of clouds, aerosols, air-sea exchanges, and their interactions over the Southern Ocean. The aircraft was stationed in Hobart, Australia, and flew 15 RFs between 15 January and 24 February 2018 [12]. HCR SOCRATES data are available at https://doi.org/10.5065/D68914PH [13]. The third major field campaign for HCR was the Organization of Tropical East Pacific Convection (OTREC) project, which aimed to determine the large-scale environmental factors that control convection over tropical oceans. It took place in the Eastern Pacific and extreme southwestern Caribbean and included 22 RFs between 7 August and 2 October 2019 [14]. OTREC HCR data are available at https://doi.org/10.26023/V9DJ-7T9J-PE0 S [15].
HCR digitizes radar observations with a resolution of~20 m in range and a temporal resolution of 10 Hz, which at typical flight speeds yields a sample volume approximately 20 m long in the direction of flight. In addition to the Doppler radar fields, such as DBZ, VEL, and LDR, reanalysis data from the ERA5 model, which have been interpolated to the HCR time-range grid via linear interpolation in four dimensions (three space and the time dimension), are also provided. For a detailed description of HCR data, data processing, and quality control, see [2]. The variables used in the melting layer detection algorithm are DBZ, LDR, VEL, spectrum width (WIDTH), and the FLAG field, which discriminates between cloud echo and other echo types [2]. The 0 • C isotherm of the interpolated ERA5 temperature field (TEMP) is used as a first estimate of the melting layer altitude and for comparison with the resulting melting layer product. The dry-bulb temperature was chosen as it is most readily available and most commonly used as alternative to melting layer measurements. Analysis of the wet-bulb temperature, which may have better correlation with the melting layer [16], is reserved for future studies.
A special challenge when using high-frequency radars, such as HCR, as opposed to S-band or C-band radars, comes from the fact that they are highly attenuated in liquid precipitation [6]. In HCR data, areas where the radar signal is completely attenuated are identified in the FLAG field. Fortunately, these areas rarely affect the melting layer detection for reasons specific to the fact, that HCR is an airborne radar. When HCR operates in nadirpointing mode at high flight altitudes, the signal only travels through frozen precipitation, which causes very little attenuation [17], until it encounters the melting layer. When the aircraft flies below the melting layer, and HCR is in zenith-pointing mode, the radar is generally in close proximity to the melting layer, as the aircraft needs to remain at a certain altitude above ground. Again, the radar signal exhibits not enough attenuation to be completely attenuated before it encounters the melting layer.

Dropsonde Data
To compare the results of the melting layer detection algorithm with dropsonde data, we use data from the Airborne Vertical Atmospheric Profiling System (AVAPS), which measures high resolution vertical profiles of ambient temperature, pressure, humidity, wind speed, and wind direction. AVAPS was also deployed on the HIAPER aircraft during all three field campaigns.  [21].

Melting Layer Detection
The melting of falling hydrometeors starts at the 0 • C isotherm and continues on their downward trajectory. Therefore, the melting layer is generally found at or below the 0 • C isotherm, hence the first step in the detection algorithm is to search for the 0 • C isotherm in the reanalysis TEMP field. It is important to recognize that there can be several 0 • C isotherms staggered on top of each other in cases of temperature inversions, and indeed several such cases are observed in the SOCRATES data set (e.g., at 00:10 UTC in Figure 1a). Therefore, in each ray, we first identify all pairs of gates where a change in the temperature sign is observed. We then work along the time dimension to connect neighboring 0 • C isotherm gates horizontally to form continuous data periods (see Figure 1a for an example). If several 0 • C isotherm data periods are staggered on top of each other, each one is processed separately in the steps that follow.
To narrow the search area, we make an initial approximation of where we expect the melting layer to be by simply applying an offset to the ERA5 0 • C isotherm, so that the approximate melting layer altitude is located just below the 0 • C isotherm. The initial offset value was~200 m for SOCRATES and CSET and~400 m for OTREC. The exact numbers are not important as they are adjusted throughout the processing steps. For more details on the offset between the model 0 • C isotherm and the melting layer, see Section 4.2.
The best HCR fields for our purpose are LDR, which shows a maximum at the melting layer (Figure 1b), and VEL, which shows a discontinuity as hydrometeors increase their fall speed when changing from the frozen to the liquid phase ( Figure 1c). These two fields stand out as prime candidates because they both show a strong signature at the melting layer and have already been used by previous studies, e.g., [9,22]. Adding the DBZ field, which sometimes shows a bright band signature at the melting layer altitude, did not yield additional benefits. Both LDR and VEL are generally more sensitive to the melting of hydrometeors, often showing characteristic signatures where no bright band is observed in DBZ. Reflectivity was therefore not used in the algorithm.  To narrow the search area, we make an initial approximation of where we expect the melting layer to be by simply applying an offset to the ERA5 0 °C isotherm, so that the approximate melting layer altitude is located just below the 0 °C isotherm. The initial offset value was ~200 m for SOCRATES and CSET and ~400 m for OTREC. The exact numbers are not important as they are adjusted throughout the processing steps. For more details on the offset between the model 0 °C isotherm and the melting layer, see Section 4.2.
The best HCR fields for our purpose are LDR, which shows a maximum at the melting layer (Figure 1b), and VEL, which shows a discontinuity as hydrometeors increase their fall speed when changing from the frozen to the liquid phase ( Figure 1c). These two fields stand out as prime candidates because they both show a strong signature at the melting layer and have already been used by previous studies, e.g., [9,22]. Adding the DBZ field, which sometimes shows a bright band signature at the melting layer altitude, did not yield additional benefits. Both LDR and VEL are generally more sensitive to the We need to be careful to remove artifacts that could lead to erroneous results, i.e., false detections. The biggest suspect in that regard is back-lobe echo, which occurs when HCR is pointing at the zenith, and the aircraft is flying near the ocean or land surface. In such cases, the back lobe of the radar is occasionally reflected off the surface and detected by the receiver. Back-lobe echo is typically characterized by a band of low reflectivity, highly variable radial velocity, and high spectrum width. It appears at a radar range equal to twice the height of the aircraft above the surface, and it recedes and approaches in range as the aircraft ascends or descends. A mechanism to detect the back-lobe echo was developed by [2], but because the LDR of the back-lobe echo strongly resembles the characteristics of the melting layer, the thresholds used by [2] were tightened to remove the echo more aggressively. We reject data where DBZ is less than −18 dBZ (as opposed to −20 dBZ in [2]) and WIDTH is greater than 1 m s −1 (1.4 m s −1 ). These new thresholds successfully identified the majority of the back-lobe echo, except for a few cases which will be discussed in Section 3.5. In addition to removing the back-lobe echo, we remove all data that is identified as something other than cloud echo in the FLAG field, such as ocean or land Remote Sens. 2021, 13, 1660 5 of 17 surface echo, radar burst echo, etc. (see [2] for details). We also set thresholds on the LDR data itself and use only data between −16 and −7 dB. Observations that are close to the aircraft also proved to be unsuitable, and we removed all data within 150 m of the aircraft in both the LDR and VEL fields.
After applying the quality checks described above, we search for a melting layer signature in the LDR field. We constrain the search to within +/− 350 m, in range, from the approximate altitude of the 0 • C isotherm. The maximum LDR values detected within this range indicate possible melting layer altitudes. Several checks of their validity are performed. First, we calculate the 30 s moving average and moving standard deviation in the time dimension. Detections for which the standard deviation exceeds 100 m are removed because we do not expect the melting layer to vary significantly in altitude over short distances. Then, we calculate the difference between the originally detected altitudes and the moving average and remove those outliers where the difference exceeds 50 m. Only those detections that pass these tests are retained as valid.
The next step is to adjust the model-approximated 0 • C isotherms up or down to match the results obtained from the LDR signature. We then search for a discontinuity in the VEL field, within 960 m of the 0 • C isotherm altitudes that were adjusted using LDR. We smooth the VEL data in range using a 10-gate moving average, to reduce the likelihood of identifying minor discontinuities. We then use the ischange MATLAB function (https://www.mathworks.com/help/matlab/ref/ischange.html, accessed 28 April 2020), which is based on [23], to identify the location of an abrupt change, i.e., a discontinuity in VEL in the range direction. We interpret the altitude of the detected discontinuity as a proxy for the melting layer altitude. The ischange function also provides the mean and variance of the velocity above and below the melting layer. We use those in later processing steps.
An example of the detected velocity discontinuity is shown as the black line in Figure 2. Once again, we verify these detections using a number of checks. A 100 s moving average is computed in the time dimension. If the difference between the unsmoothed detections and the moving average exceeds 100 m, the detection is discarded (purple dots in Figure 2b). Similarly, the results are discarded if the 30 s moving standard deviation exceeds 35 m (green dots). Results are also discarded if the difference between the LDR result and VEL result exceeds 200 m (yellow dots). We also discard times where the mean VEL above or below the melting layer indicates upward hydrometeor motion (red dots) because this upward motion is only expected in updrafts that inhibit the formation of a readily identifiable melting layer. In general, liquid hydrometeors below the melting layer exhibit larger downward motion than frozen particles above the melting layer. Therefore, it is reasonable to assume that the difference between the mean velocity above the melting layer and the mean velocity below the melting layer must exceed a certain number, and indeed it turns out that setting a threshold of 0.7 m s −1 is a good quality control step (not shown in Figure 2b because no such cases were found in this specific example).
For the final melting layer altitude, we combine the LDR-and VEL-based estimates. If the melting layer is detected in both the LDR and VEL fields, we use the LDR-based result. Remaining outliers are removed by calculating the difference between the detected altitudes and a 3 s moving average in time and excluding cases where the absolute difference exceeds 50 m. A 2 s moving average filter is applied. As a final quality control check, we exclude results that lie more than 100 m above the model-derived 0 • C isotherm as the model has proven to be very accurate (see Section 4.2). It is important to note that the melting layer altitudes detected in VEL are sometimes somewhat lower than those detected in LDR, as has been observed in previous studies, e.g., [22], but the difference was not significant, i.e., less than the overall noise, and therefore merging them into the same data field seemed justified.
it is reasonable to assume that the difference between the mean velocity above the melting layer and the mean velocity below the melting layer must exceed a certain number, and indeed it turns out that setting a threshold of 0.7 m s −1 is a good quality control step (not shown in Figure 2b because no such cases were found in this specific example). For the final melting layer altitude, we combine the LDR-and VEL-based estimates. If the melting layer is detected in both the LDR and VEL fields, we use the LDR-based result. Remaining outliers are removed by calculating the difference between the detected altitudes and a 3 s moving average in time and excluding cases where the absolute difference exceeds 50 m. A 2 s moving average filter is applied. As a final quality control check, we exclude results that lie more than 100 m above the model-derived 0 °C isotherm as the model has proven to be very accurate (see Section 4.2). It is important to note that the melting layer altitudes detected in VEL are sometimes somewhat lower than those detected in LDR, as has been observed in previous studies, e.g., [22], but the difference was not significant, i.e., less than the overall noise, and therefore merging them into the same data field seemed justified.
It is important to point out that all the constants and thresholds used in the processing steps were empirically chosen to work for HCR 10 Hz data and are not intended to translate to other radar data without adjustments. In fact, when we tried to process HCR 2 Hz It is important to point out that all the constants and thresholds used in the processing steps were empirically chosen to work for HCR 10 Hz data and are not intended to translate to other radar data without adjustments. In fact, when we tried to process HCR 2 Hz data with the same algorithm, the results were not satisfactory. Instead, we decided to run the algorithm at 10 Hz and interpolate the results to 2 Hz.

Melting Layer Interpolation in Time
It is desirable to not only provide the users with the melting layer altitude at time steps when a melting layer was identified in the LDR or VEL fields by the algorithm, but also to give an indication of its altitude when either no melting layer signature is detected or no clouds are observed at all. Of course, one may ask why we even provide melting layer data where no clouds are observed. We can envision several scenarios where an estimated melting layer might prove useful even outside of cloud areas. For example, someone may want to combine HCR data with satellite observations and the HCR melting layer estimate could be used for a satellite detected cloud in close proximity. We would also like to point out that NCAR/EOL provides a combined data set that includes both HCR and High Spectral Resolution Lidar (HSRL) data on the same grid. Having an estimate of the melting layer altitude may prove useful in areas where HSRL observes clouds but HCR does not (because of the different capabilities of the instruments).
To fill in the missing time periods, we use different methods depending on the length of the missing time interval. When a gap between two detected data points is shorter than 20 min, a simple linear interpolation between the detected data points is performed, and we retain these interpolated points as a high-quality proxy for the actual melting layer altitude. If the gap exceeds 20 min, we estimate the melting layer altitude by applying a constant offset to the 0 • C isotherm so that the estimated altitudes follow the 0 • C isotherm. To avoid jumps where the estimated altitudes join the detected altitudes, we define a 10 min transition period where a linear interpolation between the detected and estimated altitudes is applied. Even though an interpolation is applied in the transition period, it is still considered as part of the estimated data and not to be confused with the interpolated points that fill the short gaps described above. The offset used for the estimated melting layer is calculated in different ways for each field campaign. For SOCRATES and OTREC, detected melting layer points were obtained in each flight. Therefore, the offset is calculated as the mean offset over the respective flight. During CSET, cloud tops rarely exceeded the altitude of the melting layer and, therefore, the melting layer could only be detected in four flights. Calculating the offset on a flight by flight basis was therefore not possible and we chose to use the mean offset calculated from the whole project (i.e., all four flights). A detailed analysis of the offsets is presented in Section 4.2.

The Icing Level
For many applications, the actual altitude of the melting layer is not as important as the question of whether a specific point in a cloud is located above or below the melting layer, i.e., if it is part of a cold or warm cloud. Therefore, we would like to flag each 2D (time range) data point as either warm or cold. However, before we can do that, we need to decide how to handle situations where several melting layers are staggered on top of each other (see Section 3.5). We clarify the meaning of "warm" and "cold" flags by introducing the icing level. The American Meteorological Society [24] defines the icing level as the "Lowest height above sea level at which an aircraft in flight may encounter icing" which we interpret to coincide with the lowest melting layer found or estimated in HCR data.

The Melting Layer Product
The melting layer product comprises three temporal categories (Section 3.2): (1) detections: times for which a melting layer altitude was actually detected in the LDR or VEL field (dark blue in Figure 3); (2) interpolations: times with interpolated altitudes bridging shorter gaps (light blue); (3) estimates: times with altitudes estimated from the 0 • C isotherm with an appropriate offset (green). We flag each separately. We also indicate whether the melting layer coincides with the icing level (i.e., whether it is the lowest melting layer) or if it is above the icing level, and we flag the icing level itself (red). In addition, we flag the 0 • C isotherm for the users' convenience (black), again indicating whether it is above or below the icing level. Data points that are not flagged as a melting layer point are flagged as being either in the warm (below the icing level, pink in Figure 3b) or cold (above the icing level, yellow in Figure 3b) region. The actual flag values are stored in the metadata of the distributed CfRadial NetCDF files. The altitude of the icing level is provided as a 1D (time dimension) data field.

Product Quality Assessment
When developing algorithms for product creation, it is often relatively easy to write code that gives the desired results for a few selected cases. The real challenge is to develop an algorithm that is robust enough to handle at least the majority of highly diverse conditions. Tuning an algorithm for conditions as different as, for example, OTREC and SOCRATES, inevitably requires some compromise and finding the right balance between missed detections and false positives. We decided to err on the side of missing some detections because we believe that the interpolations, and, to a lesser degree, the estimates, can serve as a high-quality substitute when detections are missed. An example in which short time periods are easily assessed visually but missed by the algorithm is shown in Figure 3 at 17:28 UTC.
shorter gaps (light blue); (3) estimates: times with altitudes estimated from the 0 °C isotherm with an appropriate offset (green). We flag each separately. We also indicate whether the melting layer coincides with the icing level (i.e., whether it is the lowest melting layer) or if it is above the icing level, and we flag the icing level itself (red). In addition, we flag the 0 °C isotherm for the users' convenience (black), again indicating whether it is above or below the icing level. Data points that are not flagged as a melting layer point are flagged as being either in the warm (below the icing level, pink in Figure 3b) or cold (above the icing level, yellow in Figure 3b) region. The actual flag values are stored in the metadata of the distributed CfRadial NetCDF files. The altitude of the icing level is provided as a 1D (time dimension) data field.

Product Quality Assessment
When developing algorithms for product creation, it is often relatively easy to write code that gives the desired results for a few selected cases. The real challenge is to develop an algorithm that is robust enough to handle at least the majority of highly diverse conditions. Tuning an algorithm for conditions as different as, for example, OTREC and SOC-RATES, inevitably requires some compromise and finding the right balance between missed detections and false positives. We decided to err on the side of missing some detections because we believe that the interpolations, and, to a lesser degree, the estimates, can serve as a high-quality substitute when detections are missed. An example in which short time periods are easily assessed visually but missed by the algorithm is shown in Figure 3 at 17:28 UTC. As the melting layer is relatively easy to identify visually in the radar fields, the most straightforward way to evaluate the accuracy of the algorithm is to examine the results visually. To do this, we plotted the results in hourly increments, which were then carefully screened for both missed events and false detections. From the visual inspection, we concluded that, overall, the algorithm performs very well.
Even if the melting layer is detected in the correct gates, in principle there could still be an error in the melting layer altitude if the gate altitude is incorrect. HCR gate altitudes are derived from a global positioning system (GPS) and inertial navigation system (INS) unit, which is mounted in the nose of the radar pod. The GPS/INS altitudes are corrected for gravitational effects on the Earth's surface using the Earth Gravitation Model (EGM2008, see https://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/, accessed 8 December 2019). The altitudes can easily be verified by comparing the gate altitude of the ocean's surface with the GPS/INS altitude. Such a comparison was carried out for all three field campaigns (not shown) and shows agreement within +/− 20 m, which also happens to be the gate spacing of HCR. We therefore conclude that accurately detected melting layer altitudes have an uncertainty of~20 m.
The algorithm was originally developed for the OTREC data set, which has many suitable clouds with clear melting layer structures. Generalizing the algorithm to work well for CSET but especially SOCRATES proved more complicated. During SOCRATES, flight plans were designed to zigzag vertically through the near-surface cloud layer, which often coincided with the melting layer altitude, leading to multiple melting layer crossings by the aircraft. During these zigzag maneuvers, HCR pointing regularly switched between zenith and nadir, which led to additional difficulties for the algorithm. An example of a situation where the identification of the icing level is not straightforward is shown in Figure 4a. Another challenge came from the fact that during the southernmost parts of the SOCRATES flights, air temperature often dropped below 0 • C at the surface so the intersection of the melting layer with the surface needed to be handled appropriately (Figure 4b). During SOCRATES, we found several cases in which the back-lobe echo was not properly filtered and removed (Section 3.1), leading to false positive detections ( Figure  4c,d). Fortunately, these cases proved rare, and were unique to SOCRATES. We must also point out that the algorithm was developed for zenith and nadir-pointing operations only and is, therefore, unreliable during sea surface calibration events during which the radar is scanning 20° off nadir perpendicular to the aircraft for several minutes for reflectivity calibration monitoring purposes [2]. As sea surface calibration events are generally performed in clear air conditions, this problem has little relevance. In addition, the algorithm is designed to detect a single melting layer in the vicinity of the model 0 °C isotherm and, therefore, does not properly handle cases in which the melting layer splits into two (see Section 4.1 below). In these situations, it tends to jump from one to the other, or misses both.

Special Observations
It is expected but nonetheless interesting to observe how strongly the melting layer altitude varies, even within an individual cloud. This is especially true in regions of strong convection where areas with strong melting layer signatures are interrupted by updrafts. An example of such a region is shown in Figure 5a  During SOCRATES, we found several cases in which the back-lobe echo was not properly filtered and removed (Section 3.1), leading to false positive detections (Figure 4c,d). Fortunately, these cases proved rare, and were unique to SOCRATES. We must also point out that the algorithm was developed for zenith and nadir-pointing operations only and is, therefore, unreliable during sea surface calibration events during which the radar is scanning 20 • off nadir perpendicular to the aircraft for several minutes for reflectivity calibration monitoring purposes [2]. As sea surface calibration events are generally performed in clear air conditions, this problem has little relevance. In addition, the algorithm is designed to detect a single melting layer in the vicinity of the model 0 • C isotherm and, therefore, does not properly handle cases in which the melting layer splits into two (see Section 4.1 below). In these situations, it tends to jump from one to the other, or misses both.

Special Observations
It is expected but nonetheless interesting to observe how strongly the melting layer altitude varies, even within an individual cloud. This is especially true in regions of strong convection where areas with strong melting layer signatures are interrupted by updrafts. An example of such a region is shown in Figure 5a, where the melting layer on one side of an updraft region (~12:27 UTC) is 200 m higher than on the other side of the updraft region (~12:29 UTC). Even in regions where the melting layer is not interrupted by an updraft, we sometimes observe strong variability in the melting layer altitude (Figure 5b) or even jumps (just after 16:18 UTC, Figure 5c) as the melting layer is modified by storm internal processes.
We also observe several cases where the melting layer altitude is lower in the center of a cloud feature (01:55 UTC, Figure 6a; 13:53 UTC and 14:10 UTC, Figure 6b) in areas of increased fall speeds. We speculate that the melting layer is either pushed lower in the center of a strong downdraft as cold air is transported downward or it only appears lower as larger particles with higher fall speeds take longer to melt and therefore stay frozen to lower altitudes. An analysis of the connection between particle fall speed and vertical offset between the 0 • C isotherm and the melting layer is presented in Section 4.2.
Of special interest are regions with what we refer to as split melting layers, where two vertically staggered melting layer signatures are observed in the VEL and LDR fields in very close proximity to each other (Figure 7). Note that this phenomenon indicates a splitting of one distinct melting layer and is not to be confused with cases of several staggered melting layers caused by temperature inversions (Section 3.5). The split melting layer phenomenon was also recently observed and described by [22] and was attributed to the melting of two distinct populations of ice particles. The fact that we observed split melting layers over the Southern Ocean during the SOCRATES field campaign (Figure 7) confirms that the authors of [22] were correct to conclude that "more reports of this phenomenon are expected in the future".

Offset between the Melting Layer and the 0 °C Isotherm
In Section 4.1, we have shown that the altitude of the melting layer not only depends on the properties of the particles, such as their type or size, and external factors, such as temperature gradient and air motion, but is also determined by cloud physical and dynamical processes.
One factor that seems of particular importance in determining the vertical offset between the 0 °C isotherm and the melting layer is the fall speed of the particles. To investigate its role, the mean radial velocity of the ten gates above and ten gates below each melting layer detection is computed and compared with the vertical offset at that specific time (Figure 8). Both the velocity above (Figure 8a) and below the melting layer ( Figure  8b) show a clear correlation with the offset. A geometric mean regression [25] is performed

Offset between the Melting Layer and the 0 • C Isotherm
In Section 4.1, we have shown that the altitude of the melting layer not only depends on the properties of the particles, such as their type or size, and external factors, such as temperature gradient and air motion, but is also determined by cloud physical and dynamical processes.
One factor that seems of particular importance in determining the vertical offset between the 0 • C isotherm and the melting layer is the fall speed of the particles. To investigate its role, the mean radial velocity of the ten gates above and ten gates below each melting layer detection is computed and compared with the vertical offset at that specific time (Figure 8). Both the velocity above ( Figure 8a) and below the melting layer ( Figure 8b) show a clear correlation with the offset. A geometric mean regression [25] is performed between the offset and velocity points, and the resulting regression lines and coefficients are shown in Figure 8. The slope of the regression line is steeper for velocities above the melting layer as fall speeds are smaller in the frozen region. The correlation between velocity and offset suggests that velocity either above or below the melting layer can potentially be used to estimate the altitude of the melting layer from the 0 • C isotherm in datasets where no melting layer altitude is provided. Due to the many factors influencing the offset between the 0 °C isotherm and the melting layer, especially the fall speed of the particles, it is not surprising that the distance between the detected melting layer and the 0 °C isotherm varies between field campaigns, between individual flights of a specific field campaign, and also within each flight and indeed within each cloud. These variations are generally on the order of a few hundred meters. Mean offsets between the 0 °C isotherm, which is generally very accurate (see below), and the melting layer per flight and per field campaign are shown in Table 1.  Due to the many factors influencing the offset between the 0 • C isotherm and the melting layer, especially the fall speed of the particles, it is not surprising that the distance between the detected melting layer and the 0 • C isotherm varies between field campaigns, between individual flights of a specific field campaign, and also within each flight and indeed within each cloud. These variations are generally on the order of a few hundred meters. Mean offsets between the 0 • C isotherm, which is generally very accurate (see below), and the melting layer per flight and per field campaign are shown in Table 1. The offset in OTREC is, on average, 323 m and twice as big as the 156 m offset observed in SOCRATES. We speculate that this difference is the result of the very different types of cloud observed during these two field campaigns. The tropical convective systems sampled during OTREC likely contain larger particles with higher fall speeds than the shallow cloud layers observed in SOCRATES. Clouds sampled in CSET were also mostly shallow in nature, which may explain the smaller offset of 170 m; however, because of the limited number of detections and the large variability, we cannot draw definitive conclusions for CSET. It is remarkable that individual flights in OTREC show offsets of up to, and in one case (RF10) even above, 500 m, which underlines the importance of using the melting layer, as opposed to the model temperature, to delineate warm from cold clouds.
The offsets listed in Table 1 and the conclusions derived from them depend on the accuracy of the ERA5 temperature data, specifically the altitude of the 0 • C isotherm. To assess the quality of the 0 • C isotherm model data, we compare them with temperature data from dropsondes that were launched from HIAPER during all three field campaigns. For each dropsonde profile, we find all altitudes where a change in sign is reported, identify the one that is closest to the 0 • C isotherm in the ERA5 data, and calculate the difference between the two. It is important to use the 0 • C point that is closest to the melting layer because the dropsonde data shows even more structure in the temperature profile than the model data, as expected ( Figure 9).
Overall, the comparisons for all three field campaigns show excellent agreement between the dropsonde and ERA5 0 • C levels ( Figure 10). The mean altitude difference between the two is 78 m for CSET, 32 m for SOCRATES, and 27 m for OTREC. The better agreement for SOCRATES and OTREC is likely due to the fact that the dropsonde data was assimilated into the ERA5 model for these two field campaigns, but not for CSET. It is interesting to note that in all three field campaigns the dropsonde 0 • C isotherm was observed slightly above the ERA5 0 • C isotherm, which means that the model is generally reporting slightly lower temperatures than the dropsondes. From the slightly negative bias in the model, one would expect the actual offset between the melting layer and the 0 • C isotherm to be slightly larger than that reported in Table 1. To confirm this assumption, we calculate the offset using dropsonde data. We distinguish between offsets calculated with detected, interpolated, and estimated melting layer altitudes. Indeed, the offset calculated from dropsonde data is, at just more than 200 m, slightly larger than the 156 m calculated from ERA5 temperatures (Table 1) for SOCRATES. For OTREC, which has not only the smallest temperature bias between the dropsondes and the ERA5 data but also the largest variability in the offset calculation, there is no significant difference between the offsets calculated with dropsonde vs. ERA5 temperatures.
Remote Sens. 2021, 13, x FOR PEER REVIEW Figure 9. Example comparing melting layer (colors as in Figure 5) and ERA5 temperature dropsonde temperatures. Vertical columns show dropsonde temperatures with black dots the vertical columns indicating 0 °C data points. Background fields are (a) radial velocity a ERA5 temperatures. Times are in UTC.

Conclusions
We developed a melting layer detection algorithm for HCR data. The algorithm use the LDR and VEL fields from the radar. The final product combines confirmed detection values interpolated in time between detection events, and estimated values where no de tections are available in close time proximity. ERA5 0 °C isotherm values are provided a well. The icing level is defined as the lowest melting layer in each vertical ray, and dat points are flagged as above or below the icing level so that it is easy for the users to deter mine whether a specific location is in the warm or cold region of the atmosphere. Th algorithm was applied to observations from three fundamentally different field cam paigns; the results were carefully checked for accuracy and found to be very reliable. A the algorithm shows excellent capabilities in both tropical regions and at high latitude

Conclusions
We developed a melting layer detection algorithm for HCR data. The algorithm uses the LDR and VEL fields from the radar. The final product combines confirmed detections, values interpolated in time between detection events, and estimated values where no detections are available in close time proximity. ERA5 0 • C isotherm values are provided as well. The icing level is defined as the lowest melting layer in each vertical ray, and data points are flagged as above or below the icing level so that it is easy for the users to determine whether a specific location is in the warm or cold region of the atmosphere. The algorithm was applied to observations from three fundamentally different field campaigns; the results were carefully checked for accuracy and found to be very reliable. As the algorithm shows excellent capabilities in both tropical regions and at high latitudes, we are confident that it will translate to other regions of the world without major adjustments.
The melting layer is highly variable on very short scales and strongly modified by storm internal processes. We found cases with sudden jumps in the melting layer altitude, wave-like structures within single storms, and strong altitude differences on opposite sides of updraft regions. The melting layer sometimes has a tendency to be lower towards the middle of cloud features with strong downdrafts. In some cases, the melting layer splits into two horizontal layers.
The magnitude of the downward offset of the melting layer from the model 0 • C isotherm depends on the type of storm sampled. In tropical storms, as sampled in the OTREC field campaign, it can exceed 500 m, while in regions with shallow clouds, such as those observed during SOCRATES and CSET, it tends to be less than a few hundred meters. The offset calculations were verified with dropsonde data, which show excellent agreement with the ERA5 model. A relationship between the offset and the particle fall speed both above and below the melting layer could be established. As expected, velocities in the frozen region above the melting layer were generally smaller than those in the liquid region below the melting layer, which resulted in a steeper regression line in the frozen compared to the liquid region.