A three-dimensional variational (3DVAR) assimilation technique developed for a convective-scale NWP model—advanced regional prediction system (ARPS)—is used to analyze the 8 May 2003, Moore/Midwest City, Oklahoma tornadic supercell thunderstorm. Previous studies on this case used only one or two radars that are very close to this storm. However, three other radars observed the upper-level part of the storm. Because these three radars are located far away from the targeted storm, they were overlooked by previous studies. High-frequency intermittent 3DVAR analyses are performed using the data from five radars that together provide a more complete picture of this storm. The analyses capture a well-defined mesocyclone in the midlevels and the wind circulation associated with a hook-shaped echo. The analyses produced through this technique are used as initial conditions for a 40-minute storm-scale forecast. The impact of multiple radars on a short-term NWP forecast is most evident when compared to forecasts using data from only one and two radars. The use of all radars provides the best forecast in which a strong low-level mesocyclone develops and tracks in close proximity to the actual tornado damage path.

1. Introduction

Doppler radar observations became more widely used as an analysis tool since the 1970s, by mapping polar coordinate radial velocity data to a Cartesian grid space. Results in Brandes [1] showed the evolution of a supercell including an intensifying mesocyclone and development of the rear-flank downdraft (RFD), as observed by a dual-Doppler analysis. Additional study was done using this technique in Brandes [2] where the airflow within a tornadic thunderstorm was compared with the observed tornado damage path. Brandes [2] also found that, at the beginning of the tornadic phase of a supercell thunderstorm, an increase in radial flow near the tornado was observed in the radar data analysis. This helped determine that there was a strongly convergent flow beneath the mesocyclone during this time. Another major development in the use of Doppler radar was the discovery of a unique Tornado Vortex Signature (TVS), that coincided with the development of low-level rotation and tornadogenesis in the 1973 Union City tornado [3]. The development of this signature was closely linked to the mesocyclone structure model found in Lemon and Doswell [4].

With the frequent use of dual-Doppler analysis and increased resolution of storm scale models, data assimilation has become an important area of research that integrates observations into a numerical weather prediction (NWP) model. The development of variational data assimilation techniques for Doppler analysis has shown better results over other methods of analysis [5]. Several methods to obtain an initial state by thermodynamic retrieval using observations were originally proposed by Gal-Chen [6] and tested by Hane et al. [7] for simulated data and Gal-Chen and Kropfli [8] for radar observations of the planetary boundary layer.

Crook [9] tested a variational technique that would allow the retrieval of thermodynamic variables using radar data and a background sounding. The technique used a time tendency term to estimate the thermodynamic variables and then estimated the deviation of those variables from the background values. This technique was then tested using observational data from three gust front cases [10]. Several data interpolation methods were tested to create an initial analysis, followed by a forecast to see if the analysis could help accurately predict the propagation of several outflow boundaries. The findings from these experiments showed that the analysis and forecast were greatly improved when surface data were used in addition to radar data and upper air soundings.

Weygandt et al. [11] tested a single velocity Doppler retrieval technique on a supercell thunderstorm with the use of the Taylor frozen turbulence hypothesis as a weak constraint. The results showed that there was reasonable agreement in the analysis with the exception of the magnitude of the vertical velocities. In an additional experiment performed by Weygandt et al. [12], thermodynamic retrieval was included using three separate retrieval times for a supercell thunderstorm. The initial analysis from this study showed a supercell which exhibited pressure perturbations consistent with Rotunno and Klemp [13] linear theory. When a non-hydrostatic model was initialized using a set of retrieved fields as initial conditions, the evolution of the storm seemed reasonably well predicted [12]. Results, obtained from another experiment using a simplified version of the analysis technique which only included radial velocities, estimated mean horizontal wind, and perturbation radial divergence, were poor. One problem noted by this study was the speed of the cold pool propagation, which was faster than the observed speed. This could have been due to the lack of ice microphysics in the model [12].

A four-dimensional variational (4DVAR) analysis technique was developed by Sun and Crook [14, 15], whose goal was to accurately retrieve wind and thermodynamic fields by minimizing the difference between the background and the observations defined by a cost function. By applying this technique, all the model variables are determined at the same time. Initially, the method was tested with simulated observations which showed good results for the thermodynamic structure, and the addition of a penalty term further increased the accuracy of the model state [14]. The results from the latter studies using 4DVAR analysis of several Florida thunderstorms showed good agreement with observational data [15]. Sun and Crook [16] also used the 4DVAR technique to assimilate data from a line of strong thunderstorms. Their findings showed that the penalty term helps to decrease noise in the analysis and the background term increases the forecast accuracy (after the first 10 minutes) in areas of sparse observational data. It was also found that using a previous analysis as a background for the analysis produced better results than those obtained by using a short-term forecast as a background [16]. However, the 4DVAR analysis technique is computationally expensive.

Gao et al. [17] tested a less costly computational technique known as three-dimensional variational (3DVAR). In this method, the cost function is composed of several terms: an observational term, a background term, and several penalty terms. Simulated observations showed that this analysis technique was less sensitive to the boundary conditions than other methods that used mass continuity as a strong constraint, thereby mitigating the effects of error accumulation through explicit integration. The use of a recursive filter has been tested with the 3DVAR technique to improve the quality and efficiency of the analysis [18]. Each pass of the recursive filter employs the use of a filter in two directions, a left moving filter and a right moving filter [19]. When applied using dual-Doppler data of a supercell thunderstorm, it was shown that the use of a recursive filter produced analysis results similar to that in Gao et al. [17] but greatly improved the efficiency of the 3DVAR method.

The method developed in Gao et al. [18] was also used to run an analysis and numerical prediction cycle of a tornadic thunderstorm outbreak in north Texas [20, 21]. Two cloud analysis schemes along with the 3DVAR technique were tested to determine the sensitivity of these schemes to the analysis and subsequent forecast. The use of radial velocity in the analysis was found to play an important role in the development of vorticity in the forecasts. Also, the role of reflectivity in the analysis was shown to have a major impact on the development and maintenance of intense convection.

Recently the ensemble Kalman filter (EnKF) technique, originally developed by Evensen [22], has been widely used as a new method to assimilate observations into the model state. The EnKF technique seeks to accurately determine the flow dependent background error covariance through the use of an ensemble of nonlinear forecasts. The application of the EnKF method to storm-scale data assimilation has been examined by Snyder and Zhang [23], Zhang et al. [24], Tong and Xue [25], and Xue et al. [26] using simulated data, and by Dowell et al. [27], Snook et al. [28, 29], and other authors using observed radar data. Although EnKF avoids the linearization of the background error covariance, the use of a large ensemble is computationally expensive. Also, when decreasing the amount of members, the covariance can become underestimated.

Most of the above research used observations of severe thunderstorms from only one or two radars located close to the convection being studied, thereby limiting the information available to provide a complete picture of the storm, and which does not take full advantage of the WSR-88D radar network. By only using data from one or two nearby radars in an analysis, a large part of atmosphere located above the radar cannot be effectively observed due to the limits in the elevation of a radar volume scan, known as the “cone of silence”.

The 8 May 2003 tornado in central Oklahoma has been widely studied because of the high impact on society and the high density of weather observations. The storm spawned a F4 tornado that tracked through southern Oklahoma City causing $370 million in damage and 134 injuries. Located in the region are four WSR-88D radars: KVNX, located in northwest Oklahoma at Vance Air Force Base; KINX, located just northeast of Tulsa; KFDR, located in southwest Oklahoma near the town of Frederick; and KTLX, located southeast of Oklahoma City (Figure 1). Two additional radars are also located in close proximity to Oklahoma City: KOUN, a research dual-polarization WSR-88D radar located in Norman, and KOKC, a terminal Doppler weather radar (TDWR) used by the FAA near Will Rogers Airport in Oklahoma City. Because these two radars are too close to each other, we chose to use data from one of them, the KOKC TDWR radar, in this study. The Oklahoma Mesonet also provides surface observations every five minutes for 119 sites in all 77 of Oklahoma counties [30]. This, along with four upper-air wind profilers in the state, allows for high temporal resolution of atmospheric phenomena.

There have been several studies examining the evolution of the 8 May 2003 Oklahoma City supercell through use of radar data, assimilation techniques, and numerical models [27, 3136]. In these studies, observations from only one or two radars close to the storm were assimilated, thus omitting the upper-level structures of this supercell.

In this study, multiple radars both close to and far from the storm are used to perform an analysis using the 3DVAR technique. The use of the 3DVAR technique to study this case was chosen over the EnKF or 4DVAR techniques because the 3DVAR technique is not computationally costly and has shown good results in other related studies [18, 20, 21]. In Section 2, the advanced regional prediction system (ARPS) and its 3DVAR technique are briefly discussed. Section 2 also describes the 8 May tornadic storm case, the processing of radar data, and experiment configurations. Results about the analysis and forecast experiments are presented in Section 3. Finally, conclusions are discussed in Section 4.

2. Models, Methods, and Data

2.1. ARPS Model, ARPS 3DVAR, and Cloud Analysis Algorithm

The advanced regional prediction system (ARPS) has been developed at the Center for Analysis and Prediction of Storms (CAPS) at the University of Oklahoma over the past 15 years [3739]. The ARPS model was designed as a system suitable for explicit prediction of convective storms. It is a three-dimensional, nonhydrostatic, compressible model formulated in generalized terrain-following coordinates. The model employs advanced numerical techniques, including monotonic advection schemes for scalar transport and variance-conserving fourth-order advection for other variables. The model also includes state-of-the-art physics parameterization schemes that are important for explicit prediction of convective storms. The system has been used in real-time high-resolution prediction experiments for convective scales in the past several years over the continental United States [40, 41].

A detailed description of the ARPS 3DVAR system can be found in Gao et al. [18]. In the current version of the ARPS 3DVAR system, the spatial covariances of the background error are modeled by a recursive filter [19], and the square root of the matrix is used for preconditioning [42]. The corresponding covariance matrix is diagonal, and its diagonal elements are specified according to the estimated observation errors.

Radar data can easily be incorporated into the cost function in the observation term; the observed values represent the radial component of wind even though the data have been mapped to a Cartesian grid during preprocessing. The ARPS 3DVAR data assimilation system can ingest data from a number of different sources including surface air observations, upper air soundings, profiler data, and aircraft data. Quality control in ARPS 3DVAR is performed by determining if an innovation vector is less than a threshold value. An observation whose innovation exceeds the threshold is rejected. The quality control threshold is a function of the specified background error, in which a larger background error value corresponds to a larger threshold. With the 3DVAR analysis, the background error also affects the relative weight of the observations and background. Observations can also be rejected based on climatological error statistics of the observation stations.

One unique feature of the ARPS 3DVAR system is that multiple analysis passes can be used to analyze different data types with different spatial correlation scales to account for the variations in the observation density among different data sources. Upper-air rawinsonde data and radar observations are examples of two observing systems with very different spatial density. An initial analysis pass can be performed using only upper air observations, using a large correlation scale. A second pass is then performed with radar data with a smaller correlation scale. Such a procedure allows for the retention of multi-scale information contained in observations of vastly different spatial density. A mass divergence constraint is also employed within the 3DVAR cost function to help couple the three wind components. This constraint permits the retrieval of vertical velocity from the mostly horizontal wind components as measured by radar and also helps spread observational information to nearby grid points in the analysis [21].

After completion of the 3DVAR analysis, a cloud analysis is performed. The procedure is based on that from the local analysis and prediction system [43] and includes modifications from Zhang et al. [44] and Brewster [45]. The cloud analysis package estimates mixing ratios of water vapor, rainwater, and cloud water based on reflectivity measurements. The cloud analysis can also adjust in-cloud temperature based on reflectivity [20].

2.2. The 8 May Storm Environment

In the first ten days of May 2003, there were 434 tornadoes, more than any month in the US on record (although the record was broken in 2008). On three of the ten days, central and eastern Oklahoma were hit by significant tornado outbreaks. The 1200 UTC upper-air observations on the morning of 8 May 2003 showed a significant positively tilted trough in the Pacific Northwest characterized by 500 mb heights near 540 decameters (not shown). This trough exhibited a double-barrel structure with two regions of minimum heights, one in central Alberta and a second approaching the coast of Washington State. Several short-waves rotated around the large trough and propagated downstream. The most prominent of these short waves (shaded in red on Figure 2(a)) provided enhanced vorticity advection over Oklahoma in the hours previous to 0000 UTC 9 May which induced midtropospheric rising motion and pressure falls at the surface. These pressure falls enhanced those associated with lee troughing downstream of the Rocky Mountains, over which a strong upper-level jet existed. A poleward meridional flow then developed on the eastern side of the developing low that advected large quantities of unstable warm moist air from the Gulf of Mexico in the lower levels of the atmosphere (Figure 2(b)). An important aspect of this dynamical setting was the continued presence of the trough in the northwestern United States during the week prior to the 8 May event and the subsequent reinforcement of the lee trough and southerly winds. This continuous flow of moist air from the Gulf of Mexico helped reinforce thermodynamic instability over the central plains throughout the week.

At the surface, several outflow boundaries existed in the warm sector caused by a complex of nocturnal tornadic thunderstorms. These outflow boundaries served as a focal point for convergence and subsequent convective initiation later that afternoon. As the convection from the previous night moved east along the Red River and dissipated, skies over Oklahoma cleared, contributing to strong surface heating during the day. The low-level environment through the late morning of 8 May was characterized by temperatures near 28°C and dewpoints around 22°C. A sounding from Norman, Oklahoma, (OUN) at 0000 UTC 9 May 2003 showed an extreme level of instability within the warm sector (Figure 3). Surface convective available potential energy (CAPE) reached values near 5000 J kg−1, but convection was inhibited initially by a capping inversion. The OUN sounding also shows a well-defined “quarter circle” hodograph representative of strong directional and speed shear with associated storm relative environmental helicity (SREH) values of 317 m2 s−2 in the lowest 1 km. This strongly sheared environment with large values of CAPE was conducive to the formation of strong tornadoes.

At 3 pm local daylight time (LDT) or 2000 UTC, satellite data showed the upper-level short-wave moving eastward into western Oklahoma. Observations from the Oklahoma Mesonet indicated a dryline located just to the west of OUN, in central Oklahoma. Along this boundary dry air from the western p pmlains which converged with moist air from the Gulf (not shown). By 4 pm LDT, towering cumulus had penetrated the capping inversion and initiated strong surface-based convective thunderstorms along a bulge in the dryline. After 4 pm the storms that formed along the dryline became severe in Grady County, just southwest of Oklahoma City and by 4:30 pm LDT the most prominent of the convective cells split and deviated in storm motion by about 9 degrees to the right of the mean motion. At this time both cells showed moderate rotation at midlevels. By 5 pm LDT, the right-split storm was entering Oklahoma County and radar detected strong low-level convergence and rotation. It was at this time that sporadic wind damage was caused by the advancing RFD. Minutes later, at 5:11 pm LDT (2211 UTC) a tornado was confirmed on the south side of Moore, Oklahoma. The tornado increased in intensity as it tracked to the east-northeast just south of Interstate 240 and into Oklahoma City. At 5:38 pm LDT, after traveling 30 km, the tornado dissipated 4 km south-southeast of Choctaw, Oklahoma.

2.3. Radar Data Ingest and Preprocessing

Radar data must undergo several procedures before being used in the analysis. For this study radar data in several formats were ingested into the ARPS 3DVAR system. Level II data were obtained from the National Climate Data Center (NCDC) in binary format which included full volume scans every 5 to 6 minutes. Data from TDWR were ingested in network common data form (NetCDF) format with full volume scans being reconstructed using each tilt during the preprocessing procedure. TDWR data included 9 elevations, while WSR-88D radars had 13 elevation levels. After the volume scans are read into the program, the program checks the quality of the velocity and reflectivity data and makes corrections to the data if needed. For raw velocity data, the absolute values of velocities greater than the Nyquist velocity of the radar are folded (Figure 4(a)). This means that inbound velocities greater than the Nyquist velocity in the raw data will be displayed as if they were outbound velocities (of the opposite sign) with a magnitude less than the Nyquist velocity. To correct the velocity information, a dealiasing algorithm is used to unfold the raw velocity data by comparing the data to a background sounding or by comparing the data to velocity measurements surrounding the region being dealiased. If the velocity is determined to be folded, the deviation from the Nyquist velocity will be found and then added to the Nyquist velocity of the opposite sign. Figure 4(b) shows a sweep of velocity data after dealiasing which clearly depicts the location of the circulation. The initial radar reflectivity includes ground clutter (displayed in Figure 4(c)), anomalous propagation, and transient clear air echoes. These artifacts are removed by additional radar reflectivity quality control procedures by detecting both large changes in mean radial velocity and reflectivity over neighboring range gates at the lowest elevations.

After being cleaned, the data are projected onto a Cartesian grid with the same resolution as the analysis being performed. Radar data are then interpolated to a grid using a least square method. During the interpolation, a linear fit is used in the vertical and a quadratic polynomial fit is used in the horizontal. The domain of this method is bounded by the range of available radar data, precluding extrapolation of the data.

2.4. Experiment Design

A number of analyses are performed on the 8 May 2003 Moore/Oklahoma City supercell at 5-minute intervals starting from 2145 UTC to 2240 UTC. For all experiments, the grid resolution in the horizontal direction is 1 km, with a lower bound vertical grid resolution of 100 meters which is stretched vertically to 23 km using a cubic function with an average grid spacing of 500 meters. The domain size is grid points with the center of the domain located just east-southeast of Oklahoma City at the location where spotters reported a large tornado. The background state is obtained from a 9 km ARPS forecast using a background from the operational NCEP ETA model and upper air observations. A forecast is then performed using intermittent analyses that assimilate only sounding and profiler data starting from 18 UTC to produce the background for the 1 km analyses. The observational data that are used in the 1 km analysis consist of Oklahoma Climate Survey Mesonet observations, profiler data, and radar data. Upper air sounding data were only included in the 9 km forecast and not the 1 km analyses due to the sparse temporal resolution of the data. The radar data used in this study consist of four WSR-88D radars and one TDWR, as discussed in the introduction. Since radar observations are obtained at different times, scans closest to the time of the analysis were used.

In the control analysis experiment, all available observations are ingested into the assimilation system using four analysis passes. During the first pass the radius of influence is set at 50 km in the horizontal, 3 grid points in the vertical, and only profiler data are used. On the second analysis pass, the radius of influence is changed to 25 km in the horizontal and 2 grid points in the vertical, and Mesonet and profiler data are used. In the third pass, the radius of influence is decreased to 10 km in the horizontal and 2 grid points in the vertical, and Mesonet and radar data are used. On the fourth and final pass, the radius of influence is changed to 5 km in the horizontal and 1 grid point in the vertical, and only radar data are used.

Two other analysis experiments are performed to determine the impact of the number of radars used in the analysis at one time, 2200 UTC (Table 1). These experiments show the role of each radar on the control analysis and their impact on the divergence and vertical velocity fields. The second experiment (1RAD) is the same as the control analysis, but with the use of radar data only from KOKC, the closest radar to the storm during a large duration of the analysis period. The third experiment (2RAD) includes the two radars that only observed the low levels of the thunderstorm, KOKC and KTLX.

Finally, several forecast experiments are performed to examine the sensitivity of a forecast to the number of radars using the previous 3DVAR analyses at 2155 UTC (Table 1). These analyses are then used as initial conditions for integration of the ARPS model. To demonstrate that use of data from radars that are close to the storm may not always benefit the forecast, we perform two more forecast experiments; one with only the KOKC radar, and a second with both the KOKC and KTLX radars.

The model was run using the Lin 3-ice microphysics scheme [46]. A large time step of 2 seconds and a small time step of 5 seconds are used to perform forward integration of slow and acoustic wave modes, respectively. The drop size distribution’s intercept parameter for rain was adjusted from the default value of to to allow for expected larger raindrops associated with the supercell storm that will result in less evaporative cooling and a weaker cold pool [47].

3. Results and Discussion

3.1. The Control Analysis with all Radars

In the control analysis, we use data from all available radars except for KOUN (which is very close to KTLX and therefore has similar viewing angles), so the amount of data used in the analysis is quite large. This analysis best resolves the storm-scale features throughout the analysis period and therefore serves as a benchmark for other experiments. The evolution of the storm from the stage before tornado production until the tornado lifted is shown in Figure 5. The development of a hook feature around 2200 UTC about 10 minutes before the first touchdown and dissipation just after 2230 UTC about 10 minutes before the tornado lifts during the final stages of mesocyclone occlusion are very clear. Halfway through the life of the tornado, the storm weakens and the reflectivity decreases quite rapidly after 2230 UTC. This is due to the full occlusion of the RFD around the tornadic circulation [4] and decreased updraft intensity associated with a downward oriented pressure gradient [48, 49]. Wind vectors shown in Figure 5 depict intense inflow during the mature period of the storm and weak inflow during the decaying period.

The evolution of the low-level mesocyclone and occluding RFD is visible in the divergence field (Figure 6). The storm evolves from an elongated diffuse area of weak convergence on the south flank of the reflectivity contour at 2150 UTC (shown in Figure 6(a)) to a more concentrated area of convergence beneath the area of the mesocyclone by 2220 UTC (Figure 6(d)). Just before this at 2200 UTC (Figure 6(b)) and then at 2210 UTC (Figure 6(c)) divergence is strongest wrapping around (to the west and south of) the area of convergence associated with the RFD. The advancement of the RFD around the southwestern region of the updraft at 2200 UTC coincides with the development of the tornado. The mesocyclone tracks across an area of enhanced low-level convergence on the front flank of the storm similar to that in the Lemon and Doswell model [4]. Strong upper-level divergence is visible at the 12 km level in the analysis (Figure 7), which can be attributed to the assimilation of data from radars that are located far away from the storm. As the storm progressed to the east-northeast, the divergent nature of the outflow weakened.

The dipole vorticity structure simulated by Klemp and Wilhelmson [49] and discussed by Rotunno [50] is clearly depicted in this analysis (Figure 8). The right split storm exhibits strong positive vertical vorticity at the beginning of the period with a weaker anticyclonic counterpart on the northwest side of the mesocyclone. The left split storm also exhibits signs of the dipole vorticity structure opposite to that of the right split but with weaker magnitudes of vorticity and vertical velocities. Since the hodograph in Figure 3 favors a right split storm, the left split counterpart never evolves into a mature storm [13]. It is shown that the vertical vorticity is strongest at 2210 UTC when the tornado touched down (Figure 8(c)).

3.2. Sensitivity of the Analysis to Different Number of Radars

To investigate the contribution of individual radars to the control analysis, two more 3DVAR analyses are performed only at 2200 UTC. One uses data from KOKC only (1RAD), and another uses data from KOKC and KTLX (2RAD). When assimilating radar data over a period of time, the effects of unavailable data above the highest elevation level of the volume scan, referred to as the cone of silence, may be minimized by using a storm scale NWP model to estimate quantities in regions that are unobservable [3436]. For this experiment, however, the analysis does not make use of data assimilation cycles from a storm-scale NWP model to fill in the gaps between observations. Therefore, in regions of no radar coverage, the analysis relies heavily on the background, the spatial spreading of observation information through background error correlations, and the mass continuity constraint to fill in the gaps (conventional upper air data are usually very limited). The impact of decreased radar coverage is most clearly depicted in Figure 9.

When radar data from KOKC only (1RAD), located closest to the storm, are used, the ability to obtain cross-beam winds is very limited. Even at the 3 km level (Figure 9(e)), a significant part of the circulation in the storm is missed. At the top level (12 km), both the reflectivity and storm outflow are almost completely missed (Figure 9(h)). When data from both KOKC and KTLX (2RAD), which intersect the hook of the storm at about 45 degrees, are used, the accuracy of the analyzed wind is increased dramatically. This is most visible along the front flank of the thunderstorm at the lowest levels just east of the mesocyclone where winds are parallel to the cross-beam component of the KOKC radar (Figures 9(a), 9(b), and 9(c)). The most notable impact of the radars is in the reflectivity structure. Both experiments with one and two radars that are close to the storm miss significant parts of the reflectivity structure owing to the cone of silence (Figures 9(h), and 9(i)). This artifact is eliminated when data from three other radars (KVNX, KINX, and KFDR) far away from the storm are used, as depicted in the control analysis (Figures 9(a), 9(d), and 9(g)). The ability of the analysis scheme to resolve the structure of the updrafts and downdrafts within the storm is also improved when using all available radars. This is most evident when looking at a cross section as shown in Figure 10. The higher vertical velocities in the updraft are resolved only when at least two radars are used. In the control analysis, the additional radar data from radars located far from the storm help resolve the anvil structure above the 10 km level (Figure 10(a)) up to 20 km.

When assimilating data from only KOKC, low-level convergence is weak, and upper-level divergence is almost missing (Figures 11(b) and 11(e)). When two close radars KTLX and KOKC are used, the low-level convergence is more pronounced when compared to using a single radar (Figure 11(b) versus Figure 11(c)), but the upper-level divergence is still missing. Only when assimilating data from radars located further away from the storm are used, do areas of upper-level divergence, which are consistent with the Lemon and Doswell [4] model, appear in the analysis (Figure 11(d)).

When observations of the storm structure are incomplete (as with the analyses performed by using data from only one radar), the variational analysis technique acts to fill in the gaps with values obtained from other sources (surface and upper air), as well as information from the background state. While providing a reasonable estimate to the state of the atmosphere in synoptic and mesoscale data assimilation, sounding and profiler data, when available, are rather crude estimates of the true state of the atmosphere when deep convection is present.

In summary, the control analysis captures the structure and evolution of the supercell using information from multiple radars. Features at both the low and upper levels are well resolved. These storm scale features are most apparent when examining convergence and vorticity within the analysis that depict a developing and strengthening low-level mesocyclone. Thus, a complete set of analyses of the storm is obtained when using data from radars located both close and far away from the storm, a strategy often overlooked in previous investigations.

3.3. Sensitivity of the Forecast to Different Number of Radars

Three forecast experiments are performed using initial conditions similar to those produced with data from one, two, and five radars to examine the sensitivity of a forecast to the number of radars used. To limit the amount of noise from radar data, only radar reflectivity and radial velocity in regions of high reflectivity (greater than 50 dBZ) and close to the convection were assimilated. All three forecasts are initialized using the 3DVAR analysis at 2155 UTC.

The forecast 0.5 km vertical vorticity and horizontal wind vectors between 2200 UTC and 2220 UTC are shown in Figure 12 at ten-minute intervals with the observed tornado damage path. The KOKC forecast, assimilating data from only one radar, shows the storm moving in a north-northeasterly direction, away from the actual tornado damage path. In addition, the simulation produces only low levels of positive and negative vorticity. The 2RAD forecast (Figures 12(d)12(f)) used data from only the KTLX and KOKC radars, and although it captured more of the low-level wind structure than the KOKC experiment, their volume scans still did not observe the top portion of the storm. This forecast predicts storm movement closer to the observed damage path than the KOKC experiment, but still to the north of the observed tornado damage path. In general, both these forecasts without the observations in the upper levels of the thunderstorm deviated to the north of the tornado damage path. The control forecast (using 5 radars) shows a strong maximum in vorticity (Figures 12(g)12(i)) close to the actual storm track. The storm propagates in an east-northeasterly direction closer to the observed tornado damage path than the previous two experiments. The inclusion of radar data in the upper levels of the convection seems to greatly improve the accuracy of prediction of the simulated storms mesocyclone in relation to the actual tornado damage path.

To quantify the improvement of the predicted storm for the three experiments the equitable threat scores (ETS; [51]) of the predicted reflectivity were calculated using KTLX radar composite reflectivity. The ETS scores using a 15 dBZ threshold are shown in Figure 13(a). There is a pronounced improvement in the ETS scores from 2210 through 2235 UTC for the control experiment when compared to the experiments which only assimilated data from one or two radars. Between 2205 and 2215 UTC the control run does produce higher scores for the 45 dBZ threshold than the other two experiments but the improvement is degraded after 2215 UTC.

These results are similar to that of Xue et al. [36] for a 1 km forecast of the same thunderstorm but with the use of a cycling 3DVAR analysis. In their study, only one radar was used, KTLX, and their storm also propagated to the north of the observed tornado damage path. The cycling of the 3DVAR analysis in Xue et al. [36] seemed to have a positive impact when compared to the single-time analysis used in this study. In contrast, our study shows that the track of the storm is closer to the actual tornado damage path when assimilating all available radar data.

4. Conclusion

In this study, the ARPS 3DVAR system with a cloud analysis package is used to assimilate radar observations for the 8 May 2003 tornadic supercell thunderstorm in central Oklahoma. The analysis includes observations from five radars, four operational WSR-88D radars (KTLX, KVNX, KINX, and KFDR) and one terminal Doppler weather radar (KOKC). Results show that the whole storm can be properly analyzed by including radar observations further away from the storm in addition to those located close to the storm. This aspect has not received enough attention in previous studies of convective storm data assimilation. The radars located further away from the storm are able to observe the upper levels of the storm that cannot be captured by the radars located very close to the storm owing to the “cone of silence”, an unobservable volume centered directly above the radar due to the upper limit of the scanning elevation angles. In analyses performed using data from one or two radars located close to the storm, the vertical vorticity and upper-level divergence structure are not analyzed as well when compared to the control analyses using all five radars.

The resulting analyses are then used as initial conditions for short-term storm-scale forecasts. Forecasts initialized from analyses using data from fewer radars proved to be less accurate. These forecasts had noticeably lower values of maximum low-level vertical vorticity and the track of maximum vorticity deviated from the observations more than that from a forecast using data from five radars. The computational efficiency of the 3DVAR technique makes it very suitable for short-term operational storm-scale forecasts, as has also been demonstrated by real-time radar data assimilation experiments over the continental United States [40].