Retrieving 3 D Wind Field from Phased Array Radar Rapid Scans

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 beamwith the horizontal velocity given from the first step.Themethod 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.


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 Advances in Meteorology 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.

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.

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 (Δ  , ΔV  ) with respect to the background time-mean quasi-horizontal velocity (  , V  ), 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, (Δ  , ΔV  ,   ,  ℎ ) are estimated by minimizing the following cost function: The first term in (1) measures the difference between the predicted reflectivity  and the observed reflectivity  ob , 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,  = ( − 1)Δ 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: with  (, , ) =  ob (, , ) at the boundary of  ℎ , (0, , ) =  ob (0, , ) at the initial time, where (  , V  ) = (  + Δ  , V  + ΔV  ) is the estimated timemean 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 V rm and observed radial velocity V rob .This term is given by where  vr is the weight, V rob is the radar observed radial velocity, V rm =   sin  + V  cos  is the along-beam radial component of (  , V  ), 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 Div =   Δ  +   ΔV  and Vor =   ΔV  −   Δ  .These terms impose weak constraints on the divergence and vorticity of the incremental time-mean velocity (Δ  , ΔV  ) 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 (Δ  , ΔV  ) and   .This term is given by where   ,   , and   are the weights (given by the inverses of the background error variances associated with 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  vr 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 (Δ  , ΔV  ,   ,  ℎ ) 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: with  * (, , ) = 0 at the boundary of  ℎ , and  * (, , ) = 0 at the final time.
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 (Δ  , ΔV  ) 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).

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  2 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  vr term since the observed radial velocity V rob 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: with  (, , ) =  ob (, , ) at the boundary of   , and  (0, , ) =  ob (0, , ) at the initial time, where  2 =   cos  −   sin  is the horizontal time-mean velocity component along the azimuth of radar beam in the selected (, ) cross-section,  2 =   sin  +   cos  is the vertical time-mean velocity (including the hydrometeors' terminal velocity), and  is the slope angle of the radar beam (see Figure 1).Here, (  , V  ) is the previously defined timemean 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 (  , V  ) in the three-dimensional space.Since the vertical cross-section (, ) is along the radar beam,   is the alongbeam radial component, and hence   = V rm , V  is the cross-beam component perpendicular to the vertical crosssection 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 (  , V  ) and  ℎ are already retrieved in the first step while   ,   ,  2  , and  2   can be computed directly from the observed reflectivity field, so all their related terms, that is, +  2  ), 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.
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    =  2 −   is the time-mean vertical velocity of the air motion and   (≤ 0) is the hydrometeors' terminal velocity estimated from the observed reflectivity  ob by using the empirical formula of Kessler [16] as that in (5) of Xu et al. [15].The  V 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  2 are the weights (given by the inverses of the background error variances associated with   and  2 , resp.), (  ,  2 ) are the final control variables related to (  ,  2 ) by (  *   ,  2 *  2 ) = (  ,  2 ), (  ,  2 ) are the Gaussian correlation functions associated with the square roots of the background error covariance functions with given decorrelation lengths (  ,  2 ) for (  ,  2 ), 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 (  ,  2 ,   ) in this second step are computed from their gradients with respect to (  ,  2 ,   ) 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: with  * (, , ) = 0 at the boundary of   , * (, , ) = 0 at the final time.
The gradients of   +  V with respect to the control variables ( 2 ,   ) 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 (  ,  2 ) 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 (  ,  2 ) are set to zero, and the algorithm converges in less than 100 iterations.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 3 covered by a 81 × 81 × 32 grid, the horizontal resolution is Δ = Δ = 2 km, and the vertical resolution is Δ = 0.25 km.The background wind fields produced at 19 : 40 : 24 UTC are shown by white arrows in Figure 3(a) at  = 0.25 km and Figure 3(b) at  = 5 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.

Experiment Design.
The nested horizontal domain on the conical surface of each tilt of PAR scan is centered at the radial range of  = 20 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 2 (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 Δ = Δ = 250 m.The -axis of the vertical cross section is thus also along the 210.7 ∘ azimuth and centered at  = 20 km from the PAR, while the vertical domain size is 15 × 7.5 km 2 (as shown by the white dashed boundary lines in the left panel of Figure 2(a)) covered by a 61 × 31 grid with Δ = Δ = 250 m.The total of 26 volume scans (one volume every Δ = 34∼64 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  = 3 consecutive volume scans, so  = ( − 1)Δ ≈ 1 min.The time step is Δ = 2.5 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  = 0.5[/( + Δ)] 1/2 is the time-dependent factor for the weight   and is specified similarly to that in (3.2) of Xu and Qiu [17],   is the reflectivity observation error standard deviation and is set to 1 dBZ as in Xu et al. [15],  V 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 ,   ,  V , and   are the same as defined in (22a).
Note that the variance of the radial velocity innovation (observation minus background in the observation space), denoted by  2 vri , should be the sum of the radial-velocity background error variance and radial-velocity observation error variance, denoted by  2 vro .This implies that  V can be estimated by vro (≈1 m 2 s −2 ), while  2 vri can be estimated by the spatially averaged RMS amplitude of the radial-velocity innovation.Similarly and loosely,  2   is estimated by the spatially averaged RMS amplitude of the time derivative of  ob .With these estimates, the weights and decorrelation length scales are properly tuned to the above-specified values.

Results Obtained in the First
Step.Figures 4(a)-4(f) show the retrieved (  , V  ) 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  ≈ 0.2 km) but veered with height and gradually changed to westerly on the conical surface of  = 15.6 ∘ (around  ≈ 8km).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 stormscale 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.
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.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  ≈ 0.2 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  ≈ 0.2 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.
In the upper level (around  ≈ 8km 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 .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 upperlevel 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 2 s −1 with a standard deviation of 30 m 2 s −1 .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).

Results Obtained in the Second
Step.   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 −1 ) 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  2 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  2 also absorbs the residual errors caused in (13a) by the imperfect reflectivity observations and imperfect velocity retrievals.Figure 7 shows the  2 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  2 > 0.2 dBZ s −1 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  2 can be largely explained by the condensation and associated production of hydrometeors in the main updraft.
Similarly, the upper-level negative core of  2 < −0.2 dBZ s −1 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.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 2 s −1 , and the standard deviation is merely 0.25 m 2 s −1 .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.

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  = ( − 1)Δ is increased to 10 min (since  = 3 is unchanged).This experiment will be called Expt-5 min, while the experiment performed in Section 3 will be used as the benchmark.The spatially averaged RMS values of the two difference fields are 3.08 and 4.25 m s −1 , respectively.These RMS differences can represent roughly the RMS errors caused by coarsening the temporal resolution of the input observations from Δ = 0.5 to 5 min, because the retrievals from Expt-5 m are significantly less accurate than those from the benchmark experiment, as evaluated indirectly below.
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  vr ) in (1) that measures the RMS difference of the predicted  (or retrieved V rm ) to the observed reflectivity  ob (or radialvelocity V rob ) averaged over the retrieval time window .Note that  is nondimensional and the weight   in   (or  vr in  vr ) is inversely proportional to  as shown in (22a), so the aforementioned RMS difference is measured by   ≡  1/2  (or  vr ≡  1/2 vr ) 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  vr ) and subsequent reductions achieved through the iterative descending of  for different .To facilitate the comparison, the terms  ≡  1/2 ,   , and  vr computed from  the first step of the benchmark experiment (or Expt-5 min) are denoted by -30 s,   -30 s, and  vr -30 s (or -5 min,   -5 min and  vr -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  vr -5 min) to   -30 s (or  vr -30 s) is as large as 7.3 (or 6.9) at  = 0 and increases to 9.3 (or 10.9) at  = 100.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 11(a) shows the ( 2 ,    ) 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 −1 , 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 −1 .Thus, the retrieval in Figure 11(a) is significantly less accurate than that from the benchmark experiment in Figure 6(c). 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.

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 (  , V  ) is smooth and nearly constant on each tilt within this small domain.Because of this, the retrieved (  , V  ) fields in Figure 4 change insignificantly when the background wind field is not used, while the retrieved ( 2 ,    ) 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, (  , V  ) = 0 and thus (Δ  , ΔV  ) = (  , V  ).This experiment will be called Expt-NB.The terms ,   , and  vr computed from the first step of Expt-NB are denoted by -30 s NB,   -30 s NB, and  vr -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  vr -30 s NB) to   -30 s (or  vr -30 s) is 1.5 (or 2.2) at  = 0, decreases rapidly to 1.09 (or 1.03) at  = 10, 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 (Δ  , ΔV  ) = (  , V  ) are zero initially in Expt-NB and thus have larger errors than those in the benchmark experiment where (Δ  , ΔV  ) = (  , V  )−(  , V  ) are zero initially.Thus, using the mesoscale background velocity field can reduce the initial errors of the control variables (Δ  , ΔV  ).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.

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.

Figure 1 :
Figure1: Sketch of geometric relationship between the components of ( 2 ,  2 ) and (  , V  ), plotted by the dashed black and red arrows, respectively, for the same vector k 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.

Figures 5 (Figure 4 :
Figure 4: Quasi-horizontal (  , V  ) 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 2 (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)).

Figure 5 :
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 −1 for the labeled contour values.
Figures 6(a)-6(h) show the time series of the retrieved ( 2 ,    ) 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
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 Figures6(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 Figures6(d)-6(h).Furthermore, as shown by the white arrows in Figures6(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 −1 ) as estimated by the downward movements of the reflectivity core from the time series of the RHI displays in Figure2.Note that the vertical advection term is explicitly considered in (13a), so the source term  2 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  2 also absorbs the residual errors caused in (13a) by the imperfect reflectivity observations and imperfect velocity retrievals.Figure7shows the  2 field

Figure 6 :
Figure 6: Time series of retrieved ( 2 ,    ) 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).

Figure 7 :
Figure7: Reflectivity source field (for the  2 term on the righthand side of (13a)) produced in the vertical cross-section as a byproduct of the retrieval in Figure6(c).The color scale is shown on the bottom, and the unit is dBZ s −1 for the labeled contour values.

Figures 8 (
a) and8(b)  show the (  , V  ) 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 timemean wind fields centered at the same times (19 : 45 : 43 and 19 : 46 : 11 UTC) as the 1 min time-mean wind fields in Figures4(a) and 4(e), respectively.As we can see, the retrieved wind field in Figure8(a) (or Figure8(b)) is weaker and slightly smoother than that in Figure4(a) (or Figure4(e)) although their gross patterns are similar.The reduced intensity and spatial variability for the retrievals in Figures8(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  = 1 to 10 min) for the time-mean wind retrieval.Figures9(a) and 9(b) show the difference fields obtained by subtracting the (  , V  ) fields in Figures 4(a) and 4(e) from those in Figures 8(a) and 8(b), respectively.

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

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