Achieving accurate storm-scale analyses and reducing the spinup time of modeled convection is a primary motivation for the assimilation of radar reflectivity data. One common technique of reflectivity data assimilation is using a cloud analysis, which inserts temperature and moisture increments and hydrometeors deduced from radar reflectivity via empirical relations to induce and sustain updraft circulations. Polarimetric radar data have the ability to provide enhanced insight into the microphysical and dynamic structure of convection. Thus far, however, relatively little has been done to leverage these data for numerical weather prediction. In this study, the Advanced Regional Prediction System’s cloud analysis is modified from its original reflectivity-based formulation to provide moisture and latent heat adjustments based on the detection of differential reflectivity columns, which can serve as proxies for updrafts in deep moist convection and, subsequently, areas of saturation and latent heat release. Cycled model runs using both the original cloud analysis and above modifications are performed for two high-impact weather cases: the 19 May 2013 central Oklahoma tornadic supercells and the 25 May 2016 north-central Kansas tornadic supercell. The analyses and forecasts of convection qualitatively and quantitatively improve in both cases, including more coherent analyzed updrafts, more realistic forecast reflectivity structures, a better correspondence between forecast updraft helicity tracks and radar-derived rotation tracks, and improved frequency biases and equitable threat scores for reflectivity. Based on these encouraging results, further exploration of the assimilation of dual-polarization radar data into storm-scale models is warranted.
The assimilation of radar data into convective-scale numerical weather prediction (NWP) models has gained considerable attention in recent years with increased operational implementation and use in the warning decision process (e.g., the “Warn-on-Forecast” initiative; Stensrud et al. 2009, 2013). Weather radar is one of the only sources of data available that provides information at a temporal and spatial resolution comparable to convection-allowing NWP models. Thus, radar data assimilation’s ability to promote deep moist convection and its attendant perturbations in NWP models and to reduce the spinup time has been and continues to be explored. However, most efforts to date have been limited to measurements of reflectivity at horizontal polarization (hereafter, Z) and radial velocity.
Many methods of assimilating Z have been studied over the past two decades. One of the earliest methods examined was four-dimensional variational data assimilation (4DVAR), which uses the forecast model as a dynamical constraint during the assimilation process. While results have been encouraging (e.g., Sun and Crook 1997, 1998; Sun 2005; Sun and Zhang 2008; Wang et al. 2013b), the difficulty of developing and maintaining an adjoint model and the inherent nonlinearities of the microphysics scheme often hinder proper convergence of the cost function. As such, 4DVAR methods have not been widely used for convective-scale radar data assimilation and have typically been limited to warm rain microphysics, although recent work has begun to investigate the inclusion of some ice phases (Chang et al. 2016). The simpler and more computationally efficient three-dimensional variational data assimilation (3DVAR) method has provided positive results (e.g., Xiao et al. 2005, 2007; Gao and Stensrud 2012; Wang et al. 2013a). However, 3DVAR lacks flow-dependent error covariances, which may limit the ability to update unobserved variables, and requires limiting assumptions about the model microphysics. In recent years, the ensemble Kalman filter method (EnKF; Evensen 1994) has become increasingly popular for convective-scale radar data assimilation with very promising results for producing accurate storm-scale analyses (e.g., Dowell et al. 2004, 2011; Tong and Xue 2005; Xue et al. 2006; Aksoy et al. 2009; Yussouf and Stensrud 2010; Snook et al. 2011, 2012; Yussouf et al. 2013; Wheatley et al. 2015). However, ensemble methods are computationally expensive, may suffer from issues related to rank deficiency (e.g., filter divergence due to sampling errors; Gao et al. 2014), and have not yet seen widespread operational implementation. Hybrid methods combining the strengths of variational and ensemble methods by defining ensemble-derived flow-dependent covariances for the variational scheme are being developed (e.g., Wang et al. 2008a,b) and investigated for use with radar data assimilation at the convective scale (e.g., Gao et al. 2013, 2016; Gao and Stensrud 2014).
Other assimilation techniques assimilate state variables indirectly retrieved from Z. These methods include latent heating nudging (e.g., Jones and Macpherson 1997; Macpherson 2001; Leuenberger and Rossa 2007; Stephan et al. 2008), specific humidity nudging (Davolio and Buzzi 2004), and divergence nudging (Korsholm et al. 2015), as well as techniques that use regions of observed Z to activate a convective parameterization scheme (Rogers et al. 2000) or variationally assimilate retrieved relative humidity profiles (Caumont et al. 2010). One of the most prominent methods is the Advanced Regional Prediction System’s (ARPS) cloud analysis (hereafter “cloud analysis”; Zhang et al. 1998; Zhang 1999; Brewster 2002; Hu et al. 2006a). The cloud analysis is based on the Local Analysis Prediction System (Albers et al. 1996) and makes adjustments to the model relative humidity, hydrometeor mixing ratios, and temperature based on radar, satellite, and surface observation data. Cloud analysis techniques are conceptually straightforward, computationally efficient, and have been shown to be useful for reducing the spinup of observed storms and improving short-term convective forecasts (e.g., Xue et al. 2003, 2014; Souto et al. 2003; Dawson and Xue 2006; Hu et al. 2006a; Zhao and Xue 2009; Schenkman et al. 2011; Dawson et al. 2015; Zhuang et al. 2016). Benefits can be amplified when the cloud analysis is used in conjunction with radial velocity information. To update multiple unobserved model variables from Z alone, however, cloud analysis techniques rely on semiempirical quantitative relations (e.g., retrieving hydrometeor mixing ratios from Z) and general rules relating Z to the aforementioned variables (e.g., saturating regions within a given Z threshold). These relations and rules require simplifications that can introduce errors (e.g., Gao et al. 2009; Carlin et al. 2016).
The operational WSR-88D network in the United States has been upgraded to dual-polarization, with other countries, including Germany, Canada, and the United Kingdom, following suit. These networks now provide an unprecedented volume of polarimetric observations. In contrast to single-polarization radars, dual-polarization radars transmit and receive orthogonally polarized electromagnetic waves from which information about a target’s size, shape, orientation, and composition can be garnered (e.g., Kumjian 2013a). In addition to Z, measured variables include differential reflectivity , copolar correlation coefficient , and differential phase shift . Dual-polarization radar data have been successfully leveraged to improve attenuation correction (e.g., Bringi et al. 1990; Testud et al. 2000; Snyder et al. 2010), quantitative precipitation estimation (e.g., Zrnić and Ryzhkov 1996; Ryzhkov et al. 2005a; Tabary et al. 2011), hydrometeor classification (e.g., Park et al. 2009), tornado detection (e.g., Ryzhkov et al. 2005b; Schultz et al. 2012a,b; Bodine et al. 2013, 2014; Van Den Broeke and Jauernic 2014; Snyder and Ryzhkov 2015), and the identification (e.g., Heinselman and Ryzhkov 2006) and size discrimination (e.g., Ryzhkov et al. 2013a,b) of hail. For further review of weather radar polarimetry, see Zrnić and Ryzhkov (1999), Bringi and Chandrasekar (2001), and Kumjian (2013a,b,c).
In recent years, numerous distinct polarimetric “signatures” have been identified and tied to dynamical and microphysical processes within storms. One of the most ubiquitous polarimetric signatures observed in deep moist convection is the column. Differential reflectivity columns are vertical protrusions of positive above the environmental 0°C level and are indicative of wet ice particles and oblate, supercooled raindrops in the process of freezing being lofted by the updraft. Values of within these columns can exceed 4 dB at S band (radar wavelengths of approximately 10 cm) and can reach beyond 3 km above the 0°C level in extreme cases (Kumjian et al. 2014; Snyder et al. 2015). Because columns are associated with convective storm updrafts, they can theoretically be used as identifiers for regions of positive temperature perturbations from latent heat release due to condensation and/or freezing. Although the connection between columns and updraft location has been long known (e.g., Hall et al. 1984; Illingworth et al. 1987; Tuttle et al. 1989; Ryzhkov et al. 1994), recent work has begun to investigate the relationship between columns and updraft intensity. Simulations in Kumjian et al. (2014) showed a relationship between column depth (i.e., the distance above the 0°C level that enhanced values of extend within the column) and updraft strength. Both Picca et al. (2010) and Kumjian et al. (2014) showed a correlation between column height and hail mass at the surface at appreciable lag times, which has the potential to provide increased lead time for forecasting hail at the surface as compared to traditional metrics such as 20-dBZ echo-top height.
In addition to the aforementioned columns associated with mature updrafts, enhanced coincident with low values of Z has been observed aloft in the nascent stages of developing convection (e.g., Illingworth et al. 1987; Knight 2006) due to the rapid size sorting of initial precipitation particles (Picca et al. 2017). While technically distinct from columns, these bursts of enhanced aloft in the updraft may prove useful for identifying developing updrafts. For the sake of simplicity, both of these phenomena will be referred to as “ columns.” For a more complete review of columns, see Kumjian et al. (2014) and Snyder et al. (2015).
Despite the connection between dual-polarization radar and the microphysical and thermodynamic characteristics of deep moist convection, leveraging polarimetric data for NWP is a relatively new area of research. Predicated on the idea that a physically accurate microphysics scheme should be able to reproduce realistic polarimetric signatures, many studies have explored the use of polarimetric radar forward operators (e.g., Jung et al. 2008a, 2010a; Pfeifer et al. 2008; Ryzhkov et al. 2011) to evaluate the performance of microphysics schemes (e.g., Jung et al. 2008a, 2010a, 2012; Ryzhkov et al. 2011, 2013a; Kumjian and Ryzhkov 2012; Dawson et al. 2013, 2014; Kumjian et al. 2014; Putnam et al. 2014; Johnson et al. 2016; Snyder et al. 2017a,b). If large discrepancies are present between the model-derived polarimetric signatures and those observed in nature, it can be indicative of deficiencies in the microphysics scheme. Alternatively, if a model faithfully reproduces polarimetric signatures as they are observed in nature, the model can be used to investigate what physical processes are responsible for a given signature. Some studies have used polarimetric data to assimilate improved estimates of rainwater mixing ratio using both 3DVAR (Li and Mecikalski 2010, 2012) and EnKF (Yokota et al. 2016) methods and found positive impacts compared to experiments that assimilated mixing ratios retrieved from Z alone. However, ice phases were neglected. Wu et al. (2000) attempted to use to differentiate between liquid and ice phases for hydrometeor mass retrievals using 4DVAR but found little success attributed to inadequate model physics. Few studies have attempted to directly assimilate polarimetric data. Using an EnKF framework, simulated polarimetric data were assimilated in addition to Z to estimate state variables (Jung et al. 2008a,b) and microphysical parameters (Jung et al. 2010b), with positive impacts found in both cases. Assimilating observed polarimetric data remains difficult because of data quality concerns and uncertainties in polarimetric operators.
This study explores the impact of assimilating observed polarimetric data through a modified cloud analysis routine. The cloud analysis technique was chosen because of its proven success in reducing spinup time and ease of implementation into existing code infrastructure. Direct insertion of the retrieved temperature and moisture perturbations is currently more straightforward than assimilating the polarimetric variables using variational techniques, which require cross covariances between model state variables and the polarimetric variables that are not currently well formulated. Section 2 details the modifications made to the existing cloud analysis routine, and section 3 describes the experimental setup used in this study. Results are presented in section 4, followed by a summary and discussion in section 5.
2. Description of data assimilation routine
a. ARPS 3DVAR routine
The first step of the assimilation cycling procedure makes use of the ARPS 3DVAR routine (Gao et al. 2004; Hu et al. 2006b). ARPS 3DVAR minimizes a cost function with a recursive filter containing terms for the background and observations as well as an anelastic mass continuity term as a weak constraint to produce a more balanced three-dimensional analysis of the model state variables from multiple data sources. The system is designed to work with a number of observation types including surface and upper-air observations and radial velocity. As it was designed for use at the storm scale, the routine includes multiple analysis passes with varying scales of spatial influence to help resolve flows at different scales. The resultant analysis is then used as the background when invoking the cloud analysis routine.
b. Existing cloud analysis
In the current cloud analysis, the radar data are first quality controlled and interpolated to the model grid (Brewster et al. 2005; Brewster and Stratman 2015). The process for the polarimetric variables follows that of radial velocity and Z. Radar data at ranges between 3 and 230 km from the radar site are processed. Anomalous propagation is removed using gradients and texture fields of Z and low wind speeds, with additional filtering of nonmeteorological echoes using a user-defined threshold (default of 0.85 for S band). All variables are smoothed using a nine-point median filter. Specific differential phase shift is calculated from from a local least squares fit on smoothed data using a Z-dependent averaging window. Remapping to the model grid of all radar variables is performed using a least squares fit to a local polynomial function, which thins the data near the radar and interpolates it at distances far from the radar (Brewster et al. 2005). For the latest changes to the cloud analysis, see Brewster and Stratman (2015).
An initial cloud fraction field is diagnosed from the background relative humidity field following a similar approach as Koch et al. (1997). Subsequently, clouds are directly inserted by setting the cloud fraction to 100% above the surface-based lifted condensation level anywhere Z exceeds a threshold, set to 15 dBZ above 2 km by default. Cloud water and ice content can be determined either adiabatically or, as in this case, using the Smith–Feddes model (Haines et al. 1989) with a reduction for entrainment following Hu et al. (2006a). Next, the dominant hydrometeor species in each grid box is determined using temperature and Z thresholds, where snow (rain) is considered when the temperature is below (above) 0°C, and where hail is considered when the Z exceeds 45 dBZ (Albers et al. 1996; Pan et al. 2016). In the case of cycling, the species can also be determined by the existing species in the model background. The mixing ratios of each species are then typically retrieved using single-moment retrieval equations for rain, snow, and hail based on Smith et al. (1975) and Lin et al. (1983). Summaries of these equations can be found in Dowell et al. (2011), Carlin et al. (2016), and Pan et al. (2016). However, recent work has initialized intercept parameters (and, if needed, shape parameters) for multimoment microphysics schemes using iterative techniques (Brewster and Stratman 2015), while other studies have found positive impacts from using single-moment microphysics schemes with intercept parameters diagnosed from hydrometeor mixing ratios (e.g., Wainwright et al. 2014; Pan et al. 2016), as developed in Zhang et al. (2008).
A temperature adjustment is then made to account for latent heat release. This can be done by simply adding the latent heating associated with the added cloud water and ice content (Zhang et al. 1998) or by assuming a moist-adiabatic temperature profile from cloud base with entrainment effects included (Brewster 2002). In this study, the latter method is applied to regions with vertical velocity w −0.2 m s−1 (determined from the 3DVAR analysis) with a linear ramp from no heating to full heating between −0.2 and 0.0 m s−1. Final moisture adjustments are made by reestablishing saturation anywhere the Z threshold for clouds is exceeded (incorporating the previously made temperature adjustment) or to 95% anywhere the analyzed hydrometeor mass is less than the background hydrometeor mass to help avoid overmoistening. Further details of the cloud analysis and its latest updates can be found in Brewster and Stratman (2015) and Tong (2015).
c. Modified cloud analysis
The modifications made to the cloud analysis in this study involve the final two steps of moistening and heating in updraft areas. Many studies have shown that both temperature perturbations (e.g., Hu et al. 2006a) and the initial moisture field (e.g., Weygandt et al. 2002; Ge et al. 2013) can play primary roles in determining the accuracy of modeled convection. The insertion of too much water vapor can result in an overestimate of the intensity and areal coverage of convection, leading to a degradation of the forecast (e.g., Schenkman et al. 2011; Schenkman 2012). This issue was examined in detail in Tong (2015), who found that saturating based on a Z threshold can result in too much moisture being added and large degradations in forecast skill. Forecast skill was greatly improved when a more accurate initial moisture field was provided in an observing system simulation experiment. Because of the lack of a direct relationship between in-cloud moisture and conventional observations, Tong (2015) proposed a modification to the cloud analysis in which the relative humidity in downdraft regions, which are generally unsaturated, is reduced. Notable improvements were found for both the analysis and forecast for all state variables examined, further highlighting the importance of improving the initial moisture field for convective storm-scale modeling. Despite these encouraging results, certain issues remain. While unsaturated regions correspond well with downdrafts overall, the specific quantitative relationship between water vapor mixing ratio and w is unknown and poorly constrained. In addition, even with a perfect qυ–w relationship, the success of this method relies on an accurate model analysis of w, which is not always known and/or guaranteed, particularly when few radars are available for assimilation. As an alternative to using w, a method is proposed using columns to provide adjustments to temperature and moisture in the cloud analysis similar to the methods for assimilating lightning data at the cloud scale put forth by Fierro et al. (2012, 2014, 2015) and Marchand and Fuelberg (2014).
To investigate the validity of the proposed modifications, vertical cross sections of relative humidity, latent heating rate, , Z, and storm-relative winds from a convective storm simulated by the Hebrew University Cloud Model (HUCM) are shown in Fig. 1. The HUCM is a two-dimensional nonhydrostatic cloud model with spectral bin microphysics (Khain et al. 2004). The model is coupled to a sophisticated polarimetric radar operator (Ryzhkov et al. 2011)—at S band in this case—that has been shown to recreate realistic polarimetric signatures in deep moist convection (e.g., Ryzhkov et al. 2013a; Kumjian et al. 2014; Snyder et al. 2015). The simulation shown in Fig. 1 was initialized with a sounding from a damaging hailstorm that struck Germany in 2006 (Noppel et al. 2010) and run with a 100-m vertical and 300-m horizontal grid spacing. The initial low-level cloud condensation nuclei concentration was 3000 cm−3, representing polluted conditions. Throughout the lifetime of the storm, the columns are coincident with updrafts featuring deep plumes of saturation and with the region of latent heating directly above the columns. Notably, the area contained within the 15-dBZ contour is much more extensive than the areas that are near or at saturation, with large regions exhibiting subsaturation. It is clear that saturating everywhere within the 15-dBZ threshold would result in too much moisture being added to the system, possibly inhibiting the formation of evaporation- or sublimation-driven downdrafts. These results support the conceptual model of columns and their use as proxies for updrafts and subsequently areas of moistening and heating. It should be noted, however, that while 15 dBZ is the default threshold for saturation in the cloud analysis, it is an adjustable parameter and there is no agreed-upon Z threshold to use. Other studies have addressed overmoistening concerns by instead reducing the frequency of applications of moistening (e.g., Schenkman et al. 2011).
In this study, polarimetric data are first quality controlled and mapped to the model grid as described in section 2b. Areas of interest are limited to regions in which Z 10 dBZ and 0.85 to ensure sufficient signal-to-noise ratio and good-quality data, and to regions below the environmental −20°C level to mitigate the chance of ice crystals with enhanced causing false detections. Similar to the criteria used in Snyder et al. (2015) for their column detection algorithm, a column is defined here to exist if 1.0 dB for at least two vertically contiguous grid boxes above the environmental 0°C level. To help ensure that only legitimate columns are detected and limit the chance of noise in the field causing false detections, an additional 3 km × 3 km horizontal mode filter is incorporated in which only columns exhibiting 0.85 and 1.0 5.0 dB in at least five of the nine grid boxes within the filter are counted. A summary of these detection criteria is shown in Table 1.
As opposed to warming in areas with w −0.2 m s−1 as in the existing cloud analysis, temperature adjustments are instead made anywhere columns are detected. Adjustments are made both where columns are located and to one grid box surrounding the columns to aid in establishing wide-enough updrafts that do not mix out before becoming established. Similarly, instead of saturating based on a simple Z criterion, saturation is only applied to the model columns (within the cloud region as determined by Z) where columns have been detected. Model columns surrounding detected columns are also saturated, with the horizontal extent proportional to the detected depth of the columns (in this case, half the number of model levels in the detected columns) to prevent the added moisture from mixing out and to attempt to add more moisture for “stronger” (i.e., taller) columns. In addition to moistening and heating at observed column locations, an additional drying procedure is applied in an attempt to mitigate possible overmoistening by the microphysics scheme. At any locations satisfying Z 10 dBZ and relative humidity 80% but with no detected column, the relative humidity is reduced by half of the excess relative humidity above 80% (e.g., if the relative humidity is 90% with no column detected, the relative humidity is reduced to 85%). This is, admittedly, an arbitrary process, but remains a succinct way to provide minor drying to areas characterized by precipitation (sufficient to meet the Z 10-dBZ criterion) but that are outside of columns (where, we hypothesize, deep updrafts are less likely). Future work will attempt to examine the sensitivity to the column detection criteria and the details of the filtering and weighting procedures for moistening and drying.
An example of the differences in potential temperature and water vapor mixing ratio analysis increments between the traditional cloud analysis and the modified cloud analysis is shown in Fig. 2 for the initial 2000 UTC assimilation cycle of the 19 May 2013 Oklahoma case (discussed below). While the magnitudes of the moistening and warming are comparable, the location and extent of the increments vary between the two. The traditional cloud analysis (Fig. 2a) shows a large area of moistening with two primary areas of warming west-northwest of Oklahoma City associated with the first developing supercell, and smaller areas of moistening and warming northwest and west-southwest of Oklahoma City. In contrast, the modified cloud analysis employing detected columns (Fig. 2b) shows a smaller area of moistening and warming directly west of Oklahoma City, southwest of the area modified in the traditional cloud analysis, and with little/no moistening or warming elsewhere. The only exception is far southwest of Oklahoma City, where the modified cloud analysis shows a bit more moistening associated with developing convection than the traditional cloud analysis. While the differences in analysis increments between the traditional and modified cloud analyses vary with time, Fig. 2 provides a demonstrative example of the typical differences between the methods.
3. Experimental setup
To investigate the impact of the modified cloud analysis, two tornadic supercell events were studied: the 19 May 2013 tornado outbreak in central Oklahoma (“the OK case”) and the tornadic supercell of 25 May 2016 in north-central Kansas (“the KS case”).
a. Case descriptions
Around 2000 UTC 19 May 2013, thunderstorms initiated near a dryline just west of the Oklahoma City, Oklahoma, metropolitan area in an environment characterized by strong vertical wind shear and high potential convective instability (i.e., CAPE). These storms developed quickly, and the three supercells that emerged from the convection moved toward the east-northeast; two of the supercells produced a total of eight tornadoes, whereas the third supercell was not tornadic. The northernmost supercell produced two brief tornadoes north and northeast of Oklahoma City before producing a long-lived tornado that produced EF3 damage near Carney, Oklahoma, between 2141 and 2224 UTC that resulted in 4 injuries; the southernmost supercell spawned a tornado that produced EF4 damage near Shawnee, Oklahoma, between 2300 and 2350 UTC that resulted in 2 fatalities and 10 injuries (NWS 2017a).
In the KS case, an isolated supercell formed in north-central Kansas just north of a warm front around 2200 UTC 25 May 2016 and moved slowly east-southeastward. The storm produced a total of four tornadoes, including a long-track tornado just east-northeast of Salina, Kansas, that lasted over 1.5 h (0007–0140 UTC) (NWS 2017b).
For both cases, observed tornado tracks are retrieved from shapefiles created from damage survey reports.
b. Model setup
The model used in this study is the ARPS (Xue et al. 2000, 2001, 2003), a nonhydrostatic, compressible, numerical model designed to function at multiple scales with an emphasis on the explicit prediction of convection. Terrain data were derived from the U.S. Geological Survey 3-arc-s dataset. Subgrid-scale turbulence was parameterized using a 1.5-order TKE turbulence scheme, with the evolution of the planetary boundary layer using the formulation of Sun and Chang (1986). Cloud microphysics were parameterized using the Milbrandt–Yau double-moment scheme (Milbrandt and Yau 2005a,b), and both short- and longwave radiation were parameterized using the NASA Goddard schemes (Chou 1990, 1992). A two-layer force-restore soil model based on Noilhan and Planton (1989) was used with surface fluxes based on stability-dependent drag coefficients using surface temperature and volumetric water content. More information about the full ARPS physics suite can be found in Xue et al. (2001).
Experiments were conducted using a one-way nested grid configuration. The parent domain has a size of 1200 km × 1200 km with a horizontal grid spacing of 4 km, and the inner nest a size of 500 km × 500 km with a horizontal grid spacing of 1 km. The domains for the OK and KS cases were centered on 35.45°N, 97.25°W and 38.65°N, 97.55°W, respectively. Both nests used a stretched vertical grid containing 53 vertical levels with an average spacing of 400 m and a minimum spacing of 100 m near the surface. The model top was rigid with a Rayleigh damping layer above 12 km to absorb vertically propagating waves. Lateral boundary conditions were externally forced. The simulated Z fields were computed using the T-matrix-based algorithm of Jung et al. (2010a). The domains used for each case are shown in Fig. 3, and a summary of the model setup used for these experiments is provided in Table 2.
c. Assimilation procedures
The 12-km North American Mesoscale Forecast System (NAM) model analysis and forecast data were used to initialize the parent domain. For the OK case, the 1800 UTC 19 May 2013 NAM analysis was used, and for the KS case the 2-h forecast from the 1800 UTC 25 May 2016 NAM analysis (valid at 2000 UTC) was used. The NAM data were interpolated onto the 4-km ARPS grid, which was then integrated forward for 1 h using 3-h lateral boundary conditions derived from the NAM. This forecast was then further interpolated down to the inner nest and integrated forward another 1 h, with boundary conditions on the inner nest updated at 30-min intervals from the outer nest, for a total spinup period of 2 h. This forecast was then used as the background for all assimilation experiments performed.
Assimilation cycles were performed every 10 min following Hu and Xue (2007), who found this to be the optimal cycling frequency in their experiments. Radial velocity data were assimilated using the ARPS 3DVAR routine (Gao et al. 2004; Hu et al. 2006b), after which the cloud analysis routine was called. For the OK case, Oklahoma Mesonet data (Brock et al. 1995) were also assimilated using the 3DVAR routine. After 30 min, a separate 1-h forecast was made, with 10-min assimilation cycles continuing. Then 1-h forecasts were subsequently initiated every 30-min for 3 h after the initial analysis time. A diagram of the spinup, cycling, and assimilation process is shown in Fig. 4. For the OK case, radar data from the Twin Lakes, Oklahoma, WSR-88D (KTLX) were used, while the KS case used data from the Topeka, Kansas, WSR-88D (KTWX) (Fig. 3). For each case, two runs were performed: a control run (hereafter “Control”), in which the legacy cloud analysis is used (see section 2b), and an experimental run (hereafter “ZDRCOL”), which employed the modified polarimetric cloud analysis described in section 2c.
Specific nomenclature for each experiment will be referred to hereafter by their case and which cloud analysis method was used (i.e., “KS_ZDRCOL” refers to the 25 May 2016 KS case experiment employing the modified cloud analysis).
a. 19 May 2013 case
To investigate the performance of the column detection algorithm, a composite plot of remapped 1-km Z and analyzed column depth from the KTLX radar observations in 10-min intervals for the assimilation period is shown in Fig. 5. Each of the three swaths associated with a supercell is labeled and will be used to reference each storm in the subsequent discussions. The 15-dBZ contour is shown as that is the default threshold for saturation in the original cloud analysis routine. Distinct column tracks are evident for all three main storms, with all storms exhibiting prominent columns during their formative stages before becoming more intermittent, supporting the use of columns for spinning up storms in the model early in their life cycle. In each case, the column is found on the southwest flank of the storm where the main updraft is expected to be located. The only exception to this is for the weakening and fast-moving cell north of Oklahoma City that propagates to the left of the mean wind off to the north-northeast. The supercell that begins to the west of Oklahoma City (“Supercell 1” in Fig. 5) exhibits a large, deep column from its inception that travels toward the northeast and then turns to the east-northeast before producing the first tornado. The column then shrinks and becomes shallower near and after the first tornado dissipates; the column associated with the main updraft becomes more robust shortly after the genesis of the second tornado east-northeast of Oklahoma City, near the observed track. Two more column tracks are apparent south of Oklahoma City, with the middle track (“Supercell 2”) associated with a smaller column as it tracked northeast. The southernmost storm (“Supercell 3”) exhibited a larger and taller column that suddenly weakened and never fully reappeared. The period analyzed here ends at 2300 UTC, the approximate start time of the long-track southern tornado southeast of Oklahoma City. However, no clear column is evident here due to the close proximity of the updraft to the radar (i.e., the column likely was located within the cone of silence, which extends 12 km from the radar at a height of 4 km AGL), although additional obfuscation by hail or tornadic debris cannot be ruled out. This is a known drawback that should be taken into consideration when using any methods that use vertically integrated data or echo-top heights from a single radar.
The areas encompassed by the 15-dBZ threshold are much larger and extend farther to the north and east of the analyzed columns. Note that for this and all subsequent figures for the OK case the two easternmost tornadoes (shown in gray in Fig. 5) occurred after the period examined in this experiment.
Composite plots of the maximum analyzed w (contoured at 30 m s−1) at each grid point for both OK_Control and OK_ZDRCOL through the assimilation period (2000–2300 UTC) are shown in Fig. 6. OK_Control exhibits a rather noisy w field composed of many spurious updrafts, along with a pronounced northward bias compared to the observed tornado tracks. This northern and positive forward speed bias has been observed in many storm-scale modeling studies (e.g., Potvin et al. 2014; Xue et al. 2014; Stratman and Brewster 2015; Wheatley et al. 2015). In sharp contrast, OK_ZDRCOL features much more consolidated updraft tracks that closely follow the analyzed column paths (and observed tornado tracks). As in Fig. 5, the final analysis included is at 2300 UTC near the beginning of the long-track tornado southeast of Oklahoma City, evident with a large and strong updraft in excess of 40 m s−1 near the start of the tornado track that is not as apparent in OK_Control. The 30 m s−1 contours also appear to be larger in OK_ZDRCOL than in OK_Control, suggesting wider, stronger updrafts.
The composited 1–6 km above ground level (AGL) updraft helicity (Kain et al. 2008) swaths for three different forecast periods are shown in Fig. 7. Model output was saved every 5 min and composited over the 1-h forecast, with the maximum for the forecast period shown at each grid point. The 1–6-km updraft helicity provides a reasonable depiction of the path of mesocyclones and overall storm track. To aid in verifying the forecast updraft helicity swaths, rotation tracks derived from the Multi-Radar Multi-Sensor (MRMS; Smith et al. 2016) system, which are composited maximum values of radar-derived azimuthal shear (Smith and Elmore 2004) in a layer through a given time period, are included in Fig. 7. While the traditional azimuthal shear product uses 0–2- or 3–6-km AGL layers, the 1–6-km AGL azimuthal shear was used in this study to better correspond with the 1–6-km updraft helicity derived from the model output. The rotation tracks shown in Fig. 7 correspond to the 1-h forecast periods shown in each panel.
During the first 0–1-h forecast at 2030 UTC (Figs. 7a,b), both OK_Control and OK_ZDRCOL feature a storm track for Supercell 1 that is located too far north. The updraft helicity swath in OK_ZDRCOL, however, is more consolidated and features a smaller northward bias compared to OK_Control. Supercell 1 in OK_ZDRCOL has a slower mean storm motion, with the center of the updraft helicity swath covering approximately 10 fewer kilometers than OK_Control during the forecast period. Finally, OK_ZDRCOL features a weak updraft helicity swath associated with the second developing storm (Supercell 2, southwest of Oklahoma City) that is absent in the OK_Control run.
The improvements of OK_ZDRCOL over OK_Control are most pronounced in the forecast initiated at 2130 UTC (Figs. 7c,d), approximately 10 min before the start of the long-track tornado northeast of Oklahoma City. OK_Control features multiple updraft helicity swaths. There is no identifiable strong updraft helicity swath coincident with the observed rotation track of Supercell 1, with instead a very strong and prominent updraft helicty swath displaced far to the northeast of the observed rotation track and corresponding tornado. Moreover, there are two notable updraft helicity swaths corresponding to the weakening rotation track of Supercell 2 southeast of Oklahoma City, with no updraft helicity swath that clearly corresponds with the rotation track for Supercell 3. In stark contrast, OK_ZDRCOL captures the updraft helicty swath of Supercell 1 well, with the forecast swath nearly coincident with the observed rotation track and with only a slight bias in forward speed. It also correctly captures the updraft helicity swath associated with Supercell 2 that weakens as it moves to the northeast. Finally, the early development of strong rotation in the southernmost supercell (Supercell 3) that would go on to produce the Shawnee tornado is depicted to the south of Oklahoma City while being absent in OK_Control.
The 2230 UTC 0–1-h forecast (Figs. 7e,f) show many of the same improvements. Both Supercells 1 and 2 were nontornadic and beginning to weaken, with less pronounced updraft helicity swaths in OK_ZDRCOL. In contrast, OK_Control has strong but noisy updraft helicity swaths for these storms displaced to the northeast of their observed locations. For Supercell 3, both OK_Control and OK_ZDRCOL exhibit updraft helicity associated with the strong and broad observed rotation track south of Oklahoma City. However, the updraft helicity swath in OK_Control is primarily north and east of the observed rotation track, while OK_ZDRCOL captures the rotation (albeit with a slight north bias) and its timing well.
To further examine the improvements in the OK_ZDRCOL forecasts over OK_Control, the 1-km AGL Z is shown for the forecasts initiated at 2130 UTC in 20-min increments and compared to the observed radar fields in Fig. 8. This time period represents the duration of the northern long-track tornado northeast of Oklahoma City, which was on the ground between 2141 and 2224 UTC, as well as the leadup period to the long-track tornado produced by Supercell 3, which first touched down at 2300 UTC. For both the OK_Control and OK_ZDRCOL runs, an adjustment period is seen in the first 20 min (Figs. 8e,f) with small, yet intense, precipitation cores (Z 65 dBZ) present. These high values of Z occur within the core of the middle and northern storms (Supercells 1 and 2) in OK_Control (Fig. 8e), while in OK_ZDRCOL these high Z values are predominantly near the southern flank of the storms and/or within the hook echoes, where the columns were analyzed. Later in the forecast period (2210–2230 UTC), it is again clear that the OK_Control run features a northward and positive forward speed bias (Figs. 8h,k) compared to the observations (Figs. 8g,j). Supercell 2 fails to remain distinct, and by 2230 UTC unobserved banding features are seen in the OK_Control run (Fig. 8k). The storms are also larger than those observed, with large areal coverage of Z 45 dBZ. In contrast, the OK_ZDRCOL run is much closer to the observations. While a small northeastward bias does still exist, the forecast storms are in better agreement with the observations in terms of size and position, with three distinct storms featuring identifiable hook echoes and broad supercellular features present 1 h into the forecast (Fig. 8l).
Based on these encouraging qualitative results, the equitable threat scores (ETS) and frequency biases were computed for a quantitative look at the performance of OK_Control and OK_ZDRCOL. The ETS, also known as the Gilbert skill score (Gilbert 1884), is given by
where H is the number of hits, M is the number of misses, F is the number of false alarms, and is the number of hits expected due to random chance, given by
where N is the total number of forecast points included in the calculation. The ETS is calculated on a gridpoint basis satisfying or exceeding a defined (here Z) threshold, with a value of 1.0 indicating a perfect forecast and 0.0 indicating no forecast skill. The bias is calculated from
and provides a ratio of the number of forecast grids and the number of observed grids exceeding a threshold, normalized to zero. A bias of zero indicates no bias, while a positive bias indicates an overestimate of Z exceeding a threshold. Both the ETS and bias were calculated for the composite Z at 20-, 30-, and 40-dBZ thresholds and are shown for the OK case in Fig. 9. Note that the ETS at the analysis time is not necessarily equal to 1.0 due to both smoothing procedures and differences in how Z is calculated: the observations are being compared against simulated Z derived from a T-matrix code for the single-species hydrometeor distribution that was retrieved from the cloud analysis. The ETS for the 20-dBZ threshold (Fig. 9a) are comparable between the two experiments, but notable improvements in ETS are seen in OK_ZDRCOL over OK_Control for high Z thresholds (Figs. 9c,e). The ETS for OK_ZDRCOL remains superior for the entire 1-h duration of every forecast, showing a noteworthy positive impact of column assimilation. Both OK_Control and OK_ZDRCOL exhibit generally positive biases that increase with time at all three Z thresholds. For all forecasts at all times, however, OK_ZDRCOL features lower biases (Figs. 9b,d,f). This tendency toward lower Z biases in OK_ZDRCOL is also seen in 1-km Z for the 2130 UTC forecast (Fig. 8).
b. 25 May 2016 case
The KS case presents a somewhat more challenging forecast scenario owing to a complex evolution of the supercell and the greater distance between the supercell and the radar. After becoming mature, the main supercell began moving slowly to the southeast. A new storm developed to the southwest of the main supercell, which produced a left-moving supercell that moved off to the north-northeast before merging with the primary supercell. Additional convection also formed along, and was absorbed into, the southern flank of the forward-flank downdraft in the supercell. This storm was farther away from the radar than the storms in the OK case were (initiation occurred approximately 140 km away from the radar compared to 65 km away from the radar in the OK case), resulting in a decrease of the quality of radar data available for assimilation due to both decreased low-level coverage and increasing radar resolution volume (0.22 km3 at a 65-km range vs 1.01 km3 at a 140-km range for a 0.5° elevation angle). Additionally, in contrast to the OK case, this case lacked the assimilation of mesonet surface observations.
A long, continuous swath of columns are shown for the duration of the period analyzed in the southwest corner of the supercell where the main updraft is expected to be located (Fig. 10). The column’s width and depth increases shortly before the start of the long-track tornado northeast of Salina, and the column remains broad and deep until the end of the tornado, near the end of the assimilation period. A second column swath is seen with the left split of the supercell as it moves off to the north-northeast. This column weakens near the end of the assimilation period, and the storm weakened shortly thereafter. As in the OK case, the easternmost tornado falls outside the analyzed period for this case.
The composite plot of maximum w in the analyses for the KS case shows many of the same improvements documented in the OK case. The KS_Control case shows a more disorganized and less coherent updraft path, with many spurious updrafts to the north of the main supercell path and observed tornado tracks (Fig. 11a). Considering that the end of the assimilation period is near the ending time of the long-track tornado, the general progression of the analyzed updrafts is also too fast. In contrast, KS_ZDRCOL features a much more coherent updraft swath with a slower forward motion to the east-southeast and a path closer to the observed tornado track (Fig. 11b). KS_ZDRCOL also features less spurious convection than KS_Control in the central and southern parts of the domain.
A comparison of 1–6-km updraft helicity with MRMS-derived rotation tracks, similar to Fig. 7, is shown for the KS case in Fig. 12. The forecasts selected here were chosen to coincide with the long-track tornado. In the forecast initiated at 2300 UTC (Figs. 12a,b), both KS_Control and KS_ZDRCOL produce a developing supercell north of Salina with an unorganized updraft and an east-northeast motion. The observed rotation tracks show only slight, messy rotation during this period. Starker differences are apparent for the 0000 UTC forecast (Figs. 12c,d). KS_Control features a disorganized updraft helicity swath displaced far to the north of the observed rotation track. In contrast, KS_ZDRCOL features a consolidated updraft helicity swath through the duration of the forecast period along and just north of the observed rotation track, although a slight slow bias in forward speed is apparent. These same general patterns are also observed for the 0100 UTC forecast, with a noisy updraft helicity field too far to the northeast in KS_Control; KS_ZDRCOL exhibits a large, southeastward-directed updraft helicity swath displaced slightly southwest of the observed rotation track.
An example of observed and forecast Z for both KS_Control and KS_ZDRCOL is shown in Fig. 13 for the 0000 UTC forecast. This 1-h period begins near the start time of the primary long-track tornado and features a complex evolution involving the secondary storm to the southwest splitting and merging with the main supercell (Figs. 13a,d,g,j). As such, both KS_Control and KS_ZDRCOL struggle to accurately predict the evolution of the storm during this period. A very large and elongated forward-flank downdraft not seen in the radar observations quickly develops and extends to the east-southeast and east-northeast in KS_Control and KS_ZDRCOL, respectively. This forward-flank precipitation appears to stem from weak upper-level Z in the anvil in the observations. Despite this, KS_ZDRCOL features a more realistic supercell structure 20-min into the forecast (Fig. 13f) compared to KS_Control (Fig. 13e), with a well-defined hook echo and rear-flank downdraft near the observed tornado track. Neither KS_Control nor KS_ZDRCOL clearly capture the left-splitting supercell. An erroneous region of moderate Z (i.e., 25–35 dBZ) within the inflow region of the supercell is also seen in the KS_Control run (Fig. 13e) that is not seen in the KS_ZDRCOL run. Both KS_Control and KS_ZDRCOL generally feature Z values that are too low (by 5–10 dBZ) compared to observations outside of the forward-flank downdraft. Overall, KS_ZDRCOL features a slower and noticeably more accurate forecast track of the hook echo than KS_Control (as also seen Figs. 12c,d), as well as a more realistic looking hook echo (Figs. 13i,l vs Figs. 13h,k).
Quantitatively, KS_ZDRCOL generally exhibits improvements over KS_Control with higher ETS scores and lower biases, although the improvements in ETS scores are more mixed than in the OK case, with lower scores for the first two forecasts in the period (Fig. 14). Overall scores are lower in the KS case, in part due to the challenging nature of the forecast and in part due to the aforementioned biases in Z (e.g., Figs. 14b,d,f) and the extensive forward-flank downdrafts, which generally exceed the biases seen for the OK case, particularly for later forecasts.
5. Summary and discussion
In this study, the potential for the assimilation of polarimetric radar data observations via a cloud analysis technique to aid in the spinup and forecast of convection in storm-scale models is examined. Differential reflectivity columns are ubiquitous features of deep moist convection that are coincident with updrafts and, thus, with areas of saturation and latent heat release. Based on this premise, a column detection algorithm is developed to identify columns and 1) insert positive temperature and moisture perturbations at these locations and 2) remove modest amounts of moisture outside of these locations where Z exceeds 10 dBZ. To evaluate this method, two cases are analyzed: the 19 May 2013 tornadic supercells in central Oklahoma and the 25 May 2016 tornadic supercell in north-central Kansas. For each case, two runs were performed to gauge the impact of these changes: a “Control” run using the original cloud analysis, and a “ZDRCOL” run using the newly modified cloud analysis that incorporates dual-polarization radar data.
The column detection algorithm is shown to reliably identify columns associated with convective updrafts, and notable qualitative and quantitative improvements in the analysis and forecasts of convection are seen in both cases. The analyzed updraft tracks are more coherent and consolidated with less spurious convection in the ZDRCOL runs than in the Control runs. The short-term forecasts show a reduction in forward speed and northward position bias of the main updrafts, an issue encountered in many storm-scale modeling experiments, with updraft helicity indicating better agreement with radar-derived rotation tracks in each case. The forecast Z fields also agree better with observations in the ZDRCOL runs compared to the Control runs, in turn leading to improved quantitative verification scores and reduced frequency bias.
These experiments represent a basic proof-of-concept investigation of the potential for assimilating columns into storm-scale models and warrant further study. However, drawbacks to the method examined here exist that have yet to be addressed. First, while columns are fairly ubiquitous in deep convection and generally collocated with updrafts, they can be masked by the presence of hail or tornadic debris (as discussed in Snyder et al. 2015) in the updraft, resulting in the intermittent appearance (or complete disappearance) of columns; in some cases, columns may not be observed in deep convection at all. Future work may examine the potential of using columns in a similar manner to alleviate these issues (van Lier-Walqui et al. 2016), although the use of columns bears its own shortcomings (e.g., poor estimation of in areas of limited precipitation). The use of columns to aid spinup of precipitation in a NWP model may not be appropriate for weak convective storms and stratiform rain where columns are ill defined or may not exist. Only two cases were analyzed in this study, with each being an archetypal case of very strong convection with good radar coverage and prominent columns. The parameters both for detecting columns and for applying moisture and temperature increments were subjectively determined and should undergo further refinement. Finally, only a single radar was assimilated for each case because of lingering uncertainty in the optimal methodology for merging polarimetric data from multiple radars, although work is ongoing (J. Krause 2017, personal communication) to implement the polarimetric compositing method developed in Homeyer (2014) and Homeyer and Kumjian (2015). The polarimetric version of the newly developed Storm Labeling in Three Dimensions (SL3D) algorithm (Starzec et al. 2017), which has demonstrated success in identifying convective updrafts using a combination of weak-echo regions, columns, and columns, may also prove particularly useful going forward. It is likely that assimilating data from multiple radars would assist in the spinup and analysis of storms from radial velocity data beyond that examined in this work, as well as alleviate radar coverage concerns for detecting columns.
While being a relatively simple and efficient method for assimilating Z, cloud analysis techniques may not be optimal owing to their inherent empirical relationships that can compromise initial adjustments in the model. As temperature and moisture increments appear to play large roles in aiding the spinup of observed storms in storm-scale models, future work will seek to explore the possibility of assimilating cloud analysis-derived temperature and moisture increments based on detected columns as “pseudo observations” in a 3DVAR framework, similar to the work of Fierro et al. (2016) for lightning data assimilation. Using a variational framework to assimilate columns would allow for a more balanced analysis/solution between the kinematic and thermodynamic fields and, hence, for a smoother cycling process.
The authors thank the three anonymous reviewers for their constructive feedback, as well as Keith Brewster and Elizabeth Smith for their helpful discussions and Alexandre Fierro for providing useful feedback about the manuscript. Funding was provided by NOAA/Office of Oceanic and Atmospheric Research under NOAA–University of Oklahoma Cooperative Agreement NA11OAR4320072, U.S. Department of Commerce. Additional funding was provided by the U.S. Department of Energy Atmospheric System Research Grant DE-SC0014295 and by NSF Grant AGS-1341878. The computing for this project was performed at the OU Supercomputing Center for Education and Research (OSCER) at the University of Oklahoma (OU).