This paper presents a case study on the assimilation of observations from multiple Doppler radars of the Next Generation Weather Radar (NEXRAD) network. A squall-line case documented during the International H2O Project (IHOP_2002) is used for the study. Radar radial velocity and reflectivity observations from four NEXRADs are assimilated into a convection-permitting model using a four-dimensional variational data assimilation (4DVAR) scheme. A mesoscale analysis using a supplementary sounding, velocity–azimuth display (VAD) profiles, and surface observations from Meteorological Aerodrome Reports (METAR) are produced and used to provide a background and boundary conditions for the 4DVAR radar data assimilation. Impact of the radar data assimilation is assessed by verifying the skill of the subsequent very short-term (5 h) forecasts.
Assimilation and forecasting experiments are conducted to examine the impact of radar data assimilation on the subsequent precipitation forecasts. It is found that the 4DVAR radar data assimilation significantly reduces the model spinup required in the experiments without radar data assimilation, resulting in significantly improved 5-h forecasts. Additional experiments are conducted to study the sensitivity of the precipitation forecasts with respect to 4DVAR cycling configurations. Results from these experiments suggest that the forecasts with three 4DVAR cycles are improved over those with cold start, but the cycling impact seems to diminish with more cycles. The impact of observations from each of the individual radars is also examined by conducting a set of experiments in which data from each radar are alternately excluded. It is found that the accurate analysis of the environmental wind surrounding the convective cells is important in successfully predicting the squall line.
Quantitative precipitation forecasting (QPF) is a challenging problem (Fritsch and Carbone 2004) and its skill has been historically low. The QPF skill also varies seasonally, with the warm season skill being significantly lower than that in the cold season. The low skill and the relative lack of progress for warm season QPF can be attributed, to a large degree, to the fact that the heavy rain-producing convective processes tend to be poorly resolved by numerical weather prediction (NWP) models. Through an analysis of two warm season precipitation events, Heideman and Fritsch (1988) found that nearly all the precipitation was convective and over 80% of significant precipitation was associated with thunderstorms. To adequately resolve the convective processes and improve short-term QPF, prediction models at resolutions much higher than currently used are required. In addition, the model must be initialized by observations that can describe the mesoscale and convective-scale state of the atmosphere.
Among the number of mesoscale observation platforms that can be used for high-resolution analysis are radar networks, which are now available in several countries in the world. Active research has been carried out to investigate the feasibility of using radar observations or their derived quantities in NWP models and their impact on QPF. Doppler radar provides 3D observations of radial wind and reflectivity (associated with hydrometeors). The derived quantities include the velocity–azimuth display (VAD) wind profile, which is representative of the mean of the radar coverage area, and precipitation rate. The assimilation of the derived quantities has been studied using mesoscale models for VAD (e.g., Michelson and Seaman 2000) wind and precipitation rate (Jones and Macpherson 1997; Koizumi et al. 2005).
Recent studies have shown the feasibility of assimilating volumetric Doppler radar observations (radial velocity and reflectivity) into cloud-resolving models (Weygandt et al. 2002a, b; Sun and Crook 1997, 1998; Snyder and Zhang 2003; Dowell et al. 2004; Xiao et al. 2005, 2007; Hu et al. 2006; Sun 2005a). All of these studies have demonstrated a positive impact on very short-term QPF when a NWP model was initialized using radar observations. Such studies, although encouraging, have typically used observations from a single radar and performed in a very small domain for simulation of isolated thunderstorms. In this study, we extend previous studies to multiple Doppler radars and to larger convective systems. The purpose is to examine whether volumetric radial velocity and reflectivity observations from multiple Weather Surveillance Radar-1988 Dopplers (WSR-88Ds) can be successfully assimilated into a cloud-permitting model and whether the assimilation can improve the subsequent short-term QPF. A case study of a convective system observed during the International H2O Project (IHOP_2002) was performed as a first step in fulfilling this purpose. A similar study by Xiao and Sun (2007) was recently conducted using the Weather Research and Forecasting (WRF) model and its three-dimensional variational data assimilation (3DVAR) system.
IHOP_2002 is a multiagency and international field experiment conducted in the areas of Texas, Oklahoma, and Kansas. The overarching goal of the experiment was to obtain the characterization of the four-dimensional distribution of water vapor for improving the understanding and prediction of convection (Weckwerth et al. 2004). The IHOP_2002 field campaign took advantage of existing observing systems in the Southern Great Plains of the United States, including the well-distributed WSR-88D units in the region. There were also a large array of specially deployed instrumentation, including, but not limited to, aircraft, mobile radars, mobile soundings, dropsondes, lidar, and radiometers. Several studies have explored the role of assimilating some of the nonconventional observations obtained during IHOP on convection initiation and prediction (e.g., Childs et al. 2006; Wulfmeyer et al. 2006). In this study, our focus is on the assimilation of observations from the WSR-88D network in the IHOP_2002 domain using a four-dimensional variational data assimilation (4DVAR) system. The case being studied is a squall line that occurred during 12–13 June 2002. The benefit of the assimilation is evaluated by examining the subsequent prediction of the formation and evolution of the squall line.
This paper is organized as follows. In section 2 the IHOP_2002 squall line case and its large-scale environment are described. Section 3 gives an overview of the data assimilation system. The experimental design for the current study is described in section 4. In section 5, we present the data assimilation result, and evaluate the forecast impact. Results from a number of sensitivity experiments are also presented and discussed in this section. Summary and conclusions are given in the last section.
2. Overview of the IHOP_2002 12–13 June 2002 case
Several convective systems occurred during 12–13 June 2002 in the IHOP_2002 domain. This study focuses on the squall line that initiated east of a triple point (intersection of an outflow and a dryline; marked by “T” in Fig. 1a) near the Oklahoma–Kansas border around 2100 UTC and propagated southeastward. Figure 1 shows the convective activity at four different times by the merged reflectivity field at 0.5° elevation angle from multiple WSR-88Ds and the National Center for Atmospheric Research’s (NCAR’s) S-band dual-polarization Doppler radar (S-Pol). The merging of the radar reflectivity was done by choosing the maximum reflectivity at each grid point from data values given by different radars. A composite map of surface wind observations (white barbs) from several surface networks is overlaid. The blue wind barbs are surface observations from the Meteorological Aerodrome Reports (METAR) stations only and are used in the data assimilation. The radar locations are marked by the radars’ names in red letters. Several mesoscale boundaries, including a nearly stationary outflow boundary from the previous day’s convection (orange broken line in Fig. 1a), a dryline (pink line), and a weak cold front (brown line), are observed. A mesoscale circulation (mesoscale low) around the triple point is evidently shown by the surface observations. The warm, moist air mass to the east of the dryline and south of the outflow boundary contained CAPE near 2000 J kg−1 and little convective inhibition (CIN) around 2100 UTC. The air mass to the north of the outflow boundary is slightly cooler.
Convective cells are developing east of the triple point and are clearly visible at 2120 UTC (Fig. 1a). The thunderstorm near the triple point intensified and new storms were developed near the Oklahoma–Kansas border as the circulation associated with the mesoscale low strengthened (Fig. 1b). Damaging winds and large hail were reported in conjunction with some of the storms. The severe storm near the triple point produced golf ball–sized hail, maximum outflow wind speeds exceeding 30 m s−1, flashing flood, and at least one tornado. By 0200 UTC, a well-organized squall line had developed (Fig. 1c) and it then moved southeastward. The squall line is evidently reduced in spatial extension at 0400 UTC (Fig. 1d) and completely dissipated around 0900 UTC. The stage-IV precipitation analysis of the National Centers for Environmental Prediction reported a maximum 3-h rainfall accumulation of 89.5 mm at 0300 UTC in north Oklahoma. Our focus in the current study is not on the prediction of the initial storm development, but the prediction of the evolution of existing storms and the development of new ones from outflow interaction. Predicting storm initiation requires a higher resolution and a model that includes better physics and numerics that can resolve the finescale structure of the atmosphere.
Figure 2 presents a more detailed look at the formation and evolution of the squall line in the smaller domain indicated by the yellow rectangle in Fig. 1a. This plot only shows reflectivity above 20 dBZ to give a clearer picture of the convective system. The domain displayed in Fig. 2 is the data assimilation and model forecast domain used in this study. The reflectivity field shown in Fig. 2 is merged using the four radars (KDDC, KICT, KVNX, and KTLX), with their locations indicated in the figure. The hourly reflectivity field plotted in Fig. 2 clearly shows that the convective cells A and B (marked in Fig. 2a) moved southward, intensified, and finally merged into a well-defined squall line by 0154 UTC. The squall line then propagated southeastward. It should be noted that the merged reflectivity field shown in Fig. 2 is produced using a different method from that in Fig. 1. A technique based on the Cressman interpolation is used to merge the reflectivity and used as the first guess for 4DVAR described in the next section.
3. The radar data assimilation system
The radar data assimilation system used in this study is NCAR’s Variational Doppler Radar Analysis System (VDRAS; Sun and Crook 1997, 2001; Crook and Sun 2002; Sun 2005a). In this section, we briefly describe the system, followed in the next section by a description of the experimental design for the current study. The reader is referred to the above publications for more detailed description of the 4DVAR system and the prediction model. The dry version of VDRAS has been implemented in a number of field projects, including the World Meteorological Organization (WMO) Sydney, Australia, 2000 Olympic Project (Crook and Sun 2002), to produce real-time low-level wind analysis for use in nowcasting thunderstorms. The full VDRAS with warm rain microphysics has been used in a number of studies on the retrieval of wind, thermodynamical, and microphysical fields in convective storms (Sun and Crook 1998; Wu et al. 2000). It was also used to study the prediction of isolated thunderstorms through assimilation of observations from single Doppler radar (Warner et al. 2000; Sun 2005a).
a. Description of VDRAS
VDRAS was designed for the assimilation of high-resolution Doppler radar observations using a 4DVAR scheme. Although the main objective of VDRAS is to assimilate high-resolution remotely sensed observations from radar or lidar, it also makes use of data from radiosonde, profiler, and surface networks, wherever these data are available, through a mesoscale background analysis prior to the 4DVAR radar data assimilation.
The central process of VDRAS is the 4DVAR radar data assimilation, which includes a cloud-resolving model, the adjoint of the prediction model, a cost function, a minimization algorithm, and the specification of weighting coefficients. For details on each of these components, readers are referred to Sun and Crook (1997) and Sun (2005a). The general form of the cost function can be written by
where x represents the model state and y the observation. The subscript 0 stands for the beginning of the assimilation window and b the background forecast from previous cycle valid at the beginning of the current assimilation window. Because the model state can be uniquely determined from the initial conditions through model integration (assuming the boundary conditions are known), the control variables in (3.1) are the six prognostic variables (three velocity components, liquid water potential temperature, rainwater mixing ratio, and total water mixing ratio) at the initial time of the assimilation window. The symbol 𝗕 denotes the background covariance matrix and 𝗢 the observation covariance matrix. The symbol Jp denotes the spatial and temporal smoothness penalty term. The function of the smoothness penalty term is to ensure a smooth fit to the observations. Its exact form can be found in Sun and Crook (2001). With the recent implementation of a recursive filter (described later) in VDRAS, the Jp term is only necessary for a cold-start cycle. The symbol Jmb denotes the mesoscale constraint that ensures the analysis does not deviate too far from a mesoscale background. The mesoscale background is produced by an analysis using all of the observations other than radar and will be explained in more detail later.
When the radial velocity and reflectivity observations from radars are considered, the cost function (3.1) can be written as
We have assumed there is no spatial error correlation between observations such that the observation term [second term in (3.2)] retains only the quadratic components. The symbols ηυ and ηq represent the inverse of the observation error variance for radial velocity and rainwater, respectively. The summation is over space (denoted by σ) and time. The variable υr is the radial velocity computed from the model velocity components, υro is the observed radial velocity, qr is the rainwater mixing ratio from the model, and qro is the rainwater mixing ratio estimated from the reflectivity observation using the formula (Sun and Crook 1997)
where Z denotes the reflectivity factor and ρ the air density. This Z–qr relation is obtained by assuming a Marshall–Palmer raindrop size distribution and has been used in several previous studies (Sun and Crook 1997, 1998; Xiao et al. 2005). Note that returns from Doppler radar include echoes both from raindrops and bugs existing in the clear boundary layer. VDRAS assimilates both types of echoes for radial velocity, but only raindrop echoes for reflectivity.
A recent modification of VDRAS is the implementation of the recursive filter. The recursive filter was introduced by Hayden and Purser (1995) and is used in several variational systems, including the fifth-generation Pennsylvania State University (PSU)–NCAR Mesoscale Model (MM5) and WRF 3DVAR system (Barker et al. 2004). A variable transformation is introduced to Eq. (3.1) by
The original control variable x is replaced by the new control variable z and can be computed from z by
The first term in the cost functions (3.1) and (3.2) becomes zTz after the control variable transformation. The adjoint for the original control variable x can be derived from the adjoint for the new control variable z by
From (3.5) and (3.6), it is clearly shown that the original control variable and its adjoint can be computed from the new ones by applying the square root of the background covariance matrix. Instead of defining the values of 𝗕 and computing the matrix 𝗕1/2, which is computationally challenging, the multipass recursive filter of Hayden and Purser (1995) is used to approximate 𝗕 by a nearly Gaussian, homogenous correlation function. The decorrelation length of the forecast error background was determined based on a statistical calculation using data collected in the summer of 2005 over a domain that covered the area of Illinois and Indiana. An ensemble of VDRAS analyses and 20-min forecasts were performed over a number of convective days. The differences between the analyses and the forecasts valid at the same times were calculated and statistical properties were analyzed. Based on the results, a horizontal decorrelation length of 32 km is used for velocities and temperature, and an 8-km length is used for microphysical variables. In the current study, the recursive filter is applied only in the horizontal directions. The statistical results also indicate a typical 20-min forecast error of 1.5 m s−1 for horizontal velocity components, 0.3 m s−1 for vertical velocity, 0.7°C for liquid water potential temperature, and 0.15 g kg−1 for rainwater.
b. Mesoscale background analysis
Although VDRAS aims at the assimilation of radar observations, it also relies on other types of data to provide a background analysis. It is well recognized that the WSR-88Ds have limited coverage ranges, typically 225 km for hydrometeor returns and 100 km for clear-air returns. Large data voids exist in the area without reflectors or beyond the radar coverage range. An important question is how a data assimilation system can combine these observations optimally and smoothly in a way that the convective-scale details, as observed from radar, are preserved while larger-scale observations are used to fill the gaps between radar coverage areas.
An approach that has been used in a number of high-resolution data assimilation systems is to assimilate radar data in a second step after an analysis that assimilates all other data types is done. For example, in the Local Analysis and Prediction System (LAPS) wind analysis, radar wind is analyzed in the second pass using the same successive correction technique but a different weighting function (Albers 1995). Crook and Sun (2002) described a technique in which a two-step analysis procedure was used. First, observations from profilers, mesonet stations, and VAD analyses were analyzed using a Barnes scheme. This analysis was then used as boundary conditions and background for a subsequent 4DVAR analysis that assimilates radar observations. In this study, we use the same technique as in Crook and Sun (2002), but the data used for the mesoscale analysis include METAR surface observations, an upper-air sounding from an Atmospheric Radiation Measurement Program (ARM) supplemental sounding, and VAD profiles. Although there are a number of other mesonets during IHOP, we chose to use only METAR because of its operational availability. The ARM special sounding was launched at 2028 UTC at Vici, Oklahoma (location is marked as VCI in Fig. 1a). The mesoscale analysis for the six control variables (3D wind, temperature, humidity, and rainwater) is obtained through different procedures. These procedures are described below. First, the mesoscale analysis procedure for horizontal wind is conducted in the following steps:
A VAD wind profile is produced from each of the radars using data from the third elevation angle (2.4°). Figure 3 shows the VAD profiles from the four radars. Note how both the wind speed and direction varies between these profiles.
The winds from the upper-air sounding are then used to extend each of the four VAD profiles to the model top.
These composite soundings are analyzed to the model grid using a Barnes interpolation scheme with a radius of influence of 200 km.
The surface observations are analyzed to the model grid using the Barnes interpolation scheme with a radius of influence of 2 times of the average station spacing (∼20 km).
A local linear least squares fit is applied vertically at each grid point to combine the surface analysis and the upper-air analysis.
The mesoscale analysis of temperature and humidity is obtained by blending the profiles of these two fields from the VCI sounding with the surface observations. First, the VCI sounding is applied to every horizontal grid point, and then the surface observations are analyzed to the model grid similarly to step 4 above; finally, the values of the temperature and humidity below 1 km are modified by assuming a linear change between the radiosonde observation at 1 km and the surface observations near the surface.
We admit that the above procedures for the mesoscale background analysis are rather crude. A more sophisticated procedure that involves output from mesoscale models is being developed and tested. The impact of the mesoscale background analysis obtained through this new procedure is being evaluated. In the current study, we use a supplemental sounding launched near Vici, Oklahoma. In operational applications, the data assimilation system can only rely on the operational observations. Therefore, the impact of the special sounding used in this study versus operationally available soundings is also being examined. The results will be reported in a follow-up paper.
c. Radar data quality control and preprocessing
Radar data quality control and preprocessing were performed using NCAR’s Sorted Position Radar Interpolation (SPRINT; Mohr and Vaughan 1979; Miller et al. 1986) and Custom Editing and Display of Reduced Information in Cartesian Space (CEDRIC; Mohr et al. 1986) software. For a summary of the procedure, the reader is referred to Sun (2005a). A difference from Sun (2005a) in the current study is the data interpolation. In Sun (2005a), the quality-controlled data were interpolated to the model Cartesian grid. In the current study, the quality-controlled data were left on their original constant-elevation surfaces in the vertical while being interpolated to the model Cartesian grid in the horizontal.
4. Experimental design and verification method
The yellow rectangle in Fig. 1a shows the domain in which the data assimilation and forecast experiments are conducted. The choice of the domain size (440 km × 400 km) is limited by the available computation resource. VDRAS was run on a Dell Linux dual-processor workstation with a 3.0-GHz CPU and a 4-GB memory. VDRAS was configured such that it fits within the maximum memory with a domain covering the entire squall line during the formation and development stage. The grid spacing is 4 km in the horizontal and 500 m in the vertical. At this resolution, the model can resolve organized convection but not individual cells; hence, the simulation might be viewed as “convection permitting” but not “convection resolving.”
VDRAS was set up to run with a cycling 4DVAR procedure. Figure 4 shows a schematic diagram illustrating the assimilation cycles. The typical volumetric scan frequency for the radars is 5–6 min, depending on the scanning pattern used for each radar. Each 4DVAR cycle has a preset temporal window of 12 min, followed by a 10-min short forecast cycle to produce a forecast background for the next assimilation cycle. After three assimilation cycles, a 5-h forecast is conducted. The short assimilation window of 12 min allows for a closer fit to observations and enables the analysis to capture the features that have shorter predictability length. However, the downside of a short assimilation window is that it is not favorable for producing larger-scale balances that are important for longer-term forecasts. However, we expect that the cycling with the short window can compensate this weakness to a certain degree. At the start time of each data assimilation cycle (i.e., 2158, 2220, and 2242 UTC, as indicated in Fig. 4), a new mesoscale analysis is performed using updated surface observations and VAD analysis.
A series of experiments were conducted to investigate the impact of observations from multiple WSR-88Ds on forecasting the squall line of IHOP_2002 12–13 June 2002. Table 1 summarizes all of the experiments. The first two experiments are conducted without 4DVAR radar data assimilation. In the NO4DV1 experiment, a mesoscale analysis is performed according to the procedure described in the last section. The forecast starts from the mesoscale analysis and there is no cloud and rain initialization. The NO4DV2 experiment is similar to NO4DV1 except that, in addition to the mesoscale analysis, the rainwater mixing ratio field derived from the composite radar reflectivity shown in Fig. 2 using Eq. (3.3) is employed to initialize the rainwater by direct insertion. (Note that the derived rainwater mixing ratio will be referred to as observed rainwater mixing ratio hereafter.) The next two experiments, CYC3_2254 and CYC3_2355, assimilate radar volumetric data with three 4DVAR cycles. In all of the experiments, the number after the underscore indicates the forecast starting time (UTC). Figure 3 shows the design of the cycles for CYC3_2254. Each 4DVAR cycle has an assimilation window of 12 min, followed by a 10-min forecast cycle. The CYC3_2355 cycles are the same as CYC3_2254, except that the cold start begins 1 h and 1 min later. These two experiments are designed to examine the impact of the radar data assimilation and also are used as controls for the next two groups of experiments. The third group of experiments was conducted to test the sensitivity of the assimilation and forecast with respect to different settings of assimilation cycles. There are five experiments in this group. The CYC1_2254 experiment is the same as CYC3_2254, except that only one cold-start cycle is run. CYC3_2355 also has three continuous cycles, but the free forecast starts at 2355 UTC. CYC6_2355 is the same as CYC3_2355, but six assimilation cycles are employed. The last experiment in this group (FCYC2_2254) extends the short forecast length from 10 min in the above experiments to 30 min. Consequently, only two assimilation cycles are run within the period from 2158 to 2254 UTC. The last group of experiments evaluates the relative importance of each radar by excluding data from each radar, respectively. Each of the NOKVNX, NOKDDC, NOKICT, and NOKTLX experiments withdraws the KVNX, KDDC, KICT, and KTLX radars, respectively. The YSKVNX experiment assimilates only the data from KVNX.
For all the experiments, 5-h forecasts are performed after the data assimilation. Longer-range forecasts are not conducted because of the consideration of the boundary influences. Because the main purpose of this study is to investigate the impact of assimilating radar observations on QPF, we evaluate the results of the experiments by computing the rainwater correlation coefficient between the model forecast rainwater and that derived from the observed reflectivity, although in the first three experiments we also compute the threat score of 5-h-accumulated rainfall and radial velocity root-mean-square (RMS) difference between the model-produced and the observed radial velocity. The radar data used to compute the rainwater correlation and the radial velocity RMS error during the free forecast period were subject to the same quality control and preprocessing procedures (see section 3c) as the data used during the assimilation period. Because the preprocessed data are on radar conical surfaces, the same observation operator that is used in the cost function to interpolate the model data from the model grid points to the verification data grid points is applied in order to compute the rainwater correlation and the radial velocity RMS error. The model produces a forecast output every 15 min; hence, the verification quantities are computed with a 15-min interval centered at each model output time. At any verification time t, all radar volumes (including all tilts) that fall into the time window (t − 7.5 min, t + 7.5 min) are used for the computation of the rainwater correlation and the radial velocity RMS error. The time differences between the radar volumes are disregarded.
a. Experiments without 4DVAR radar data assimilation
The experiments without 4DVAR radar data assimilation (NO4DV1 and NO4DV2) were conducted to serve as benchmarks for the rest of the experiments. The initial fields at 2254 UTC for these two experiments are from the mesoscale analysis that was obtained following the procedures described in section 4. Although these two experiments do not assimilate radar volumetric data, they do use radar-derived VAD profiles in the analysis. Figure 5a shows the results from the mesoscale analysis. The grayscale depicts the vertical velocity field at the height of 2.25 km. The dashed contours show the temperature deviation from the mean profile at the lowest model level (0.25 km). The wind vectors are also from the lowest model level. The vertical velocity field clearly shows that the maximum vertical motion is located nearby the north Oklahoma border and west of the KVNX radar. The location of this vertical velocity maximum corresponds to the intersection of the dryline and the outflow boundary (see Fig. 1a). The dryline and the outflow boundary are shown as broad ascending zones in the vertical velocity field. The air is colder in the northeast area of the domain, which is the result of the cold air outflow from previous convections.
Figure 6 shows the 3- and 5-h forecasts starting from the initial conditions of the mesoscale analysis (Fig. 5a). Because there is no radar data assimilation, the hydrometeor fields (cloud water and rainwater) are initially zero. By 3 h into the forecast (0154 UTC), convection begins to develop in the center of the forecast domain near the Oklahoma border, as shown in Fig. 6a by the forecasted reflectivity [shading; converted from the model output of the rainwater mixing ratio using (3.3)]. The storm initiated in NO4DV1 seems to correspond to storm A (see Fig. 2a), observed at 2254 UTC. Storm B, which also plays an important role in developing the squall line, is completely missed in the model forecast. On the other hand, the observed reflectivity (plotted by contours in Fig. 6a) shows that a squall line has formed in north Oklahoma at 0154 UTC. The model-forecasted storm expands in the west–east direction in the next 2 h, while the observed squall-line system moves to the southeast corner of the domain with a southwest-to-northeast orientation.
The other experiment (NO4DV2) without 4DVAR radar data assimilation is similar to NO4DV1 except that the observed rainwater mixing ratio field is inserted at the beginning (2158, 2220, and 2242 UTC; see Fig. 2) of each assimilation cycle. Although it was anticipated that the rain would fall to the ground without the dynamical adjustment to maintain the inserted rainwater, this experiment was run to examine whether such a process will produce a negative impact on the forecasts of precipitation and wind fields. Figures 7a,b,c show 1-, 3-, and 5-h forecasts from NO4DV2, respectively. The rainwater inserted at the initial time falls out of the domain completely before the first forecast hour, which results in a cold-air outflow centered at northwest of the radar KVNX near the Oklahoma and Kansas border (Fig. 7a). The wind disturbance south of the border associated with the outflow and the southerly prestorm environmental wind formed a convergence boundary that is responsible for the subsequent initiation of convection (Fig. 7b). However, comparing Fig. 7c with Fig. 6b, it is seen that the location and orientation of the convective system at the fifth hour is improved over that in NO4DV1.
To quantitatively compare the forecasts, rainwater correlation between forecast and observation, domain maximum vertical velocity, and radial velocity RMS error from NO4DV1, NO4DV2, and CYC3_2254 (will be described in the next subsection) are plotted in Fig. 8 with respect to forecast time. Note that the computation method of the rainwater correlation and radial velocity RMS error was described in section 4. The results from the experiment that has the 4DVAR radar data assimilation are also shown (solid line). The improvement resulting from the initialization of the rainwater field in NO4DV2 is shown by the slight increase of the rainwater correlation over NO4DV1 throughout the 5-h forecast. The low correlation at the initial time (2254 UTC; forecast time 0) of the forecast period indicates that the rain has mostly fallen to the ground 12 min after the last insertion of the rainwater at 2242 UTC. In Table 1, the averaged rainwater correlation coefficient over the 5-h forecast period is shown for each experiment. NO4DV2 has an averaged value of 0.034, while NO4DV1 has a value of 0.008.
The threat score of the 5-h-accumulated precipitation is plotted in Fig. 8b with respect to the precipitation amount. The verifying precipitation is estimated from radar reflectivity observations based on Eq. (3.3). Figure 8b shows that the threat score for the low precipitation amount (2 mm) is only 0.1 when radar data are not used. The scores increase significantly in the two experiments that include radar data in particular for the CYC3_2254 experiment that assimilates both radial velocity and reflectivity. The threat score in all of the three experiments decreases as the precipitation amount increases, suggesting that the high precipitation is less predictable and the traditional verification technique is unable to evaluate the forecast skills of the convective storms.
The domain maximum velocity is plotted in Fig. 8c to compare the model spinup. It is clearly shown that the continuous insertion of radar-derived rainwater in NO4DV2 reduces the model spinup period by about 50 min, compared with NO4DV1. These results indicate that the inclusion of the radar-derived rainwater can still have a small, positive impact on QPF although there is no dynamical adjustment to match the inserted rainwater field during the analysis. To examine the response of the wind forecast to the insertion of the rainwater field, the radial velocity RMS error during the forecast period is shown in Fig. 8d. It reveals that the NO4DV2 experiment has a greater wind error during the forecast period between 85 and 240 min. This finding of positive impact on QPF and negative impact on wind is consistent with that from a recent study on nudging radar-derived rainwater to MM5 model (M. Xu 2007, personal communication).
b. Experiments with 4DVAR radar data assimilation
Two experiments that assimilate radar full-volume observations are conducted at different stages of the squall line. The CYC3_2254 experiment has a cold start at 2158 UTC and a free forecast starting at 2254 UTC. The CYC3_2355 experiment has a cold start at 2259 UTC and a free forecast starting at 2355 UTC. The purpose for conducting the two experiments at two different times of the convection development is to examine the robustness and consistency of the data assimilation system.
The analyzed wind vector and perturbation temperature fields (line contours) at the height of 0.25 km at 2254 UTC are shown in Fig. 5b. The perturbation temperature is computed by subtracting the horizontal mean from the total temperature. The vertical velocity field at the height of 2.25 km is also shown in Fig. 5b by the shaded contours. In comparison with Fig. 5a, it is clear that the analysis captures convective-scale variations both in the wind and temperature fields. Two cold air maxima, with values lower than −8° and −5°C corresponding to storms A and B, respectively, are analyzed. The maximum vertical velocity in Fig. 5b reaches a value of near 9 m s−1 on the plotted level and is dramatically increased from the analysis without radar data (Fig. 5a). The domain maximum vertical velocity from CYC3_2254 for the entire forecast period is shown in Fig. 8c by the solid line. In comparison with the curves for NO4DV1 and NO4DV2, it is evident that the strong updraft resulting from the resolved convective-scale balance significantly reduces model spinup as required in the experiments without radar data assimilation.
To examine the vertical distribution of the analysis from CYC3_2254, the cross section along the north–south line indicated in Fig. 5b is shown in Figs. 9a–c for vertical velocity, perturbation temperature, and relative humidity, respectively. The model-generated reflectivity field is also shown in all the three panels by shaded contours. The analysis shows an updraft centered at the midlevel with a maximum value exceeding 15 m s−1 and a downdraft in the lower levels with a maximum value greater than 5 m s−1. The air in the updraft region is saturated or near saturated. The low-level cold pool, produced by evaporative cooling of rainwater, has a depth of about 3 km. The latent heating of condensation of saturated air plays an important role in producing the positive temperature anomaly above 3 km.
Next, we examine the evolution of the forecasted squall line and its agreement with observations. Figure 10 shows the model’s first level of reflectivity (shaded contour) from CYC3_2254 and observed reflectivity (line contour) at 0, 1, 3, and 5 h into the forecast. The forecasted squall line overlaps with the observed one, resulting in significant increase of the rainwater correlation coefficient (Fig. 8a). The forecasted wind field from CYC3_2254 also improves over these two experiments, as shown in Fig. 8d. To further examine the forecasts, the perturbation temperature at 0, 1, 3, and 5 h at the height of 0.25 km are plotted in Fig. 11, with the vertical velocity contours overlaid. The cold pool associated with thunderstorm A located at the Oklahoma–Kansas border expands with time and merges with the earlier cold pool in the northeast of the domain. To qualitatively verify the cold pool analysis and forecasts, surface observations from a composite surface network are shown in Fig. 12. The temperature contours of 22° and 26°C, corresponding to −8° and −4°C of the perturbation temperature from VDRAS, are drawn in the figure. Comparing Fig. 12 with Fig. 11, we found that the surface network is insufficient to capture the strong cold pool shown in the VDRAS analysis at 2254 UTC because of insufficient spatial resolution. As the cold pool expands in later times, however, it is well observed by the surface network. The general patterns and locations of the cold pool agree quite well between the forecast and the observation at all forecast times shown in Fig. 11. Both the observation and the forecast show that the west portion of the cold pool advances faster prior to 0054 UTC. Nevertheless, the observation at 0154 UTC indicates that the faster advancement is shifted to the east while the model does not forecast this shifting. This may explain the overprediction of the convection in west of Oklahoma at 0354 UTC (Fig. 10d).
The CYC3_2355 experiment is the same as CYC3_2254, but started 1 h and 1 min later. Figure 13 shows the model’s first level of reflectivity (shaded contour) from CYC3_2355 and observed reflectivity (line contour) at 0, 2, 3, and 4 h into the forecast. The forecasts well capture the location and the extent of the major convective band. However, spurious convections are seen at all the three forecast times plotted in Fig. 13, especially west of Oklahoma. Close examination revealed that they were initiated by the disturbance created by the falling of precipitation particles in the convective cell near the southwest corner of the domain. This convective cell is located on the dryline that is not represented very well by the VDRAS analysis; a possible cause is the small domain used in this study. As mentioned before, the choice of the domain size is limited by the available computation resource. A new version of VDRAS that has parallelization capability is being tested, and a study on the sensitivity of storm forecast with respect to domain size and boundary conditions are underway. [The rainwater correlation plotted in Fig. 15b (solid line) clearly shows that the high correlation in the analysis (forecast time 0 min) drops rapidly because of the falling of some precipitating particles, suggesting that the analysis overfits to the rainwater observations and some initial adjustment occurred during the forecast.]
c. Sensitivity to cycling configuration
Three experiments are designed to examine the sensitivity of the analysis and forecast with respect to VDRAS cycling configuration. The configurations of the three experiments—CYC1_2254, CYC6_2355, and FCYC2_2254—along with the two control experiments—CYC3_2254 and CYC3_2355—are illustrated in Fig. 14. The sensitivity of these experiments is evaluated by the rainwater correlation between the forecast and the observation. The rainwater correlations from CYC1_2254 and FCYC2_2254 are compared with that from CYC3_2254 in Fig. 15a. It is shown that three continuous cycles lower the correlation at the analysis time but increase the correlation almost at every output forecast times. The experiment with increased short forecast length (FCYC2_2254) shows a decreased skill initially but comparable skill after 100 min into the forecast. This is an encouraging result as far as real-time applications are concerned. For real-time applications, we are often forced to use a longer forecast period to keep up with the wall-clock time because the forecast takes less time than 4DVAR analysis for the same forecast length.
To examine the impact of using more than three assimilation cycles, the rainwater correlation from CYC6_2355 is compared with that from CYC3_2355 in Fig. 15b. Results of the rainwater correlations indicate that the experiment with six cycles has a lower skill in the first 3 h but slightly higher skill in the last 2 h. The model’s first level of reflectivity (shaded contour) from CYC6_2355 and observed reflectivity (line contour) at 0, 2, 3, and 4 h into the forecast are shown in Fig. 16. Comparing the forecasts with those from CYC3_2355 (Fig. 13), we can see that the area coverage of the forecast convective system from CYC6_2355 is reduced from that in CYC3_2355. The results also show that the number and strength of spurious convections are reduced, which is a benefit resulting from a longer assimilation period. Another point worth mentioning from the above sensitivity experiments is that the 4DVAR scheme is quite robust with regard to the change of the cycling configuration, which is another way to confirm the capability of the 4DVAR technique.
d. Data impact from each radar
The final group of sensitivity experiments is designed to test the impact of each radar to the analysis and subsequent forecast. These experiments are similar to CYC3_2254, except that radar data from KICT, KDDC, KVNX, and KTLX are withdrawn (including the VAD profile from that radar in the mesoscale analysis) for the NOKICT, NOKDDC, NOKVNX, and NOKTLX experiments, respectively. The YSKVNX experiment uses data from KVNX, excluding data from all of the other three radars. The purpose of this set of experiments is to examine the sensitivity of the analysis and the prediction of the convective system with respect to the amount and location of the data. The last experiment is to show how much the impact of a single radar, located in the center of the convective system, can have on the analysis and forecast. Figure 17 shows the horizontal wind vector difference between the control experiment CYC3_2254 and each of the five experiments—NOKICT (Fig. 17a), NOKDDC (Fig. 17b), NOKVNX (Fig. 17c), NOKTLX (Fig. 17d), and YSKVNX (Fig. 17e)—at the lowest model level z = 0.25 km at the analysis time of 2254 UTC. The horizontal wind vector at the same level from CYC3_2254 is shown in Fig. 17f for reference. The reflectivity (shading) at z = 0.25 km and the vertical velocity (contour) at z = 2.25 km are also shown in all the six panels. As expected, the greater wind differences generally occur around the radar that is excluded and exceed 10 m s−1 in magnitude (see Figs. 17 a,b,c). The only exception is for NOKTLX. For this experiment, the difference wind vector is much smaller in magnitude and is more uniformly distributed. (Note that the scale of the wind vector in Fig. 17d is nearly 4 times smaller than that in the other panels). The wind difference in YSKVNX is more widespread and reaches a maximum of 19.4 m s−1 west of KVNX around storm A. The vertical velocity patterns are not significantly different from the CYC3_2254 control experiment, except for YSKVNX. The band pattern is broken into a cellular pattern in YSKVNX. When the vertical distribution of the analyzed vertical velocity is plotted (Fig. 18a), it is found that the NOKICT, NOKDDC, and NOKTLX experiments result in greater velocity magnitude in the mid- and high levels, while the other two experiments significantly underestimate the velocity. However, all of these experiments produced precipitation forecasts with lower skills than the control experiment, CYC3_2254 (see Fig. 18b and Table 1), as measured by the rainwater correlation coefficient. Not surprisingly, the YSKVNX experiment results in the lowest skill. It is interesting to note from Fig. 18b that the rainwater correlation coefficient from this experiment is highest at the analysis time, but drops rapidly to near zero by 100 min into the forecast because of the lack of the correct wind analysis around the convective system. In addition, YSKVNX and NOKDDC have the lowest overall skill score, suggesting that the inflow in northwest of the storm that is underestimated (reflected by the large vector difference in Fig. 17b in the northwest of the domain) has a significant impact on the forecast. When KVNX is excluded from the data assimilation, a much weaker vertical velocity is obtained due to the fact that KVNX is located closest to the convective cells and, as a result, observes the detailed structure of these cells. However, the negative effect associated with the exclusion of KVNX data is rather temporary. In Fig. 18b, it is clearly shown that rainwater correlation for NOKVNX drops sharply in the first 39 min or so, but rebounds to the level comparable with the control experiment, indicating that the lack of the southerly wind immediately south of the convective cells (as shown in Fig. 17c) only has a temporal effect on the subsequent forecasts. These experiments suggest that accurate analysis of the mesoscale flow structure plays a more important role than the analysis of the detailed structure of the initial convective cells in predicting the development of the convective system.
These five experiments were repeated, but with the VAD profiles used in the background analysis. Similar behaviors were observed although the forecast skills for all of the experiments were improved over the corresponding experiments presented above that not only exclude the volumetric data from the specific radar but also the VAD profile from that radar.
6. Summary and conclusions
In this paper, we present a case study to investigate the feasibility of using observations from multiple Doppler radars for initialization of a convective system and the impact on its forecast. The squall-line case of 12–13 June 2002 observed during IHOP_2002 was used in the study. A 4DVAR system Doppler radar assimilation system was employed for the initialization of a cloud model. The system applies a two-step procedure in which the nonradar observations are first analyzed to produce a mesoscale background for the subsequent 4DVAR radar data assimilation. A cycling procedure with a 12-min assimilation window, followed by a 10-min forecast window, is applied to take advantage of the frequent radar observations. The impact of the initialization on the subsequent short-term forecast was examined by comparing experiments both with and without radar observations. The results were evaluated by computing the rainwater correlation coefficient between the forecast and observation. A set of experiments was performed to test the sensitivity of the analysis and forecast with respect to data assimilation cycling configuration. Another set of experiments was conducted to examine the relative importance of observations from each radar.
The major conclusions from this case study is summarized below:
The assimilation of multiple Doppler radar radial velocity and reflectivity significantly reduces the model spinup and results in much improved forecasts for at least 5 h.
The initialization of rainwater mixing ratio by inserting radar reflectivity to the initial conditions results in slightly improved precipitation forecast compared with an experiment without rainwater initialization, but the process degrades the forecast accuracy of the wind.
Forecasts from the analysis that used three continuous assimilation cycles are improved over those from the analysis that employs only one cold-start cycle, but no further benefit is found with more cycles.
A longer forecast period of 30 min in between two 4DVAR cycles results in forecasts that are only slightly degraded in comparison with that of 10 min.
Accurate analysis of the mesoscale flow structure from radar observations plays a more important role than the analysis of the detailed structure of the existing convective cells in predicting the squall-line development.
The current study, though from only one case, has shown the potential in using multiple Doppler radar observations for initializing a convection-permitting model. Needless to say, more case studies and real-time demonstrations are needed to confirm the value of Doppler radar data on QPF. One important question that is not answered through this study is the sensitivity of the convective-scale analysis and prediction with respect to the accuracy of the mesoscale and large-scale background analysis. This topic will be studied in our future work. The model used in the current study is a simple cloud model; a more realistic model, such as WRF, will be used in the future to evaluate the applicability of the 4DVAR methodology in full-blown mesoscale models.
The authors thank Jay Miller and Sherrie Frederic for their help in preprocessing and performing quality control of the radar data. This research was supported by the U.S. Weather Research Program (USWRP). We appreciate reviews of the manuscript by Qingnong Xiao and Morris Weisman.
* The National Center for Atmospheric Research is sponsored by the National Science Foundation
Corresponding author address: Dr. Juanzhen Sun, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307. Email: email@example.com