A Preliminary Study of Retrieving Vortex Circulation Using A Lagrangian Coordinate and Single Doppler Radar Data

A moving frame of reference technique is proposed so that the retrieval of cross-beam wind components is feasible from observations by single Doppler Radar without sacrificing data resolution. In this algorithm, the reflectivity fields observed by consecutive radar scans are used to find a moving Lagrangian coordinate on which the reflectivity measurements are as stationary as possible. After interpolating all the observational data to this optimal moving frame, and assuming the wind field is steady state within several radar scans, one can formulate a cost function which contains the following weak constraints: (1) reflectivity is preserved; (2) a geometric relationship between radial velocity Vr and the Cartesian components u, v, w ; (3) mass conservation ; and ( 4) weak vertical vorticity. By minimizing this cost function, the unobserved cross-beam wind component can be re­ covered. Applying this algorithm to a simulated, idealized Rankine vortex cir­ culation, it is found that with an appropriate scan strategy, that is, a suffi­ cient number of radar scans (>2), separated by a short enough period of time (3 minutes), low level circulation can be successfully retrieved using data detected by single Doppler radar. This technique is robust in compen­ sating for observational errors, and, in a qualitative sense, works well when only reflectivity measurements are available.


INTRODUCTION
The law of Lagrangian conservation for reflectivity can be expressed by:

321
(1) wherel] is the radar reflectivity, u, v, and ware the wind velocities along x, y, and z coordi nates respectively, and t stands for time. Equation (1) is widely used by many single Doppler wind retrieval techniques (e.g., Tuttle and Foote, 1990;Sun et al. 1991 ;Qiu and Xu, 1992;Laroche and Zawadzki, 1994;Sun and Crook, 1994;Xu et al., 1994;Shapiro et al. 1995, etc.). However, this conservation law is valid only for some special cases, such as when alumi num chaff is used a passive tracer. In reality, ( 1) is merely an approximation.
Due to the discrete character of radar measurements, estimates of local derivatives of radar data by finite differencing between two observations from consecutive volume scans may lead to large errors. Gal-Chen (1982) suggested that this error can be significantly re duced by introducing a moving frame of reference in which the local derivatives are as small as possible. This is equivalent to finding a set of optimal velocities by which the reflectivity field is advected while (1) is least violated. Based on this concept, Zhang and Gal-Chen (1996) developed an advection retrieval scheme that retrieves three-dimensional wind fields from multiple time levels of single Doppler radar data. However, the wind fields retrieved by their scheme are not guaranteed to satisfy, even approximately, the continuity equation, because the constraint of mass conservation is not included. This may cause "shocks" when the retrieved wind fields are used to serve as initial conditions, or assimilated into a numerical model.
In this study, the moving frame technique is improved by implementing the continuity equation, as well as the weak vertical vorticity constraints, into the formulation. An idealized Rankine circulation is utilized to explore the applicability of this modified scheme.
This paper is organized as follows. In the next section, the method of estimating a set of optimal advection velocities is briefly reviewed. Section 3 introduces the method to retrieve perturbation velocities. Section 4 reveals the experimental designs. The retrieval results are presented in Section 5, followed by conclusions and discussion of future works in Section 6.

ESTINATION OF OPTIMAL ADVECTION VELOCITY
In order to satisfy (1) as nearly as possible, Gal-Chen (1982) suggested that a set of opti mal advection velocities U, V, W can be obtained by minimizing the local temporal derivative of a scalar field (or reflectivity in the radar observation) in a moving frame. If the formulation is rewritten with a fixed coordinate with respect to the ground, one gets: (2) W,represents the terminal velocity of the precipitation particles, and can be estimated empiri cally from the reflectivity data (e.g. Kessler, 1969). 0 covers the whole area where the re trieval is intended, and extends in the time domain from the first to the last radar scans. a is a weighting coefficient for each grid point, and should be selected to balance the magnitude of each penalty term in the cost function. Similar to Zhang and Gal-Chen (1996), we have: (3 ) Here the overbar denotes a temporal average over a time period of several volume scans dur ing which the retrieval is conducted. Note that U, V, and Ware not included in (3) because they are usually not available a priori. Unlike (2), the integration is performed within a small re gion ( ro) centered at each grid point. For this paper the range of mis chosen, based on trial and error, to be 7 dxX 7 dyX 1 dz, where dx, L\y, L\z are the grid distances. Setting the derivatives of J with respect to U, V, W to vanish, one obtains: where the coefficients are: The solution of (4) only provides a first guess for the optimal advection winds. Neverthe less, a new reference frame can be established by: Here t}s the reference time. All the radar data are redefined in this moving frame. Then, ( 4) -(6) are solved again for this newly defined moving frame in order to get a second adjustment of the advection wind. This procedure is repeated until the solutions for U, V, W converge.
Note that only reflectivity data are involved in this calculation.

Formulation of the Cost Function
The true wind field, assumed to be steady, is decomposed into: where U and V are the optimal advection velocities described above, and LI ,V' ,w" are the unknown perturbation velocities. A cost function is formulated as: where , r( is the reflectivities defined for the moving frame, p; , p� , andp; are the Cartesian coor dinates of the radar for the moving frame, while r' is the distance between the radar and the grid point. They are defined by: The first constraint J1in equation (8) simply states that when moving along with the Lagrangian coordinate, the retrieved Cartesian perturbation wind components LI ,v ' and w ' should make the reflectivities obey the conservation law as much as possible. In addition, constraint 12 requires the satisfaction of the well-known geometric relationship between Cartesian winds LI ,v ' , w' and the radar detected radial velocity v;. 13 is implemented to reduce the diver gence, while 14, the weak vertical vorticity condition, is added to serve as a smoothness con straint to prevent ill-conditioning and suppress the noise (Sasaki, 1970;Qiu and Xu, 1996;Xu et al. 1995). The importance of J3and 14 are examined below. Equation (8) is a functional with very large number of dimensionalities. In order to mini mize it, a limited memory quasi-Newton conjugate-gradient algorithm named V Al SAD (Liu and Nocedal, 1988) is employed through which the minimum is sought by iteration. However, the searching process requires an estimation of the derivatives of J with respect to the un knowns (U', v', w ' ) at each grid' point during each iteration. To measure these quantities, the variation of functional J is taken, and one obtains: :� == JT a 4� � dt � , :::: :: fr -a 3 d�d t (13) Using an initial guess for the wind field, and the gradients of J computed by (10)-(16), the minimization algorithm finds a set of optimal wind perturbations at each grid point which makes the functional J reach a steady state near the minimum. Finally, the perturbations ( u' , v', w') and the advection wind fields ( U, V) are added together to get the total wind field, as shown in (7).

Determination of the Weighting Coefficients
The weighting coefficients control the relative magnitudes among the penalty terms in a cost function. In this study, a1, a2 and a3 are selected so that the first three terms in (8) have the same order of magnitude. To achieve this, a1 is estimated by a formula similar to (3), except that the data defined for the moving frame are now utilized for evaluation, and a time variation is allowed. That is: On the other hand, a2 in each radar scan is measured by: To specify a3, the radar observed radial velocities are projected along the x and y directions to get an estimation of the u and v fields. These partially measured wind fields Cu'', v") are then applied to assess the divergence, and an estimation of a3 is given by: Here the overbar denotes a time average. As mentioned earlier, the vertical vorticity term is introduced in (8) as an operator to eliminate the noise. To prevent the results becoming over smooth, the magnitude of this term, controlled by a4, should be the smallest among all the penalty terms. Experiments indicate that letting I a4 X141 be approximately I% of the previous three terms would be a suitable choice. Therefore, a4 is determined by: where the vertical vorticity is also calculated using the partially estimated ii' and v" fields from the radial velocities. As shown below, although the magnitude of this 41h constraint is small, its existence can improve the retrieval results.

Verification
In order to verify the retrieval results, several indices are evaluated. If A is a variable, such as for u, v, or the cross-beam wind component "4i, and the subscripts "t" and "r" stand for the "true" or "retrieved" wind field respectively, then the root-mean-squared error (rmse) is ob tained by: (21) where N is the total number of grid points. The spatial correlation coefficient (SCC) score of two variables is defined by: where the overbar is the spatial average in the retrieval domain, or -1 "" A=-L-A N In addition, the directional difference is evaluated by the following formula:

Data Sets
Only artificial data are utilized in this study to explore the performance of this method. An idealized circulation is generated by a Rankine vortex in which the wind field is decided as: where vr}lmax 40mls, VRmax= -10 mis, R is the distance from the vortex center, and � = 15 km. e�and e r denote unit vectors along azimuthal and radial direction respectively, 11,rf and /LR equal to 1 for 0 < R < l\nax , and -1 for R > l\nax. Equation (25) is superimposed onto a 10 mis mean flow from northwest, or 315°. Figure 1 shows the complete wind field associated with the idealized Rankine vortex. The reflectivity field is advected by the wind field described above according to the fol lowing governing equation:

J7J
J7J J7J J7J The last term on the right hand side of the equation (S) denotes a source/sink term, and can be viewed as a measure of the extent to which the reflectivity is not conserved. In this experiment, its magnitude is set to be 10% of the advection term.
Equation (26) is integrated for 25 minutes, with horizontal and vertical resolution to be 1000 and 500 meters respectively. The radial wind and reflectivity data are collected at the 1 ()lh , 13th, 16th, 19th and 22th minute by an imaginary radar located about 20 km southwest of the lower left comer of the domain, to imitate five volume scans with three minutes of separation. With this geometric configuration, the direction of the mean flow is nearly perpendicular to the radar beam, thus cannot be directly measured. The contributions from the vertical and terminal velocity to the radar-observed radial velocity are considered to be negligible. This is equivalent to assuming that the radar scans at low elevation angles. Figure 2 shows the reflectivity field at the 16th minute of the simulation, along with radial velocity detected by the radar in Figure 3. The un, observed azimuth velocity is presented in Figure 4.  field. Assumed to be the "true" field .

RETRIEVAL RESULTS
A reference run is first conducted for the purpose of comparison. In this case, three sets of radial wind and reflectivity fields, with 3 minutes of separation, measured by the imaginary Doppler radar at the 13th , 16th and 19th minute of the simulation, are utilized as input data. Figure 5 plots the retrieved total wind vectors, and shows that the circulation structure is suc cessfully recovered, although the retrieved maximum wind speed is reduced to approximately 80% of the true wind. It should be pointed out that, unlike the traditional single Doppler radar analysis technique such as volume velocity processing (VVP, Lhermitte and Atlas, 1961;Waldteufel and Corbin, 1979), the resolutions of the retrieved results remain the same as the input radar data. The contours of the retrieved radial and azimuthal components are displayed in Figures 6 and 7 respectively. The magnitudes of wind extremes in both retrieved fields are much smaller than the true ones. Although this may be a limit to the single retrieval scheme, as suggested by Yang and Xu ( 1996), it may also be attributed to the formulation of (8) in which the information of the radial wind is designed to be distributed to three Cartesian wind compo nents u', v' and w'. A possible alternative is to rewrite (8) on a radar spherical coordinate so that the shape and strength of the radial wind can be preserved optimally.
The remaining retrieval experiments are organized into three groups to investigate the property of this technique, and the results are discussed as follows:

Group A
Experiment A.1 explores the retrieval performance with only two radar scans, separated by the same time interval, as in the reference run (i.e. 3 minutes). Therefore, only data col-    lected at the 13th and 16th minute are adopted. Experiment A.2, on the other hand, also uses three sets of radar data, but with a longer time interval (i.e., 6 minutes) between each scan. Thus, radial wind and reflectivity observed at the 1Qlh ,16th and 22th minute during the simula tion are employed. The retrieved total wind fields for A. 1 and A.2 are shown in Figures 8 and  9, while the resulting statistics, including these from the reference run, are listed in Table 1. Apparently, when the number of radar scans is too sparse, or the time period between scans is too long, the results degrade significantly. Therefore, a suitable observational sequence for this single radar retrieval technique is a sufficient number of scans in a short period of time.

X(KM)
The former is required to make (8) an over-determined problem, and to prevent the ill-condi tioning from happening. The latter, on the other hand, is needed so that the assumption of quasi-steady wind field is not violated too seriously during the retrieval period. Nevertheless, the scan strategy proposed by the reference ruri, that is, at least three scans, separated by three minutes, is achievable by today's Doppler radar capacity. When a typhoon is approaching, and determining the structure of the low levels circulation is desired, then it should be justifiable to conduct intensive sector scans at low elevation angles in order to produce appropriate data sets for this single Doppler retrieval scheme.

GroupB
Experiments in group B are designed to investigate the role played by the 3rd and 4ih constraints in (8), namely, the weak constraint on mass conservation and small vertical vortic ity. In Experiment B.1, the constraint of mass conservation is eliminated, while in Experiment B.2, the retrieval is performed without the weak vertical vorticity constraint. Retrieval results are exhibited in Figure 10 and 11 respectively, as well as in Table 2. They clearly demonstrate   the importance of these two constraints. Suppose there are N (N > 2) data sets, using cubic spline function to estimate the time derivative terms, and the complete form of (8) with a total of four terms, one obtains an over-determined problem which contains N X 4 equations for three unknowns (u', v ' , w') at each grid point. By contrast, without the 3rd and 41h penalty terms, the possibility of getting ill-conditioning increases. Consequently, the retrieved results deteriorate significantly. In Experiment B.3, the 3 rct and 4u. constraints are retained, but a random error field, with 10 mis of maximum amplitude, is introduced to the radial velocity to imitate observational errors. From Figure 12 and Table 2, it is found that due to the existence of the 3rd and 4ih constraints, the retrieval scheme is quite robust in compensating for observational errors. This property is encouraging since eventually it will be desirable to apply this method to real Dop pler radar data sets in which errors will be inevitable.

Group C
This group contains one experiment to study the impact on the retrieval results when only reflectivity data are available (e.g., from a conventional radar). Throughout the whole retrieval process, it is found that information of the radial velocities is only needed in (8) when calculat ing 12, and the weighting coefficients a 3 , CX4 in (19) and (20). By setting a2 at zero, one eliminates contributions from the radial velocity constraints. Unfortunately, without the radial wind, the measurements of the divergence and vertical vorticity, even approximately, are im possible. Consequently, the magnitudes of a3, a4 cannot be determined either. Under this condition, we select a rough estimate for the divergence, and the weighting coefficients are defined as: . (28) The retrieved wind field is displayed in Figure 13. It is shown that, at least qualitatively, the principal cyclone-type flow feature is still recognizable. For an operator of a conventional    radar, this extra piece of information should be very helpful in terms of getting a general picture of the complete wind structure.

CONCLUSIONS AND FUTURE . WORKS
The concept of performing the single Doppler wind retrieval on a moving frame of refer ence is further improved by implementing the incompressibility and weak vertical vorticity constraints into the formulation. The resulting scheme is able to construct the complete wind field without sacrificing data resolution. A simulated Rankine vortex circulation is utilized to RETRIEVED WIND (EXPT. B3) +00

+25
.405E+02 test this modified method. Several conclusions, as well as the plan for future works, are dis cussed in the following: a. In order to apply this method, a general rule for the scan strategy is: conduct a sufficient number of radar scans within a short period of time in the retrieval area. Experiments shown in this study imply that at least three scans, with three minutes of separation between each scan, are required. This strategy is achievable by modern Doppler radar. When a severe weather event, such as a typhoon, occurs, it is justifiable to execute a special scan sequence for the purpose of using this technique.
b. This modified method is robust in compensating for observational errors. This property is desired because the ultimate goal is to apply this method to real Doppler radar data in which errors are inevitable.
c. At least in the qualitative sense, this modified method performs well in retrieving the total wind field when only reflectivity data sets are available. In other words, the application of this method to a conventional radar is feasible. d. For future works, application of this method to a real case study -Typhoon Herb (1996), is in process. In order to make the best use of the radial wind observations, another version of this scheme, formulated on a radar spherical coordinate, is currently under development. It is believed that, by using this coordinate, one can quantitatively increase the accuracy of the retrieved wind.