The U.S. Weather Surveillance Radar-1988 Doppler (WSR-88D) network has provided meteorologists and hydrologists with quantitative precipitation observations at an unprecedented high spatial–temporal resolution since its deployment in the mid-1990s. Since each single radar can only cover a maximum range of 460 km, a mosaic of multiple-radar observations is needed to generate any national-scale products. The Multi-Radar Multi-Sensor (MRMS) system utilizes a physically based two-dimensional mosaicking algorithm of the WSR-88D data to generate seamless national quantitative precipitation estimation (QPE) products. For areas covered by multiple radars, the mosaicking scheme first determines if precipitation is present by checking the lowest-altitude observation. If the lowest observed radar data indicate no precipitation, then the mosaicked value is set to no precipitation. Otherwise, a weighted mean of multiple-radar observations is taken as the mosaicked value. The weighting function is based on multiple factors, including the distance from the radar and the height of the observation with respect to the melting layer. The mosaic algorithm uses the physically lowest radar observations with no/little blockage while maintaining a spatial continuity in the mosaicked field. The performance of the MRMS seamless radar mosaic algorithm was examined for various precipitation events of different characteristics. The results of these case evaluations are presented in this paper.
The Weather Surveillance Radar-1988 Doppler (WSR-88D) network has been currently considered as a unique tool capable of providing precipitation fields at a high spatiotemporal resolution over the conterminous United States (CONUS; Crum and Alberty 1993; Klazura and Imy 1993; Howard et al. 1997; Serafin and Wilson 2000; Zhang et al. 2016). The WSR-88D network is composed of 146 radars over CONUS as of 2015 (Fig. 1). To mitigate the single-radar restrictions of finite coverage, as well as spatial and temporal gaps in radar data coverage (Zhang et al. 2005; Langston et al. 2007), the current study plans to combine the multiple-radar observations over the CONUS using a physically based methodology.
As well known, observations from the WSR-88D network have been utilized by government agencies as well as private companies for many meteorological, hydrological, and aviation applications. One of the main applications is generating quantitative precipitation estimation (QPE). Accurate precipitation estimates are pivotal for numerous applications related to climate studies, mesoscale model validation, flash flood detection (Yao et al. 2012; Yao et al. 2014), and other hydrometeorological applications (Wan et al. 2015; K. Zhang et al. 2007, 2011; Zhang and Gao 2016; Zhang et al. 2014). These applications range from the simplest use of single radar to mosaicked multiple-radar observations with sophisticated algorithms. The earliest basic use of radar observations is in a single-radar field.
A single WSR-88D radar can provide a three-dimensional (3D) structure of the storm and bring meteorologists critical information. Numerous techniques have been developed for gridding WSR-88D and other ground- or airborne-based radar observations in Cartesian space. Commonly used interpolation techniques include the nearest-neighbor methodology (e.g., Jorgensen et al. 1983); linear interpolation (e.g., Mohr and Vaughn 1979; Miller et al. 1986; Fulton 1998); the Cressman weighting scheme (e.g., Weygandt et al. 2002); and the Barnes, or exponential, weighting scheme (e.g., Askelson et al. 2000; Shapiro et al. 2003). Zhang et al. (2005) proposed several remapping schemes to interpolate and smooth the radar data field with simplified logic for applying radar observations to a Cartesian grid.
The aforementioned interpolation techniques are limited to a single-radar field, and the inclusion of multiple radars is a natural improvement upon single-radar applications. A single radar has limitations for providing accurate data over a particular area. These restrictions include limited coverage area (Maddox et al. 2002), spatial and temporal gaps in radar data coverage (Zhang et al. 2012a, 2016), and sensitivities to vertical temperature profiles in the atmosphere (Zhang and Qi 2010; Qi and Zhang 2013; Qi et al. 2013a,b). Mosaicking multiple-radar observation fields onto two- or three-dimensional 3D Cartesian grids can mitigate some of these challenges. To provide more physically and scientifically sound depictions of meteorological phenomena through multiple observing systems, Zhang et al. (2005) tested three kinds of 3D mosaicking techniques: nearest neighbor, maximum value, and weighted mean. The nearest-neighbor scheme would select the observation from the closest radar at each grid point and does not impose any smoothing when generating a mosaicked field from multiple radars. This would result in discontinuities midway between adjacent radars. The maximum value scheme would just use the highest reflectivity intensities from multiple radars that cover the same grid cell, which would result in overestimation biases from the vertical profile of reflectivity (VPR) issue, and radars may produce artificially high reflectivities due to the calibration problems. The weighted mean scheme weighs multiple-radar observations covering a grid cell based on the distance between the grid cell and each corresponding radar site location and takes the weighted average as the mosaicked value. It creates a smooth transition between neighboring radars and does not bias toward “hot” radars when generating the mosaicked reflectivity field. Langston et al. (2007) proposed a four-dimensional dynamic radar mosaic (4DDRM) scheme that applies weighted mean techniques similar to Zhang et al. (2005) but accounts for the different ages of the observations covering the same grid cell. This scheme gives more weight to newer observations and less weight to older ones, and reduces position error of fast-moving storms.
However, in order to generate accurate radar QPE, the current study is focused on generating a 2D mosaicked reflectivity field for deriving surface precipitation rate. Tabary (2007) proposed a weighted linear mosaicking scheme, and the weight was dependent on a quality index associated with ground-clutter identification, reflectivity-to-rain-rate conversion, partial beam blocking correction, VPR correction, and synchronization. The quality index reflects roughly the quality of the surface rainfall rate. J. Zhang et al. (2011, hereafter Z11) proposed a 2D mosaicking scheme of hybrid scan reflectivity (HSR) in the National Mosaic and Multi-Sensor QPE (NMQ) system and the same technique is adapted in the early version of Multi-Radar Multi-Sensor (MRMS) system (Zhang et al. 2016). The HSR mosaic scheme is based on two weighting factors, one in horizontal and one in vertical (Z11), and it generates QPE fields with better spatial continuity than the aforementioned nearest-neighbor or maximum value schemes. However, the vertical weighting factor in Z11 allowed significant impact from data of higher altitudes and was one of the factors that caused MRMS QPE underestimation, where the observations from the neighboring radars are at far ranges.
To generate an accurate surface precipitation intensity estimation in a regional or national domain, the main goal of this study is to optimally merge radar observations with two considerations: mitigate nonprecipitation echoes’ impact, and mitigate/reduce the impacts of the observations from the melting layer and ice region. Multiple radars can scan different levels of the same storm structure at the given WSR-88D network density (Fig. 2). To obtain the accurate ground rainfall observation, the lowest radar data with little blockage should be used at any given location to avoid potential QPE errors associated with discrepancies between the radar observation height and the ground. Such discrepancies include evaporation, condensational and coalescence growth, drifting, and orographic forcing. However, a simple “taking the lowest” approach could result in discontinuities midway between neighboring radars, as shown in the nearest-neighbor results in Zhang et al. (2005). To assure a spatially continuous precipitation field, a weighted mean of multiple-radar observations is desired for the radar mosaic. Which and how many overlapping radar observations covering a given location should be combined to create the Cartesian 2D reflectivity mosaic, and what weighting function should be used? This paper tries to answer these questions by studying a number of different radar coverages and precipitation distribution scenarios. Three steps are involved in the 2D reflectivity mosaic for the surface precipitation estimation: 1) generating a seamless hybrid scan reflectivity (SHSR) from a single-radar volume scan in the polar grid, 2) remapping the polar SHSR onto a Cartesian grid, and 3) mosaicking multiple-radar SHSR fields to create a regional/national SHSR map. Section 2 provides a detailed description of the three steps. Case study results representing various geographical regions and different seasons in the United States are presented in section 3. A summary and discussion of future work follows in section 4.
2. SHSR mosaic methodology
a. Single-radar SHSR
The first step is to generate SHSR from a single-radar volume scan in the polar grid. The lowest-altitude radar reflectivity bins from the HSR (Fulton 1998; Z11; Zhang et al. 2016) are derived from quality-controlled single-radar reflectivity data (Tang et al. 2014) based on the radar beam blockage data (Zhang et al. 2012a, 2016). Radar beam blockage data are calculated based on the radar beam propagation under standard atmospheric conditions and high-resolution topographic data generated from the NASA Shuttle Radar Topography Mission (SRTM; Farr et al. 2007). The HSR does not account for nonstandard blockages due to objects above the ground, such as high buildings, man-made towers, and trees, which often result in gaps in the reflectivity field (Tang et al. 2013; Zhang et al. 2016). Such gaps are identified from precipitation accumulation maps, and a linear cross-azimuth interpolation is applied to fill in small gaps. Gaps larger than 5° azimuth are filled with data from the unblocked upper tilts. This process is called nonstandard blockage mitigation (NSBM; Tang et al. 2013), and the resultant HSR field is called “seamless” HSR (i.e., SHSR; Fig. 3). At the same time, the SHSR field has been corrected with apparent vertical profile of reflectivity (AVPR; Zhang and Qi 2010; Qi and Zhang 2013; Qi et al. 2013a,b) to account for variations of reflectivity between the height of SHSR observations and the ground.
The second step is to remap the polar SHSR onto a Cartesian grid. The single-radar SHSR field is interpolated from the native spherical coordinate onto a 2D Cartesian grid that is centered at each radar site and in a cylindrical equidistant map projection (Snyder 1987, 1993). For WSR-88D, the Cartesian grid covers a 460-km range for coastal radars and a 300-km range for inland radars in the MRMS system, and it has a horizontal resolution of 0.01° (approximately 1 km × 1 km).
b. Multiple-radar SHSR mosaic
The current mosaic scheme was developed based on Z11, which created a national mosaic of SHSR using a weighted mean scheme. In Z11, the weight given to each single-radar SHSR was a function of the distance di from the radar (Fig. 4) and the height hi of the SHSR above the mean sea level (MSL). The weight is described by the following equations:
Here, MR means the mosaicked reflectivity, i is the radar index (where i = “radar 1,” “radar 2,” etc.), and R is the single-radar reflectivity field. There are two components in the weighting function, one for the horizontal wL and another for the vertical wH. The parameters L and H are scale factors of the distance and height weighting components, respectively. Any SHSR data at 5 km MSL or higher were not included in the mosaic unless there were no other lower-altitude SHSR data at the same location. The weighting function reflects the level of confidence given to each radar observation when it is used to produce an estimation of the ground precipitation intensity. Because of the beam spreading and height increase with range, the radar-derived rainfall estimates are usually less accurate at far ranges than at close ranges.
Because of the unconditionally weighted mean of all SHSR data below 5 km MSL in Z11, false light precipitation was often found in areas of cloud anvil and virga. In addition, some SHSR data from ice regions or areas where the radar beam partially overshot cloud tops were allowed in the mosaic and caused underestimation in the Z11 QPE. To mitigate these over- and underestimation issues, a new physically based mosaic scheme (QZ16 hereafter) with additional constraints was developed in the current study. The new constraints are based on the radar QPE quality index (RQI; Zhang et al. 2012b) that is a function of the radar beam blockage, beam height, and the freezing-level height in the atmosphere.
RQI values range from zero to one, with lower values associated with significant blockage, the radar beam positioned well above the melting layer, or a combination of these two. Values closer to one would mean minimally obstructed (or nonobstructed) beam sampling below the melting layer. Lower RQI values more often occur in the mountainous terrain due to blockage and in winter months when the melting layer is closer to the ground. If the RQI = 1.0, the radar bins are in rain without any blockage. If the RQI > 0.6 and <1.0, the radar bins are in the melting layer–affected area, or in rain with a beam blockage of 10%–50%. If the RQI ≤ 0.6, the radar bins are in the ice region above the melting layer or within rain or the melting layer but have a beam blockage of 10%–50%. Note that no data with more than 50% blockage are allowed in the single-radar SHSR field.
In a relatively dense weather radar network such as the U.S. Next Generation Weather Radar (NEXRAD), a given grid point can be monitored by multiple radars; thus, varying precipitation types, including those found in the rain area, melting layer, and ice region, can be observed at the different altitudes (Fig. 5). At such a grid point with multiple-radar coverage, the QZ16 mosaicking scheme goes through a decision tree to determine what the value is (Fig. 6). QZ16 first determines if precipitation is present by checking the lowest-altitude observation. If the lowest observation indicates no precipitation, then the mosaicked value is set to no precipitation. Otherwise, the precipitation area is divided into transitional or nontransitional zones. The transition zone is defined as
where RQIlowest means the RQI value comes from the lowest elevation and RQIlowest+1 means the RQI value comes from the next elevation above the lowest elevation. In the transitional zone, a weighted mean of three radar observations with the highest RQI values is taken as the mosaicked value. The weighting function is the same as in Eqs. (1)–(3), except that the height is above the ground level (AGL) instead of above mean sea level. In nontransitional areas, the mosaic value was set to the lowest-altitude SHSR value if it also has the highest RQI value among all the SHSRs. Otherwise, a weighted mean of the lowest-altitude SHSR and the SHSR with the maximum RQI is assigned to the mosaic value. This usually occurs when the lowest-altitude SHSR is partially blocked.
The new SHSR mosaic scheme was examined for a large number of precipitation events of different characteristics, and example results from several representative events are presented in the current section.
a. Virga and cloud anvil removal
Abnormal propagation (AP) radial, cell discontinuities, and large cloud anvils can exist prior to the application of the QZ16 scheme (Fig. 7a). The large cloud anvils were observed by KICT (Wichita, Kansas) [KAMA (Amarillo, Texas)] at the distance of 200 km (150 km) and in the lowest beam height of 4.5–6 km (4–10 km) MSL (not shown); however, KTLX (Oklahoma City, Oklahoma), which is 20 km away, observed no rain in this area in the lowest beam (height of 0.2–2 km MSL). Since Z11 included all SHSR data below 5 km MSL, the cloud anvil echoes (white circle in Fig. 7a) got into the mosaic and resulted in false light rain marked in Fig. 7b. The QZ16 scheme, on the other hand, correctly identified the area as no precipitation (Fig. 7c) by checking the lowest-altitude SHSR from KTLX radar. The same rule in QZ16 also helped to avoid a false storm echo presented in Z11 at ~110 km north-northwest of KFDR (Frederick/Altus Air Force Base, Oklahoma) in Fig. 7a. The false echo was due to the large beam size of KAMA (~230 km away) and KICT (~290 km away) in the area that resulted in smeared SHSR observations (Figs. 8a,b) of a strong storm. The new mosaic logic in QZ16 produced a more physically realistic depiction of the storm by eliminating echoes outside the storm areas observed by the closest radar (KFDR, Fig. 8c). The cloud anvil echoes and the storm false echoes have been successfully removed in the new mosaicked reflectivity (Fig. 7d).
b. Dampening of SHSR intensity
The weighting of precipitation throughout the vertical column can dampen the reflectivity values and the subsequent precipitation rates (e.g., Fig. 9). The intensity of SHSR field mosaicked with Z11 has been largely reduced (Fig. 9a), especially in the highlighted region, where there is a dampening of ~5 dBZ compared to the QZ16 SHSR (Fig. 9b). A detailed investigation at point A (Fig. 9a) indicated that the Z11 SHSR was obtained from a weighted mean of single-radar SHSRs (not shown) from KEOX (Ft. Rucker, Alabama; 54.5 dBZ, 1.15 km MSL), KFFC (Atlanta, Georgia; 33.5 dBZ, 4.33 km MSL), KMXX (Maxwell Air Force Base, Alabama; 46.5 dBZ, 2.8 km MSL), and KEVX (Northwest Florida/Eglin Air Force Base, Florida; 42.5 dBZ, 3.43 km MSL). The QZ16 SHSR was a value only obtained from KEOX. Since the freezing-level height was ~3.5 km, KFFC SHSR (4.33 km MSL) was largely from the ice region and thus a relatively low reflectivity. The inclusion of this high-altitude, low-intensity SHSR resulted in a lower value in the Z11 SHSR mosaic (Fig. 9a) than in the QZ16 (Fig. 9b). The dampening of the mosaicked SHSR field from Z11 partially contributed to an underestimation of the subsequent 1-h QPE accumulation, which was ~32% lower than the gauge observations on a domain average (Fig. 9c). The QZ16 QPE, by mosaicking only selected SHSRs below the melting layer, matched significantly better with gauges (Fig. 9d) with a 17% underestimation. The new physically based mosaic also reduced the root-mean-square error (RMSE; from 3.73 to 2.91 mm) and improved the correlation coefficient (CC; from 0.90 to 0.93; Figs. 9e,f).
c. Brightband contamination
Bright band (BB) is a region of artificially inflated reflectivity around the melting layer in radar observations. The inflated power reflected back to the radar is related to the type of the droplets (ice or liquid water) and its dielectric constant. The dielectric constant for ice is 0.197 and for liquid water is 0.93. Since snowflakes are usually larger than raindrops, the sheer size of the melted snowflakes can contribute to greater reflectivity values, since the reflectivity factor is proportional to the sixth power of the particle size. Therefore, the radar beams striking the melting layer of snow cause overestimation of precipitation (Fabry and Zawadzki 1995; Zhang and Qi 2010). An event on 5 January 2015 with a freezing-level height of approximately 2 km MSL shows the BB appearing at a range of 70–150 km and at tilt 0.50° (Figs. 10b,e) across three radars with light beam-broadening affect. If no VPR correction is applied, the BB will cause huge overestimation of radar QPE. In a common region covered by the three radars, the SHSR fields from different radars intersected the melting layer (ML) at different ranges (Figs. 10a–c). KLGX (Langley Hill, Washington) SHSR was at tilt 0.5° and intersected the melting layer at a ~150 km range (Fig. 10a). KATX (Seattle, Washington) SHSR was at tilt 0.5° and intersected the ML at a much closer range (~60 km) with the bright band extended to a ~180 km range (Fig. 10b). KRTX (Portland, Oregon) SHSR was at tilt 0.5° and intersected the ML at around 100–170 km (Fig. 10c). All the SHSR fields are below 5 km MSL in the area (Figs. 10d–f), and Z11 SHSR mosaic (Fig. 11a) combined the SHSRs from all three radars and showed a large area of 35–45 dBZ due to the impact from the KATX and KRTX SHSRs (Figs. 10b,e and 10c,f). QZ16 SHSR mosaic (Fig. 11b) shows the mitigation of the BB impacts by using the reflectivity data from the lowest available radar data from KLGX (Figs. 10a,d). The reflectivity values were on average reduced by 7 dB.
When dealing with BB contamination, one challenge that has to be considered is inevitable multiphase combination of radar observations for multiradar mosaicking. Recall that the MRMS RQI product can provide a proxy for the quality of the radar observation based on beam blockage and the location of the beam with respect to the BB. The freezing height used in RQI is derived from the model sounding at each radar site, and that value is used for the entire radar coverage area. So, the current RQI product does not fully represent the precipitation phase distribution, especially when there are large variations of the freezing height over an area. In the example event shown in Fig. 12, the KBBX (Beale Air Force Base, California) radar had considerable beam blockage to the north and east while the KRGX (Reno, Nevada) radar was unobstructed. More importantly, the freezing-level heights varied at each site. The freezing-level height was approximately 2.1 km MSL at KBBX and 3.1 km MSL for KRGX. This is likely a by-product of the radars being sited on the opposite side of the Sierra Nevada.
The SHSR values shown suggest that the melting layer observed by KBBX may be wider because the freezing-level height becomes higher from the bottom to the top of the Sierra Nevada (Fig. 12). Only ice regions or partial melting layers are observed by radar KRGX because the freezing height becomes lower from radar KRGX site to the bottom of the Sierra Nevada. Figure 13a shows that the mosaicked reflectivity field between KBBX and KRGX as generated by the physically based mosaicking scheme without multiphase combination (i.e., no transition zone is set) only has precipitation information from two radar sites (Fig. 13a). Radar KBBX (68 m MSL) is located in the valley and radar KRGX (2560 m MSL) is on the top of the mountain. Radar KBBX can provide the precipitation information from rain area and melting layer, where radar KRGX can only provide the precipitation information from the ice region. This leads to discontinuities between radar KBBX and KRGX since only one observation existed for each precipitation phase. The removal of the discontinuities by setting a transition zone based on the logic shown in Eq. (4) allows for multiphase observations to be combined and the subsequent removal of the discontinuity (Fig. 13b).
4. Summary and future work
A physically based mosaicking scheme (QZ16) was developed to combine radar data in a dense radar network for accurate quantitative precipitation estimation. A radar beam can go through the melting layer and different phases of precipitation as the range increases. The new mosaicking scheme takes into account different phases of precipitation particles and largely avoids the multiphase contamination (overestimation from observations in the melting layer and underestimation from observations in the ice region) when combining multiradar observations. Further, the new scheme greatly reduced false light precipitation from overhanging precipitation that was present in the Z11 scheme by checking the precipitation/no precipitation condition in the lowest-altitude SHSR. Case study results showed that the new weighted mean mosaic conditioned on RQI values provided improved precipitation estimates over the Z11 scheme by limiting the impact of bright band and overshooting radar observations in the mosaic.
While the new mosaicking scheme provides improvements over the previous Z11 scheme, discontinuities may occur in the cold season when strong surface temperature gradients exist near the freezing point. The discontinuities were a result of a deficiency in the RQI calculation (Zhang et al. 2012b) where a single freezing-level height at the radar site was used for the whole radar domain. As a result, the RQI values at the midpoint of two neighboring radars, one on the warm side and another on the cold side of the large temperature gradients, are substantially different even though the radar beam heights and the freezing level are the same at the midpoint. The mismatch of RQI fields can cause undesired selection of SHSR data in the mosaic where a higher-altitude SHSR from the radar on the warm side of the temperature gradients is applied right next to a lower-altitude SHSR from the radar on the cold side. Work is ongoing to refine the RQI product using spatially varying 2D freezing-level height over the radar domain. Polarimetric radar fields, especially the correlation coefficient ρHV and differential reflectivity ZDR, will be used to better delineate the bottom of the melting layer (Fig. 14). The new RQI product is planned for operational implementation in MRMS by the end of 2017 and the QZ16 scheme will be further improved.
Funding was provided by NOAA/Office of Oceanic and Atmospheric Research under NOAA–University of Oklahoma Cooperative Agreement NA11OAR4320072, U.S. Department of Commerce. The authors also thank Steven Martinaitis [Cooperative Institute for Mesoscale Meteorological Studies (CIMMS)/University of Oklahoma (OU)–National Severe Storms Laboratory (NSSL)] for his thoughts and review of the manuscript.