Radar-Based Automatic Identification and Quantification of Weak Echo Regions for Hail Nowcasting

Abstract

The identification of some radar reflectivity signatures plays a vital role in severe thunderstorm nowcasting. A weak echo region is one of the signatures that could indicate updraft, which is a fundamental condition for hail production. However, this signature is underutilized in automatic forecasting systems due to the lack of a reliable detection method and the uncertain relationships between different weak echo regions and hail-producing thunderstorms. In this paper, three algorithms related to weak echo regions are proposed. The first is a quasi-real-time weak echo region morphology identification algorithm using the radar echo bottom height image. The second is an automatic vertical cross-section-making algorithm. It provides a convenient tool for automatically determining the location of a vertical cross-section that exhibits a visible weak echo region to help forecasters assess the vertical structures of thunderstorms with less time consumption. The last is a weak echo region quantification algorithm mainly used for hail nowcasting. It could generate a parameter describing the scale of a weak echo region to distinguish hail and no-hail thunderstorms. Evaluation with real data of the Tianjin radar indicates that the critical success index of the weak echo region identification algorithm is 0.61. Statistics on these data also show that when the weak echo region parameters generated by the quantification algorithm are in a particular range, more than 85% of the convective cells produced hail.

1. Introduction

Severe thunderstorms are the most frequent weather threats to life and property in North China [1,2]. Since thunderstorms have small spatial scales and rapidly evolve, their monitoring and nowcasting commonly rely on weather surveillance radars. Weather radars can capture three-dimensional thunderstorm structure information by volume scanning with a high time and space resolution to determine their severity and development stage.

Upward-moving air in a thunderstorm, which is known as the updraft, is a critical factor in assessing thunderstorm potential [3]. In radar-based severe thunderstorm nowcasting, the identification of updrafts can help forecast some high impact events [4,5,6]. Since the vertical component of radar radial velocity is very small, it is impossible to measure vertical drafts directly. Therefore, the observation of updrafts in single-polarization radars relies on specific morphology signatures of vertical reflectivity, such as the echo overhang, weak echo region, and bounded weak echo region. For dual-polarization radars, differential reflectivity and differential phase are also good indicators of updrafts [7,8,9].

The weak echo region (WER) is a region with low radar reflectivity that is bounded on one side and above by strong echoes [10,11]. It is usually located on the low-altitude inflow side of a thunderstorm. The WER is caused by a strong updraft that carries precipitation particles to the middle levels before they grow to a radar-detectable size [12]. The strong tilted echo above WER is called overhang echo. In Figure 1, the schematic diagram of WER formation is shown. When both sides of a WER are bounded, it is a bounded weak echo region (BWER, also known as an echo-free vault). A BWER is considered to be a sign of a particularly severe thunderstorm.

The formation and growth of hailstones is also closely related to updrafts [14,15,16]. When updrafts carry hydrometeors up above the melting layer, small frozen particles become hail embryos and grow as additional supercooled liquid water freezes onto them. On the condition that updraft strength is sufficient to support large hailstones, and not too strong to eject them out of the growth region, hailstones would reside in the updraft for a long time and become sizeable [10,15,16,17,18,19]. Previous studies also found that broader updrafts in severe thunderstorms are conducive to the growth of hail [19,20,21].

Although updraft-related signatures, such as WERs and BWERs, are helpful to identify hail-producing thunderstorms, they are rarely used in operational forecasting systems due to the lack of a reliable identification algorithm. Most existing identification algorithms aim at BWERs, such as the work of [22,23,24,25]. As part of their storm-labeling system to identify convective updrafts, the authors in [9] proposed a WER identification algorithm. This algorithm searches for WERs by finding if there are at least six of the eight horizontally-adjacent grid volumes with horizontal reflectivity gradients greater than or equal to 8 $\mathrm{dB}Z·\mathrm{k}{\mathrm{m}}^{-1}$ and column maximum reflectivity greater than or equal to 40 $\mathrm{dB}Z$. Obviously, it only considers the interface of the WER and the overhang echo, while ignoring whether there are tilted wall echoes that may lead to false alarms.

Another obstacle to utilize WERs for automatic hail forecasting is uncertainties in the presence of WERs and hailfall. The presence of a WER in a thunderstorm does not necessarily indicate that the updraft has sufficient strength for hail growth [19,26]. In manual forecasting, forecasters can make decisions through comprehensive analysis of WERs based on their experience. In automatic forecasting, categorical values or quantitative values are needed for machine decision making. There are various radar parameters for hail detection. Some of the parameters can be used independently by setting a warning threshold, such as the Waldvogel parameter [27], vertically-integrated liquid water content [28], vertically-integrated liquid water density [29], severe hail index, and probability of severe hail [30].

Some studies proposed various multiparameter approaches. In [31], a simple combine criterion was defined by requiring that at least one of the three parameters exceeds the threshold. In [32], several forecasting variables were combined into one parameter using linear discriminant analysis. In [32,33], logistic regression was used to develop a relationship between radar parameters and hail occurrence. In [34], new combined parameters were defined with the help of principal component analysis. In [35], the performance of mesocyclonic hailstorm detection was improved by combining mesocyclone detection results with radar hail parameters. In [36], hail and no-hail thunderstorms were distinguished by using life-cycle analysis of multiple new radar parameters. In [37,38], semi-supervised clustering and a support vector machine were applied to classify polarimetric radar hydrometeors, respectively. A parameter used in multiparameter models is not necessarily capable of independently distinguishing all samples. However, it may improve overall performance by better distinguishing certain samples.

As a 3D signature, WERs cannot be seen in plan views, but can be exhibited in vertical cross-sections. In retrospective studies and operational hail-suppression programs, vertical cross-sections could help people obtain more information about thunderstorms and make better decisions [39]. However, manually determining the position of a cross-section is time consuming since the information of multiple radar elevations should be considered at the same time.

Based on the above background, we conducted a series of studies around WERs to make better use of them in severe thunderstorm forecasting. The first is the WER detection algorithm, which utilizes their morphological features. It is also the basis of the other algorithms. The second is the automatic vertical cross-section-making algorithm. It can be used as a convenient visualization tool for forecasters in diverse applications. The last is the WER quantification algorithm used for hail detection, which can generate a parameter for each WER. The generated parameter cannot be used to independently detect hail, but it can make better decisions in some cases, so it may improve overall performance in multiparameter models.

The paper is structured as follows: The datasets and study area are specified in Section 2. Section 3 gives a detailed description of the proposed methods, including the WER detection algorithm, the automatic cross-section-making algorithm, and the weak echo region quantification algorithm. The results of validations using real weather radar data and two case examples are provided in Section 4. Finally, Section 5 briefly draws the most important conclusions.

2. Study Area and Data

Our study area was focused on the Beijing-Tianjin-Hebei region in North China. The Beijing-Tianjin-Hebei region is located in the mid-latitude coastal zone, with mountains in the north and west and the Bohai Sea in the east. The terrain is high in the northwest and low in the southeast. This particular terrain leads to frequent thunderstorms [40,41]. Thunderstorm season in this area begins in May and goes until September [42,43]. A general view of the study area is shown in Figure 2a.

Data used in this paper included severe weather event observational data, Doppler weather radar data, and radiosounding data. All radar data were from the single-polarization S-band Tianjin radar. The interval of the radar volume scans was 6 min. Each volume scan included 9 elevations, and the lowest and highest elevation angles were ${0.5}^{\circ }$ and ${19.5}^{\circ }$, respectively. The altitude of the radar was 70.3 m, and its scan range was 230 km. Therefore, at the far end of the radar beam, the minimum detectable height was 2.077 km. All algorithms in this paper were based on the radar plan position indicator (PPI) image, which was generated by bilinear interpolation with a resolution of 1 × 1 km. Radar data were quality controlled, including eliminations of ground clutter [44] and superrefraction [45].

Hail case data were extracted from the severe weather event observational database provided by the China Meteorological Administration. In the database, hail cases are from hail reports of 175 manual observation stations in the study area, which are based on human eye-observations of hailstones of any hail size. The location of the Tianjin radar and the geographical distribution of the manual observation stations are shown in Figure 2b. Data were quality controlled, and only hail records with detailed start and end times were used in this study.

Some analyses needed the heights of the melting layer and the −20 degree layer even though they were not used in the algorithms. Heights were obtained from the products of the Meteorological Information Combine Analysis and Process System (MICAPS) provided by the China Meteorological Administration. The raw files were grid-point field data with a resolution of ${2.5}^{\circ }\times {2.5}^{\circ }$, which was estimated from radiosounding data. The heights of a specific location were calculated by bilinear interpolation in our system.

In this paper, the algorithms were performed on each convective cell. The convective cell identification algorithm used in this paper was based on Storm Cell Identification and Tracking (SCIT) algorithm [46], but some modifications were made to identify cells in PPI images instead of radial images. SCIT identifies 1D storm segments in radial reflectivity data and then combines individual storm segments into 2D storm components based on spatial proximity. We substituted these two steps by a border-following algorithm proposed by [47]. After binarizing the image with different reflectivity thresholds, the border-following algorithm can directly extract all the 2D connected components, which is similar to the 2D storm components of SCIT. Finally, 3D convective cells were obtained by using the same method as SCIT. Since we mainly focused on severe thunderstorms, in order to relieve computing pressure, the cells with a maximum reflectivity of less than 40 $\mathrm{dB}Z$ were excluded from the study.

The dataset of convective cells for evaluation was divided into two subsets, one for the detection algorithm and the other for the quantification algorithm. The subset for the detection algorithm contained the cells from May–September 2015, which were manually labeled when there were WERs inside. The method of labeling is described in Section 4.1. The subset for the quantification algorithm involved cells from May–September of 2011–2015. For this subset, we used the following methods to label cells as hail-producing convective cells (positive samples) according to hail case data. First, the cell nearest to the record location within 10 km at the record time from the hail reports was manually labeled. Then, we tracked backward from this cell to find all related convective cells in the same hail case because report time was commonly lagging and there were also WERs in the cumulus stage. There are, however, automatic tracking algorithms such as thunderstorm identification, tracking, analysis, and nowcasting (TITAN) [48], SCIT [46], optical flow [49], and some other advanced methods [50,51]. Errors were unavoidable for automatic algorithms, especially on the condition that there was cell merging and splitting. Therefore, cells were manually tracked to ensure that hail samples were reliable. Finally, 823 cells from 93 hail cases were labeled as hail samples.

Besides the positive set, a negative set that contained no-hail convective cells was prepared for comparison. However, hail cases were easy to miss due to the observation means and the large variability of a hail-swath [52]. Since the probability of hail occurrence was low, the numbers of positive and negative samples were also hugely imbalanced. Consequently, treating all unlabeled cells as negative samples may enormously influence the results. Therefore, after obtaining all identified cells in the dataset, we performed the following filter work to make the negative sample set as clean as possible. On the one hand, all cells that were close to any hail record in time and space were excluded. The specific rule was that when an observation station within 50 km from a cell had a hail record of within one hour, we removed this cell from the negative sample set. On the other hand, if a cell was not located above any observation station, it was filtered out. Most of these samples were located in the Bohai Sea, east of the radar. Through this filtering, although the negative set could not be completely clean, the impact of false-negative samples on the results was reduced.

3. Methodology

In this section, we present detailed descriptions of our methods. The first two subsections are about the WER morphology detection algorithm: Section 3.1 describes the method to find candidate WER areas using their morphological features, and Section 3.2 gives the rules to filter out candidate areas that only look like WERs. In Section 3.3, the automatic vertical cross-section-making algorithm is described. Lastly, Section 3.4 presents the WER quantification method.

3.1. Weak Echo Region Morphology Detection

As the weak echo region is a 3D structure, its identification needs to view multiple elevation scans. However, constructing a 3D image or simultaneously processing multiple images consumes a large amount of computing and storage resources. Considering the real-time detection capabilities, we designed the following method utilizing the 2D echo bottom image of convective cells.

According to the definition, a weak echo region is bounded by at least two strong echo parts, one part above and a tilted wall echo part on one side. The feature of the strong echo part above is that its echo bottom height is high, while echo bottom heights of the tilted wall echo part change rapidly in the same direction. Therefore, the basic idea of our WER identification algorithm was to find a pair of the region with a high echo bottom height and the region with high echo bottom gradients with the same direction. Figure 3 shows the flowchart of the WER detection algorithm. For convenience, we refer to the region with a high echo bottom height as HEBHR and to the region with a high echo bottom gradient region as HEBGR.

To improve the precision of echo bottom images, echo bottom heights were estimated with a linear interpolation method similar to [53]. The ${\mathrm{Z}}_T$$\mathrm{dB}Z$ echo bottom height $H_{eb}$ of a pixel is calculated as follows:

• Find the minimum elevation where reflectivity exceeds ${\mathrm{Z}}_T$. Call the angle of this elevation ${\theta }_u$ and reflectivity value ${\mathrm{Z}}_u$.
• If ${\theta }_u$ is the lowest elevation of the radar, set $H_{eb}=0$; otherwise, skip to the next step.
• Obtain the next lower elevation angle ${\theta }_d$ and its reflectivity value ${\mathrm{Z}}_d$. Then, echo bottom height $H_{eb}$ is given by:

  $$H_{eb}=d·tan\left[\left(Z_T-Z_d\right)\frac{{\theta }_u-{\theta }_d}{Z_u-Z_d}+{\theta }_d\right],$$

  (1)

  where d is the horizontal distance from this point to the radar location.

Figure 4 shows the schematic diagram of echo bottom estimation.

To use the echo bottom gradient information more conveniently, an echo bottom gradient direction image and an echo bottom gradient magnitude image were generated. The gradient of a non-boundary point is given by:

$$\nabla f(x,y)=\left[\begin{array}{c}{g}_x(x,y)\hfill \\ g_y(x,y)\hfill \end{array}\right]=\left[\begin{array}{c}{H}_{eb}(x+1,y)-H_{eb}(x-1,y)\hfill \\ H_{eb}(x,y+1)-H_{eb}(x,y-1)\hfill \end{array}\right],$$

(2)

where $H_{eb}(x,y)$ is the echo bottom height of the point with coordinates $(x,y)$ and $g_x(x,y)$ and $g_y(x,y)$ are components of the gradient in the direction of the two coordinate axes. The gradient direction of point $(x,y)$ can be calculated by the formula:

$$\theta (x,y)=arctan\frac{g_y(x,y)}{g_x(x,y)},$$

(3)

and the magnitude is given by:

$$∥G(x,y)∥=\sqrt{g_x^2(x,y)+g_y^2(x,y)}.$$

(4)

After obtaining the needed images, the first step was to find the HEBGRs in which the echo bottom gradient magnitude of each point was large and the gradient direction was similar. They were the possible positions of tilted wall echoes. The way to find HEBGRs is as follows.

• Select a gradient magnitude threshold $G_{min}$. Threshold the gradient magnitude image using $G_{min}$, that is a pixel in the image was set to zero if its value was less than $G_{min}$. Obtain the other connected components using the edge-following algorithm proposed by [47].
• Select an area threshold $S_{gmin}$. Filter out connected components whose area is less than $S_{gmin}$.
• Estimate the standard deviation of the gradient directions in each connected component. As a circular quantity, the standard deviation of gradient directions is given by [54]:

  $${\sigma }_{\theta }=\sqrt{-2ln\left(\overline{R}\right)}.$$

  (5)

  $\overline{R}$ is given by the expression:

  $$\overline{R}=\frac{1}{n}\left[{\left(\sum _{i=1}^{n}cos{\theta }_i\right)}^2+{\left(\sum _{i=1}^{n}sin{\theta }_i\right)}^2\right],$$

  (6)

  where n is the number of samples.
• Select a gradient direction radian standard deviation threshold ${\sigma }_{\theta max}$. Filter out connected components whose standard deviation is greater than ${\sigma }_{\theta max}$.

Then, for each HEBGR, determine whether there are paired HEBHRs. This can be done by the following steps:

• In the HEBGR, find the point with maximum echo bottom height $h_{g\max }$ and the point with minimum echo bottom height $h_{g\min }$. Calculate the height threshold by:

  $$h_H=0.5\left(h_{g\max }-h_{g\min }\right)+h_{g\min },$$

  (7)

  where the constant 0.5 is used to make sure that the heights of HEBHR are more than half of the heights of the paired HEBGR.
• Threshold the echo bottom image by $h_H$. Obtain the other connected components.
• Select area threshold $S_{hmin}$, and filter connected components whose area is less than $S_{hmin}$.
• If an HEBHR has overlapping parts with HEBGR and is not totally included in it, pair them as one candidate location of a weak echo region. Note that an HEBGR may be paired with multiple HEBHRs.

The schematic diagrams of the algorithm are shown in Figure 5.

The thresholds used in the detection algorithm included the minimum echo bottom gradient of HEBGR $G_{min}$, the minimum area of HEBGR $S_{gmin}$, the minimum area of HEBHR $S_{hmin}$, and maximum gradient standard deviation ${\sigma }_{\theta \max }$. $S_{gmin}$ and $S_{hmin}$ were used to eliminate noise. We set their values to 10 $\mathrm{k}{\mathrm{m}}^2$ since a typical WER is larger than that. To acquire the appropriate values of $G_{min}$ and ${\sigma }_{\theta \max }$, we conducted preliminary research by making statistics on several real samples. According to the results, we set $G_{min}$ to 1 and ${\sigma }_{\theta \max }$ to 1.2 (the units of both the horizontal and vertical distance are km).

3.2. Result Filter

The detection algorithm presented in the previous subsection identifies WERs by their morphological definition. However, in practice, there are many WER lookalikes that are not caused by updrafts. On the one hand, signatures such as overspreading anvil and three-body scatter spike [55] are morphologically similar to WER. On the other hand, some causes such as fast-moving storms [23] and echo-connecting may make echoes look like WERs. These echoes should be filtered out as much as possible. Therefore, we designed two filters to screen candidate WER locations. Some examples of WER lookalikes are shown in Figure 6.

As WERs are located within thunderstorms, the reflectivity of bounded echoes is strong. In contrast, reflectivity in an anvil or three-body scatter spike is relatively weak. Hence, identification results with only weak reflectivity should be filtered out. Considering the variety of thunderstorms, we used six reflectivity thresholds to run the detection algorithm and kept a result if at least three of them in the same position were identified. The reflectivity thresholds used in this paper were 20, 25, 30, 35, 40, and 45 $\mathrm{dB}Z$.

When updrafts flow into convective cells, reflectivity gradients are higher on the low-level inflow side. Therefore, one kind of solution is to recognize tight gradient regions of cells in low-elevation PPI images and discard detection results far away from these locations. However, in our preliminary attempt, we found that the actual situations were too complicated to identify high reflectivity gradient regions and determine whether a high gradient region was the inflow side. In order to improve stability and applicability, we designed a filter method by comparing the gradients of eight different cell directions. The method is presented as follows:

• Build the composite reflectivity image of the lowest three elevations. Calculate the reflectivity gradient of every grid point on this image. The gradient equation is similar to Equation (2), but echo bottom heights are replaced with reflectivity values in units of $\mathrm{dB}Z$.
• For each cell, find the centroid of the second maximum reflectivity region. Centering on this point, divide the cell image into eight equal angles, as shown in Figure 7a.
• Calculate the average gradient magnitude of each part. To eliminate noise points, the points whose gradient directions do not point near the center are excluded from calculating.
• Find the part with the maximum average gradient magnitude, and mark it and its two adjacent parts as the possible inflow region. Filter out detection results that have no overlapping area with the possible inflow region, as shown in Figure 7b.

3.3. Automatic Vertical Cross-Section Making

The radar reflectivity vertical cross-section is an important visualization tool to observe the internal structure of thunderstorms. In some applications, such as hail-suppression operations, rapidly and accurately generating vertical cross-sections can help people make decisions [56]. Vertical cross-sections are also able to exhibit the structures of WERs and BWERs. However, manually determining the location of such vertical cross-sections is time consuming and fallible because it needs to observe and analyze multiple elevations at the same time. Consequently, based on the automatic WER detection algorithm, we designed an algorithm to determine the location of a vertical cross-section automatically that can show a clear WER structure.

After a WER is detected, its vertical cross-section line could be determined by a point and a direction as follows. On the one hand, in order to display both the tilted wall echo part and the overhang echo part, the line should pass through the two regions at the same time, that is the line should pass through a certain point in the overlapping area of HEBGR and HEBHR. Thus, we chose the centroid of the maximum reflectivity region in the overlapping area of HEBGR and HEBHR as the point through which the line should pass. Using the centroid of the maximum reflectivity region is to show the most severe part of the thunderstorm. On the other hand, in order to show a clear WER structure, the direction of the straight line is preferably in the direction of the echo tilt. Therefore, we used the average direction of HEBGR $\overline{R}$ given by Equation (6) as the direction of the vertical cross-section line. Examples of vertical cross-sections obtained automatically by this method are shown in Section 4.3.

3.4. Weak Echo Region Quantification

The presence of a WER or BWER means there are updrafts inside thunderstorms, and a strong updraft is a prerequisite for hail formation. Unlike BWERs that commonly indicate particularly severe thunderstorms, the presence of WERs does not always mean the updraft strength is sufficient for hail growth. It would be helpful for hail nowcasting if we could find which WERs are related to hail-producing thunderstorms. Based on this motivation, we attempted to design a WER quantification method that could generate a parameter that reflects a correlation between WERs and hail-producing thunderstorms.

Considering the diversity and complexity of a WER structure, it is difficult to measure directly and represent it with one quantity. An alternative solution is to construct a quantity that is a combination of multiple components. Each component is a simple measurable value of one part of the WER. Several points were considered while selecting the components. First, previous studies found that broader updrafts are favored for enhanced hail growth [17,19,20]. The volume of the WER could reflect updraft width. Second, the echo overhanging the WER contains supercooled liquid water on which hail growth is dependent [26]. Thus, the intensity and volume of overhang echo may be relevant. Third, multiple reflectivity thresholds can be used to adapt to different thunderstorm events. Therefore, the following two kinds of volumes with different reflectivity thresholds $R_T$ were used as the components:

• Volume between the undersurface of the overhang echo and plane of height $h_T$, denoted by $V_{\mathrm{we},R_T}$.
• Volume of the overhang echo whose undersurface is higher than $h_T$, denoted by $V_{\mathrm{oh},R_T}$.

The formulas for estimating $V_{\mathrm{we},R_T}$ and $V_{\mathrm{oh},R_T}$ are:

$$V_{\mathrm{we},R_T}=\sum _{p_i\in S_{\mathrm{we}}}\left(h_{EB,R_T,p_i}-h_T\right),$$

(8)

and:

$$V_{\mathrm{oh},R_T}=\sum _{p_i\in S_{\mathrm{we}}}\left(h_{ET,R_T,p_i}-h_{EB,R_T,p_i}\right),$$

(9)

where $S_{\mathrm{we}}$ is the horizontal grid-point set of the area where the WER is located, which is obtained by the detection algorithm, $p_i$ is the $i^{\mathrm{th}}$ grid point of $S_{\mathrm{we}}$, $h_{ET,R_T,p_i}$ and $h_{EB,R_T,p_i}$ are the echo top height and echo bottom height at point $p_i$, respectively, and $h_T$ is a user-defined height threshold, discussed later. The schematic diagram describing the two parts of volumes is shown in Figure 8.

The setting of height threshold $h_T$ is to eliminate the influence of different minimum detectable heights on volume estimation. As shown in Figure 9a, if the same thunderstorm is detected in different radar ranges, the volumes between the undersurface of the overhang echo and the plane where the lowest radar beams are located are unequal. Therefore, we can make volumes stay constant by specifying the minimum height, as shown in Figure 9b. $h_T$ should be carefully selected according to radar information. When using multiple radars, thresholds should be consistent to ensure that volumes estimated by each radar are of the same standard. If radar altitudes vary, a solution is to limit the effective detection range of some radars. In this paper, $h_T=2\mathrm{km}$, since the lowest detectable height of the far end of the radar beams is 2.077 km.

After obtaining the components, a combination method should be used to construct the quantity. Since the quantity is mainly used to distinguish hail from non-hail, a supervised dimensionality reduction method, such as linear discriminant analysis (LDA), is an ideal choice [57]. LDA is a data-driven method that finds a linear combination of features that characterizes or separates two or more classes of objects. It uses features and class labels to find a projection hyperplane that minimizes interclass variance and maximizes the distance between the projected means of the classes. Compared with unsupervised dimensionality reduction methods, such as principal component analysis, LDA can use class information to obtain a projection with more interclass discrimination. In this task, we applied linear discriminant analysis to a historical dataset that included both hail and non-hail samples and then used the obtained coefficients for a linear combination to construct the quantity (hereinafter referred to as the WER parameter). The equation for WER parameter $V_W$ is:

$$V_W=\sum _{n=1}^N{C}_n{V}_{\mathrm{we},R_n}+\sum _{n=1}^N{C}_{N+n}{V}_{\mathrm{oh},R_n}$$

(10)

where N is the number of reflectivity thresholds, $R_n$ is the $n^{\mathrm{th}}$ reflectivity threshold value, $V_{\mathrm{we},R_n}$ and $V_{oh,R_n}$ are the volumes estimated using reflectivity threshold $R_n$, and $C_1,\cdots ,C_{2N}$ are coefficients obtained by linear discriminant analysis. The unit of $V_W$ is $\mathrm{k}{\mathrm{m}}^3$.

4. Results and Discussion

4.1. Detection Evaluation

We used manually-identified samples to evaluate the performance of the WER detection algorithm. As WERs are easy to miss and some samples are ambiguous, three experienced forecasters from the Tianjin Meteorological Observatory participated in this work. After independent identifications, results were cross-checked, and only samples identified at least twice were confirmed. The data for evaluation were Tianjin radar base data from May–September 2015, which included 32,747 valid volume scans and 16,170 identified severe convective cells. Finally, 294 WER samples were identified.

The contingency table shown in

Table 1. Contingency table for WER detection.

		Detections
		Yes	No
Events	Yes	$a=229$	$b=65$
	No	$c=81$	$d=$ 15,795

contains the evaluation results. Detection errors caused by radar data quality were not counted; where, in the contingency table, a is the number of events that were correctly detected, b is the number of events that failed to be detected, c is the number of false detected nonevents, often referred to as false alarms, and d is the number of undetected nonevents. The probability of detection (POD) and false alarm rate (FAR) are, respectively, given by:

$$\mathrm{POD}=\frac{a}{a+b}=0.78,$$

(11)

$$\mathrm{FAR}=\frac{c}{a+c}=0.26,$$

(12)

and the critical success index (CSI) [58] is:

$$\mathrm{CSI}=\frac{a}{a+b+c}=0.61.$$

(13)

From Table 1, it can be seen that both misdetection and false detection were inevitable since the method only used the reflectivity morphological features. By analyzing misses, we found that almost all of them were filtered by the low-elevation reflectivity gradient because the gradient was not obvious at the beginning of the inflow. However, with no filtering procedure, there were more false detections.

Then, we performed the automatic cross-section-making algorithm on all detected samples and manually checked the cross-sections generated by the algorithm. Results showed that all generated cross-sections can exhibited clear WER structures. Examples of cross-sections generated automatically by the algorithm are shown in Figures 14 and 16.

Our experiment system was written in Python with the OpenCV package [59]. The maximum running time of the WER identification and the cross-section making of each cell was 30 ms on a computer with an Intel i7-8700k processor. Therefore, these methods can meet the real-time requirements in most applications.

One thing to note is that the minimum detectable height had a significant impact on the identification of WERs. In this evaluation, we did not limit the range of detections because the true value was manually identified, and WERs occurring below the lowest radar beams are also invisible to people. However, in operational forecasting, effective detection range should be carefully considered according to the altitude and lowest elevation angle of the radar. Otherwise, the algorithm may miss some WERs that could provide valuable information for forecasting. When using multiple radars for detection, a better choice is to set the optimal detection range of each radar.

4.2. Quantification Evaluation

To evaluate the WER quantification method for hail nowcasting, we compared the WER parameters of WERs within hail-producing convective cells with ones in other cells. The dataset used in this evaluation is described in Section 2.

We ran the WER detection algorithm to obtain all cells with WERs and computed the values of $V_{\mathrm{we}}$ and $V_{\mathrm{oh}}$ with reflectivity thresholds of 30, 35, 40, 45, and 50 $\mathrm{dB}Z$ on these cells. The values of the cells without WERs were all set to zero. Figure 10 shows the box plot of the results. From this figure, it can be seen that the median of all components was zero since most cells do not possess WERs. Except for the 50 $\mathrm{dB}Z$, the distributions of $V_{\mathrm{we}}$ and $V_{\mathrm{oh}}$ with other reflectivity thresholds were significantly different between hail and no-hail convective cells. Rare no-hail cells also have WERs with larger reflectivity than 40 $\mathrm{dB}Z$. Therefore, we used $V_{\mathrm{we}}$ and $V_{\mathrm{oh}}$ with reflectivity thresholds of 30, 35, 40, and 45 $\mathrm{dB}Z$ as the components of linear discriminant analysis.

Normalized coefficients determined by linear discriminant analysis are shown in Figure 11. It can be seen that $V_{\mathrm{we},45}$ and $V_{\mathrm{oh},45}$ had the first two large coefficients. Generally, a component with a larger absolute coefficient value is a better predictor in linear discriminant analysis. It is worth noting that $V_{\mathrm{oh},30}$ had a negative coefficient. As its absolute value was small, this does not mean that $V_{\mathrm{oh},30}$ was negatively related to hail, but their correlation was not significant and the number of no-hail samples larger.

Figure 12a shows the probability density function of the WER parameters of cells with WERs, which is estimated by Gaussian kernel density estimation with automatic bandwidth selection [60,61]. For no-hail samples, most of their WER parameters were around zero, and very few samples exceeded 150. For hail samples, although most of their values were also near zero, and there were many samples with a WER parameter larger than 150. To compare the proportion of samples with different WER parameters, percentage changes are shown in Figure 12b. From this figure, it can be seen that, when the WER parameter was less than 150, no-hail samples accounted for the majority, while, when the value was in the range of 200–450, hail samples were in the majority. Detailed statistics for samples with values between 200 and 450 found that more than 85% were hail samples.

According to the above results, it is clear that the WER parameter cannot be used to forecast hail-producing thunderstorms alone. However, when values are in a certain range, more than 85% of the samples produced hail. Therefore, it can give a more confident forecast in some cases, and could be used as a parameter in multiparameter models to improve the overall performance of hail nowcasting.

4.3. Case Examples

Two hail cases from the Tianjin radar were used as examples to show the performance of the WER detection algorithm and the automatic vertical cross-section-making algorithm. Figure 13 shows the motion trajectories of cell centroids. Figure 14 shows the composite reflectivity images and the vertical cross-section images of the cells at each time step of Case 1. The line chart in Figure 15 displays the change of $V_{\mathrm{we},R_T}$ over time steps of Case 1. The images and line chart of Case 2 are shown in Figure 16 and Figure 17, respectively. Since the duration of Case 2 was long, only part of the time steps is shown. The vertical cross-sections of reflectivity in these figures were obtained automatically by the proposed algorithm.

For Case 1, the thunderstorm initiated in the south of Beijing and moved toward the south after 14:24 local time (LT = UTC + 8), 10 June 2014. Observation stations on the moving path started to report hail from 15:42 LT. At 14:54 LT, the algorithm started to detect WERs in the south of the cell. As convective cells evolved, WERs became larger. From Figure 15, it can be seen that the WER became largest at 07:54 LT and then began to decay. The cross-section images from 14:54–16:12 LT in Figure 14 show that WERs were bounded on both sides (BWERs).

For Case 2, the strong convective cell formed in the west of Tianjin at 14:30 LT, 27 July 2015. Then, it moved toward the northeast. The earliest hail record was at 16:48 LT. In the life cycle of this case, there were two significant WER processes. First, at 15:18 LT, the algorithm identified a WER in the northeast of the cell, but it lasted for a short time and did not grow large. A more significant WER was detected from 16:30 LT, and it became largest at about 17:00 LT. Figure 16 clearly shows the development of the two WER processes.

5. Conclusions

This paper proposed three algorithms related to weak echo regions, including a WER morphology detection algorithm, an automatic vertical cross-section-making algorithm, and a weak echo region quantification method.

The WER morphology identification algorithm finds candidate WER regions by pairing regions with a high echo bottom height and regions with a high echo bottom gradient, then filters candidate regions that are not caused by updrafts by using low-level reflectivity gradients. Validation results on real data of the Tianjin radar indicate that the identification algorithm had a CSI score of 0.61. Based on the identification algorithm, the automatic vertical cross-section-making algorithm determines the cross-section line by finding a point in the overlapping area of a high echo bottom height region and a high echo bottom gradient region, calculating the average tilt direction of the wall echoes. Two case examples showed that each vertical cross-section generated by the algorithm could exhibit obvious WERs. These two algorithms provide convenient tools for severe thunderstorm analysis that can be used in many applications, such as retrospective studies and operational hail suppression. Compared with manual cross-section making, it can save much time and be more accurate.

The WER quantification algorithm quantifies a WER with the linear combination of several components. The components are volumes of different WER parts, and the coefficients are generated by linear discriminant analysis trained on historical data. Results on Tianjin data showed that, when the WER parameter was between 200 and 450, more than 85% of the samples produced hail. It can be used as an auxiliary parameter in multiparameter forecasting approaches such as logistic regression [33] to improve the performance of hail nowcasting.

Although we used several reflectivity features to filter out WER lookalikes that were not caused by updrafts, detection performance was still limited because the single-polarization radar cannot provide more information. Therefore, applying the algorithm on dual-polarization radars and introducing more parameters may improve performance. As dual-polarization radars are not widely used in some regions, our algorithm is still useful.

In this paper, we only analyzed the relationships between WER parameters and hail-producing thunderstorms. In a future study, we plan to use the WER parameter as a feature in a machine learning classifier [62,63,64] with other classic hail parameters to test how much it can improve overall performance.