Advances in Meteorology

Volume 2013 (2013), Article ID 792631, 16 pages

http://dx.doi.org/10.1155/2013/792631

Research Article

Retrieving 3D Wind Field from Phased Array Radar Rapid Scans

Xiaobin Qiu,1,2 Qin Xu,³ Chongjian Qiu,² Kang Nai,⁴ and Pengfei Zhang3,4

¹Tianjin Institute of Meteorological Science, Tianjin 300074, China

²College of Atmospheric Sciences, Lanzhou University, Lanzhou 730000, China

³NOAA/OAR National Severe Storms Laboratory, David L. Boren Boulevard, Norman, OK 73072, USA

⁴Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, Norman, OK 73072, USA

Received 11 May 2013; Revised 1 July 2013; Accepted 19 July 2013

Academic Editor: Jidong Gao

Copyright © 2013 Xiaobin Qiu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The previous two-dimensional simple adjoint method for retrieving horizontal wind field from a time sequence of single-Doppler scans of reflectivity and/or radial velocity is further developed into a new method to retrieve both horizontal and vertical winds at high temporal and spatial resolutions. This new method performs two steps. First, the horizontal wind field is retrieved on the conical surface at each tilt (elevation angle) of radar scan. Second, the vertical velocity field is retrieved in a vertical cross-section along the radar beam with the horizontal velocity given from the first step. The method is applied to phased array radar (PAR) rapid scans of the storm winds and reflectivity in a strong microburst event and is shown to be able to retrieve the three-dimensional wind field around a targeted downdraft within the storm that subsequently produced a damaging microburst. The method is computationally very efficient and can be used for real-time applications with PAR rapid scans.

1. Introduction

Updrafts and downdrafts are the essential components of storms. Their strengths often determine the type and evolution stage of storms. Quickly detecting updrafts and downdrafts and estimating their strengths in storm wind fields will make timely and accurate assessments of hazardous weather conditions. It is thus desirable to develop an efficient method to retrieve both the horizontal and vertical winds, including updrafts and downdrafts, in real time from phased array radar (PAR) rapid scans of storms. A key advantage of PAR over Weather Surveillance Radar 1988-Doppler (WSR-88D) is the capability to rapidly and adaptively scan storms. With its agile electronic beam steering, the PAR scan strategy can be optimized on particular weather phenomena with the volume scan time reduced from minutes to seconds (Zrnic et al. [1], Torres et al. [2]). High spatial and temporal resolution volumetric radar data are often necessary to resolve very fine echo structures, their transient developments, and movements inside storms (Heinselman et al. [3]). Previous research also indicates that the retrieval errors can be reduced if the reflectivity and radial-velocity fields are sampled more frequently (Qiu and Xu [4], Shapiro et al. [5]).

Since Doppler radar observations are limited mainly to reflectivity and radial-component velocity (along the radar beam) and there is no direct measurement of the remaining two wind components perpendicular to the radar beam, a two-dimensional simple adjoint (2D-SA) method was developed by Qiu and Xu [6] to retrieve the horizontal wind field from radar scans at low-elevation angles. In this 2D-SA method, a simplified reflectivity advection equation which is used to predict the reflectivity and the time-mean velocity that advects the reflectivity field in this equation is estimated by minimizing the difference between the radar observed and predicted reflectivity fields. The method was then refined and successfully tested with many real radar observations (Xu et al. [7, 8]). Built on the above success, a three-dimensional simple adjoint (3D-SA) method was developed by Xu et al. [9]. By using the full momentum equations and the mass continuity equation as weak constraints, this 3D-SA method can retrieve the three-dimensional wind field and the perturbation potential temperature field, similarly to the four-dimensional adjoint (4DA) method with a complete system of dynamic and thermodynamic equations (Sun and Crook [10]). Computationally, this 3D-SA method is more efficient than the 4DA method but still too expensive to apply to real-time observations from PAR rapid scans. As the control variable dimension in the 3D-SA method is much larger than that in the 2D-SA method, the 3D-SA method is not only computationally more expensive but also less flexible to adapt to incomplete data coverage than the 2D-SA method. In view of the above limitations, the 3D-SA method has not been applied to PAR data. Instead, the 2D-SA method can be further developed into a two-step SA method in which the vertical velocity is retrieved in a selected vertical cross-section along the radar beam in the second step after the horizontal velocity is retrieved in the first step on the conical surface at each tilt of PAR scans in a targeted domain. The basic idea and formulations of this two-step SA method are described in the next section. The method is applied to PAR observations for a selected case in Section 3. The benefits of rapid scans and the usefulness of mesoscale background wind field are examined in Section 4. Conclusions follow in Section 5.

2. Description of the Method

2.1. Basic Idea

To reduce the computational cost, the horizontal and vertical winds will be retrieved separately in two steps. In the first step, the 2D-SA method is used to retrieve the horizontal winds on each conical surface of the radar scans in a targeted domain of convective scale, while the mesoscale background horizontal wind field is provided by the existing radar wind analysis system that was developed based on the statistic interpolation for real-time applications with the operational WSR-88D radars (Xu et al. [11]). The vertical velocity is then retrieved in the second step in the along-beam vertical cross-section that cuts through the concerned feature at the center of the targeted domain. Since the horizontal winds are retrieved in the first step and their related terms are known in the forecast equation in the second step, the control variable dimension is reduced in the second step.

2.2. Horizontal Wind Retrieval in the First Step

The 2D-SA method is used in the first step to estimate the incremental time-mean quasi-horizontal velocity with respect to the background time-mean quasi-horizontal velocity , the time-mean source term (that includes the effect of vertical advection), and the horizontal turbulent diffusivity coefficient in the equation of quasi-horizontal advection of reflectivity over the time period of consecutive volume scans. Here, the quasi-horizontal velocity is defined as the nearly horizontal component of the three-dimensional vector velocity projected onto the conical surface of radar scan, and the quasi-horizontal advection is the advection produced by the quasi-horizontal velocity on the conical surface of radar scan. In this first step, are estimated by minimizing the following cost function:

The first term in (1) measures the difference between the predicted reflectivity and the observed reflectivity , and this term is given by

where is the horizontal area on the conical surface of radar scan in the targeted retrieval domain, is the weight, is the time period covering the sequential volume scans, and is the time elapsed for each volume scan. Here, is predicted by the following quasi-horizontal advection equation:

where is the estimated time-mean quasi-horizontal velocity (including the projection of hydrometeors’ terminal velocity onto the radar beam).

The second term in (1) measures the difference between the estimated time-mean radial velocity and observed radial velocity . This term is given by

where is the weight, is the radar observed radial velocity, is the along-beam radial component of , and is the azimuthal angle of the radar beam (positive for clockwise rotation from the -coordinate pointing to the north) at the observation point in (on the conical surface).

The third and fourth terms in (1) are given by

respectively, where and . These terms impose weak constraints on the divergence and vorticity of the incremental time-mean velocity to suppress spurious divergence and vorticity caused by data noises in the same way as in Xu et al. [7]. Their associated differential operators enhance the background error correlations, according to Xu [12], in addition to the Gaussian correlations used for in (7).

The last term in (1) is the background term for and . This term is given by

where , , and are the weights (given by the inverses of the background error variances associated with , , and , resp.), are the control variables related to by , and are the stream function and velocity potential defined by and , are the background error correlation functions with given decorrelation lengths for , and denotes the spatial convolution between the two functions (on the two sides of ) over . Numerically, these convolutions are computed by a recursive filter (Purser et al. [13, 14]). Their discrete formulations (in matrix forms) are similar to those in (8a), (8b), (8c), and (8d) of Xu et al. [15]. The time integrations for in (2) and in (4) are computed by summing their respective integrands at each time step over the time period , where is the time step used for the numerical integration of (3a) with (3b)-(3c). The spatial integrations for the cost-function terms in (2) and (4)–(7) are computed by summing their respective integrands at each grid point over the horizontal area .

The gradients of the first cost-function term in (2) with respect to the control variables ) are computed from their gradients with respect to by using the recursive filer and related transformations, while the latter gradients are given by

where is the adjoint variable obtained by integrating (backward in time) the following adjoint equation:

The derivations of (8a)–(9c) follow (2.5)–(2.13) of Qiu and Xu [6].

The gradients of the subsequent three cost-function terms in (4)–(6) with respect to the control variables are obviously zero. Their gradients with respect to the control variables are computed from their gradients with respect to by using the recursive filer and related transformations, while the latter gradients are given by

The gradient of the background term in (7) with respect to is zero. Its gradients with respect to the control variables are derived directly from (7) in the following forms:

The standard conjugate-gradient descending algorithm is used with the above computed gradients to minimize the cost-function in (1).

2.3. Vertical Wind Retrieval in the Second Step

The second step retrieves the time-mean vertical crossbeam velocity component (that includes the projection of hydrometeors’ terminal velocity but with zero background vertical velocity) in the selected vertical cross-section along the radar beam in the targeted domain. In addition to , the control variables also include a time-mean source term and the vertical turbulent diffusivity coefficient in the reflectivity advection equation. The cost function is formulated similarly to that in (1), but there is no term since the observed radial velocity has zero projection on the crossbeam velocity component . Thus, the cost function consists of only four terms, that is,

These four terms have similar forms as those in (2) and (4)–(7), but the spatial integrations in these terms are over the vertical cross-section , instead of in the retrieval domain. In addition to this difference, there are several other differences as described below.

The first cost-function term in (12) has the same form as that in (2) except that the spatial integration is over and is predicted by the following advection equation:

where is the horizontal time-mean velocity component along the azimuth of radar beam in the selected cross-section, is the vertical time-mean velocity (including the hydrometeors’ terminal velocity), and is the slope angle of the radar beam (see Figure 1). Here, is the previously defined time-mean quasi-horizontal velocity vector (projected onto the conical surface of radar scan), while is the time-mean quasi-vertical velocity component that is perpendicular to in the three-dimensional space. Since the vertical cross-section is along the radar beam, is the along-beam radial component, and hence , is the cross-beam component perpendicular to the vertical cross-section and hence is exactly horizontal. The -component is also perpendicular to the radar beam but is confined within (rather than perpendicular to) the vertical cross-section. Note that and are already retrieved in the first step while , , , and can be computed directly from the observed reflectivity field, so all their related terms, that is, , can be treated as known and thus moved to the right-hand side of (13a). This treatment reduces the control variable dimension and increases the computational efficiency significantly.

Figure 1: Sketch of geometric relationship between the components of and , plotted by the dashed black and red arrows, respectively, for the same vector plotted by the solid blue arrow. The radar beam is plotted by the gray dashed line, and is the slope angle of the radar beam.

The term has the same form as that in (5) except that the spatial integration is over , and Div is replaced by the mass continuity constraint:

where is the time-mean vertical velocity of the air motion and (≤ 0) is the hydrometeors’ terminal velocity estimated from the observed reflectivity by using the empirical formula of Kessler [16] as that in of Xu et al. [15]. The term has the same form as that in (6) except that that the spatial integration is over , and Vor is replaced by

The background term is now changed into the following form:

where and are the weights (given by the inverses of the background error variances associated with and , resp.), are the final control variables related to by , are the Gaussian correlation functions associated with the square roots of the background error covariance functions with given decorrelation lengths for , and denotes the spatial convolution between the two functions (on the two sides of ) over . Again, like the constraints in (5)-(6), the constraints in (14)-(15) can also suppress the spurious divergence and vorticity caused by data noises, and their associated differential operators also enhance the background error correlation in addition to the Gaussian correlation used above for .

The gradients of the first cost-function term with respect to the control variables in this second step are computed from their gradients with respect to by using the recursive filer as described in the first step, while the latter gradients are given by

where is the adjoint variable obtained by integrating the following adjoint equation:

The derivations of (17a)–(18c) are similar to those that lead to (8a)–(9c).

The gradients of with respect to the control variables are zero. Their gradient with respect to the control variable is computed from their gradient with respect to by using the recursive filer as described in the first step, while the latter is given by

The gradient of the background term with respect to is zero. Its gradients with respect to the control variables are derived directly from (16) in the following forms:

The standard conjugate-gradient descending algorithm is used again to minimize the cost-function in (12) in this second step. The initial guesses of are set to zero, and the algorithm converges in less than 100 iterations.

3. Application to PAR Data

3.1. PAR Observations and Mesoscale Background Wind Field

The multimission PAR, located at the National Weather Radar Testbed (NWRT), Norman, Oklahoma, is a research radar using a converted U.S. Navy SPY-1A phased array antenna. This PAR has essentially the same wavelength (9.4 cm in S band) and range resolution (250 m) as the WSR-88D radars, and it can mimic WSR-88D volume coverage patterns (VCPs) and collect data at similar pulse repetition intervals. The most significant difference between the PAR and the WSR-88D is that the phased array antenna forms each beam electronically by controlling the phases of transmit-receive elements and thus can scan storms rapidly and adaptively. With its rapid scan capability, the PAR captured the rapid evolution of a severe storm produced microburst at the 20 km radial range to the south-southwest during the early evening of July 10, 2006. The PAR applied a beam multiplexing scanning strategy to volume coverage pattern 12 (VCP12) for a 90° sector scan, so each volume scan contained up to 53 tilts (with the elevation angles from 0.51° to 41°) and is completed in just 34 seconds. Since the retrieval domain is small and not very close to the PAR, only 14 tilts (from 0.51° to 19.5°) are intercepted by the retrieval domain and will be used for the retrievals in this paper. The storm and its produced microburst were fully sampled in time and space by the PAR. In particular, the PAR detected a reflectivity core aloft between 19 : 40 : 21 and 19 : 42 : 20 UTC at about 6 km above ground level. This reflectivity core produced a strong downdraft maximized around 19 : 49 : 07 UTC as shown in Figure 2, and this downdraft generated a strong outflow near the surface sampled by the PAR from 19 : 50 : 15 to 19 : 58 : 07 UTC (Heinselman et al. [3]).

Figure 2: Reflectivity images from 20 consecutive PAR 90° sector scans of a severe microburst event from 19 : 41 : 12 to 19 : 53 : 19 UTC on July 10, 2006. Each panel consists of two subpanels: (i) the range height indicator (RHI) display along the 210.7° azimuth on the left side and (ii) the plan position indicator (PPI) display at 0.51° elevation on the right side. For the RHI display, the radial range is leftward from to 35 km and the height range is upward from to 10 km. The horizontal (or vertical) domain used in the first (or second) step of retrieval is shown by the white dashed boundary lines in the right (or left) subpanel of panel (a). The white letter C in each left subpanel marks the reflectivity core. As shown from panel (a) to (t), the reflectivity core falls from  km to 1.8 km in 667 s, so the estimated falling speed of the reflectivity core is  m s⁻¹, where (≤ 0) is the hydrometeors’ terminal velocity as defined in (14).

As mentioned in Section 2.1, the radar wind analysis system (Xu et al. [11]) is used to produce the mesoscale background horizontal wind field. For the storm case shown in Figure 2, this system is applied to radial-velocity data from the operational KTLX radar at the Oklahoma City Twin Lakes in combination with surface wind data from the Oklahoma Mesonet. The analysis domain is centered at the KTLX radar site. The domain size is 160 × 160 × 8 km³ covered by a 81 × 81 × 32 grid, the horizontal resolution is  km, and the vertical resolution is  km. The background wind fields produced at 19 : 40 : 24 UTC are shown by white arrows in Figure 3(a) at  km and Figure 3(b) at  km (above the KTLX radar site) superimposed on the KTLX reflectivity images at 0.5° and 3.3° tilts, respectively. The small red square in panel (b) of Figure 3 marks the nested domain used by the two-step SA method. As shown in Figure 3(a), the low-level background flow is divergent (or convergent) on the mesoscale to the southwest (or northeast) on the upstream (or downstream) side of the nested domain. This mesoscale divergence (or convergence) is consistent with the weakened (or enhanced) upper-level reflectivity in Figure 3(b), relative to the lower-level reflectivity in Figure 3(a), in the area to the southwest (or northeast) on the upstream (or downstream) side of the nested domain, because the weakened (or enhanced) upper-level reflectivity indicates that its associated convective-cell cluster had started to decay (or grow) and thus produced lower-level divergence (or convergence). In the nested domain, the background wind field is smooth and nearly constant at each vertical level. Clearly, the mesoscale background wind field cannot resolve the small-scale flow structures associated with the convective cell in the nested domain. As we can see from Figure 3, the lower-level and upper-level reflectivity fields have about the same intensity in the nested domain, and this suggests that the convective cell in the nested domain was fully developed around 19 : 40 : 00 UTC. The small-scale flow structure associated the downdraft produced by the convective cell in the nested domain can be retrieved by the two-step SA method as shown in the next two subsections.

(a)

(b)

Figure 3: Background wind fields (shown by white arrows) produced by the radar wind analysis system at 19 : 40 : 42 UTC at (a)  km and (b)  km superimposed on the KTLX reflectivity images at 0.5° and 3.3° tilts, respectively. The domain size is 160 × 160 × 8 km³ covered by an 81 × 81 × 32 grid centered at the KTLX radar site (marked by the blue dot). The blue + sign marks the PAR site. In panel (b), the red square box marks the nested domain used by the two-step SA method, and the yellow dashed line shows the radar beam direction from PAR along the 210.7° azimuth.

3.2. Experiment Design

The nested horizontal domain on the conical surface of each tilt of PAR scan is centered at the radial range of  km from the PAR (marked by the blue + sign in Figure 3) with the -axis along the PAR beam at the 210.7° azimuth (as shown by the yellow dashed line in Figure 3). The domain size is 16 × 16 km² (as shown by the white dashed boundary lines in the right panel of Figure 2(a)) covered by a 65 × 65 grid with a horizontal resolution of  m. The -axis of the vertical cross section is thus also along the 210.7° azimuth and centered at  km from the PAR, while the vertical domain size is 15 × 7.5 km² (as shown by the white dashed boundary lines in the left panel of Figure 2(a)) covered by a 61 × 31 grid with  m. The total of 26 volume scans (one volume every  s) for the time period from 19 : 40 : 00 to 20 : 00 : 00 UTC are used by the two-step SA method to retrieve the wind fields over 24 time windows. Each time window contains consecutive volume scans, so  min. The time step is  s for the forward integrations of (3a), (3b), and (3c) and (13a), (13b), and (13c) and the backward integrations of (9a), (9b), and (9c) and (18a), (18b), and (18c).

The weights and decorrelation length scales used by the cost function in (1) are specified as follows:

where is the time-dependent factor for the weight and is specified similarly to that in of Xu and Qiu [17], is the reflectivity observation error standard deviation and is set to 1 dBZ as in Xu et al. [15], is the error standard deviation for the background wind field, and is the error standard deviation for the (zero) background reflectivity source field.

The weights and decorrelation length scales used by the cost function in (12) are specified as follows:

where , , , and are the same as defined in . Note that the variance of the radial velocity innovation (observation minus background in the observation space), denoted by , should be the sum of the radial-velocity background error variance and radial-velocity observation error variance, denoted by . This implies that can be estimated by if (≈1 m2 s⁻²), while can be estimated by the spatially averaged RMS amplitude of the radial-velocity innovation. Similarly and loosely, is estimated by the spatially averaged RMS amplitude of the time derivative of . With these estimates, the weights and decorrelation length scales are properly tuned to the above-specified values.

3.3. Results Obtained in the First Step

Figures 4(a)–4(f) show the retrieved fields in the nested domain on the conical surfaces of six tilts (selected from the total 14 tilts from 0.51° to 19.5°) from the PAR 90° sector scans over the time period of 19 : 45 : 43~19 : 46 : 13 UTC. Note that the -coordinate is along the 210.7° azimuth (as shown by the dashed yellow line in Figure 3(b)). Thus, as shown in Figure 4, the prevailing winds in the nested domain were southwesterly on the conical surface of = 0.51° (around  km) but veered with height and gradually changed to westerly on the conical surface of = 15.6° (around  km). This feature is consistent with the vertical variation of the background winds in Figure 3. In addition to this feature, the retrieved quasi-horizontal wind field exhibits strong storm-scale variations inside and around the main reflectivity core area (with reflectivity >40 dBZ) on each conical surface. The detailed storm-scale flow structures and associated horizontal divergence/convergence patterns are examined below.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Quasi-horizontal () velocity fields (plotted by black arrows) retrieved in the nested domain on the conical surfaces of (a) = 0.51°, (b) 1.3°, (c) 2.4°, (d) 5.1°, (e) 12.5°, and (f) 15.6° from PAR 90° sector scans over the time period from 19 : 44 : 43 to 19 : 46 : 13 UTC on July 10, 2006. In each panel, the colored field shows the PAR-observed reflectivity (gridded and smoothed by the spatial interpolation on each tilt), and the black contours (with solid for positive and dashed for zero and negative) plot the field of horizontal divergence. The reflectivity of color scale is shown at the bottom of the figure, and the velocity vector scale is shown at the low-right corner of the figure. The nested domain size is 16 × 16 km² (as marked by the red square box in Figure 3(b)). The -coordinate is originated from the PAR site and directed along the 210.7° azimuth (as shown by the dashed yellow line in Figure 3(b)).

From Figures 4(a) and 4(b) we can see that before the lower-level inflow (from the left boundary of the domain) reached the reflectivity core area, the flow was not only deflected rightward but also became strongly convergent (shown by the dashed negative contours of horizontal divergence) in a narrow elongated area to the left of the reflectivity core area. Inside the reflectivity core area, the lower-level flow was strongly divergent, as shown by the solid positive contours of the horizontal divergence in the reflectivity core area in Figure 4(a). This strong lower-level divergence was caused by and tied up with a strong downdraft aloft, while this strong downdraft is not only revealed by the downward movement of the reflectivity core as exhibited by the time series of RHI display in Figures 2(a)–2(t) but also well retrieved in the second step (as shown later in Figure 6(c)).

On the other hand, as we can see from Figures 4(e) and 4(f), the upper flow was divergent in a narrow elongated area to the left of the reflectivity core and became convergent inside the reflectivity core. From Figures 4(c) and 4(d), we can see that the middle-level flow was weakly convergent or divergent, and the convergence-divergence pattern for the middle-level flow is intermediate between the two nearly opposite patterns at the lower and upper levels. Thus, at least qualitatively, the retrieved quasi-horizontal wind fields in the first step captured the storm scale convergence-divergence patterns associated with the strong downdraft generated by the storm.

Figures 5(a) and 5(b) show the reflectivity source fields (see the term on the right-hand side of the reflectivity advection equation in (3a)) obtained as by-products of the retrievals in Figures 4(a) and 4(e) on the conical surfaces of = 0.51° and 12.5°, respectively. As shown in Figure 5(a), the retrieved reflectivity source field on the 0.51° tilt is characterized by small-scale fluctuations around zero, and the RMS amplitude of these small-scale fluctuations is smaller by several times than the RMS amplitudes of the local time derivative and advection terms (not shown) in the reflectivity advection equation. Note that the vertical advection term is implicitly absorbed into the source term on the right-hand side of (3a), and this vertical advection term includes the effect of hydrometeors’ downward terminal velocity, whereas the vertical gradient of reflectivity was small and more or less positive in the lower level (around  km on the 0.51° tilt as shown by the RHI display in Figure 2(d)). Thus, the contribution of the vertical advection to is weakly positive in the lower level. This positive contribution could offset the negative contribution to produced by precipitation evaporation (since the evaporation could be also weak as implied by the small vertical gradient of reflectivity around  km). This may explain the smallness of the retrieved in Figure 5(a), although the source term also absorbs residual errors of the reflectivity advection equation in (3a) especially those caused by the imperfect reflectivity observations and imperfect velocity retrievals.

(a)

(b)

Figure 5: Reflectivity source fields (for the term on the right hand-side of (3a)) produced as by-products of the retrievals in Figures 4(a) and 4(e) on the conical surfaces of (a) = 0.51° and (b) = 12.5°. The color scale is shown on the bottom of the figure, and the unit is dBZ s⁻¹ for the labeled contour values.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 6: Time series of retrieved fields for air motions (plotted by black arrows) in the vertical cross-section (along the -coordinate in Figure 4) superimposed on the PAR-observed reflectivity fields (gridded and smoothed by the spatial interpolation in the vertical cross-section) from 19 : 42 : 03 to 19 : 55 : 43 UTC on July 10, 2006. The white arrows plot the vector velocities of hydrometeors (i.e., air velocity plus the hydrometeors’ downward terminal velocity).

In the upper level (around  km on the 12.5° tilt), as shown in Figure 5(b), the retrieved source term is maximized in the narrow curved area along and inside the left boundary of the reflectivity core (as seen by overlapping Figure 5(b) onto Figure 4(e)). In this curved upper-level area, the upward vertical velocity of air largely cancelled the downward terminal velocity for hydrometeors (as seen later from the retrieved vertical velocities in Figure 6(c)), and the vertical gradient of reflectivity was very small (on the 12.5° tilt as shown by the RHI display in Figure 2(d)). The vertical advection of reflectivity was thus very small, so only the condensation caused by the upward air motion may explain why the source term is maximized in the curved upper-level area as we have seen from Figure 5(b).

As another by-product, the horizontal turbulent diffusivity coefficient is also retrieved, and the retrieved value varies slightly with the scan elevation and time window. The mean value of the retrieved is 274 m² s⁻¹ with a standard deviation of 30 m² s⁻¹. These retrieved values of are in the same range as those obtained or used in the previous studies of the 2D-SA method (see the second paragraph in Section 2 of Xu et al. [8] and reference cited there).

3.4. Results Obtained in the Second Step

Figures 6(a)–6(h) show the time series of the retrieved fields for air motions (plotted by black arrows) in the vertical cross-section along the -coordinate in the nested domain superimposed on the PAR-observed reflectivity fields from 19 : 42 : 03 to 19 : 55 : 43 UTC on July 10, 2006. As we can see from Figures 6(a)–6(c), the air motions were upward within the reflectivity core, especially in the upper levels, and their produced condensation may largely explain the rapid intensification of the reflectivity core. However, as the reflectivity core moved down to the middle levels and finally reached the ground, the upward air motions inside the reflectivity core became weak and even changed to slightly downward motions (as shown in Figures 6(d)–6(h) and thus ceased to produce condensation during the later times. This may largely explain why the reflectivity core ceased to grow and even became slightly weak as it fell into and below the middle levels during the later time period from 19 : 48 : 50 to 19 : 55 : 43 UTC, as shown in Figures 6(d)–6(h). Furthermore, as shown by the white arrows in Figures 6(a)–6(h), the vector velocities of hydrometeors were overwhelmingly downward especially inside the reflectivity core, and their vertical-component velocities are around the value (≈6.3 m s⁻¹) as estimated by the downward movements of the reflectivity core from the time series of the RHI displays in Figure 2.

Note that the vertical advection term is explicitly considered in (13a), so the source term on the right-hand side of (13a) should be related primarily to the production/reduction of hydrometeors caused by condensation/evaporation, although the source term also absorbs the residual errors caused in (13a) by the imperfect reflectivity observations and imperfect velocity retrievals. Figure 7 shows the field produced in the vertical cross-section as a by-product of the retrieval in Figure 6(c). By overlapping Figure 7 onto Figure 6(c), it is easy to see that the red-colored positive core of  dBZ s⁻¹ in the middle levels is inside the main updraft (where the upward air motions are maximized) within the reflectivity core, so this middle-level positive core of can be largely explained by the condensation and associated production of hydrometeors in the main updraft. Similarly, the upper-level negative core of  dBZ s⁻¹ between 20 km > > 18 km along the top boundary of Figure 7 can be largely explained by the evaporation and associated reduction of hydrometeors due to the downward advection of hydrometeors into a relatively dry area (with reduced reflectivity) as shown in Figure 6(c). However, all other relatively weak positive and negative cores in Figure 7 do not seem to be physically meaningful and they may merely represent the residual errors caused in (13a) by the imperfect reflectivity observations and imperfect velocity retrievals.

Figure 7: Reflectivity source field (for the term on the right-hand side of (13a)) produced in the vertical cross-section as a by-product of the retrieval in Figure 6(c). The color scale is shown on the bottom, and the unit is dBZ s⁻¹ for the labeled contour values.

As another by-product, the vertical turbulent diffusivity coefficient is also retrieved. The retrieved values for the different time windows are almost the same. The mean value of the retrieved is 201 m² s⁻¹, and the standard deviation is merely 0.25 m² s⁻¹. Note that all values of are retrieved for the same along-beam vertical domain over a short time period, and the bulk effect of the vertical turbulent diffusion could remain nearly the same over this short time period. This may partially explain why the standard deviation is small. The retrieved values of are roughly within the same range as those of retrieved in the first step.

4. Benefits of PAR Rapid Scans and Usefulness of Mesoscale Background Wind Field

4.1. Benefits of PAR Rapid Scans

As mentioned in the introduction, the retrieval errors can be reduced if the reflectivity and radial velocity fields are scanned more rapidly than the operational WSR-88D radar scans (Qiu and Xu [4], Shapiro et al. [5]). A similar benefit is expected for the two-step SA method developed in this paper and applied to PAR rapid scans. To demonstrate this, an additional experiment is performed in this subsection with the input PAR sector scan data used every 5 min instead of every 30 s, so the temporal resolution of the input observations is reduced by 10 times and becomes about the same as that of the operational WSR-88D radar scans. The parameter settings in this experiment are the same as those described for the experiment in Section 3.2 except that is increased to 5 min and thus is increased to 10 min (since is unchanged). This experiment will be called Expt-5 min, while the experiment performed in Section 3 will be used as the benchmark.

Figures 8(a) and 8(b) show the fields retrieved from the first step of Expt-5 min on the conical surfaces of = 0.51° and 12.50°, respectively, and they are the 10 min time-mean wind fields centered at the same times (19 : 45 : 43 and 19 : 46 : 11 UTC) as the 1 min time-mean wind fields in Figures 4(a) and 4(e), respectively. As we can see, the retrieved wind field in Figure 8(a) (or Figure 8(b)) is weaker and slightly smoother than that in Figure 4(a) (or Figure 4(e)) although their gross patterns are similar. The reduced intensity and spatial variability for the retrievals in Figures 8(a) and 8(b) can be attributed largely to the reduced temporal resolution of the input observations and its caused increase of the time window (from to 10 min) for the time-mean wind retrieval. Figures 9(a) and 9(b) show the difference fields obtained by subtracting the fields in Figures 4(a) and 4(e) from those in Figures 8(a) and 8(b), respectively. The spatially averaged RMS values of the two difference fields are 3.08 and 4.25 m s⁻¹, respectively. These RMS differences can represent roughly the RMS errors caused by coarsening the temporal resolution of the input observations from to 5 min, because the retrievals from Expt-5 m are significantly less accurate than those from the benchmark experiment, as evaluated indirectly below.

(a)

(b)

Figure 8: As in Figure 4 but the retrieved fields are from the first step of Expt-5 min on the conical surfaces of (a) = 0.51° and (b) 12.50° at the same times (19 : 45 : 43 and 19 : 46 : 11 UTC) as those in Figures 4(a) and 4(e), respectively.

(a)

(b)

Figure 9: The difference fields obtained by subtracting the fields in Figures 4(a) and 4(e) from those in Figures 8(a) and 8(b), respectively.

The 1 min time-mean wind fields retrieved in the benchmark experiment are less smeared in time and therefore expected to be more accurate than the 10 min time-mean wind fields retrieved in Expt-5 min. However, since the true fields are not known, it is difficult to directly evaluate whether and how much the retrieval accuracy is improved by rapid scans. To overcome this difficulty, we resort to the square root of the first (or second) cost-function term (or ) in (1) that measures the RMS difference of the predicted (or retrieved ) to the observed reflectivity (or radial-velocity ) averaged over the retrieval time window . Note that is nondimensional and the weight in (or in ) is inversely proportional to as shown in , so the aforementioned RMS difference is measured by (or ) in the same way for different . Thus, the relative accuracies of retrievals obtained with different can be evaluated indirectly by comparing the initial values of (or ) and subsequent reductions achieved through the iterative descending of for different . To facilitate the comparison, the terms , , and computed from the first step of the benchmark experiment (or Expt-5 min) are denoted by -30 s, -30 s, and -30 s (or -5 min, -5 min and -5 min), respectively. These terms are plotted in Figure 10 as functions of iteration number by the dashed curves for the retrieval in Figure 8(a) versus those plotted by the solid curves for the retrieval in Figure 4(a). By comparing the red (or blue) solid curve with the red (or blue) dashed curve in Figure 10, we can see that the ratio of -5 min (or -5 min) to -30 s (or -30 s) is as large as 7.3 (or 6.9) at and increases to 9.3 (or 10.9) at . Similar large ratios are seen (not shown) from the retrieval in Figure 8(b) versus that in Figure 4(e). These large ratios imply that the first-step retrievals from Expt-5 m are significantly less accurate than those from the benchmark experiment.

Figure 10: -30 s, -30 s, and -30 s for the retrieval in Figure 4(a) plotted by the gray, red, and blue solid curves, respectively, as functions of iteration number . -5 min, -5 min, and -5 min for the retrieval in Figure 8(a) plotted by the gray, red, and blue dashed curves, respectively, as functions of .

Figure 11(a) shows the field for air motions (plotted by black arrows) retrieved from the second step of Expt-5 min around the same time as that (19 : 46 : 00 UTC) in Figure 6(c). As shown, the retrieved vertical-velocity field in Figure 11(a) has much stronger variations than that in Figure 6(c). In particular, the hydrometeors’ downward vertical velocities (shown by the white arrows) inside the reflectivity core in Figure 11(a) are much larger than the value of 6.3 m s⁻¹, estimated from the downward movements of the reflectivity core in Figure 2 (see figure’s caption for the estimated value), whereas the downward vertical velocities of hydrometeors inside the reflectivity core in Figure 6(c) are very close to and around the estimated value of 6.3 m s⁻¹. Thus, the retrieval in Figure 11(a) is significantly less accurate than that from the benchmark experiment in Figure 6(c). Figure 11(b) shows the difference field obtained by subtracting the field in Figure 11(a) from that in Figure 6(c). The spatially averaged RMS value of this difference field is 13.31 m s⁻¹. This large RMS difference is caused mostly by the RMS error of the retrieval in Figure 11(a), and the latter is caused both directly by the reduced temporal resolution of the input observations in the second step and indirectly by the reduced accuracies in the horizontal velocities retrieved in the first step and used in the second step. This large increase of retrieval error is also reflected by the large ratio of -5 min/-30 s (≈15), where -5 min (or -30 s) is the square root of the first cost-function term in (12) computed from Expt-5 min (or benchmark experiment). Thus, the PAR rapid scans can improve not only the horizontal-wind retrieval in the first step but also the vertical-velocity retrieval in the second step, and the improvement is more significant in the second step.

(a)

(b)

Figure 11: (a) As in Figure 6(c) but the retrieved field is from the second step of Expt-5 min (around 19 : 46 : 00 UTC). (b) The difference field obtained by subtracting the field in Figure 6(c) from that in (a).

4.2. Usefulness of Mesoscale Background Wind Field

For the particular case considered in this paper, the targeted retrieval domain is small and the mesoscale background velocity is smooth and nearly constant on each tilt within this small domain. Because of this, the retrieved fields in Figure 4 change insignificantly when the background wind field is not used, while the retrieved fields in Figure 4 change even less due to the fact that the background vertical velocity is zero. These properties are verified by another additional experiment performed with no background wind, that is, and thus . This experiment will be called Expt-NB.

The terms , , and computed from the first step of Expt-NB are denoted by -30 s NB, -30 s NB, and -30 s NB, respectively. These terms are plotted in Figure 12 as functions of the iteration number by the dashed curves versus those plotted by the solid curves for the retrieval in Figure 4(a). By comparing the red (or blue) solid curve with the red (or blue) dashed curve in Figure 12, we can see that the ratio of -30 s NB (or -30 s NB) to -30 s (or -30 s) is 1.5 (or 2.2) at , decreases rapidly to 1.09 (or 1.03) at , and then gradually approaches 1.000 (or 1.000) as further increases toward 100. These ratios are large initially, and this is because the unknown control variables are zero initially in Expt-NB and thus have larger errors than those in the benchmark experiment where are zero initially. Thus, using the mesoscale background velocity field can reduce the initial errors of the control variables . The mesoscale background velocity should be more useful than just reducing the initial errors of the control variables for the two-step SA method if the retrieval domain is not too small (as that in this paper) to resolve some mesoscale variations from the background field. This speculation needs to be verified in future applications of the two-step SA method.

Figure 12: As in Figure 10 but the gray, red, and blue solid curves plot -30 s NB, -30 s NB, and -30 s NB, respectively, for the retrieval from the first step of Expt-NB (instead of Expt-5 min).

5. Conclusions

In this paper, the simple adjoint method (Qiu and Xu [6], Xu et al. [7, 8, 17]) is further developed into a two-step method to retrieve high-resolution horizontal and vertical wind fields from the PAR rapid scans of convective storm cells in a targeted domain. The first step retrieves the horizontal vector wind field on the conical surface of radar scan at each elevation angle, and the second step retrieves the vertical velocity in an along-beam vertical cross-section. As the horizontal winds can be retrieved in parallel on different elevations in the first step and only the vertical velocity needs to be retrieved in the second step, this two-step method is computationally efficient. In particular, on a workstation with Intel(R) Xeon(R) CPU X7550 (2.00 GHz), the computer time for the first step is about 20 s (with the retrievals performed in parallel for different tilts) and the computer time for the second step is merely 10 s. Thus, the method can be applied very efficiently to real-time PAR rapid scans. The performance and expected capability of the method are demonstrated by a real-data example in Section 3 where the method is applied to PAR rapid 90° sector scans of a severe storm that produced a strong downdraft and subsequent damaging microburst during the early evening of July 10, 2006.

The above severe storm was scanned not only by the PAR but also by the operational KTLX radar (about every 5 minutes per volume) and the terminal Doppler weather radar (also about every 5 minutes per volume). These two operational radars, however, are all located in the same northeast quadrant as the PAR relative to the storm, so real dual-Doppler observations are not available for quantitative verifications of the single-Doppler retrievals obtained in this paper. Nevertheless, the method used in the first step is a refinement of the previous 2D-SA method, and the previous method has been tested for many real cases with the retrieved wind fields well verified by dual-Doppler analyses (Xu et al. [7, 8]). Thus, the method in the first step should remain at least as reliable as the previous 2D-SA method, and this speculation is supported by the detailed analysis of the retrievals obtained in the first step (see Section 3.3). The method developed for the second step is new. This new method performs well with the PAR rapid scans as suggested by the detailed analysis of the retrievals obtained in the second step (see Section 3.4), but the performance deteriorates significantly when the temporal resolution of the input observations is coarsened from 0.5 to 5 min to mimic the operational WSR-88D radar scans (see Section 4.1).

Note that the 2D-SA can be also formulated in the polar coordinates centered at the radar (instead of the local Cartesian coordinates as shown in Figure 2(a)) with the retrieval performed in a targeted sector area on each tilt of radar scan essentially in the same way as that in Xu et al. [18] and Xu and Qiu [17]. After these sector-area retrievals are performed in parallel for different tilts in the first step, the second step can be performed efficiently in parallel for different along-beam vertical cross-sections. In this way, the method can retrieve high-resolution three-dimensional wind fields in real time from PAR rapid scans over various adaptively-nested sector areas threatened by damaging winds generated by severe storms. This real-time capability will be developed and tested in future studies.

Acknowledgments

The authors are thankful to Dr. Alan Shapiro and anonymous reviewers for their comments and suggestions that improved the representation of the results and to GuangYu Wang for his help in editing Figure 4. The research was supported by ONR Grant N000141010778 to the University of Oklahoma (OU), NOAA-OU Cooperative Agreement NA17RJ1227, and Tianjin Municipal Meteorological Bureau Grant 201316.

References

D. S. Zrnic, J. F. Kimpel, D. E. Forsyth et al., “Agile-beam phased array radar for weather observations,” Bulletin of the American Meteorological Society, vol. 88, no. 11, pp. 1753–1766, 2007.

S. Torres, P. Heinselman, R. Adams et al., “ADAPTS implementation: can we exploit phased-array radar's electronic beam steering capabilities to reduce update times?” in Proceedings of the 28th Conference on Interactive Information and Processing Systems (IIPS) for Meteorology, Oceanography, and Hydrology, American Meteorological Society, New Orleans, La, USA, 2012.

P. L. Heinselman, D. L. Priegnitz, K. L. Manross, T. M. Smith, and R. W. Adams, “Rapid sampling of severe storms by the national weather radar testbed phased array radar,” Weather and Forecasting, vol. 23, no. 5, pp. 808–824, 2008.

C. Qiu and Q. Xu, “Least squares retrieval of microburst winds from single-Doppler radar data,” Monthly Weather Review, vol. 124, no. 6, pp. 1132–1144, 1996.

A. Shapiro, P. Robinson, J. Wurman, and J. Gao, “Single-Doppler velocity retrieval with rapid-scan radar,” Journal of Atmospheric and Oceanic Technology, vol. 20, no. 12, pp. 1758–1775, 2003.

C. Qiu and Q. Xu, “A simple adjoint method of wind analysis for single-Doppler data,” Journal of Atmospheric and Oceanic Technology, vol. 9, no. 5, pp. 588–598, 1992.

Q. Xu, C. Qiu, H. Gu, and J. Yu, “Simple adjoint retrievals of microburst winds from single-Doppler radar data,” Monthly Weather Review, vol. 123, no. 6, pp. 1822–1833, 1995.

Q. Xu, H. Gu, and C. Qiu, “Simple adjoint retrievals of wet-microburst winds and gust-front winds from single-Doppler radar data,” Journal of Applied Meteorology, vol. 40, no. 8, pp. 1485–1499, 2001.

Q. Xu, H. Gu, and S. Yang, “Simple adjoint method for three-dimensional wind retrievals from single-Doppler radar,” Quarterly Journal of the Royal Meteorological Society, vol. 127, no. 573, pp. 1053–1067, 2001.

J. Sun and A. Crook, “Wind and thermodynamic retrieval from single-Doppler measurements of a gust front observed during Phoenix II,” Monthly Weather Review, vol. 122, no. 6, pp. 1075–1091, 1994.

Q. Xu, K. Nai, L. Wei, P. Zhang, Q. Zhao, and P. R. Harasti, “A real-time radar wind data quality control and analysis system for nowcast application,” in Proceedings of the International Symposium on Nowcasting and Very Short Range Forecasting (WSN '09), WMO, Whistler, Canada, 2009, http://www.nowcasting2009.ca/images/stories/Presentations/WSN09_Presentations.htm.

Q. Xu, “Representations of inverse covariances by differential operators,” Advances in Atmospheric Sciences, vol. 22, no. 2, pp. 181–198, 2005.

R. J. Purser, W. S. Wu, D. F. Parrish, and N. M. Roberts, “Numerical aspects of the application of recursive filters to variational statistical analysis—part I: spatially homogeneous and isotropic Gaussian covariances,” Monthly Weather Review, vol. 131, no. 8, pp. 1524–1535, 2003.

R. J. Purser, W. Wu, D. F. Parrish, and N. M. Roberts, “Numerical aspects of the application of recursive filters to variational statistical analysis—part II: spatially inhomogeneous and anisotropic general covariances,” Monthly Weather Review, vol. 131, no. 8, pp. 1536–1548, 2003.

Q. Xu, L. Wei, W. Gu, J. Gong, and Q. Zhao, “A 3.5-dimensional variational method for Doppler radar data assimilation and its application to phased-array radar observations,” Advances in Meteorology, vol. 2010, Article ID 797265, 14 pages, 2010.

E. Kessler, “On the distribution and continuity of water substance in atmospheric circulation,” Meteorological Monographs, vol. 10, no. 32, p. 84, 1969.

Q. Xu and C. Qiu, “Adjoint-method retrievals of low-altitude wind fields from single-Doppler reflectivity and radial-wind data,” Journal of Atmospheric and Oceanic Technology, vol. 12, no. 5, pp. 1111–1119, 1995.

Q. Xu, C. Qiu, and J. Yu, “Adjoint-method retrievals of low-altitude wind fields from single-Doppler reflectivity measured during Phoenix II,” Journal of Atmospheric and Oceanic Technology, vol. 11, no. 2, pp. 275–288, 1994.