Retrieval of three-dimensional wind fields from Doppler radar data using an efficient two-step approach

In this work, we describe an efficient approach for wind retrieval from dual Doppler radar data. The approach produces a gridded field that not only satisfies the observations, but also satisfies the anelastic mass continuity equation. The method is based on the so-called three-dimensional variational approach to the retrieval of wind fields from radar data. The novelty consists in separating the task into steps that reduce the amount of data processed by the global minimization algorithm, while keeping the most relevant information from the radar observations. The method is flexible enough to incorporate observations from several radars, accommodate complex sampling geometries, and readily include dropsonde or sounding observations in the analysis. We demonstrate the usefulness of our method by analyzing a real case with data collected during the TPARC/TCS-08 field campaign.


Introduction
The retrieval of the wind field from Doppler radar observations is critical for the study of fluxes and time tendencies associated with convective systems (Marks and Houze, 1987;Kingsmill and Houze, Jr., 1999).However, Doppler radars are only able to measure velocities toward or away from their antenna.Hence, at least 3 independent views of the same volume of the atmosphere are required to infer the full velocity Correspondence to: C. López Carrillo (clopez@kestrel.nmt.edu)vector of the wind at the region of observation.Furthermore, most radars used in the study of convective systems operate within the centimeter wavelength range, so they can only detect precipitation-size particles.Therefore, an assumption must be made about the relation between wind and particle velocities.The standard assumption is that precipitating particles move horizontally with the wind, while falling relative to it with a particle terminal velocity.An estimate of the terminal velocity is generally obtained from the reflectivity measured by the radar, which is a strong function of the particle diameter.The fact that most meteorological radars can only sense precipitation size particles precludes the observation of wind motions at high elevations and hence the observation of divergence patterns at those levels.Nevertheless, when a measurement is made, it is typically very accurate (within 1 m s −1 ).Since a measured Doppler velocity is actually an average over the particles within the observation volume, most Doppler radars also provide some indicators of the quality of this velocity, for instance the spectral width and the signal-to-noise ratio associated with every measurement.
When observations of the same volume of the atmosphere can be obtained from three or more independent directions, the wind field can in principle be recovered unambiguously.The synthesis of the wind field in this case is called multiple Doppler synthesis.Data sets amenable for this kind of synthesis are rare, however.A more typical situation is that only two radars are available and therefore only two independent measurements, which may be amenable to the so called dual-Doppler synthesis.A single Doppler radar mounted on a moving platform can be made to scan so that it produces data close to the dual-Doppler situation (Marks and Houze, 1987;Raymond et al., 1998).The sampling strategy in this case is Published by Copernicus Publications on behalf of the European Geosciences Union.called pseudo-dual Doppler, and the techniques to synthesize the wind field are essentially the same as in the real dual-Doppler case.Despite the fact that many advances have been made in the retrieval of wind fields, the sampling strategy of the wind field is what ultimately determines what information can be retrieved from the measurements.Among the problems that have been studied are the lack of sampling at low and high elevation angles, clutter contamination, evolution and propagation of the system, see for instance Ray et al. (1980), Shapiro and Mewes (1999), and references therein.
In this work we focus on Dual-Doppler techniques, which have their origins in Armijo (1969) and Miller and Strauch (1974).Lacking enough information to produce a threedimensional wind field, dual-Doppler techniques have to apply an extra constraint on the field in order to retrieve the full velocity vector at the observation point.Typically, the constraint is the anelastic mass continuity equation.There are, in general, two approaches to deal with this constraint: (1) in the traditional method a direct integration of the mass continuity equation is performed to obtain the vertical velocity; (2) in the variational approach, the residual of the mass continuity equation is incorporated into the cost function, which measures the misfit between the observations and the synthesized field.In the variational approach, the anelastic mass continuity equation can be incorporated as a weak or strong constraint.Originally it was incorporated as a weak constraint by Sasaki (1970) and illustrated with model and real data by Gao et al. (1999).More recently the usefulness of this approach as an analysis tool for mesoscale fields has been demonstrated by Reasor et al. (2009).The anelastic mass continuity as a strong constraint has been implemented using adjoint techniques by Laroche and Zawadzki (1994) and Protat and Zawadzki (1999).Lagrange multiplier techniques have also been used for this purpose, see for instance Ray et al. (1980).
Several authors have pointed out that, in general, the main problem with traditional approaches is the accumulation of error during the direct integration of the mass continuity equation.In particular for methods that first obtain the horizontal velocities on a Cartesian grid, due to the loss of information during the synthesis of the horizontal velocities from the radar observations.In order to address this problem, traditional methods has been refined to incorporate an iterative approach to the integration of the mass continuity, Ray et al. (1980), Hildebrand and Mueller (1993).This refinement improves the retrievals, but unfortunately this approach becomes unstable for large elevation angles Dowell and Shapiro (2003).
By avoiding a direct integration, the variational approach produces better estimates of the three-dimensional wind field, Matejka and Bartels (1998).Gao et al. (1999) have also pointed out the potential of the variational approach for including various kinds of observations and a background field to further constrain the wind field.Furthermore, the variational approach is also able to include other constraints.For instance, some authors have also explored the inclusion of the anelastic vorticity equation with the aim at improving the estimates of the vertical velocity.The problem with the vertical velocities arises mainly because for some sampling geometries the radar is unable to make observations at low and high elevations angles.Shapiro et al. (2009) have recently used the voritcity and mass continuity as weak constraints and found a substantial improvement in their retrievals when the vorticity equation is included.On the other hand, Protat and Zawadski (2000) who have also used the vorticity equation as a weak constraint, but in combination with the mass continuity as a strong constraint did not find much improvement in their wind synthesis, although they reported an improvement in their thermodynamical retrievals.Because the vorticity equation is a time dependent equation, its usage comes with additional sampling problems and in order to take advantage of it one must be able to compute time derivatives of the vorticity field.We will not explore its usage in this work, but note that we will use the anelastic mass continuity as a strong constraint.
Given its advantages, it is no surprise that the variational approach is now the tool of choice for the synthesis of the wind field from dual-Doppler observations.The main drawback of this approach may be that it is computationally very intensive, due to the enormous amount of data that have to be processed by the minimization algorithm.This problem could be avoided were the data reduced before it enters the minimization process.However, such reduction is generally avoided because its use with traditional methods that use Cartesian grids has resulted in the loss of valuable information, in particular information associated with large elevation angles.However, if one does not limit oneself to a Cartesian geometry, the so called COPLANE approach, pioneered by Armijo (1969), illustrates the possibility to reduce radial velocities without loss of information.This method is able to recover the component of the velocity vector that lies on the plane defined by the two observing directions.It has been shown (Testud and Chong, 1983;Chong and Testud, 1996) that the reduction of the radial velocities by the traditional COPLANE approach results in no loss of information, at least in principle.The problem with Armijo's procedure is not the reduction of radial velocities, but the integration of the mass continuity equation (in particular with its boundary conditions).
In this paper, we present an implementation of the variational approach that reduces the amount of data processed by the minimization algorithm, while keeping the most relevant information from radial observations.We do this in two steps: first, we develop a technique to produce a reduced data set with minimal loss of information; second, we use the reduced data set to build the cost function used during the variational step.Our method is able to estimate the variances associated with the reduced set of observations.These variances can then be used to properly weight the reduced set of observations in the cost function.We also show how sounding measurements can be incorporated into the reduced set of observations to further constrain the velocity field.
Our approach to the variational method is presented in Sect. 2. In Sect.3, we present an application of the method to real data.Section 4 offers a discussion of the method.

Methodology
In this section, we detail the two processing steps required by our method to produce a three-dimensional wind field.Because the steps are based on the weighted-least squares technique, which is not robust to the presence of outliers in the data, the steps depend on the data being outlier-free.Therefore a preprocessing step is always necessary to clean up the data as much as possible, in order to mitigate the outliers' effects on the results.

Preprocessing step
Outliers in Doppler radar measurements come from several sources and can be difficult to remove.Examples of sources of outliers are ground and sea clutter contamination and folded velocities due to limited Nyquist interval.Typically, many hours of manual processing are needed to obtain research-quality data.While this processing preserves most of the data, it is highly time consuming and quite subjective.
Our approach to outliers is to apply conservative thresholds that eliminate most of the questionable data.The approach tends to discard more data than is perhaps necessary, but it is quick, objective, and in our experience still preserves enough data to do the retrievals.
The same conservative threshold approach has also proven effective in eliminating bad data after the first step of our analysis, the gridding step.This is because the gridding step also generates various fields associated with the quality of the gridded velocities, so we can easily eliminate gridded velocities that do not satisfy our quality control.Bad gridded velocity components may occur due to outliers in data, but can also be the result of a poor sampling geometry or large time differences among the observations contributing to a grid point.
It should be noted that elimination of outliers is a projectdependent task, and the approach to their elimination has to be decided on a case-by-case basis.In Sect.3, we present a specific approach for the case of a developing typhoon sampled using a pseudo-dual Doppler technique.

Data gridding step
In this step, data are reduced to a regular three-dimensional grid.The grid may be in Cartesian or spherical coordinates.Grid steps are fixed, but may be different in each direction.In Sect.3, we illustrate our technique using spherical coordinates where a point above the surface of the Earth is identified by its longitude, latitude and altitude.
The reduction of the data around each grid point proceeds essentially as a local fit using a weighted-least squares technique.The weights for the fitting are defined in terms of the distance between the observation and the grid point.Only observations within a volume of influence are considered for the local fit, and the volume of influence is defined here as the eight cells surrounding a grid point.
It is important to emphasize that our gridding step is not a simple interpolation of the radial velocities to the grid points, but it is a search for the particle velocity vector at each grid point that best fits all the observations inside the volume of influence.It is also important to keep in mind that observations within the volume of influence may come from one or several radars.
The target function used in the local fit of the Doppler velocities is where v p is the particle velocity vector, v ri is the observed radial velocity, n i is a unit vector pointing radially away from the radar antenna, and N is the total number of radial velocity observations influencing the grid point.The observations are assumed to be uncorrelated and the variance associated with each observation is defined as σ 2 i = σ 2 0 /w i , where σ 0 is the error in the measurement, and w i is a weighting factor that depends on the distance from the observation to the grid point.Observations closer to the grid point are given larger weights: where the normalization factor Q is adjusted so that i w i = 1.All observations used are within one grid step of the grid point.The points (x i ,y i ,z i ) and (x g ,y g ,z g ) represent the observation and grid point, respectively, and δx,δy,δz are the grid steps.It should be noted that the coordinates (x,y,z) are not necessarily Cartesian coordinates, but may be latitude, longitude, and elevation.Minimization of (1) with respect to the components of v p leads to a system of equations called the normal equations, see for instance Rice (1988), which can be written as is the vector of observations; W is an N × N diagonal matrix whose entries are 1/σ 2 i ; and A is an N ×3 matrix, whose rows are the direction cosines of the radar rays, n i .Here the direction cosines are taken with respect to a Cartesian coordinate system, whose horizontal plane is parallel to the surface of the Earth at the location of the radar and with horizontal axes aligned in the usual east-west, north-south directions.We call this the standard Cartesian coordinate system.Note that by solving Eq. ( 3), we obtain the components of v p in these coordinates.
The matrix S = A T WA on the left-hand side of Eq. ( 3) is referred to as the system matrix.Since in general S is nondiagonal, the components of v p are not independent and the system of Eq. ( 3) must be solved simultaneously.Therefore if there is not enough information to obtain the three components, then all of them are uncertain.This situation can be avoided by noting that we can find a rotated reference frame where the components of v p are independent.The rotated frame can always be found because the system matrix is real and symmetric and therefore there is an orthogonal transformation R such that D = R T SR is a 3 × 3 diagonal matrix whose elements, a α , are the eigenvalues of S. Furthermore, for some vector U in the rotated reference frame, we can write v p = RU .So, in terms of U , the normal equations can then be written as Left-multiplying both sides of Eq. ( 4) by R T , we obtain a diagonal system of equations for the particle velocity in the rotated reference frame At this point the set of equations represented by Eq. ( 5) are independent of each other.Calling b α the elements of the 3×1-column vector that results on right-hand side of Eq. ( 5), we can write the individual equations as follows where no sum over α is implied.This form is very useful for actual calculations, because once we have the eigenvalues of S, it immediately tells us what components of U are not well defined, so we can discard them without affecting the other components.
On the other hand, the formal solution to Eq. ( 5) can help in elucidating the nature of the covariance matrix associated with gridded particle velocities in the rotated reference frame.The formal solution can be written as Rice (1988).
Since we are assuming that the radial velocities are uncorrelated, and their covariance matrix is given by the inverse of W, the covariance matrix for each gridded velocity can be written as that is, the covariance matrix of the gridded velocities is also diagonal, indicating that the errors associated with the components of U are uncorrelated in the rotated reference frame.Hence, the radial velocities around a given grid point are reduced to the single velocity U , and the variance of each of its components is given by the reciprocal of the associated eigenvalue.That is, the error associated with each component of the eigen-velocity is given by This is in agreement with the case where an eigenvalue is zero and the associated component of the particle velocity in the rotated frame is undefined.
One of the advantages of working in the rotated system is that the eigenvalues of the system matrix immediately tell us how the observing geometry affects the relative quality of the observations of the particle velocity.For instance, if we only have measurements along one direction (single Doppler), the eigenvalue associated with that direction will be the biggest of the three and the others will be close to zero.On the other hand, in a multiple Doppler situation where three or more observing directions are available, all the eigenvalues are different from zero.For the radar scanning strategy of interest here (dual or pseudo-dual Doppler), the observing directions of the data being reduced at each grid point loosely define a plane.In this case, our technique results in one eigenvalue being much smaller than the other two.The eigenvectors associated with the two large eigenvalues define the plane that best fit the observations.The eigenvector associated with the small eigenvalue indicates the direction where there is almost no information.Thus, the synthesized particle velocity vectors are obtained on planes where they best fit the contributing data.These planes will generally not coincide with the normal Cartesian planes.For an observing system using a mobile platform, the orientation of the planes will change as the platform moves and turns.
Since for the dual or pseudo-dual Doppler case, our technique uses data from observing directions that only approximately define a plane, it can be considered as an extension of the COPLANE technique.

Reflectivity
The gridding of the radar reflectivity is important because of its relationship to the particle terminal velocity.A local weighted least squares fit is used for reflectivity data.We use the same volume of influence and the same form of σ i as in the case of the velocity fit, but σ 0 now represents the error measurements on the reflectivity.The target function for the reflectivity synthesis is where R i and R g are the reflectivities at the observation and grid points, respectively.Minimization of Eq. ( 10) provides the best fit to the data inside the volume of influence.If there are no reflectivity data inside the volume of influence, the grid point is assigned a missing value; otherwise it gets the result of the fit.

Sounding or dropsonde velocities
In addition to basic thermodynamic information like temperature, humidity, and pressure, sounding or dropsonde data contain information about the horizontal wind velocity.All the fields are typically sampled at about 1 mb intervals in the vertical.However, the horizontal sampling is very sparse, at least when compared to radar data.In order to obtain the value of the horizontal components to be assigned to the grid points, averages of the zonal and meridional components inside each grid box are obtained.The results are then assigned to the nearest corner of the grid cell.For the purpose of incorporating gridded sounding velocities in our variational analysis, we assign eigenvectors pointing in the east and north directions to the zonal and meridional components, respectively.Also the eigenvalue associated with each component is set to 1.This value corresponds to the maximum quality that a gridded velocity from radar data can have, for example when the eigenvector is along the observing direction of the radar.We have chosen this value because the zonal and meridional components from sounding measurements are of better quality that their radar counterparts; typically the error along the ray of the radar is about 1 m s −1 , while for soundings is about 0.5 m s −1 .

Variational step
After gridding the data and eliminating bad gridded velocities, we proceed to generate a three-dimensional wind field.Following the weighted least squares methodology, we seek a field that minimizes the misfit between the resultant field and the gridded data.However, in this case the minimization problem is constrained by the anelastic mass continuity equation.Also a smoothing constraint is imposed on the resultant wind field.Hence, according to the penalty-function methodology for solving constrained minimization problems (Nash and Sofer, 1996), we introduce the constraints as penalty terms in the target function.The resultant target function is given by where G indicates the sum is over the whole grid; v is the resultant wind field; v t is the particle terminal velocity, and U α are the gridded data components computed from Eq. ( 7).The first term is a measure of the misfit with the gridded particle velocities, the second a measure of the smoothness of the solution, and the third is the residual of the mass continuity.
The misfit term is a sum over the components of the gridded particle velocity U α .Each term of this sum is weighted by the associated eigenvalue, see Eq. ( 8).Note that the direction of each of these components is defined by the associated eigenvector, e α .Furthermore, the misfit term requires an estimate of the particle terminal velocity, which is obtained from the gridded reflectivity field using the empirical relationships of Joss and Waldvogel (1970).For this estimate, we assume that only rain exists below the height h r and only snow above the height h s .For the in-between layer, a mixed-phase model is built by a simple linear combination as follows: where h is the height of the grid point and v tr and v ts are the terminal particle velocities of rain and snow, respectively.The smoothing term is defined in terms of operators that act independently in each of the grid directions.The horizontal smoothing is weighted equally in both directions, while the vertical smoothing has its own weighting coefficient.Furthermore, note that the smoothing is only applied to the horizontal components of the velocity.At this point the smoothing operators, P h1,2 and P v , are quite general, but in the appendix, we present their form as currently implemented in our system.
Note that in an effort to counteract the effect of the reduction in density with altitude, we have decided to weight the residual of the mass continuity by 1/ρ(z).In the current implementation of our system, the anelastic mass continuity is imposed as a strong constraint in that its penalty parameter, W m , is systematically increased until its residual is lowered to a desired level.On the other hand, following Sasaki (1970), the penalty associated with the smoothing operators is imposed as a weak constraint.That is, the penalty coefficients W hs and W vs are kept constant during the minimization.
No lateral conditions are imposed during the minimization procedure.However, the values of the vertical velocities at the lower and upper boundaries are kept constant by setting to zero the partial derivatives of the functional, with respect to the vertical velocities, at these boundaries.Vertical velocities at these boundaries are set to zero on the initial solution.
Once the smoothing penalty parameters have been selected, the minimization of Eq. ( 11) proceeds as follows: 1. Set wind field to zero everywhere.
2. Find the initial solution for the wind field by fitting only the gridded data.(penalty terms for smoothing and mass continuity are set to zero).
3. Restore penalty terms: fix parameters for the smoothing operators, and initialize the penalty parameter for the residual of the mass continuity.4. Starting at the previous solution, search for a new wind field that minimizes the full target function for the current penalty parameters.
5. Check if the solution at step 4 satisfies the prescribed tolerance for the residual of the mass continuity.If it does, stop; otherwise continue, provided we have not used the maximum number of iterations.
6. Update the solution with the result of step 4 and increase the penalty associated with the residual of the mass continuity.Then return to step 4.
The six-step algorithm given above describes our solver for the constrained minimization problem at hand.It cycles from step 4 to 6 until the wind field is mass-balanced within the given tolerance or until a maximum number of iterations have been performed.Since we are working with real data, we find it sufficient to accept a solution for which the residual of the mass continuity equation at any grid point is less than 10 −3 kg m −3 ks −1 .We chose this criteria because it is easy to implement and insures that every column in the domain will be mass balanced within 10 −2 kg m −2 s −1 .Therefore if a vertical integral of the detrained mass flux were used to calculate the vertical velocity at the top of the domain, the accumulated error on those velocities will be on the order of 10 −1 ms −1 .
The global search for the wind field that minimizes the target function is done using the conjugate-gradient implementation of Shanno and Phua (1980) for unconstrained nonlinear minimization.This implementation requires the user to provide the code to calculate the target function and its gradient.In the appendix, we show the expressions used in our current code.3 Wind retrieval from real data

Data sources
As an example of the application of our method to real data, we present the retrieval of mesoscale winds from radar and dropsonde data corresponding to the tropical depression stage of typhoon Nuri (2008).The data were collected during the field campaign TPARC-TCS08, which is described in detail by Elsberry and Harr (2008).
The Doppler radar data were collected from an altitude of 2-4 km using NCAR's (National Center for Atmospheric Research) ELDORA Doppler radar, carried by the P-3 aircraft of the US Naval Research Laboratory (NRL).The configuration of the ELDORA radar during TPARC-TCS08 is shown in Table 1.Dropsondes were also launched periodically by the P-3.In addition, the aircraft was equipped with a Doppler wind lidar, and ancillary instrumentation to measure in-situ winds and thermodynamic variables.Additional dropsondes were launched by a WC-130J aircraft of the US Air Force Reserve 53rd Weather Reconnaissance Squadron, typically from a height of 9 km.
Of primary importance to us are the navigational fields from the P-3 aircraft, the radar data and the dropsonde data from both aircraft.From the radar data, we pay special attention to the radial velocities (v r ), reflectivity (R), and the normalized coherent power (NCP).From the dropsondes we only analyze the horizontal wind velocities.This figure shows plan views of the error associated with the component of the gridded velocity along the eigen-direction with the second largest eigenvalue (see Eq. 9).The marks represent grid points that had more than 50 radar observations associated with them.Marks in red represent errors larger than five ms −1 ; errors between five and four ms −1 are shown in orange, between four and three ms −1 in green, between three and two in blue, and between two and one in magenta (no points had errors smaller than 1).The data for this study were quality controlled and are maintained by the Earth Observing Laboratory at the National Center for Atmospheric Research (NCAR).NCAR is sponsored by the National Science Foundation (NSF).The data, as well as a documentation file including information on quality control procedures and subsequent findings are available from the TPARC project data page (http: //data.eol.ucar.edu/masterlist/?project=T-PARC), under Upper Air and Radar.For the radar data, NCAR also generated an unfolded radial velocity field from the raw radial velocity and corrections for the antenna angles, (Bell et al., 2009).In addition, they provided us with the software to apply these corrections.We use the unfolded version of the radial velocity in our analysis and will refer to it simply as v r .
Radar observations started on 16 August at 22:20:49 UTC.At this time, the storm velocity was (−9.0,0.8)ms −1 and it was centered at (145.5 • E,15.8 • N).The P-3 was flying at a nominal altitude of 2.4 km.We have chosen for our analysis the radar observations taken during a period of 12 ks (3 h and 20 min) starting at 23:15:00 UTC.We used nine dropsondes from the P-3 and thirty from the WC-130J.The collection interval for the P-3 dropsondes is about 4 h starting at 22:43:00 UTC, while for the WC-130J it is ∼7 h starting at 20:34:56 UTC.

Data analysis
Pre-processing of the unfolded radial velocity consisted in discarding gates that satisfy any of the following criteria: (a) their NCP was less than 0.01; (b) they were located within 200 m of the radar; (c) they had an altitude lower than 500 m.No pre-processing was applied to the reflectivity field.
Individual dropsondes were interpolated in the vertical to regular intervals of 150 m.The pressure-to-height conversion assumed a lapse rate of 5.5 Kkm −1 , and that surface temperature and pressure were 300 K and 1020 hPa, respectively.
Although we show the earth-relative winds in our results, the retrieved winds were obtained in a reference frame moving with the storm.Therefore, during the gridding step, the radar data were renavigated to the location they would have had in the comoving frame at a specified reference time -the reference time was taken to be the middle of the collection time interval of the radar data, 00:55:00 UTC on 17 August 2008.This procedure assumes that the storm does not evolve significantly during the collection time.A horizontal plan view of the grid chosen for our analysis is shown in Fig. 1, along with the aircraft track flown during the observations and the winds from the dropsondes.Since this figure shows a snapshot of the co-moving frame at the reference time, the aircraft track and dropsonde locations have also been re-navigated to locations in the co-moving frame at the reference time.The grid step is 0.25 • in both longitude and latitude.The grid extends vertically from the surface to 20 km and has a grid step of 0.625 km.
After gridding the data, we examined histograms of their measures of goodness, which are produced during the gridding step.For instance, we examined the χ 2 value of the local fit, the number of contributing observations, the eigenvalues, and the collection times, at the grid points.Based on that information, we applied conservative thresholds to eliminate questionable gridded velocities.For the current example, it was sufficient to eliminate gridded velocities that met any of the following criteria: (a) their second largest eigenvalue is smaller than 0.03; (b) the number of contributing measurements at each grid point was less than 50.It should be noted that our criteria eliminates points that only had enough information to obtain a single component of the gridded velocity (single Doppler).However, it does not preclude cases where there is enough information to obtain the three components of the velocity (triple Doppler.)Even in cases where the smallest eigenvalue is to small and the associated component of the velocity has a large error.Nevertheless, as Fig. 2 shows there are many instances where even the smallest eigenvalue is larger than the threshold value.The horizontal distributions of the errors associated with these eigenvalues (see Eq. 9) are shown on Figs. 3, 4, and 5, for some selected levels.Figure 6 shows distributions of these errors as functions of their relative eigen-velocities.The relative velocity is used in this figure to eliminate the contribution of the storm motion.
Finally the gridded data were passed to the second step, which requires the smoothing weights.For the present case, we have used W hs = 0.3 and W vs = 0.1 for the horizontal and vertical smoothing, respectively.

Results
Since the variational step incorporates the anelastic mass continuity equation as a strong constraint, our algorithm cycles through steps 4 to 6 until it satisfies this constraint within the required tolerance.For the results presented here, Fig. 7 shows how the residual of the mass continuity is reduced as its penalty parameter is increased.
The results of the retrieval using only radar data are shown in Figs.A comparison between the results using only radar data and those obtained when dropsonde data are included in the retrieval is shown in Fig. 14.This figure also shows the gridded dropsonde velocities.The reflectivity field is also included to show the area where radar data were available.
Although gridded dropsonde velocities are certainly affected by the convective environment, while our retrievals of the mesoscale wind field tend to smooth out convective contributions, Figs. 15 and 16 offer a comparison between gridded dropsonde data and the retrieved velocities from radar only and radar plus dropsonde data, respectively.

Summary and discussion
We have presented a dual-Doppler analysis technique of the variational type for the retrieval of wind velocity fields.The technique is efficient, has the ability to quickly incorporate dropsonde observations of the wind velocity, and it produces a mass-balanced wind field.The efficiency of the technique stems from its capacity to reduce Doppler observations to data over a regular grid while minimizing the loss of relevant information, in particular information collected at high elevation angles.Furthermore, the reduction of the radar data preserves the diagonal structure of the covariance matrix of the raw data and generates estimates for the variances of the reduced velocities on the analysis grid.These variances reflect how the observing geometry affects the estimate of the gridded velocity.This is an important piece of information, especially when the observations are collected by a moving platform.For instance, they allow us to properly weight the results obtained even if the platform is not moving in a straight line.The mass balanced wind, results from the implementation of the anelastic mass continuity equation as an strong constraint.It is worth nothing that our implementation uses the penalty-function methodology, which does not requires the solution of either the adjoint or the Euler-Lagrange equations.
We have illustrated how observations from soundings can be incorporated into the variational analysis to further constrain the resultant wind field.Basically, once the sounding velocities are gridded and weighted, the variational step treats them as if they were radar gridded velocities.Figure 14 shows that the inclusion of the soundings in the analysis has the most impact on regions void of radar data, where the wind field would otherwise be constrained only by the mass continuity equation and the smoother.Although dropsonde data may not be representative of the mesoscale wind at a given grid cell, the comparison shown in Fig. 15 between the gridded dropsonde data and the retrieved wind from radar only data is reasonable, especially for grid points that have radar observations (indicated by blue boxes in the figure).As expected Fig. 16 shows that when dropsondes are added to the synthesis, the agreement between the retrieved wind and the dropsondes is improved.So, it is worth nothing that other types of wind observations can be incorporated into the analysis as well, so long they have estimates for their variances.
We have also shown that when the mass continuity equation is used as a strong constraint during the variational analysis, the resultant wind field is mass-balanced.This characteristic of the resultant wind field is particularly important when such a field is to be used for initializing model runs.
The reduction of the data produces a robust mean value at a grid point and we are able to perform analyses on large domain sizes in reasonable times.For instance, the grid size used for this work has 33 × 33 × 33 grid points on which the analysis takes a few minutes on a laptop computer.Our tests on larger grids (257×257×33) indicate that the analysis can be performed on a PC in about 24 h.Furthermore, the analysis can be done in Cartesian or spherical coordinates, depending on the size of the domain.Cartesian coordinates are simpler, but may be appropriate only for small domains.For large domains, however, spherical coordinates are generally necessary.
In regards to quality, the technique offers plausible massbalanced winds that are consistent with the observations.However, the reader should keep in mind that, like any technique based on least-squares, ours is sensitive to outliers.Therefore, a clean data set is critical to its performance.In this sense, a two-step approach has the advantage of reducing the propagation of the outliers' impact over one region to the global domain.This is because the approach offers the possibility of examining the goodness of the local fit.This possibility allows, for instance, to detect whether the evolution of the system is having a bad impact on the results.
Finally, although not illustrated here, we note that, during the gridding step, the technique allows us to combine observations from multiple Doppler radars.Thus, we not only recover better estimates of the gridded velocities, but also reduce the amount of data processed in the variational step in this situation.
Table A1.Constants to be used with our finite-difference approximations to the derivatives.For inner points, central differences are used.For the first point in the domain, forward differences are used, while for the last point backward differences are employed.the residual of the mass continuity, the smoothing operators, and the data misfit.The gradients of these functions will be presented as updating formulas for the gradient components, which means that the values given by the formulas have to be added to the current value of the component.Updating formulas are preferred over the full form of the gradients, because they are computationally more efficient in iterative processes.We have implemented our method on a regular grid in both spherical and Cartesian coordinates.Here, we present the implementation that uses spherical coordinates.In these coordinates the radial distance is measured from the center of the Earth, so we let r = R E + z, where R E is the radius of the Earth and z is the height above the surface of the Earth.The radial direction is our vertical direction.The zonal direction will be represented by the angle φ, and the meridional direction by the angle λ.They are considered positive towards the east and north, respectively.The number of points in each of these directions is N z , N φ and N λ , respectively.The finitedifference approximations for the derivatives are defined as follows where φ i = i φ , λ j = j λ , z k = k z , and φ , λ , z are the grid steps.The rest of the constants are specified in Table A1.

A1 Anelastic mass continuity
The velocity vector is given as v(φ,λ,z) = v φ (φ,λ,z)e φ + v λ (φ,λ,z)e λ + v z (φ,λ,z)e r (A4) where e φ , e λ and e r are unit vectors pointing in the zonal, meridional and radial directions, respectively.Therefore, the expression for the mass continuity, in spherical coordinates, can be written as Making use of the finite-difference approximations given by Eqs.(A1-A3), we can write a discrete version for the residual of the anelastic mass continuity equation as a linear combination of the velocities in the grid: D(φ i ,λ j ,z k ) = A φ (φ i ,λ j ,z k )v φ (φ i + a φ ,λ j ,z k ) −A φ (φ i ,λ j ,z k )v φ (φ i + b φ ,λ j ,z k ) +A λ (φ i ,λ j ,z k )cos(λ j + a λ ) v λ (φ i ,λ j + a λ ,z k ) −A λ (φ i ,λ j ,z k )cos(λ j + b λ ) v λ (φ i ,λ j + b λ ,z k ) +A z (φ i ,λ j ,z k )(R E + z k + a z ) 2 ρ(z k + a z )v z (φ i ,λ j ,z k + a z ) where Hence, the contribution from the residual of the mass continuity to our target function is given by where the index G indicates that the sum is over all of the grid points, ρ(z k ) is the air density, and W m is the weighting factor for this contribution.

A1.1 Gradient
The components of the gradient of T r can be obtained by considering its variation with respect to an arbitrary component of the velocity, v ψ (φ l ,λ m ,z n ), at point (φ l ,λ m ,z n ): ∂T r ∂v ψ (φ l ,λ m ,z n ) = W m g D(φ i ,λ j ,z k ) ρ 2 (z k ) ∂D(φ i ,λ j ,z k ) ∂v ψ (φ l ,λ m ,z n ) (A11) Since D is just a linear combination of velocities in the grid (Eq.A6), its partial derivative with respect to the given component of the velocity at a given point is given by

Fig. 1 .
Fig. 1.Horizontal plan view of analysis grid.The grid is anchored to the storm, which is moving with velocity (−9.0,0.8)ms −1 .The axes are longitude and latitude at the reference time, 00:55:00 UTC on 17 August 2008.The grid extends from the surface to 20 km with a vertical grid step of 0.625 km.The horizontal grid steps are 0.25 • .The aircraft track flown during the radar collection time is shown, along with the winds from dropsondes, The aircraft and dropsonde positions have been re-navigated to the co-moving grid.The scales for the zonal and meridional winds are shown at the bottom of the plot.Regions where reflectivity exists at an altitude of 2.5 km are shown in yellow.

Fig. 2 .
Fig.2.The lower panel shows the eigenvalues for each grid point that has at least 50 radar observations associated with it.For each point, the largest eigenvalue is shown in black, the second largest in blue, and the smallest in red.The horizontal axis is just the grid point index obtained by traversing the grid in C-style order(the last index varying the fastest).The horizontal line in magenta indicates the threshold value used to accept the associated velocity.The upper panel is a zoom-in of the lower panel.

Fig. 5 .
Fig.5.Same as in Fig.3, but for the error associated with the component of the gridded velocity along the eigen-direction with the largest eigenvalue.This plot shows that the error associated with the largest eigenvalue is always between one and two ms −1 .

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .
Fig.6.This figure shows the overall distributions of the errors for gridded velocities that have passed our quality control.The scatter plots show the behavior of the error as a function of its associated relative eigen-velocity (u r = u − v s , where v s is the storm velocity and u is the gridded eigen-velocity).The lower right panel shows histograms of the errors (solid lines) and their corresponding cumulative distribution functions (dashed lines).

Fig. 10 .Fig. 11 .
Fig. 10.Results of the synthesis using radar data only.A vertical cross section at latitude 14.5 • is displayed.The arrows show the projection of the Earth-relative wind field, while colors and contours show the reflectivity field.The scales for the zonal and vertical winds are indicated at the bottom of the plot.Contours are shown every 5 dbz.

Fig. 12 .Fig. 13 .
Fig. 12. Results of the synthesis using radar data only.A vertical cross section at longitude 140.5 • is displayed.The arrows show the projection of the Earth-relative wind field, while colors and contours show the reflectivity field.The scales for the meridional and vertical winds are indicated at the bottom of the plot.Contours are shown every 5 dbz.

Fig. 14 .
Fig. 14.Horizontal plan view, at height 1.875 km, of two estimates of the wind field along with winds measured by dropsondes.All velocities are Earth-relative.Dropsonde winds are shown in green; the wind retrieved, using radar data only, is in blue; the wind retrieved, using radar plus dropsondes, is shown in red.The scale of the meridional and zonal winds is shown at the bottom of the plot.Areas where reflectivity exist are shown in yellow.Contours of the reflectivity are also shown every 5 dbz.The thick black line represents the aircraft track.

Fig. 15 .
Fig.15.The upper panel shows the difference between the meridional component of the retrieved velocity (using radar data only) and the gridded sounding data.Red marks show this difference for all grid points where there are dropsonde observations, while blue boxes show which of those grid points have radar observations that constrained the retrieved velocity.The horizontal axis is just the grid point index obtained by traversing the grid in C-style order(the last index varying the fastest) The lower panel shows the same, but for the zonal component of the velocity.

Fig. 16 .
Fig.16.The upper panel shows the difference between the meridional component of the retrieved velocity (using radar and dropsonde data) and the gridded sounding data.Red marks show this difference for all grid points where there are dropsonde observations, while blue boxes show which of those grid points have radar observations that constrained the retrieved velocity.The horizontal axis is just the grid point index obtained by traversing the grid in C-style order(the last index varying the fastest) The lower panel shows the same, but for the zonal component of the velocity.