Abstract

Using radar observations, the performances of the ensemble square root filter (EnSRF) and an indirect three-dimensional variational (3DVar) data assimilation method were compared for a mesoscale convective system (MCS) that occurred in the Front Range of the Rocky Mountains, Colorado (USA). The results showed that the root mean square innovations (RMSIs) of EnSRF were lower than 3DVar for radar reflectivity and radial velocity and that the spread of EnSRF was generally consistent with its RMSIs. EnSRF substantially improved the analysis of the MCS compared with an experiment without radar data assimilation, and it produced a slight but noticeable improvement over 3DVar in terms of both coverage and intensity. Forecast results initiated from the final analysis revealed that EnSRF generally produced the best prediction of the MCS, with improved quantitative reflectivity and precipitation forecast skills. EnSRF also demonstrated better performance than 3DVar in the prediction of neighborhood probability for reflectivity at thresholds of 20 and 35 dBZ, which better matched the observed radar reflectivity in terms of both shape and extension. Additionally, the humidity, temperature, and wind fields were also improved by EnSRF; the largest error reduction was found in the water vapor field near the surface and at upper levels.

1. Introduction

Mesoscale convective systems (MCSs) are severe convective storms that can cause injury and damage property. However, accurate prediction of MCSs remains a challenge because of their rapid development and evolution. Doppler radars are platforms capable of observing MCS structure with high resolution and frequency, and such observations have been used for convective-scale data assimilation to improve the initial conditions of numerical weather prediction (NWP) models.

The three-dimensional variational (3DVar) data assimilation (DA) system, which includes mass continuity equations and other model equations as weak constraints, is efficient for storm-scale DA. Xiao et al. [1, 2] developed a Weather Research and Forecasting DA (WRFDA) 3DVar system to assimilate radial velocity and radar reflectivity by considering the total water vapor mixing ratio as a control variable, which improved quantitative precipitation forecasts for a hurricane. However, their system is limited to warm rain microphysics. Gao and Stensrud [3] used the background temperature to classify hydrometers to include ice and snow microphysics. Wang et al. [4] developed an approach for indirect radar reflectivity assimilation that assimilates retrieved rainwater and in-cloud water vapor estimated from radar reflectivity. An advantage of this new approach is that it avoids linearization errors attributable to the nonlinear relationship between reflectivity and microphysical variables. It has been used in several recent studies and it has demonstrated capability in improving short-term forecasts of convective storms [57]. However, it cannot overcome the inherent weakness of the 3DVar method that results from the neglect of the flow-dependent nature of the background error covariance (BEC). This problem is most severe in storm-scale DA because few state variables are observed directly and large-scale balance relationships are invalid.

Another advanced DA method for convective-scale NWP is the ensemble Kalman filter (EnKF), which estimates BEC using an ensemble of forecasts. EnKF is able to evolve the BEC dynamically throughout multiple assimilation cycles. Cross covariance produced by EnKF is especially important for convective-scale DA because state variables that are not observed directly can be updated. Snyder and Zhang [8] and Zhang et al. [9] first demonstrated the capability of EnKF for convective-scale DA using radar radial winds acquired during observing system simulation experiments. Subsequent studies by Xue et al. [10, 11] showed the performance of EnKF could be improved further when radar reflectivity was also assimilated into a compressible model with complex ice microphysics. The application of EnKF for assimilating real radar data has also produced successful results for several convective cases [12, 13]. It has been found that EnKF can handle complex and highly nonlinear processes involved in DA, which makes it attractive for convective-scale DA [1416].

Several studies have demonstrated that the performance of EnKF was better than 3DVar in global- to mesoscale DA [1720]. Whitaker et al. [18] found that EnKF produced an improvement in analysis and forecasting relative to 3DVar when using the National Centers for Environmental Prediction (NCEP) Global DA system, especially in data-sparse regions. Meng and Zhang [20] showed that EnKF generally outperformed 3DVar by assimilating observations from either individual or multiple data platforms (e.g., soundings, and surface and wind profilers) for a mesoscale convective vortex event. However, few studies have compared the performances of 3DVar and EnKF for real convective storms. Johnson et al. [21] compared EnKF and 3DVar at multiple scales, and they found that the method of radar assimilation using EnKF could maintain the storm features throughout the entire forecast period, whereas a 3DVar forecast produced some deficits after the first hour.

In the context of convective-scale DA, the present study systematically compared the ensemble square root filter (EnSRF) and a newly developed indirect 3DVar method. The indirect 3DVar method used was that proposed by Wang et al. [4], which includes a procedure to assimilate an estimated in-cloud humidity field. A case study was conducted to investigate the impact of using different radar DA techniques on both the analysis and the subsequent prediction of an MCS that occurred over the Front Range of the Rocky Mountains in Colorado (USA) on 30 July 2014, which demonstrated the promising results of EnSRF over 3DVar in MCS analysis and forecasting. The remainder of this paper is arranged as follows. Section 2 presents the indirect 3DVar and EnSRF radar DA methods. The case description and the experimental setups are outlined in Section 3. Section 4 describes the analysis and forecast results, and the summary is presented in Section 5.

2. Method

2.1. WRFDA 3DVar Indirect Radar Data Assimilation Method

The indirect reflectivity assimilation approach of the WRFDA 3DVar system for convective-scale radar DA [4] was used in this study. The cost function in an updated system can be expressed as follows:where the subscript represents the background and the subscript represents the observation terms. To avoid the linearization error of the reflectivity–rainwater equation, two additional observation terms, corresponding to the rainwater mixing ratio and in-cloud water vapor (both estimated from reflectivity), are added in the cost function. Here, we extended the system by including snow and graupel in and , where

Here, represents the hydrometeor mixing ratios (rainwater, snow, graupel–hail) of the atmospheric state, are their background first guesses, are their observations retrieved from radar reflectivity, and and are their related observational variances and background error matrix, respectively; , , and are the water vapor mixing ratios, their observations retrieved from radar reflectivity, and their related observational variances, respectively. The radar radial velocity assimilation followed the procedure in Xiao et al. [1]. In the assimilation of radial velocity, the operator includes the rainwater terminal velocity, and the terminal velocity can be calculated according to the rainwater mixing ratio.

2.2. EnSRF

For the standard EnKF, the update equations can be formulated as follows:where is the ensemble state vector, is the Kalman gain matrix, is the observation operator that projects state variables to observed quantities, and is the linearized version of . The covariance matrices for the observation and background errors are and , respectively, denotes the NWP model, and is the current analysis time. Here, subscript denotes the ensemble member order ranging from 1 to , and is the ensemble size. The superscripts , , and denote the background, analysis, and observation, respectively, and the overbar represents an ensemble mean.

To avoid introducing additional sampling errors to the ensemble, we used the EnSRF algorithm [22], which replaces (4) with the following equation:where is a factor in the square root algorithm and is the covariance inflation factor [23], which is used to deal with sampling and model errors. is an identity matrix. The prime symbol (′) in the equations represents the perturbation of the ensemble members.

3. Case Description and Experiment Design

An MCS that occurred over the Front Range of the Rocky Mountains in Colorado (USA) on 29–30 July 2014 was selected to investigate the differences between 3DVar and EnKF in storm-scale DA. The MCS originated in an environment associated with a cold front and a cyclone. The low-level precipitable water and wind field at 850 hPa 0000 UTC 30 July 2014 are shown in Figure 1. Convergence between the northerly cold dry air and southerly warm moist air caused intense convection over northern Colorado, where the precipitable water was >40 kg·m−2 (Figure 1(a)) and the maximum wind speed was up to 15 m·s−1 (Figure 1(b)). The system persisted for 13 h and it caused heavy rain, strong winds, and frequent lightning.

Figure 2 shows the evolution of this MCS in terms of radar composite reflectivity at different times. At 2100 UTC 29 July 2014 (Figure 2(a)), convective cells had become organized into a convective line over eastern Colorado, with maximum reflectivity of 55 dBZ. There were also some isolated convective cells in southern Wyoming and western Colorado. At 0000 UTC 30 July (Figure 2(b)), the MCS developed and became substantial. The northern part developed rapidly with a large area of reflectivity of 30–45 dBZ. At 0400 UTC 30 July (Figure 2(c)), the MCS was in the mature stage and its southern part became larger and intensified further. During the southeastward movement, some isolated convective cells became incorporated into the moving convective system. By 0900 UTC 30 July (Figure 2(d)), the area of intense radar echoes had contracted, weakened, and moved into Kansas.

The Advanced Research WRF (ARW) version 3.9 was adopted as the forecast model used in this study. The model domain was set using two-way nesting between three nested grids, which had horizontal grid spacings of 15, 3, and 1 km that generated 211 × 161, 321 × 281, and 802 × 751 horizontal grid points, respectively. All domains had 51 vertical grid levels from the surface up to 50 hPa.

The model domains and the locations of the radar stations are shown in Figure 3. The Kain–Fritsch (KF) cumulus parameterization scheme [24] was used only in domain 1. Other model physics schemes adopted included the Thompson microphysics scheme [25], Mellor–Yamada–Janjić (MYJ) planetary boundary layer model [26], Noah land surface model [27], and Rapid Radiative Transfer Model for GCMs (RRTMG) for longwave and shortwave radiation [28].

As shown in Figure 3, nine high-resolution radars of the Weather Surveillance Radars-1988 Doppler network were used in the DA. Quality control that included despeckling, removal of ground clutter for reflectivity, and velocity de-aliasing (unfolding) [29] was conducted before the radar data were assimilated in the 3 km domain. For the overlapping radar observations, a data location check is conducted in both 3DVar and EnSRF DA systems to avoid counting the same radar observations. When the distance between radar data and radar site is longer than another data, the radar observation is not assimilated. Observations used for rainfall verification were quantitative precipitation estimates from the Multi-Radar/Multi-Sensor System developed by the National Severe Storms Laboratory [30].

Three experiments comprising NoDA, 3DVar, and EnSRF (Table 1) were designed to examine the impact of different radar DA methods on both the analysis and the prediction of the MCS. The NoDA experiment did not assimilate any data, and the initial and lateral boundary conditions were obtained from the 0.5° × 0.5° NCEP operational Global Forecast System (GFS). For the 3DVar experiment, a 12 h spin-up forecast was performed starting from 1200 UTC 29 July 2014. The radar data were assimilated every 5 min for 1 h in the 3 km domain, and the National Meteorological Center (NMC) method [31] was used to estimate background error covariance for assimilation of radar observations. Then, a 5 h deterministic forecast was conducted from 0100 to 0600 UTC 30 July 2014 in the 1 km domain.

For the EnSRF experiment, as shown in Figure 4, 40 members were generated using the random-CV [32] method after an 11 h spin-up forecast initialized from the NCEP GFS analysis at 1200 UTC 29 July 2014. The amplitudes of the perturbations for the horizontal wind (), potential temperature (), and water vapor mixing ratio () were approximately 2 m·s−1, 2 K, and 1 g·kg−1, respectively. Then, a 1 h ensemble forecast was conducted to develop the BEC for sampling the mesoscale environmental uncertainties. Radar data were assimilated every 5 min until 0100 UTC 30 July 2014 in the 3 km domain. The horizontal and vertical covariance localization radii were 12 and 4 km, respectively. A multiplicative covariance inflation following Xue et al. [11] with a factor of 1.25 was used to help maintain the ensemble spread. Additive noise [33] was also added to the analysis members with standard deviations of 0.5 m·s−1 for and , 0.5 K for , and 0.1 g·kg−1 for . Finally, 5 h deterministic and ensemble forecasts were produced at 0100 UTC 30 July 2014 in the 1 km domain from the interpolated ensemble mean and the analysis members, respectively.

4. Results

4.1. Analysis Result

To evaluate the performances of the 3DVar and the EnSRF analyses quantitatively, the root mean square innovations (RMSIs) and the ensemble spread were calculated for radar reflectivity and radial velocity during the 1 h assimilation period (Figure 5). The RMSIs provide a measure of the overall fit of the model state to the observations, and the ensemble spread can be used to examine analysis uncertainty. The calculation was limited to regions where reflectivity was >15 dBZ. Generally, the RMSIs of radar reflectivity and radial velocity for both 3DVar and EnSRF tended to decrease with time. EnSRF showed lower RMSIs than 3DVar during all analysis and forecast cycles, especially for the radial velocity at t = 10–30 min. This suggests EnSRF had smaller analysis error compared with 3DVar. The spread of EnSRF was slightly lower than the RMSIs for radar reflectivity and radial velocity in the first 20 min, suggesting underdispersion of the ensemble. Such underdispersion is a common problem in real radar DA at the convective scale [12, 34]. However, the ensemble spreads of radar reflectivity and radial velocity were consistent with the RMSI values in the following cycles, indicating the forecast error was representative of the ensemble spread.

Figure 6 shows the analyzed composite radar reflectivity of the NoDA, 3DVar, and EnSRF experiments against observations at 0100 UTC 30 July 2014. The observed reflectivity at this time showed a strong MCS covering eastern Colorado and southeastern Wyoming, with a convective line at the Colorado–Wyoming border. Some weak convective cells also formed behind the MCS during its southeastward movement (Figure 6(a)). The NoDA experiment clearly overestimated a northward shift of the northern MCS, and it produced only disorganized convection in the southern MCS (Figure 6(b)). Conversely, both 3DVar and EnSRF successfully captured the distribution and location of the MCS, reflecting the positive effect of radar DA. However, 3DVar overpredicted the intensity of convection with reflectivity of 30–40 dBZ (Figure 6(c)), while EnSRF reduced these overestimations, producing results much closer to the observations (Figure 6(d)). The surrounding weaker convection with reflectivity of 20–30 dBZ was also analyzed better by EnSRF, particularly west of Kansas.

Surface temperature and wind field differences among the experiments at 0100 UTC were also pronounced (Figure 7). EnSRF presented the strongest cold pool with a broad area of low temperatures ranging from 10 to 18°C in the convective region, resulting from precipitation evaporation. Note that the solid black line represents the contour of reflectivity > 20 dBZ. The temperatures in the NoDA and 3DVar experiments were 12–20°C and 16–22°C, respectively. The different features of the cold pool between the 3DVar and EnSRF experiments might be due to the lack of cross-variable correlations in the static BEC for hydrometeors in 3DVar [21]. Moreover, EnSRF adjusted the wind field better by producing the strongest wind, which plays an important role in maintaining convection.

4.2. Forecast Result

To assess the overall forecast performances quantitatively, four verification metrics were adopted: the fractions skill score (FSS) [35], equitable threat score (ETS), BIAS, and false alarm ratio (FAR) [36]. The neighborhood-based FSS is defined aswhere and are the forecast and observed fractional coverage of an elementary area by reflectivity or rainfall larger than a given threshold value, respectively, and represents the number of grid points in the verification domain. The FSS was calculated at each forecast time and each grid based on neighborhood sizes of 6 km. It was then averaged over the domain every 10 min.

The ETS calculates the fraction of observed events predicted correctly, while the BIAS and FAR represent the bias and false alarms of the reflectivity or precipitation forecast, respectively. They are defined as follows:where , , and are the numbers of hits, false alarms, and missing grids, respectively, and is the number of hits expected through random chance. BIAS is greater (less) than 1.0 when the frequency of precipitation events, for a given threshold, is over-forecasted (under-forecasted). A perfect forecast would be that the FSS, ETS, and BIAS are equal to 1.0 and FAR is equal to 0.0.

The results of FSS, ETS, BIAS, and FAR for the forecast reflectivity from the NoDA, 3DVar, and EnSRF experiments are compared in Figure 8 for reflectivity thresholds of 5, 15, and 25 dBZ. NoDA had the lowest FSS and ETS for all thresholds during the entire forecast period. 3DVar had substantially higher FSS and ETS than NoDA, but the scores decreased more quickly than EnSRF for all thresholds, except for FSS during the first 10 min and during 160–200 min for 25 dBZ (Figures 8(a)8(f)). The BIAS shown in Figures 8(g)8(i) above or below 1 indicates reflectivity higher or lower than the observation. NoDA generally overestimated reflectivity above 5 dBZ, whereas it was reduced by 3DVar and EnSRF, particularly for the period 160–300 min. This was also true for 15 dBZ at t = 100–300 min and for 25 dBZ at t = 210–300 min. It was also found that EnSRF reduced the overestimation by 3DVar to some extent for all thresholds. 3DVar and EnSRF also showed consistent improvements over NoDA in terms of smaller FAR, while EnSRF showed slightly lower FAR than 3DVar (Figures 8(g)8(l)).

A neighborhood probability (NP; Schwartz et al. [37]) with a radius of 6 km was also calculated to verify the forecast radar reflectivity of the three experiments. Figure 9 shows the NP of reflectivity (>20 dBZ) at 0230, 0300, and 0330 UTC 30 July 2014 for the NoDA, 3DVar, and EnSRF experiments. The NP forecast of reflectivity (>20 dBZ) of NoDA, for which no radar data were assimilated, was less skillful than 3DVar and EnSRF. It produced a large area of probabilities of 0.0 in the southern MCS, indicating convection was not predicted (Figures 9(a)9(c)). The reflectivity forecast was generally improved in the 3DVar experiment; however, it contained some underpredictions in the southern MCS (Figures 9(d)9(f)). EnSRF improved the forecast of the entire MCS compared with the other two experiments. It predicted broad areas of high probability that closely matched the observed reflectivity >20 dBZ in terms of shape, extent, and position (Figures 9(g)9(i)). In addition, the spurious convection along the Oklahoma panhandle, seen in both NoDA and 3DVar, was suppressed by EnSRF. The higher forecast skill of EnSRF was also evident for the threshold of 35 dBZ, which is associated with stronger convection (Figure 10). Although NoDA and 3DVar exhibited greater forecast error than for reflectivity > 20 dBZ, EnSRF had large overlap of high probabilities in the observed reflectivity (>35 dBZ) during the 1 h forecast, particularly at 0330 UTC.

Figure 11 shows the observed composite radar reflectivity as well as the forecast results for the NoDA, 3DVar, and EnSRF experiments at 0230 UTC 30 July 2014. By this time, the observed MCS had developed into a well-defined structure with an extensive area of stratiform precipitation. The convective center of the northern MCS in NoDA was highly overestimated with a northward displacement (Figure 11(b)). 3DVar produced a more realistic spatial pattern of the MCS, but the intensity of convection was overpredicted in the observed convective region (Figure 11(c)). EnSRF continued to present some improvements, including better representation of the stratiform precipitation and weaker reflectivity of 20–30 dBZ. The false storms produced over northwest Colorado by NoDA and 3DVar were also corrected (Figure 11(d)). At 0330 UTC (Figure 12), NoDA (panel b) presented a weaker MCS with a reduced area of high reflectivity (>55 dBZ); however, the main system remained displaced to the north. 3DVar (panel c) successfully predicted two parts of the MCS, and the degree of overestimation was reduced to some extent compared with 0230 UTC. The moderate stratiform precipitation was simulated better by EnSRF (panel d), especially for the southern MCS, although with a larger areal coverage than observed. Unlike the analysis results, the cold pool in 3DVar was stronger than NoDA at this time, while EnSRF still had the strongest cold pool and a near-surface wind that helped maintain a more realistic MCS structure (Figure 13).

Precipitation is another important variable in convective weather. Figure 14 compares the FSS, ETS, BIAS, and FAR for the 5 h forecast from all three experiments. It can be seen that the assimilation of radar data in 3DVar (blue curve) and EnSRF (red curve) resulted in much higher FSS and ETS than NoDA for the thresholds of 1.0, 2.5, and 5.0 mm (1 h)−1 during the entire forecast period. EnSRF showed relatively higher FSS and ETS than 3DVar, and the difference was larger in terms of FSS for the threshold of 5.0 mm (1 h)−1 at t = 60–180 min. EnSRF reduced the dry bias from 3DVar and NoDA for the threshold of 1.0 mm (1 h)−1, and the BIAS of EnSRF (which was close to 1) was between 3DVar and NoDA for the thresholds of 2.5 and 5.0 mm (1 h)−1. NoDA presented the largest FAR, followed in descending order by 3DVar and EnSRF.

Figure 15 compares the forecast 1 h accumulated precipitation of the NoDA, 3DVar, and EnSRF experiments. 3DVar and EnSRF were noticeably superior to NoDA in terms of spatial pattern and amount but with different strengths. 3DVar overpredicted the MCS and it produced a larger area of heavy precipitation (>25.6 mm); however, it also missed some moderate (1.6–12.8 mm) and weaker precipitation (0.2–1.6 mm). In contrast, EnSRF predicted a broader area of precipitation than the other experiments, which agreed better with the observations, particularly for precipitation in the range 0.2–1.6 mm. For the 3 h accumulated precipitation (Figure not shown), the observations produced a wider rainband with a larger area of heavy rain. 3DVar again produced too much heavy rain, whereas the EnSRF simulated precipitation had weaker amplitude but with wider coverage, as in the observations.

Figure 16 shows the vertical distribution of domain-averaged root mean square errors (RMSEs) of the horizontal wind components, temperature, and water vapor forecasts against all 10 radiosondes from the 1 km domain at 0300 UTC 30 July 2014. In comparison with NoDA and 3DVar, EnSRF (red curve) had slightly lower RMSEs for horizontal wind components and temperature below 70 hPa, and much smaller RMSEs for water vapor below 500 hPa. Figure 17 shows the RMSEs of horizontal wind components, temperature, and water vapor analyses against the data from all 781 surface METAR stations within the 1 km domain at 0300 UTC 30 July 2014. Evidently, EnSRF had the smallest RMSEs, while the errors of NoDA were the largest and significantly reduced by 3DVar and EnSRF for all variables.

5. Summary

In this study, the EnSRF and a recently developed indirect 3DVar method were used to assimilate radar data for an MCS that occurred over the Front Range of the Rocky Mountains in Colorado (USA) on 29–30 July 2014. Three experiments (NoDA, 3DVar, and EnSRF) were compared to investigate the impact of radar DA and different assimilation methods on both the analysis and the subsequent reflectivity and precipitation forecasts. Both reflectivity and radial velocity data were assimilated from nine radars of the Weather Surveillance Radars-1988 Doppler network in the 3 km domain.

The analysis results showed that EnSRF reduced the RMSIs compared with 3DVar and that its ensemble spread and RMSI values were of comparable magnitude for both reflectivity and radial velocity in the analysis–forecast cycles. Assimilation of radar data resulted in considerable improvement of the analysis, and the analysis produced by EnSRF was more realistic than 3DVar in terms of intensity. EnSRF showed a stronger surface cold pool and wind outflow than 3DVar in the convective regions.

The positive impact of assimilating radar data was maintained for the forecast period, and EnSRF improved the quantitative reflectivity forecast skills measured by FSS, ETS, BIAS, and FAR over 3DVar. The NP of the reflectivity forecasts indicated that EnSRF produced the most skillful probabilistic forecast with high probabilities of the observed reflectivity with different thresholds that included moderate and strong convection. Although 3DVar improved the prediction of the extent of the MCS compared with NoDA, it was no better than EnSRF, especially for the southern MCS. The MCS structure was improved by both 3DVar and EnSRF in terms of location, intensity, and areal coverage, with EnSRF being the best. For precipitation, EnSRF also increased FSS and ETS and reduced FAR compared with NoDA and 3DVar, and it consistently produced the best BIAS for the different thresholds. The precipitation patterns and locations were represented well by 3DVar and EnSRF, although EnSRF reduced the overestimation of the heavy rainfall center compared with 3DVar. Verification against radiosonde sounding data in the 1 km domain suggested that the EnSRF reduced the RMSE of wind, temperature, and water vapor mixing ratio in comparison with NoDA and 3DVar.

The present study demonstrated the encouraging results of EnSRF over 3DVar in MCS analysis and forecasting; however, further research involving additional cases and longer forecast periods will be necessary to reach conclusions that are more reliable. Further studies on radar DA using more advanced methods such as 4DVar and hybrid DA are needed.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This research was supported by the National Key Research and Development Program of China (Grant no. 2017YFC1502102) and the Research Innovation Program for College Graduates of Jiangsu Province (Grant no. KYLX_0829). The author Shibo Gao is supported by scholarship from the China Scholarship Council.