Three-Dimensional Wind Field Retrieved from Dual-Doppler Radar Based on a Variational Method: Refinement of Vertical Velocity Estimates

Chenbin XUE, Zhiying DING, Xinyong SHEN, and Xian CHEN

1Key Laboratory of Meteorological Disaster, Ministry of Education/Joint International Research Laboratory of Climate and Environment Change/Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters,Nanjing University of Information Science and Technology, Nanjing 210044, China

2Jiangxi Meteorological Observatory, Nanchang 330096, China

3Guangdong Province Key Laboratory of Regional Numerical Weather Prediction,Institute of Tropical and Marine Meteorology, CMA, Guangzhou 510080, China

4Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai 519082, China

5Jiangxi Institute of Land and Space Survey and Planning/Jiangxi Geomatics Center, Nanchang 330000, China

(Received 25 January 2021; revised 21 June 2021; accepted 30 June 2021)

ABSTRACT In this paper, a scheme of dual-Doppler radar wind analysis based on a three-dimensional variational method is proposed and performed in two steps. First, the horizontal wind field is simultaneously recovered through minimizing a cost function defined as a radial observation term with the standard conjugate gradient method, avoiding a weighting parameter specification step. Compared with conventional dual-Doppler wind synthesis approaches, this variational method minimizes errors caused by interpolation from radar observation to analysis grid in the iterative solution process, which is one of the main sources of errors. Then, through the accelerated Liebmann method, the vertical velocity is further reestimated as an extra step by solving the Poisson equation with impermeable conditions imposed at the ground and near the tropopause. The Poisson equation defined by the second derivative of the vertical velocity is derived from the mass continuity equation. Compared with the method proposed by O'Brien, this method is less sensitive to the uncertainty of the boundary conditions and has better stability and reliability. Furthermore, the method proposed in this paper is applied to Doppler radar observation of a squall line process. It is shown that the retrieved vertical wind profile agrees well with the vertical profile obtained with the velocity-azimuth display (VAD) method, and the retrieved radial velocity as well as the analyzed positive and negative velocity centers and horizontal wind shear of the squall line are in accord with radar observations. There is a good correspondence between the divergence field of the derived wind field and the vertical velocity. And, the horizontal and vertical circulations within and around the squall line, as well as strong updrafts, the associated downdrafts, and associated rear inflow of the bow echo, are analyzed well. It is worth mentioning that the variational method in this paper can be applied to simultaneously synthesize the three-dimensional wind field from multiple-Doppler radar observations.

Key words: dual-Doppler radar, three-dimensional wind, a variational method, vertical velocity, wind synthesis

1. Introduction

Compared with conventional observations, Doppler radar data can reflect the characteristics of convective-scale systems well with high temporal and spatial resolution. Therefore, it is one of the most important tools for studying mesoscale convective systems. In order to study the dynamic mechanism of the mesoscale systems, it is desirable to gain the three-dimensional wind field of the atmosphere. For this reason, meteorologists have proposed a variety of Doppler radar wind field retrieval techniques. Based on the number of radars required, radar wind synthesis technology can be classified into single radar retrieval and dual- or multipleradar retrieval.

Early wind field analysis from single-Doppler radar data such as the Velocity Azimuth Display (VAD) method(Lhermitte and Atlas, 1961; Caton, 1963; Boucher et al.,1965) and Volume Velocity Processing (VVP) method(Waldteufel and Corbin, 1979; Koscielny et al., 1982),mostly based on the assumption that the local wind field is linearly distributed, is feasibly used to obtain the vertical wind profile through spatial constraints and least squares minimization methods. However, it is less accurate to estimate the wind field when these assumptions cannot be satisfied. With the development of data assimilation methods, a simple adjoint method is developed for retrieving the wind field from a number of consecutive single-Doppler measurements (Qiu and Xu, 1992). It is found that using data over multiple time levels not only provides more information and increases the accuracy of the retrieval but also makes the method less sensitive to errors in the observed reflectivity.In a cost function, the conservation equations of reflectivity or radial velocity can be applied as weak constraints(Laroche and Zawadzki, 1994; Qiu and Xu, 1996) or strong constraints (Qiu and Xu, 1992). This method has been widely used in many cases (Xu et al., 1994, 1995; Xu and Qiu, 1995).

When two or more Doppler radar observations are available, it is possible to retrieve a higher-precision three-dimensional wind field. Armijo (Armijo, 1969) first theoretically deduced the equations for the joint retrieval of the threedimensional wind field in the Cartesian coordinate system from two or three radar observations with the anelastic mass continuity equation. Since the 1960s, with the improvement of radar network construction (Ray et al., 1979), great progress has been made in the joint detection of mesoscale weather systems with multiple Doppler radars. Ray et al.(Ray et al., 1978, 1980) proposed a variational scheme for retrieving the wind field by combining multiple radar and ground observations to improve the error caused by insufficient sampling of the near-ground divergence field, which has greatly promoted the research progress of Doppler radar three-dimensional wind synthesis. Chong and Campos(1996) proposed the EODD (extended ODD) technology with a variational method based on the ODD (overdetermined dual-Doppler, Ray et al., 1980) technology to overcome geometrical limitations and provide solutions to mitigate the inevitable contribution of errors in estimating the vertical velocity. Subsequently, Bousquet and Chong (1998)presented the multiple-Doppler synthesis and continuity adjustment technique (MUSCAT) derived from the EODD variational formalism through a noniterative solution of a dual- or multiple-equation system and an anelastic mass continuity equation, which could remove potential drawbacks of an iterative solution of Cartesian dual-Doppler analysis techniques (Ray et al., 1980; Ray and Sangren, 1983;Chong and Campos, 1996). Gao et al. (1999) proposed a variational dual-Doppler synthesis method for the analysis of three-dimensional winds by imposing the mass continuity equation as a weak constraint and using the radar data directly at observation locations. They found that their method offered much more flexibility in its use of observational data and various constraints, and their results were not very sensitive to the specifications of the boundary conditions.Chong and Bousquet (2001) and Liou and Chang (2009) discussed the recovery of the wind field along the radar baseline. In Liou et al. (2019), a variational method was presented to incorporate the wind observations obtained by wind profilers and radiosondes. In addition, in Shapiro et al.(2009) and Liou and Chang (2009), the vertical vorticity equation was directly built into the cost function to improve the quality of the multiple-radar retrieved wind fields when radar observations were lacking near the ground. Furthermore, Chong and Cosma (2000) and Liou et al. (2012)developed an extension of a variational-based multiple-Doppler radar synthesis technique to construct the three-dimensional wind field over complex topography. It has the advantage of eliminating the need to explicitly analyze the vertical velocity associated with the slope wind near the surface.

It is important to improve the accuracy of vertical velocity when retrieving the three-dimensional wind field from Doppler radar data. In order to minimize the errors caused by direct integration of the horizontal divergence to a certain extent, most of the retrieval methods mentioned above make use of the mass continuity equation as a weak constraint to solve the vertical velocity. Unfortunately, the vertical velocity analyzed through the variational dual-Doppler analysis techniques may still contain large errors if the horizontal divergence field within the data void region is not accurately determined by these constraints alone. In view of the limitations of the above methods for the calculation of vertical velocity, a new method is proposed for calculating vertical velocity based on the mass continuity equation in Cartesian coordinates. The vertical velocity can be directly derived through the accelerated Liebmann method with the Poisson equation defined by the second derivative of the vertical velocity. Because this scheme puts the mass continuity equation in a second-order form, error accumulation caused by direct vertical integration of the horizontal divergence in the iterative solution process is minimized. As a consequence, it is less sensitive to the vertical boundary uncertainties.

This paper is organized as follows. In section 2, the scheme of dual-Doppler radar wind retrieval based on a three-dimensional variational method is introduced, followed in section 3 by the refinement of vertical velocity estimates. In section 4, the performance of the variational method is verified by means of a set of idealized data sampled from a simulated supercell storm. Furthermore, the evaluation of this scheme is carried out on a squall line that occurred in Jiangxi on 11 May 2017, and the dynamic and fine structure of the squall line based on the derived wind field is analyzed in section 5. Finally, the summary and concluding remarks are given in section 6.

2. Description of variational methodology

As shown in Fig. 1, wind field retrieval is performed from dual- or multiple-Doppler radar observations in the Cartesian coordinate system. The radial velocityVrat a given point (x1,y1,z1) can be expressed as

Fig. 1. A schematic view of the geometry of a ground-based radar in the Cartesian coordinate system (r is the distance between the target and the radar, θ is the azimuth angle, φ is the elevation angle, and u, v, and w′ are the velocity components in the x, y, and z directions, respectively).

where α=(x1-x0)/r, β=(y1-y0)/r, γ=(z1-z0)/r, and define the distance between the target and the radar located at (x0,y0,z0).u,v,andware the Cartesian velocity components to be determined, andwtis the terminal velocity of precipitation particles that can be estimated from an empirical relationship with the observed radar reflectivity (Foote and Toit,1969).

whereZis the observed reflectivity, ρ is the empirical formula of air density, ρ0is the air density near the ground,gis the acceleration of gravity,Nis the atmospheric molecular weight,his the height above the ground,Ris the ideal gas constant, andTis the absolute temperature.

According to Eq. (1), the cost function in our dual-Doppler radar analysis may be written as follows:

where subscriptqdefines theqth measurement of a total numbernq(≥2) that can be observed from thepth radar which falls inside an ellipsoid with an influence radius in the horizontal twice as large as in the vertical centered on the grid point.Vis the observed radial velocity defined in observation space, and ω is the Cressman weight function (Cressman, 1959) depending on the distance between the observed radial velocities within the ellipsoid and grid point.Vqis pre-interpolated by using the Cressman weighting function from the observation position along the radar beam onto the Cartesian grid of interest to be included into the data fit. In order to obtain one radial velocity from each radar at each grid point, Eq. (5) is defined to overcome the limitation resulting from the radar data interpolation. In a conventional interpolation step (Ray and Sangren, 1983; Chong and Campos,1996), interpolated grid values are defined as distanceweighted averages of all measured radial velocities falling inside an ellipsoid with a prescribed influence radius,thereby forcing averaging radial vectors pointing in different directions and hence causes a sampling problem. Since these vectors are different due to the changes in elevation and azimuth angles, they should not be simply averaged. Considering that Eq. (5) is defined according to the cumulative weight of radial velocity from each radar, it is therefore more reasonable than the data fitting term in the MUSCAT developed by Bousquet and Chong [see reference (Bousquet and Chong, 1998) Eq. (6)]. In Eq. (5), the interpolation operator which transfers radial velocityVqfrom the observation space to the Cartesian space is included in the pre-interpolated step. This is one of the main differences from the scheme used in Gao et al. (1999) which applies the radar data directly at observation locations, avoiding an interpolation step. In Gao et al. (1999), a background field, the mass continuity equation, and smoothness constraints were incorporated into a cost function, allowing a flexible application of radar data in combination with other sources of wind observation (e.g., soundings or wind profilers). The pre-calculated values of the cost function weights determine the relative influence of each constraint in the analysis. Each weight must be explored when determining the optimal tuning parameter values. In previous studies (Xu et al., 1995; Gao et al., 1999, 2001; Liu et al., 2005), the weights were chosen empirically based on past experience or batch trials. Actually, a series of numerical experiments should be conducted to find the optimal weighting parameters for each experiment type by evaluating the sensitivity of each weight to each constraint. However, difficulties exist in specifying the optimal weights for each constraint to ensure that the retrieved wind fields are not sensitive to this set of weighting coefficients. In this study, our wind analysis technique is performed through minimizing the cost function Eq. (4)defined as a radial observation term, therefore avoiding a weighting parameter specification step. Considering the minimization with respect tou,v, andwin Eq. (4), the solution foru,v, andwcomponents at the grid point is derived according to a simultaneous resolution of

The final form of the variational analysis ofu,v, andwis given by the following matrix equation

such that

where nine terms for each individual grid point are defined as follows:

Assuming anm×nCartesian grid on the horizontal (mgrid points alongx-axis,ngrid points alongy-axis), the coefficient matrixCis a (3 mn) order, positive definite symmetric matrix,Vis a column vector composed of three successive vectors containingu,v, andwgrid elements, andEis a(3 mn) column vector containing radar observation information. Finally, the matrix in Eq. (8) is simultaneously solved with the standard conjugate gradient method. This is a major advantage over the conventional dual-Doppler wind synthesis approaches (Ray et al., 1980; Ray and Sangren,1983; Chong and Campos, 1996). As in the conventional dual-Doppler analysis, an iterative process is required to determine the horizontal wind components from observations unambiguously, and the vertical velocity is deduced from the mass continuity equation. In consequence, it requires a step-by-step correction of the vertical velocity contribution to the observed radial velocityVr. Since these averaged radial velocities with less accurate estimate are needed multiple times depending on the number of iterations during the iterative procedure, exacerbating the errors in estimating the horizontal wind vectors, we suspect that this type of error is one of the main sources of error. Thus, Eq. (8) in this paper is put forward to overcome this limitation. It should be clarified that the pre-interpolation step as shown in Eq. (5) is needed only once in our method. In other words, it is not required to perform the pre-interpolation step in each iteration during the minimization procedure.Therefore, the three wind components can be retrieved simultaneously. However, it is worth mentioning that this technique requires high memory storage in a three-dimensional space if too many grid points are included, which would be problematic in practice.

Considering that the contribution of vertical velocity to the radial velocityVris much smaller than that from the horizontal wind components due to quasi horizontal radar scanning geometries as shown in Fig. 1, the vertical velocity is the most poorly sampled wind component. This set ofu,v,andwobtained at each grid point from the solution of Eq. (8)only satisfies the radar observations in a least squares sense and thus the resultantwfield from such solution can be very problematic. Therefore, the retrieved vertical velocities from the solution of Eq. (8) are not reliable enough, and they are simply discarded in the following computation. It is worth mentioning that our wind analysis technique is performed in two steps. First, the horizontal wind field is directly recovered through solving Eq. (8) with the standard conjugate gradient method. Then, a refinement of the vertical velocity estimates is performed as an extra step by solving the Poisson equation derived from the mass continuity equation which is imposed as a strong constraint of the wind vectors (u,v, andw) according to the following approach proposed in section 3.

3. Estimation of vertical velocity

Once the horizontal wind components are found, the vertical component is usually obtained by integrating the mass continuity equation upward (or downward) according to the boundary condition at the bottom boundary (or top boundary). The anelastic mass continuity equation can be written as (Ogura and Phillips, 1962)

wherek=-∂lnρ(z)/∂zand ρ is the mean air density at a given horizontal level defined in Eq. (3). Then the shallow convective continuity equation can be obtained with the air density-related termkwomitted.

A wind adjustment (O'Brien, 1970) is typically applied to ensure thatw=0 at both the upper and lower boundaries.It can be rewritten as

Here, a novel method to estimate the vertical velocity is proposed. First, we take the partial derivative of Eq. (10)with respect tozand obtain

whereD=∂u/∂x+∂v/∂yis the horizontal divergence term related with the horizontal wind components derived above.Eq. (13) is a one-dimensional Poisson's equation with unknown vertical velocity (w) on the left and known quantities (ρ andD) on the right. The problem is now second-order inw, and both the bottom and top boundaries can be applied(Miller and Strauch, 1974). Considering that mesoscale convective systems should satisfy rigid boundary conditions as dynamic constraints, Eq. (13) is solved by imposing the impermeability conditionw= 0 both at the ground and near the tropopause. The vertical velocity can be obtained through solving Eq. (13) with a relaxation method, such as the Richardson method (Richardson, 1911) or Liebmann method (Liebmann, 1918), allowing a simultaneous solution ofwto ensure thatw= 0 at both the upper and lower boundaries.[At the initial step,w0= 0; then it is the solution of (13).Note that two boundary conditions are satisfied.] In this study, the vertical velocity is directly solved through the accelerated Liebmann method, which is based on the Liebmann method multiplied by a coefficient to accelerate convergence. We call this method the Poisson method in the following sections.

The potential advantages and limitations of the Poisson method solution are discussed with respect to two aspects.On one hand, as discussed in O'Brien (1970), imposing two boundary conditions onwwill often introduce errors smaller than those resulting from using only a lower boundary condition. The O’Brien method is an effective solution used to reducewerrors that accumulate during integration of the mass continuity equation. The weighted correction forwis almost quadratic in Eq. (12), with little correction in the lower part of the atmosphere and a large correction in the upper atmosphere. The O’Brien method and Poisson method proposed herein both use the same mass continuity equation and make the problem second-order inw, thus it seems that there should be no difference. It should be pointed out that in Eq. (12), the solution of O’Brien method is obtained under special adjustment criteria, considering that the weighting parameter is a linear function ofkor pressure and performs a direct correction tow. However, the error is not inwbut in the estimates of the horizontal divergence field. By contrast, the Poisson method directly solves the second-order equation with respect tow, and the secondorder adjustment scheme implied in the Poisson method herein is equivalent to adjusting the divergence at each level by a constant, as Lateef (Lateef, 1967) points out. Hence, considering the reasons mentioned above, the Poisson method has advantages over the O’Brien method. On the other hand, in reality, one of the main problems is the data void region caused by radar detection since the weather radar scans can neither reach the ground nor the echo top, and a data void cone exists directly above the radar at a scan elevation of 19.5°. This will inevitably introduce errors to the retrieved wind fields, even though the variational analysis techniques are employed (Bousquet and Chong, 1998; Gao et al., 1999). The solution for the wind field within the data void region can be obtained from the background wind field, the vertical vorticity equation, and spatial smoothness terms imposed as weak constraints (Shapiro et al., 2009;Liou et al., 2012). And in both the O’Brien and Poisson methods, the adjustment in retrievedwwill be severely limited because of lacking horizontal divergence information. From Eq. (12), it is shown that the adjustment of the O’Brien method strongly depends on the vertical velocity at the top level. And, ifwat the top level is missing or invalid, adjustment is impossible, particularly when the storm is intensifying or decaying. This is the main difference between the O’Brien and Poisson methods as the latter directly depends on the horizontal divergence field in Eq. (13). In spite of that, if a domain covering more radar data is selected for wind retrieval, the above problems can be avoided. It should be pointed out that in Shapiro et al. (2009) and Liou and Chang (2009), the vertical vorticity equation has been used to improve retrievals ofwin common situations where lowlevel data are lacking and the impermeability condition cannot be reliably imposed aloft.

4. Tests with a simulated supercell case

4.1. Model setup and simulation data

In order to evaluate the performance of our three-dimensional variational method proposed in the previous section,we utilize a set of assumed dual-Doppler data from a simulated supercell storm. The Advanced Research Weather Research and Forecasting Model (WRF, Skamarock et al.,2008) is used here to perform a 3-h simulation of a supercell storm. As shown in Fig. 2, the WRF model domain consists of 83 × 83 × 51 grid points with a uniform grid interval of 1 km in the horizontal and model top of 18.5 km in the vertical. The WRF model makes use of the third-order Runge-Kutta explicit time integration scheme with an integration time step of 6 s and turns off the cumulus convection parameterization scheme. Kessler warm cloud microphysical parameterization is used together with a 1.5-order turbulent kinetic energy subgrid parameterization. Open boundary conditions are utilized at the lateral boundaries while a high-level Rayleigh damping layer is introduced to reduce wave reflection from the top of the model.

Except for subtle differences in water mixing ratio above the troposphere, the sounding data used in this case is basically similar to those introduced by Weisman and Klemp (Weisman and Klemp, 1982) (using a surface potential temperature of 300 K and a surface water vapor mixing ratio of 14 g kg-1). For details on the morphology and evolution of supercell storms, readers can refer to Ray (Ray et al.,1981), Klemp et al. (Klemp and Rotunno, 1983), and Schenkman et al. (Schenkman et al., 2011).

The analysis domain shown in Fig. 2 consists of 41 ×41 × 25 grid points with a uniform grid interval of 1 km in the horizontal and 0.5 km in the vertical, in which a supercell is formed and developed. As shown in Fig. 3a, a mature supercell has been formed by 30 min into the simulation. A strong rotating updraft (with maximum vertical velocities exceeding 14 m s-1) near the center of the hook-echo pattern is evident at 2.5 km and is associated with a low-level maximum downdraft of 5.9 m s-1in the strong echo area.The vertical cross section (Fig. 3b) exhibits a gravitational oscillation pattern, implying that downstream of the overshooting updraft (with maximum vertical velocities reaching 35 m s-1at 6.5 km) at the tropopause level are downward returning flows (with maximum vertical velocitiesexceeding 8 m s-1). The structure of simulated supercell storm is similar to that described by Ray et al (Ray et al.,1975), Klemp et al (Klemp and Rotunno, 1983), and Wilhelmson et al (Wilhelmson and Klemp, 1981).

The model-simulated three-dimensional wind field at 30 min is sampled by two assumed radars located at (10, 0)and (80, 0) at ground level in Fig. 2 working in a scanning mode similar to VCP11 with a maximum range of 230 km,gate spacing of 1 km, azimuth resolution of 1°, and elevation angles of 0.5°, 1.0°, 1.5°, 2.0°, 2.5°, 3.0°, 4.0°, 5.0°,6.0°, 8.0°, 10.0°, 12.0°, 15.0°, and 18.0°, which are selected based on most operational radar configurations. The elapsed times for the volume scans of two pseudo-radars are neglected, and thus the radial winds are observed simultaneously. The simulated wind components are interpolated using a Cressman weighting function (Cressman, 1959)from the grid points to the sampling positions along the radar beams with horizontal and vertical influence radii of 2.4 km and 1.2 km, respectively, and then are synthesized to obtain radial velocities according to Eq. (1). In order to simulate the actual observational errors of radars, random errors with a standard deviation of 1 m s-1are added to these radial velocities [i.e.,is taken as the final observation in the retrieval procedure, where ε represents normally distributed random data with an expectation of 0 and a standard deviation of 1]. The wind retrieval is performed in the regions covered by both radars.

To give a quantitative assessment of the accuracy of the wind field derived from dual-Doppler radar data, we define the following indices between the retrieved variables and the true counterparts for verification. The mean absolute error (MAE) and correlation coefficient (CC) of wind field are defined as:

the root-mean-square error (RMS) and relative RMS error(RRE) of the horizontal velocity are defined as:and the RMS error and relative RMS error (RRE) of the vertical velocity are defined as:whereNis the total number of grid points,Fstands for any of the wind components,u,v, andw, in thex,y, andzdirections, and the subscripts cal and ref represent the retrieved variables and model-simulated variables, respectively.

4.2. Verification of the horizontal velocity

In this section, we first evaluate the quality of the retrieved horizontal wind field. It is clearly shown that the distribution of horizontal wind at 2.5 km in Fig. 3c is similar to the simulated field in Fig. 3a. By comparing the retrieved wind fields with the true ones from the WRF model, all the features in the horizontal wind field, including the cyclonic convergence and rotating updrafts around the supercell, are derived well, even though the retrieved wind components are a little stronger than the simulated ones. The vertical cross section (Fig. 3d) reflects the dynamic characteristics of the supercell well, with the retrieved maximum ascending motion of 34 m s-1, which is very close to the simulated one. The downstream oscillations near the top of the main updraft due to gravity waves are also evident in Fig. 3b. As a result, the retrieved three-dimensional wind field agrees well with the true one in Fig. 3. The quality of the retrieved winds can be further quantitatively exhibited by the MAE,CC, RMS, and RRE, comparing the retrieved winds with those from the model. The error statistics of the differences between retrieved and simulated horizontal wind components at each level are given in Table 1. The RMS_V,RRE_V, and MAE errors are all small while the CC is high below 6 km, where most mesoscale convective systems take place, because radars usually work at low elevation angles and the interpolation errors are relatively small with plenty of observation samples. At altitudes above 6 km, the RMS_V error first increases and then decreases with height.Although the RMS_V error is up to 2.5 m s-1at 9 km and winds at this altitude are rarely used, it is still within an acceptable range. Overall, the general features of the horizontal wind field at all levels are retrieved well, with relatively small RMS error (1.718 m s-1) and high CC (0.916). It is noted that due to radar detection, the coverage of the data decreases with height, which may affect the subsequent estimate of the vertical velocity.

4.3. Verification of the vertical velocity

4.3.1.Experiment design

To examine the reliability of the retrieved vertical velocity in more detail, we design the following five experiments listed in Table 2.

Table 1. Error statistics for the retrieved horizontal winds compared to the simulated winds. The bold values at the bottom of the table indicate the averages of each indice.

Table 2. List of experiments for vertical velocity estimation.

On the vertical velocity retrieval, experiments SO (SP)and DO (DP) are designed to compare the differences between the shallow and deep convective continuity equations with the O'Brien (Poisson) method, and experiments DP and DPE are used to evaluate the sensitivity of different air density distributions. Since the solution of the numerical model satisfies the basic equations of atmospheric motion,we take the vertical velocities from the WRF model as the true values for evaluation. Firstly, the simulated horizontal wind fields are used to estimate the vertical velocity, and the errors of each experiment are verified.

4.3.2.Vertical velocity derived from model-simulated horizontal wind

The advantage of using a simulated horizontal wind field is that the data at each grid point is valid, which can eliminate the uncertainties caused by insufficient data in the retrievals. Height profiles of the statistical deviation of the differences between the vertical velocities retrieved from the simulated horizontal wind components by each scheme and those obtained directly from the WRF model are shown in Fig. 4. On the whole, the errors of retrieved-simulated vertical velocity comparisons derived with the deep convective continuity equation (experiments DO, DP, and DPE) are reasonably smaller than those derived with the shallow convective continuity equation (experiments SO and SP). The RMS errors of experiments SO and SP are larger than the other three schemes at each layer in Fig. 4c. At altitudes below 11 km, the correlation coefficients in experiments DO, DP,and DPE remain relatively high, varying from 0.95 to 0.98(Fig. 4b), while the RMS error in Fig. 4c and the RRE in Fig. 4d remains smaller than 1 m s-1and 0.3 at each layer,respectively. As a result, the vertical wind field considering the density change with height by means of the deep convect-ive continuity equation is closer to the true state (i.e., the air density should not be ignored since the supercell storm develops in a deep and strong convective environment). In addition, there seems to be little difference in the vertical velocity error between the simulated air density experiment(DP) and the empirical air density experiment (DPE). The above results lay the experimental foundation for the threedimensional wind field derived from dual- or multiple-Doppler radar observations for a real case, as presented in the next section.

4.3.3.Vertical velocity obtained from retrieved horizontal wind

Invalid data inevitably exist at upper levels due to the noncontinuous radar elevation angles. This problem may have a certain impact on the wind field retrieval, especially on the estimation of the vertical velocity. As shown in Fig. 5,the errors of two experiments (SO and DO) performed with the O’Brien method increase faster with height than the other three experiments above 5.5 km. One of the main reasons is that the data at the highest level is used to adjust the vertical velocity from top to bottom in the O'Brien method[see Eq. (12)], thus it will be impossible to modify the data at a lower layer if the upper data is missing. The data coverages at 12 km and 12.5 km are only 86% and 64% respectively, which is common in practice, so the O’Brien method has certain defects. In contrast, the errors of the three experiments (SP, DP, and DPE) performed with the Poisson method remain relatively smaller and more stable. The RMS error reaches its maximum (about 2.5 m s-1) near 9 km and then decreases with height. The errors of experiments DP and DPE are smaller than those of experiment SP at each layer (Fig. 5), indicating that the vertical velocities derived with the deep convective continuity equation are more accurate than those derived with the shallow convective continuity equation. These features are basically consistent with the results outlined in the previous paragraph.Although the errors of experiment DPE at the upper level(＞7 km) are slightly larger than those of experiment DP, the differences are negligible. Approximate error distributions of experiments DP and DPE are also presented in Fig. 5,demonstrating that the true air density can be replaced by empirical air density. This conclusion has important practical significance, because the true air density of the atmosphere is difficult to obtain. Moreover, the Poisson method is rather stable and robust, as concluded from the results of experiments DP and DPE.

Fig. 4. Error distribution of vertical velocity obtained from model-simulated horizontal wind by scheme SO, DO, SP,DP, and DPE compared with simulated vertical velocity at each level. (a) Mean absolute error MAE (m s-1). (b)Correlation coefficient CC. (c) Root-mean-square error RMS_W (m s-1). (d) Relative root-mean-square error RRE_W.

Fig. 5. As in Fig. 4, but for vertical velocity derived from two virtual Doppler radars denoted in Fig. 2.

The mean error statistics of retrieved velocities at all levels for the above five experiments are listed in Table 3.As shown, both the RMS error (＜1.78 m s-1) and the RRE(0.609) remain reasonably small, and the CC is as high as 0.78 in experiments DP and DPE, indicating that the vertical wind is recovered with plausible accuracy. The above results are surprisingly similar to the conclusions of Gao et al.(Gao et al., 1999) (refer to experiment CNTL in the literature, with RMS error of 1.937 m s-1, RRE of 0.609, and CC of 0.825), proving that the method described in this paper is basically consistent with the variational solution with the mass continuity equation imposed as a weak constraint in the cost function (Gao et al., 1999, 2001) for the calculation of vertical velocity.

Table 3. List of experiments with mean error statistics of retrieved vertical wind compared with simulated wind.

5. Application to a real squall line case

In the previous section, we use a set of assumed dual-Doppler data from a model-simulated supercell storm to examine the performance of our three-dimensional variational method for both horizontal and vertical wind fields.In order to prove the effectiveness of the variational method for real radar data, we apply it to a squall line that developed in Jiangxi on 11 May 2017.

5.1. Data preprocessing

The locations of two radar stations are shown in Fig. 6,each with nine elevation angles (0.5°, 1.45°, 2.4°, 3.4°, 4.3°,6.0°, 9.9°, 14.6°, and 19.5°). The raw radar data needs to be preprocessed in advance of wind synthesis, which includes quality control for reflectivity data and radial velocity dealiasing (Zhang and Wang, 2006). Then, the Cressman weighting function is utilized to interpolate the observed wind components from the sampling positions along the radar beams to grid points with horizontal and vertical influence radii of 2.4 km and 1.2 km, respectively. Our technology of dual-Doppler radar wind synthesis is performed in the solid box shown in Fig. 6, which consists of 81 × 81 × 25 grid points with a uniform grid interval of 1 km in the horizontal and 0.5 km in the vertical.

5.2. Results of retrieved winds

In this section, we make comparisons with the VAD wind profile (VWP) and the observed radial velocity, respectively, to verify the reliability of the wind field retrieved from the two radars.

5.2.1.Comparison with VAD wind profile (VWP)

VWP is one of the most popular radar products used in operation, owing to its frequent intervals (updated every 6 min) and high resolution (up to 300 m) in the vertical direction. It can be seen from the hodographs of the two wind profiles shown in Fig. 7 that the wind speed retrieved from dual-Doppler radar observations is generally in accord with the VAD wind speed of the Nanchang radar, except for a small difference in the wind direction. Moreover, they both unanimously analyzed the vertical 1-3.5 km wind shear pointing to the northeast (the black and gray dashed lines in Fig. 7).

5.2.2.Comparison with the observed radial velocity

The accuracy of the retrieved three-dimensional wind fields can also be strictly verified with the observed radial velocities by synthesizing retrieved winds to the sampling positions along the radar beams according to Eq. (1). This comparison is rigorous, because any difference in value or location will cause a large deviation in mesoscale convective systems such as squall lines, supercell storms, etc. From a detailed point of view (Fig. 8), the retrieved radial velocities retain the observed positive and negative velocity centers and horizontal wind shear of the squall line, except that the value is slightly smaller than observed. Our variational dual-Doppler wind analysis plays an important role in filling data gaps in the retrieval due to the Cressman interpolation algorithm employed. Overall, the retrieved radial velocities match the observed velocities well.

Fig. 6. Locations of Nanchang and Fuzhou radar stations (solid five-pointed star) that observed the squall line on 11 May 2017 and the dual-Doppler analysis domain (solid box)with terrain heights (shaded and contoured every 250 m). The dotted circles denote the 150 km range of measurement while it is suitable to perform wind synthesis in the nonoverlapping domain of two solid cycles. The grid point marked with an asterisk represents the position 47 km away from the Nanchang radar. The red dots indicate the mosaic of composite radar reflectivity greater than or equal to 45 dBZ at 1242 UTC 11 May 2017.

Fig. 7. Hodographs showing the wind profile (black solid line)retrieved from dual-Doppler radar observations and the VAD wind profile (grey solid line) of the Nanchang radar at 1242 UTC 11 May 2017. The black dashed and gray dashed lines represent the vertical 1-3.5 km wind shear based on the retrieved wind profile and the VAD wind profile, respectively.The position of the retrieved wind profile is shown in Fig. 6.

Figure 9 shows the comparisons between the observed and retrieved radial velocity from all elevation angles of the Nanchang and Fuzhou radars. The RMS error is a little larger than that of the simulated tests (see Table 1) in section 4, which may be related to random errors from radar observations. In spite of that, the statistical MAE, RMS error, and RRE of the radial velocities from the two radars remain reasonably small. The observed and retrieved radial velocities from the two radars are concentrated in the range of -20 m s-1to -10 m s-1. From the least squares fitting formula of the scatters (Fig. 9), the slopes of the two lines are both ＜ 1 with intercept ＜ 0, indicating that the absolute value of the retrieved wind field is slightly smaller than the observed one, which is consistent with the conclusions outlined in the previous section. It is worth pointing out that the detection of vertical velocity at low elevation angles from Doppler radar is highly limited. In other words, the observed radial velocityVrfield does not contain enough information about the vertical movement of the atmosphere. Therefore, the vertical velocity obtained from the solution of Eqs. (7) and (8)is less accurate and should be re-estimated by using the Poisson method proposed above, although the horizontal velocity is reasonably accurate, as shown in Fig. 9. Furthermore,from an observation point of view, there is currently no effective detection equipment to directly measure vertical velocity over a large area, and verification of vertical velocity is challenging.

5.2.3.Dynamic verification

The results of our analysis for the squall line that occurred from 1224 UTC to 1300 UTC 11 May 2017 are illustrated in Fig. 10. The observed squall line moved towards the east-northeast at a nearly constant speed of 15 m s-1, calculated from the composite reflectivity of the squall line in two consecutive radar scans. The storm-relative winds are obtained by subtracting the mean squall line propagation speeds averaged during this period from the ground-relative winds. The convergent (divergent) centers in the horizontal divergence field correspond to the updrafts(downdrafts) shown in Figs. 10b, e, and h, indicating that the divergence field and the vertical velocity are resolved well. The vertical velocity increases rapidly at 1242 UTC(Fig. 10e) and substantially weakens at 1300 UTC (Fig. 10h).

Many studies have shown that the rear inflow is regarded as a main contributor to the bowing structure of a squall line (Fujita, 1979; Weisman, 1992, 1993; Meng et al.,2012). As shown in Figs. 10a, d, and g, a large bowing structure exists within the squall line with associated strong storm-relative rear inflow of 20-24 m s-1behind the leading edge of active convection at 2.5 km, following the conceptual model for the formation and evolution of squall lines proposed by Fujita (Fujita, 1978). A characteristic cyclonic vortex is resolved well and can be seen in the storm-relative wind field at the head of the squall line (Figs. 10a, d,and g). At 1224 UTC, the rear inflow grows larger and more intense (Fig. 10a), contributing to further enhancement and maintenance of the squall line. During the period with the strongest downdrafts (Fig. 10f), the echo (Fig. 10d) takes the shape of a spearhead pointing toward the heading of the squall line. The squall line starts to dissipate while the rear inflow and the vertical velocity gradually weaken over the next hour, after 1300 UTC (Figs. 10g-i). The vertical cross sections along line A-B in Figs. 10c, f, and i also indicate that there is a strong downward-rear inflow behind the bow echo below 5 km. This makes the convective storm tilt forward, forming a strong echo pendency and limitary weak echo region on the front side of the bow echo and strong vertical velocities up to 15 m s-1on the southern part of the bowing structure at its developing and formation stage. The accelerated outflow above 7 km causes strong echoes (greater than 45 dBZ) to reach a height of 6.5 km. The enhanced downdrafts behind the squall line accounting for an observed surface wind gust of 18.5 m s-1(not shown)together with the inflow from the front to the back of the squall line converge on the front side of the squall line, forcing the squall line to be maintained for a long time. These results suggest that the kinematic structures of this squall line are clearly illustrated, indicating that the retrieved horizontal and vertical winds are dynamically reasonable.

Fig. 8. The observed radial velocity (a, c, shaded, m s-1, positive away from and negative toward the radar) at 1.45°elevation angle from the Nanchang (a, b) and Fuzhou (c, d) radars and the retrieved radial velocity (b, d, shaded, m s-1)based on the derived wind from dual-Doppler radar data according to Eq. (1) at 1242 UTC 11 May 2017. The locations of the Nanchang and Fuzhou radars are shown in Fig. 6. The black arrow represents the heading of the squall line while the red ellipse indicates the positions of the positive and negative velocity centers and the horizontal wind shear. The negative values of x- and y-axis represent the west and south sides of the radar.

Fig. 9. Scattergraph of the observed-retrieved radial velocity comparisons from all elevation angles of the Nanchang(a) and Fuzhou (b) radars, with the concentration (0-1) represented by different colors of the scatters.

Fig. 10. Evolution of the squall line together with radar reflectivity, horizontal storm-relative wind field,horizontal divergence, and vertical velocity retrieved from dual-Doppler radar observations at (a-c) 1224 UTC 11 May, (d-f) 1242 UTC 11 May, and (g-i) 1300 UTC 11 May. (a), (d), (g) Horizontal storm-relative wind vectors(arrow, m s-1; the scale is in the bottom-right corner) and radar reflectivity (shaded, dBZ) at 2.5 km AGL. The blue solid lines describe horizontal storm-relative winds of 20 m s-1 and 24 m s-1. (b), (e), (h) Vertical velocity(shaded, m s-1; positive for upward and negative for downward) and horizontal divergence (contours, 10-4 s-1;dashed line for convergence, solid line for divergence) at 2.5 km AGL. (c), (f), (i) Vertical cross section of the synthesized wind field (storm-relative wind vectors projected onto the cross section, arrow, m s-1), vertical velocity (contours, solid line for upward and dashed line for downward; every 5 m s-1), and radar reflectivity(shaded, dBZ) along line A-B in (a), respectively. The thick gray solid line represents the updraft and downdraft.

Overall, the dynamic structure of the squall line reflected by the three-dimensional wind field recovered with reasonable accuracy from dual-Doppler radars in this study is in agreement with the well-documented conclusions of previous studies. Therefore, the retrieved wind field based on our three-dimensional variational method is fairly reliable.

6. Summary and concluding remarks

This paper proposes and evaluates a variational wind synthesis scheme which is capable of retrieving a three-dimensional wind field from dual-Doppler radar observations of convective storms. This scheme is performed in two steps.First, a cost function constructed by a radial observation term is minimized with the standard conjugate gradient method to directly obtain three-dimensional wind fields,avoiding a weighting parameter specification step. Then, the mass continuity equation is used as a strong constraint to reestimate the vertical velocity by solving the Poisson equation defined by the second derivative of the vertical velocity with impermeable conditions imposed at the ground and near the tropopause. It is shown that this solution has notable advantages in improving the accuracy of wind fields.The main conclusions can be summarized as follows.

(1) Compared with the true field from the WRF model,the horizontal wind field is recovered with reasonable accuracy based on our variational method and reflects the mesoscale dynamic structure of the simulated supercell storm well.

(2) The retrieved vertical wind field considering the air density change with height is much closer to the true state,which should not be ignored for strong convective storms.And, sensitivity experiments show that the true air density can be replaced by the empirical air density defined by Eq. (3).Furthermore, compared with the method proposed by O'Brien, this method outlined above is less sensitive to the uncertainty of the boundary conditions and therefore has better stability and reliability.

(3) The retrieved vertical wind profile during a severe squall line is consistent with the VAD wind profile, and the retrieved radial velocity together with the analyzed positive and negative velocity centers and horizontal wind shear of the squall line are in accord with radar observations.Moreover, the horizontal and vertical circulations within and around the squall line, as well as strong updrafts, the associated downdrafts, and associated rear inflow of the bow echo, are analyzed well, indicating that the dynamic structure of the squall line is reasonably retrieved.

It is worth mentioning that our variational method can be applied to synthesize the three-dimensional wind field from multiple-Doppler radar observations simultaneously.Despite the satisfactory results achieved in this study, there are still a lot of improvements worth making to the retrieval algorithm. For instance, will the accuracy of the wind field be greatly improved if three or more radars are used? How can we introduce additional wind data sources such as wind profilers, radiosondes, and surface station observations as well as additional dynamic constraints (e.g., the vertical vorticity equation, a Laplacian smoothing filter) into the cost function directly so that the wind field can be better resolved in common cases where low-level data are lacking? How can we develop a thermodynamic retrieval technique to derive the three-dimensional thermodynamic field from multiple Doppler radars, so as to better analyze the kinematic and thermodynamic structures of the squall line? Further research is being conducted to answer these questions and will be presented in a future paper.

Acknowledgements. This research was supported by the National Key Research and Development Program of China (Grant No. 2019YFC1510400), the National Natural Science Foundation of China (Grant Nos. 41975054 and 41930967), and the Special Fund for Forecasters of China Meteorological Administration(Grant No. CMAYBY2018-040). We would like to thank the anonymous reviewers for their insightful and detailed comments that appreciably improved the quality of this manuscript.