A diagnostic pressure equation constraint has been incorporated into a storm-scale three-dimensional variational (3DVAR) data assimilation system. This diagnostic pressure equation constraint (DPEC) is aimed to improve dynamic consistency among different model variables so as to produce better data assimilation results and improve the subsequent forecasts. Ge et al. (2012) described the development of DPEC and testing of it with idealized experiments. DPEC was also applied to a real supercell case, but only radial velocity was assimilated. In this paper, DPEC is further applied to two real tornadic supercell thunderstorm cases, where both radial velocity and radar reflectivity data are assimilated. The impact of DPEC on radar data assimilation is examined mainly based on the storm forecasts. It is found that the experiments using DPEC generally predict higher low-level vertical vorticity than the experiments not using DPEC near the time of observed tornadoes. Therefore, it is concluded that the use of DPEC improves the forecast of mesocyclone rotation within supercell thunderstorms. The experiments using different weighting coefficients generate similar results. This suggests that DPEC is not very sensitive to the weighting coefficients.

1. Introduction

A dynamic consistent initial condition is very important for making a quality storm-scale numerical weather prediction (NWP) forecast. For this purpose, a large number of studies have been focused on utilizing high-resolution radar data to provide better storm-scale initial conditions (e.g., [17]). Since radars primarily observe the radial velocity and reflectivity, most state variables have to be “retrieved” in the data assimilation (DA) process. This makes the assimilation of radar data a very challenging problem.

Three-dimensional variational (3DVAR), four-dimensional variational (4DVAR), and ensemble Kalman filter (EnKF) methods have been applied to the previously mentioned radar DA problem. The 4DVAR method uses a NWP model as a strong constraint and hence naturally produces a dynamically consistent analysis. Sun and Crook [8, 9] and Sun [10] have shown encouraging results using a 4DVAR cloud model. However, it is very difficult to develop and maintain complex adjoint codes for NWP models. Complex ice microphysics, which are important for storm-scale applications but contain discontinuities and strong nonlinearities, introduce more difficulties in this situation. All of these difficulties limit the adoption of the 4DVAR method in storm-scale NWP operations. The EnKF technique is expected to generate similar analysis quality as 4DVAR but avoids the coding of an adjoint model. Many EnKF-based radar DA studies have been carried out in recent years (e.g., [1124]). These studies have shown great potential for the EnKF method. However, EnKF is not as mature as the variational methods and so far successful applications to real data assimilation problems are still limited. Computationally it has similar cost as the 4DVAR approach.

The 3DVAR method is more efficient than 4DVAR and EnKF in terms of computation cost. For this reason, 3DVAR has been applied in many real-time applications. Past studies (e.g., [46, 25]) have used 3DVAR to assimilate radar data for the prediction of tornadic supercell storms. The advanced regional prediction system (ARPS, [2628]) 3DVAR system, and its cloud analysis package have been used to make continental-US-scale real-time weather predictions at up to 1 km resolution [29, 30]. However, the 3DVAR method is still challenged by its theoretical suboptimality due to the use of a static diagonal background error covariance matrix and the lack of suitable balances among model variables. Research has been done to address this problem. For example, hybrid approaches (e.g., [3135]) have been proposed to provide flow-dependent background covariances obtained from a forecast ensemble, for the variational system. Because of the use of an ensemble, the method is still computationally expensive and requires further research.

Suitable weak constraints can also be built into the 3DVAR cost function to improve the balance among model variables and help spread observational information to state variables which are not directly observed. Gao et al. [3638], Hu et al. [4, 5], and Hu and Xue [39] incorporated an anelastic mass continuity equation into the ARPS 3DVAR system in the form of a weak constraint to couple three wind components. Ge et al. [7] further incorporated the diagnostic pressure equation into the 3DVAR cost function in addition to the previously mentioned mass continuity equation constraint (MCEC). The main goal of this diagnostic pressure equation constraint (DPEC) is to improve the consistency between dynamic and thermodynamic fields. Ge et al. [7] demonstrated, using observational system simulation experiments (OSSEs), that DPEC can improve the data assimilation results slightly after a given period of data assimilation. The application of DPEC to a real tornadic supercell thunderstorm case showed that DPEC can improve the forecast in terms of the general evolution of storm cells and mesocyclone rotation near the time of the observed tornado.

For the real case study by Ge et al. [7], only radial velocity data was assimilated. In this paper, we will further examine the impact of DPEC by applying it to two real tornadic supercell thunderstorm cases by assimilating the radial velocity and radar reflectivity data together. The rest of the paper is organized as follows. Section 2 will brief the updated ARPS 3DVAR system with the implementation of DPEC. Section 3 applies the updated system to the 5 May 2007 Greensburg, Kansas tornadic supercell thunderstorm case, while Section 4 applies the system to the 8 May 2003 Oklahoma City tornadic supercell thunderstorm case. The summary and future work will be presented in Section 5.

2. The Scheme for the 3DVAR System

A 3DVAR system within the ARPS model framework  [2628] has been developed and applied to the assimilation of weather radar and other data [46, 29, 36, 38, 39]. The system consists of two components: one is the 3DVAR subsystem, whose purpose is to assimilate radar radial velocity data as well as other conventional observations in a three-dimensional variational framework; the other is the cloud analysis subsystem whose purpose is to assimilate the radar reflectivity data based mainly on semiempirical rules. The cloud analysis system not only updates the hydrometeor fields, but also can adjust the in-cloud temperature and water vapor fields.

2.1. The 3DVAR Subsystem

In the 3DVAR subsystem, the cost function, , is written as the sum of the background () and observational () terms plus a penalty or equation constraint term ():

Following the standard notion of Ide et al. [40], and are the analysis and background state vectors, and is the observation vector. and are the background and observation error covariance matrices, respectively. is the nonlinear observation operator. To improve the conditioning of the minimization problem and avoid the need for the inverse of , a new control variable is introduced, which is related to the analysis increment according to where is the square root of ; that is, .

In terms of , the background term becomes Consequently, the minimization is performed in the space of . The recursive filter proposed by Purser et al. [41, 42] is used to model the effect of the background error covariance, or more precisely, the square root of . Currently, in our 3DVAR system, the background state vector, , can be provided by a sounding profile, a previous ARPS model forecast, or a forecast from another model. The analysis vector contains the three wind components (, , and ), potential temperature (), pressure (), and water vapor mixing ratio (). The observations include Doppler radar radial velocity, single-level (such as surface observations), and multiple-level conventional observations (such as those of rawinsondes and wind profilers). For the study in this paper, only Doppler radial velocity data is used by the 3DVAR subsystem and a 2 m s−1 observation error variance is assumed for the data.

Term in (1) includes any penalty or equation constraint terms. Currently, it includes two terms as defined in the following:

The first term is intended to minimize the 3D anelastic mass divergence so as to provide the key coupling among the three wind components. The definition and impact of this constraint have been investigated by Gao et al. [36, 38] and Hu et al. [5].

The second term is the DPEC term defined as follows: where The vector is the forcing term of the vector Euclidian momentum equation. The includes hydrometeor mixing ratios. The , , and are unit vectors in the , , and directions. The overbar represents base state and the primed variables are perturbations from a base state, is the acoustic wave speed, and is the ratio of the gas constants for dry air and water vapor. The Coriolis coefficients are and , where is the angular velocity of the earth and is latitude. The terms , , and contain the subgrid scale turbulence and computational mixing terms in the , , and directions, respectively. When the mass continuity equation is applied, (6) becomes where represents the right-hand side (R.H.S) of (6).

Equation (6) is derived by applying the divergence operator to the three momentum equations of the ARPS model [26]:

The in (5) is the error covariance matrix associated with the DPEC constraint, which is assumed to be diagonal with empirically defined constant diagonal elements as the variances. The inverse of is called the weighting coefficient and determines the relative importance of the DPEC constraint and its optimal value can be determined through numerical experiments, similar to the way to determine certain weights in cloud-scale variational data assimilation systems (e.g., [3]). Usually, the constraint terms with their weights should be of similar orders of magnitude as other terms in for them to be effective.

2.2. The Cloud Analysis Subsystem

The cloud analysis subsystem is based on the local analysis and prediction system (LAPS, [43]) with significant modifications by Zhang et al. [44], Brewster [45], and Hu et al. [4]. It is used to assimilate radar reflectivity data into the model. It updates the following model fields: rain water mixing ratio, snow mixing ratio, hail mixing ratio, cloud ice mixing ratio, cloud water mixing ratio, water vapor mixing ratio, and temperature. The general procedure is as follows. For each model grid point inside an area with observed reflectivity, a precipitation type (rain, snow, or hail) is first determined according to the reflectivity observation and the background state. After the precipitation type is determined, its mixing ratio is computed using reflectivity equations that link the precipitation species with reflectivity data. The cloud water and cloud ice mixing ratios are estimated by assuming that air parcels ascend moist-adiabatically from cloud base to cloud top. The water vapor mixing ratio is adjusted so that the air is saturated in precipitation area. The temperature field is also changed in order to make the in-cloud temperature consistent with the changed cloud fields. Readers can refer to Hu et al. [4] for more details.

2.3. Connection between the Two Subsystems

Under the context of ingesting radar data alone (radial wind and reflectivity), the analysis variables in the 3DVAR subsystem are the three wind components , , , and the analysis variables in the cloud analysis subsystem can be potential temperature , water vapor mixing ratio , rain water mixing ratio , snow water mixing ratio , hail mixing ratio , cloud water mixing ratio , and ice mixing ratio . Currently, the cloud analysis subsystem is a follow-up step after the 3DVAR subsystem finished running. These two subsystems are separate from each other and there is no suitable coupling between the wind fields and the thermodynamic fields. Therefore, there may be inconsistencies between the different model variables in the data analysis. These inconsistencies may harm the quality of subsequent data assimilation cycles and the ensuing forecast.

To alleviate these kinds of inconsistencies, we propose that the cloud analysis subsystem is done first when it is used in the assimilation runs. The results from the cloud analysis package will then be treated as pseudoobservations and be ingested, as well as the radar radial velocity data, by the 3DVAR subsystem. DPEC will then act to couple all model variables so as to help improve the balance between the dynamic and thermodynamic fields. In this way, it is expected that a more dynamically consistent analysis will be achieved and the following data assimilation cycles and the subsequent forecast will be improved. In practice, the hydrometeors are updated directly by the cloud analysis subsystem. The temperature and water vapor mixing ratio fields are updated in the 3DVAR subsystem by treating the corresponding results obtained from the cloud analysis subsystem as pseudoobservations. The assumed error variances for the pseudoobservations and are 0.5 K and  g (kg)−1, respectively.

3. The 5 May 2007 Greensburg Tornadic Supercell Storm Case

The 5 May 2007 Greensburg, Kansas (KS) tornadic thunderstorm complex produced 18 tornadoes in the Dodge City area and additional 47 tornado reports in Kansas, Nebraska, and Missouri. One tornado was the strongest observed in recent years. This tornado started moving through Greensburg at 0245 UTC 5 May 2007 (2145 CDT 4 May) and destroyed over 90% of the town. The tornado damage was rated at EF5—the highest rating on the Enhanced Fujita scale [46]. A detailed description of the supercell that spawned this tornado and its environmental setting can be found in Bluestein [47] and Stensrud and Gao [6].

For this real data case, we used a 3 km grid spacing with 200 × 200 grid points in the horizontal. The ARPS model domain is shown in Figure 1. The domain was selected with sufficient coverage to contain the principal features of interest while maintaining some distance between the primary storms and the lateral boundaries. The model uses 47 terrain following vertical layers with nonlinear vertical stretching via a hyperbolic tangent function that yields a spacing of 100 m at the ground and expands to approximately 800 m at the top of the domain. The Lin three-ice microphysical scheme [48] was used together with a 1.5-order turbulent kinetic energy subgrid parameterization. A wave radiation condition was applied at the top boundary and rigid-wall conditions were applied to the bottom boundary.

The impact of the DPEC will be discussed in terms of the quality of ensuing forecasts instead of the analysis because no truth or high-resolution observation is available for verification of the analysis. Four experiments were conducted for this case (Table 1). The first experiment did not include DPEC in and will be referred to as experiment NoDP1. The second experiment used DPEC with the DP weighting coefficient of 1.0E8 and is referred to as experiment DP1. The third and fourth experiments were similar to DP1 except that the DP weighting coefficients were multiplied and divided by 5, respectively. They are referred to as experiments DP1m5 and DP1d5, respectively. In all of the previously four experiments the mass continuity equation constraint was used with the MC weighting coefficient of 1.0E8.

For all the previously mentioned four experiments, data from five radars at Dodge City, Kansas (KDDC), Vance Air Force Base, Oklahoma (KVNX), Wichita, Kansas (KICT), Oklahoma City, Oklahoma (KTLX), and Amarillo, Texas (KAMA), was used (Figure 1). A quality control procedure was applied before the use of the radar data, which included clutter removal and velocity dealiasing using SOLOII software from the National Center for Atmospheric Research (NCAR). The initial analysis background and the boundary conditions came from the mean of a mesoscale ensemble assimilation system run at 30 km grid spacing [6]. While Stensrud and Gao [6] performed a 3DVAR analysis at only one time before the launch of the forecast, the present study performed cycled 3DVAR analyses with a 1 h long assimilation period before the forecast. A five-minute ARPS forecast followed each analysis, and this process was repeated until the end of the 1 h assimilation period. From the final analysis, a 1 h forecast was launched. In this way, each experiment consisted of a 1 h assimilation period (from 0130 and 0230 UTC) and a 1 h forecast period (0230–0330 UTC).

We now focus the discussion on the dominant supercell thunderstorm at the southernmost portion of the storm complex, which produced the EF-5 tornado that hit the Greensburg area between 0245 UTC–0305 UTC. A hook echo signature was evident at 0230 UTC. As the storm reached Greensburg, the hook echo signature became less prominent due to the tornado moving in a more northerly direction and toward the storm’s main core. During this period, the radar velocity observations indicated strong cyclonic rotation associated with the violent tornado. The entire storm complex moved gradually toward the northeast. After passing the town of Greensburg, a second tornado, rated EF-3, developed at the end of the Greensburg tornado’s path just northeast of the town [46]. A radar reflectivity mosaic was created from the aforementioned five WSR-88D radars by interpolating reflectivity data from all radars onto model grid points and keeping the largest reflectivity value for each grid point. The reflectivity mosaic was then used for forecast verification. The evolution of the storm as indicated by the radar reflectivity mosaic at 2 km MSL is shown in Figure 2 from 0230 to 0330 UTC every twenty minutes. Note that the hook echo is not evident in these figures due to the use of 3 km resolution and a smoothing procedure applied in the mosaic generating process.

To demonstrate the impact of DPEC, we investigated these data assimilation experiments ingesting both the radial velocity data and reflectivity data. Figures 2(e)–2(l) show the reflectivity, horizontal wind vector, and vertical vorticity at  km MSL from 0230 UTC to 0330 UTC every 20 minutes for the NoDP1 and DP1 experiments. It is shown that after 1 hour of data assimilation (Figures 2(e) and 2(i)), the storm had already spun up in terms of the reflectivity pattern. The reflectivity pattern, strength, and location agree well with the observed values (Figure 2(a)). A rotating circulation and a strong vertical vorticity column are collocated at the observed hook-echo region. The storm then moves gradually toward the northeast, which also agrees with the observations. After 0300 UTC, the predicted storm moves faster than what was observed. In spite of this, both NODP1 and DP1 still made reasonable forecasts in terms of the general evolution of the major storm. DP1d5 and DP1m5 produced very similar forecasts as DP1 and are therefore not shown in Figure 2.

In Figure 2, in terms of the reflectivity pattern, there is no significant difference in the general evolution of the storm between the NODP1 and DP1 experiments. The computed forecast scores (not presented here) also show little difference, consistent with the previous result. However, there is some difference in the predicted low-level mesocyclone rotation as indicated by larger maximum vertical vorticity in Figures 2(j), 2(k), and 2(l) than that in Figures 2(f), 2(g), and 2(h). As a further demonstration, Figure 3 shows the time series of the maximum vertical vorticity below two kilometers every minute from 0230 UTC to 0330 UTC for all four experiments. It is illustrated in Figure 3 that beginning at 0245 UTC and through the end of the forecast, the low-level maximum vertical vorticity from the experiments applying DPEC (the red, blue and green lines) is much larger than that from the NODP1 experiment (the black line). Our detailed examinations show that larger low-level vertical vorticity corresponds to a better-defined mesocyclone vortex, which is stronger and deeper than those with smaller values of vertical vorticity. This kind of behavior is very similar to findings in Ge et al. [7]. As an example, Figure 4 presents the vertical vorticity at the vertical cross-section through the center of the major storm at  km at 0250 UTC 5 May 2007. It is noticeable that the experiments using DPEC (Figures 4(b), 4(c), and 4(d)) predicted stronger and deeper rotation than the “NODP1” experiment (Figure 4(a)). Therefore, it can be concluded that, for the experiments here, although the use of DPEC does not evidently improve the forecast of the general evolution of the major storm in terms of the reflectivity pattern, it does help improve the forecast of the mesocyclone rotation associated with the observed Greensburg tornado.

4. The 8 May 2003 Oklahoma City Tornadic Supercell Storm Case

During the late afternoon on 8 May 2003, a major tornado hit the southern Oklahoma City metropolitan area (Figure 5). It first touched down at Moore, a suburban city close to and south of Oklahoma City, then traveled east north-east through south of Oklahoma City to Choctaw. The life span of the tornado was about 28 minutes from 2210 UTC to 2238 UTC. It caused up to F4 (Fujita scale) damages but no deaths. The tornado is hereafter referred to as the OKC tornado and the parent storm as the OKC tornadic thunderstorm.

The synoptic environment on 8 May 2003 over Oklahoma was very favorable for the development of supercell storms and even tornadoes, as discussed by Hu and Xue [39] and Romine et al. [49]. The low-level flow was southerly over Oklahoma all day. A meridionally oriented dryline moved eastward approaching Moore, Oklahoma. A large amount of potential instability with 4004 J kg−1 convective available potential energy (CAPE), 1 J kg−1 convective inhibition (CIN) and about 25 m s−1 vertical shear over the lowest 6 km was present in the 1800 UTC 8 May Norman, Oklahoma (OUN) sounding. All of these conditions indicated that there was a high possibility for tornadic supercell thunderstorms to develop.

At about 2030 UTC, the first sign of the OKC tornadic storm showed up as a weak echo in the KTLX radar reflectivity field. By 2101 UTC, the storm had developed into a strong cell. In the following hour, the storm grew rapidly and moved northeastward. By 2201 UTC, the storm displayed an obvious hook echo signature at its southwestern end. The hook echo at this time was located just northwest of Moore, only several miles away from the center of the city. The pronounced hook echo signature was present until at least 2235 UTC while the parent supercell storm propagated east northeastward. The storm began weakening at 2240 UTC and dissipated by 0020 UTC 9 May. In addition to the OKC tornadic thunderstorm, there were three other short-lived storms (not shown). Here, we will just focus on the dominant thunderstorm. Figures 6(a), 6(b), 6(c), and 6(d) show the general evolution of the major thunderstorm every twenty minutes from 2200 UTC to 2300 UTC represented by the radar reflectivity mosaic at 2 km MSL.

All experiments were conducted with a horizontal resolution of 3 km. There were 195 grid points in both and directions. In the vertical direction, a stretched grid scheme was used. It contained 53 layers with an average grid spacing of 400 m, stretching from about 20 m at the surface to 770 m at the model top. The model domain is shown in Figure 7. It covers nearly the entire state of Oklahoma. The evolution of the 8 May 2003 Oklahoma City tornadic supercell thunderstorm was roughly at the center of the domain. The four WSR-88D radars KTLX, KVNX, KINX, and KFDR and their associated coverage region are also shown in Figure 7. The outline near the KTLX radar is the damage path of the 8 May 2003 OKC tornado.

The ARPS system was used as the prediction model. The parameterization schemes and vertical boundary conditions used in Section 3 were adopted for the present experiments. The initial first guess and the lateral boundaries were provided by a 9 km data assimilation experiment. This 9 km experiment was done in the same way as in Hu and Xue [39]. It assimilates rawinsonde data and wind profiler data every hour for a total of six hours. The Eta model analysis and forecast provide the background and lateral boundaries for the 9 km experiment.

Similar to before, four experiments (Table 1), that is, NODP2, DP2, DP2d5, and DP2m5, were conducted in order to examine the impact of DPEC and the sensitivity of DPEC to different weighting coefficients. The assimilation experiments start at 2100 UTC and assimilate radial velocity data and radar reflectivity data every 5 minutes in a cycled manner similar to the procedure described in Section 3. The assimilation window was 1 h long and the final analysis was at 2200 UTC. From the final analysis, a 1 h free forecast (2200 UTC–2300 UTC) was made.

Figures 6(e)–6(l) show the reflectivity, horizontal wind vector, and vertical vorticity at  km MSL from 2200 UTC to 2300 UTC every 20 minutes for the NODP2 and DP2 experiments. After 1 h of data assimilation, the storm had been successfully spun up (Figures 6(e) and 6(i) versus Figure 6(a)). The area of strong vertical vorticity is located where a hook echo is observed. The storm then moves east-northeastward. The direction and speed of the predicted storm are very close to what was observed (Figures 6(f), 6(g), 6(h) and 6(j), 6(k), and 6(l) versus 6(b), 6(c), and 6(d)). Therefore, both NODP2 and DP2 performed well in predicting the general evolution of the storm. DP2d5 and DP2m5 made similar forecasts to DP2 (not shown).

Comparing Figures 6(j), 6(k), and 6(l) (for DP2) with Figures 6(f), 6(g), and 6(h) (for NODP2), we can see that DP2 predicted larger low-level vertical vorticity than NODP2. This is further confirmed by examining the evolution of low-level vertical vorticity. Figure 8 shows the time series of the maximum vertical vorticity below two kilometers every one minute from 2200 UTC to 2300 UTC for all four experiments. It is illustrated in Figure 8 that after 22:17 UTC and until the end of the forecast, the low-level maximum vertical vorticity from the experiments applying DPEC (the red, blue, and green lines) was generally larger than that from the NODP2 experiment (the black line). As mentioned before, larger low-level vertical vorticity corresponds to a better-defined mesocyclone vortex, which is stronger and deeper. As an example, Figure 9 presents the vertical vorticity at the vertical cross-section through the center of the storm at  km at 2220 UTC, which is during the tornado touchdown period. The experiments using DPEC (Figures 9(b), 9(c), and 9(d)) predict a deeper column of high vertical vorticity (>0.008 s−1), extending from as low as 1.0 kilometers to as high as 8.5 kilometers. The region of high vertical vorticity (>0.008 s−1) predicted by the NODP2 experiment is mainly in the middle part of the atmosphere, roughly from 2.4 kilometers to 6.0 kilometers. Therefore, it can be concluded that the use of DPEC helps make a better forecast of low-level mesocyclone rotation.

5. Summary and Conclusions

A diagnostic pressure equation was added into the ARPS 3DVAR system as a weak constraint with the goal of coupling the dynamic and thermodynamic variables so as to improve the analysis of convective storms and their subsequent forecast. The updated ARPS 3DVAR system was tested using OSSEs in Ge et al’s. [7] and applied to a real tornadic supercell case where only radial velocity has been assimilated. This study further applied the newly updated ARPS 3DVAR system to tornadic supercell thunderstorm studies by assimilating both the radial velocity and radar reflectivity data.

For both the 5 May 2007 Greensburg tornadic supercell storm case and the 8 May 2003 Oklahoma City tornadic supercell storm case, four data assimilation experiments were conducted with three of them using different DP weighting coefficients and the other one without DPEC. The four experiments assimilated the same amount of observations from multiple Doppler radars and imposed the mass continuity equation constraint.

After 1 h of intermittent data assimilation, it was found that DPEC did affect the final analysis. However, since there is no reliable high-resolution analysis of the storm, it is not easy to tell directly which analysis is better. The evaluation of the benefit of DPEC to radar data assimilation in these real cases is examined mainly based on the ensuing forecasts.

It was demonstrated that the experiments using DPEC generally predict larger low-level vertical vorticity than the experiments not using DPEC. Therefore, it is concluded that the use of DPEC improves the forecast of supercell mesocyclone rotation of the major thunderstorm. The experiments using different weighting coefficients generated similar results. This suggests that DPEC is not very sensitive to the weighting coefficients, although very small values should still be avoided as found in Ge et al’s. [7].

Overall, the addition of DPEC in the ARPS 3DVAR system had a positive impact on storm-scale 3DVAR data assimilation of Doppler radar data and on the subsequent forecast. In the future, the system needs to be tested with more real data cases, including tornadic supercell thunderstorms and other storm-scale phenomena, to further demonstrate the robustness of these conclusions.


This work was supported by NSF Grants ATM-0331756, ATM-0738370, NSF ATM-0530814, NSF AGS-0802888, and NOAA’s Warn-on-Forecast Project. The third author was also supported by NSF OCI-0905040, AGS-0941491, and AGS-1046171. Robin Tanamachi kindly provided manually dealiased radial velocity dataset of the Dodge City, Kansas (KDDC) radar. Computations were performed at the Pittsburgh Supercomputing Center (PSC) and Oklahoma Supercomputing Center for Research and Education (OSCER). The authors appreciate two anonymous reviewers and the Editor Louis Wicker for their very thoughtful comments which helped significantly to improve the quality of this paper from its original form. The authors also thank Jacob Carlin and Edward Natenberg for proofreading the paper.