A 3 . 5-Dimensional Variational Method for Doppler Radar Data Assimilation and Its Application to Phased-Array Radar Observations

1 National Severe Storms Laboratory, Norman, OK 73072, USA 2 Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, Norman, OK 73072, USA 3 Global Modeling and Assimilation Office, NASA Goddard Space Flight Center, Science Application International Corporation, Greenbelt, MD 20771, USA 4 National Meteorological Center, China Meteorological Administration, Beijing 100081, China 5 Marine Meteorology Division, Naval Research Laboratory, Monterey, CA 93943-5502, USA


Introduction
Because the radar network provides only single-Doppler scanning over most areas in the U.S., research efforts have been undertaken to develop various methods for meteorological parameter retrievals from single-Doppler observations (Sun et al. [1]; Kapitza [2]; Qiu and Xu [3]; Sun and Crook [4]; Xu et al. [5,6]; Laroche and Zawadzki [7]; Shapiro et al [8]; Zhang and Gal-Chen [9]; Liou [10]; Gao et al. [11]; Weygandt et al. [12]).These previous efforts, however, were focused mainly on retrievals with zero background information.By utilizing the background information provided by model predictions, some of the previous retrieval methods can be upgraded for radar data assimilation (Sun and Crook [13]; Xu et al. [14]; Gu et al. [15] Gao et al. [16]; Hu et al. [17]).This paper presents a 3.5-dimensinal variational method for radar data assimilation developed by upgrading and combining the previous retrieval methods (Qiu and Xu [18]; Xu et al. [19]; Gal-Chen [20]; Hane and Scott [21]).This method uses simplified dynamical and thermodynamical equations of a numerical weather prediction (NWP) model (in this study we use the Coupled Ocean/Atmosphere Mesoscale Prediction System or COAMPS (COAMPS is a registered trademark of the Naval Research Laboratory.), Hodur [22]) as constraints while the background information is provided by the model forecasts.The concept and basic design of the method are described in the next section in connection with the general variational formulation derived from the estimation theory.The detailed formulations are presented for the three steps of the method

Estimate Theory and Variational Methods
Equations in an NWP model can be expressed symbolically in the following vector form: where x k is the state vector representing the discrete fields of the prognostic variables, f the forward model operator, q k the model error, and the subscript k denotes the kth time level of the model integration.A prior estimate of the state vector x 0 at the initial time (k = 0) is given by the model prediction from the previous data assimilation cycle.The state vector x 0 at the initial time (k = 0) can be expressed by where ( ) denotes the probabilistic mean (expectation) of ( ), and b the random part of x 0 .A prior estimate of x 0 is given by the model prediction from the previous data assimilation cycle.This background estimate is assumed to be unbiased.The observations can be expressed by where y m is the vector of observations, h the observation operator, and r m the observation error, and the subscript m denotes the mth time level of the observations.In general, we assume that the data assimilation period covers M observation time levels (from m = 1 to M) and K model time levels (from k = 1 to K), and the M observation time levels are a subset of the K model time levels.Assume that q k is a white Gaussian sequence and x k generated by (1) is Markov, which would be the case with linear dynamics.Assume also that b and r m are Gaussian random, unbiased (with zero mean), not cross-correlated and not correlated with q k .Then, according to Sherman and Bayes theorems (see Jazwinski [23,Chapter 5]), given observations y m (m = 1, 2, . . .M), the conditional maximum likelihood estimate of x k is also the conditional mean of x k (k = 0, 1, . . .K), and is given by the minimizer of the following cost function: where B = bb T , R m = r m r T m and Q k = q k q T k are the error covariance matrices for the background, observations, and model equations, respectively, m denotes the summation over m (from 1 to M), k denotes the summation over k (from 1 to K), ( ) T denotes the transpose of ( ), and λ k is the vector Lagrangian multiplier for the vector constraint given by (1).
The cost function in (4) presents a general formulation for four-dimensional variational (4dVar) data assimilation.The related methods of solution (Jazwinski [23]; Bennett [24]), however, are computationally too expensive for radar data assimilation.By setting q k = 0 in (4), the above 4dVar reduces to the perfect-model 4dVar (P4dVar), so the cost function minimum can be searched effectively by the adjoint method especially if B −1 and R −1  m are simply prescribed (Lewis and Derber [25]; Le Dimet and Talagrand [26]).Note that a volume scan from a single Doppler radar can contain up to 10 6 observations.To assimilate high-resolution radar observations over a mesoscale area, the model state vector x k often needs to contain at least as many as 10 6 gridded variables.In this case, the P4dVar is still too expensive and unpractical for real-time radar data assimilation.
A further simplification can be made by setting not only q k = 0 but also M = 1, K = 0 ( k = 0) and λ = 0 in (4).In this case, the observations are assimilated at a single time level without using any model equation constraint, so the 4dVar reduces to a three-dimensional variational (3dVar) formulation (Parrish and Derber [27]; Gao et al. [16]).Such a 3dVar is computationally efficient for real-time operational radar data assimilation, but it assimilates radar observed reflectivity and radial velocity (along the radar beam) at only a single time level.The tangential-velocity component (perpendicular to the radar beam) is not observed by the radar but often critical for model initialization and prediction.The same is true for the thermodynamic variables such as pressure and temperature which are not observed by the radar.These unobserved variables can be partially retrieved from four-dimensional radar observations (at two or more consecutive time levels) by using a P4dVar (Sun and Crook [4]) or a simple-model 4dVar (Xu et al. [14]).The required computations, however, are not practical for realtime applications unless the model size is sufficiently small (such as the one used in Xu et al. [14] with a 20×20×14 grid for a 8 × 8 × 2.8 km 3 domain).To solve the above problems, the simple-model 4dVar (Xu et al. [14]) needs to be further simplified to achieve the required computational efficiency but not overly reduced to a 3dVar, so that four-dimensional radar observations can still be used to extract information for the unobserved variables.Such simplifications are made in three steps based on the previous retrieval methods (Qiu and Xu [18]; Gal-Chen [20]; Hane and Scott [21]), which leads to the method presented in this paper.
As shown by the flowchart in Figure 1, the method assimilates four-dimensional radar observations (from three consecutive time levels of radar volume scans) in three steps.In step 1, raw level II radial velocity data (from a single or multiple radars) are quality-controlled and analyzed on the model grids by using a 3dVar.The cost function of this 3dVar is a simplification of (4) by setting M = 1 and K = 0 (with k = 0 and λ = 0) to retain only the background and observation terms (see ( 6)-( 10)).The above analysis is performed for each of three consecutive volume scans at t 0 − τ t 0 Time Step-1 Step-2 Step-3 Step 1 controls the quality of radial velocity observations v o r and produces radial velocity analyses v a r (on model grid) by minimizing J 1 in (6) at t 0 − τ, t 0 and t 0 + τ, respectively, for three consecutive radar volume scans.Step 2 produces the vector velocity increments v i by minimizing J 2 in (11) at t 0 − τ/2 and t 0 + τ/2, respectively.Step 3 produces the perturbation pressure increment π i by minimizing J π in ( 16) first, and then produces the perturbation potential temperature increment θ i by minimizing J θ in (18) at t 0 , and finally adjusts the water vapor mixing ratio q v and the rain (or snow) mixing ratio q r (or snow mixing ratio q s ) below (or above) the freezing level.their respective central time levels: t 0 − τ, t 0 and t 0 + τ (see Figure 1), where τ is the duration of one volume scan (about 1.5 minutes for the phased-array radar data used in this paper).The analyzed radial velocity fields are then used as gridded observations in step 2. The cost function used in step 2 is a simplification of (4) by setting M = K = 1 and λ = 0.This cost function retains the mass continuity equation and radial momentum equation (i.e., the radial component of the vector momentum equation) in q k in addition to the background and observation terms, but the two retained equations are simplified into incremental forms (see (13) and ( 15)).By fitting the radial momentum equation to the analyzed radial velocity observations (from step 1) at the time levels t 0 − τ and t 0 (or t 0 and t 0 + τ) in four-dimensional space, the vector velocity increment can be estimated at t 0 − τ/2 (or t 0 + τ/2).Once the vector velocity fields are updated at t 0 −τ/2 and t 0 + τ/2 in step 2, they are used as input data in step 3 to produce incremental analyses for the perturbation pressure and potential temperature at t 0 .The cost function in step 3 is a simplification of (4) by setting M = 0 (with m = 0), K = 1 and λ = 0, but only the two horizontal momentum equations are retained and simplified into incremental forms for the perturbation pressure analysis (see (16)) and only the vertical momentum equation is retained and simplified into an incremental form for the perturbation potential temperature analysis (see (17)).At the end of step 3, the water vapor and hydrometeor mixing ratios are also adjusted based on the radar observed reflectivity.In the above step 2 and step 3, the input data are fitted at two time levels in four-dimensional space, but the output incremental analyses are three-dimensional.In this sense, the method is called a three-and-half-dimensional variational (3.5dVar) method.The detailed formulations and computation procedures are presented for each step in the following three sections.

Data Quality Control and Radial
Velocity Analysis ).Project the interpolated v f onto the radial direction (along the radar beam) to obtain the background radial velocity v f r = r o • v f , where r o is the unit vector tangential to the radar beam at the observation point and is computed by considering the effects of the standard atmospheric refraction and Earth curvature (see Doviak and Zrnić [29, (2.28) and (9.9)]).Calculate the innovation v f r at each observation point, where v o r denotes the observation.
(ii) Take the integer n nearest to v d r /(2v N ) (where v N is the Nyquist velocity) and adjust the observational radial velocity to v o r − 2nv N so that the adjusted innovation (i.e., the adjusted observation minus the background value) is between ±v N .
(iii) If the absolute value of the adjusted innovation is less than 0.5v N , then the adjusted observation is accepted as a "good" one and used to replace the original v o r .Otherwise, perform the continuity check (buddy check) which is similar to the third step described in Gong et al. [28,Section 2d] but the threshold value is set to 0.25v N for the acceptable absolute difference between the adjusted observation and the averaged value of "good" neighborhood observations.For the second type of correction, the data quality control removes the error caused by the precipitation terminal velocity.The terminal velocity, w T , is estimated empirically (Kessler [30]) by where Z is the observed reflectivity in dBZ.The projection of the terminal velocity on the radial direction is w T sin α where α is the local tilt angle of the radar beam.Thus, with the above two types of correction, the originally observed v o r is replaced by v o r − 2nv N − w T sin α at each observation point where n / = 0 or w T / = 0.For simplicity, the corrected radial velocity observation will still be denoted by v o r in the remaining part of this paper.

Radial Velocity Analysis.
The radial velocity analysis is performed by minimizing the following cost function: where J bk is the background term and J ob the observation term.This cost function is a 3dVar formulation derived from (4) by setting M = 1 and K = 0, as explained in Section 2.
The control variable for J 1 is the analysis increment defined by where v a and v b denote the analysis and background velocities, respectively.Since the analyzed radial velocity fields obtained in step 1 will be used as observed "tracer" fields to retrieve the vector wind field in step 2 (see Figure 1), the analysis in step 1 should "interpolate" as accurate as possible the observed radial velocity fields onto the model grid.Because of this, the analysis is performed with zero background (v b = 0) but the background error covariance is nonzero in the background term J bk , and the radio between the background and observation error variance is set to σ 2 vr /σ 2 ob = 10, where σ 2 vr and σ 2 ob denote the radial velocity background and observation error variances, respectively.
To facilitate the formulation and computation of the background term, v a or, equivalently, v i (since v b = 0 in step 1) is treated as an intermediate control variable with its two horizontal components (u i , v i ) expressed in terms of stream function and velocity potential by The grid field of (u i , v i ) can be computed from the grid field of (ψ i , χ i ) by using the standard central finite difference scheme (in the horizontal).Denote by (ψ i , χ i , w i ) the state vectors of the grid fields of (ψ i , χ i , w i ).These three state vectors can be related to three new intermediate vectors (ψ, χ, w) by where B 1 , B 2 , and B 3 are the univariate error covariance matrices for the background fields ψ f , χ f , and w f , respectively.The background term is then given by Here, it is assumed that the background errors are not crosscorrelated between ψ f , χ f , and w f .With this assumption, ψ i , χ i , and w i are decoupled in the background term as shown by (9), so the computation is simplified.The background error covariances for (ψ b , χ b , w b ) are modeled by horizontally isotropic Gaussian functions.The decorrelation length and depth are the parameters that need to be predefined empirically.Since the analyzed radial velocity fields obtained in this step will be used as observed "tracer" fields to retrieve the vector wind field in step 2, a relatively small horizontal decorrelation length should be used in step 1 to prevent the analyzed "tracer" field from becoming overly smoothed.On the other hand, the decorrelation length should not be too small to filter spurious 2Δx waves.Based on these considerations, we set the horizontal decorrelation length to L ψ = L χ = 3 √ 2Δx (which is 1/2 of that used in step 2 as shown Section 4) and the vertical decorrelation depth to and D w = Δz for w b according to (A.8) and (A.9), where Δx (=6 km in this study) is the horizontal grid spacing and Δz is the vertical grid spacing (from 0.02 to 0.25 km in the boundary layer, and about 1 km in the middle troposphere).The error variance for (ψ b , χ f ) and w b is then related to σ 2 vr by , respectively, (see (A. 19) and (A.20)).The covariance matrices in (8) are too large to compute explicitly, but their represented linear transformations in ( 8) can be simulated efficiently by either recursive filters (Purser et al. [31]) or differential operators (Xu [32]).Each horizontal field of (ψ, χ, w) is further expressed by an expansion of quadratic B-spline basis functions on coarse finite element meshes (Xu et al. [19]), so the final control variables are the B-spline coefficients of (ψ, χ, w).
The observation term has the following general form (for step 1 and step 2) where {{ }} denotes the summation over all observation points, is the incremental radial velocity at each observational point, and r o is the unit vector tangential to the radar beam at the observation point (computed by considering the effects of the standard atmospheric refraction and Earth curvature as in Doviak and Zrnić [29, (2.28) and (9.9)]).As explained earlier, since In (10), the observation errors are assumed to be not correlated between different observation points, so R m (with m = M = 1) reduces to Iσ 2 ob in (4).This assumption should be a valid approximation because radial velocity observation errors are correlated only between neighboring observation points (Xu et al. [33]).By setting σ 2 vr /σ 2 ob = 10 as explained earlier, there is no need to estimate σ 2 ob and σ 2 vr separately for the analysis in step 1.As in (9), v i and v i r in (10) are also converted to (ψ, χ, w) by using (7) and (8).Each horizontal field of (ψ, χ, w) is then further expressed by an expansion of quadratic B-spline basis functions on coarse finite element meshes (Xu et al. [19]), so the final control variables are the B-spline coefficients of (ψ, χ, w).The B-spline representation enhances the filtering and reduces the dimension of the analysis space.By using the standard conjugate gradient algorithm, the cost function in ( 6) is minimized efficiently in the reduced space spanned by the final control variables-the B-spline coefficients.As v i is linearly related to the B-spline coefficients, v i is readily estimated once the cost function is minimized.Since only one time level of radial velocity observations is used, the estimated v i is not accurate, but its radial component r in step 1) should be accurate and can be used as gridded observations for the vector velocity analysis in step 2 in the next section.

Vector Velocity Analysis
The vector velocity analysis is performed by minimizing the following cost function: This cost function is derived from ( 4) by setting M = K = 1 and λ = 0 to retain only the mass continuity equation and the radial momentum equation in q k , so the equation term in (4) reduces to the sum of the mass continuity equation term J ms and radial momentum equation term J rm in (11).The background term in (11) has the same form as in (9) except that the background velocity is no longer set to zero (as in step 1) and is given by v b = v f , where v f is the model forecast.The error covariance functions for the model background winds were estimated from a time series of the radar innovations for the same case considered in the paper.As shown in Xu et al. [34,Figure 3], the main part of the error correlation function for the horizontal background winds can be modeled approximately by a Gaussian function over the range of 40 km, but the tail part oscillates around zero and reaches a secondary peak (≈0.1) at the separation distance of 125 km.The decorrelation length estimated from this error correlation function was L = 43 km, but this estimated value was affected by the second peak in the tail part.For the main Gaussian part, the estimated decorrelation length should be reduced at least to 36 km (= 6Δx).Thus, the decorrelation lengths for (ψ f , χ f ) and w f can be estimated by ) and (A.10)).The analysis time level for J bk in ( 11) is at t 0 − τ/2 (or t 0 + τ/2)-the middle time level between t 0 − τ and t 0 (or t 0 and t 0 + τ), so the incremental field v i (= v a − v f ) and background velocity field v f are also all at t 0 − τ/2 (or t 0 + τ/2).
The observation term J ob for J 2 in (11) has the same form as that in (10) for J 1 , but the background velocity is given by v b = v f (instead of v b = 0 as in step 1) and the observations are binned from the two time levels at t 0 − τ and t 0 (or t 0 and t 0 + τ) to the middle-time level at t 0 − τ/2 (or t 0 + τ/2).Since the background and observation errors are uncorrelated (Xu and Wei [35]), the sum of their error variances is given by the innovation variance σ d − σ 2 ob with σ ob given by a pre-estimated value.For the phased-array radar data used in this paper, σ ob = 2 ms −1 is pre-estimated according to the statistical analysis of the phased-array radar innovation data (Xu et al. [34]).
The third term in (11) is the mass continuity equation term given by where [[ ]] denotes summation over all grid points, E ms is the mass continuity equation constraint for v i , and σ 2 ms is the equation error variance or, specifically, the variance of the residual error of E ms .The equation errors are assumed to be not correlated between different grid points and different equations.In COAMPS, the mass continuity equation is embedded in the pressure tendency equation (see Hodur [22, (5) and ( 20)]).From the COAMPS pressure tendency equation, the extracted mass continuity equation has the following incremental form: where , H and z s are the heights of the model top and bottom terrain boundaries, respectively, Θ v the basic state virtual potential temperature, and ρ the basic state density.Here, the equation errors are assumed to be not spatially correlated, so Q k (with k = K = 1) reduces to Iσ 2 ms and the equation term in (4) reduces to J ms .The first two terms on the right-hand side of (13) are the horizontal divergence and the associated RMS error can be estimated by σ div /L div = √ 2σ vr /L χ (see (A. 16) and (A.18)).The last term in (13) is smaller than the first two terms, and its RMS error should be relatively small.Based on this consideration, σ ms = 2σ vr /L χ is estimated.
The last term in (11) is given by where E rm is the radial momentum equation in an incremental form, and σ 2 rm is the variance of the residual error of E rm .The radial momentum equation is obtained by projecting the model vector momentum equation (Hodur [22, (2)-( 4)]) onto the radial direction along the radar beam, and then converted into the following incremental form: where v a = v f + v i , f is the Coriolis parameter, r is the radial distance from the radar, ∇ = (∂ x , ∂ y , ∂ σ ) is the gradient operator in the σ-coordinates, D is the diffusion operator Advances in Meteorology for subgrid-scale turbulent mixing as in COAMPS, and k is the unit vector in the vertical direction.Note that the thermodynamic variables are not updated yet, so they have zero increment at this step, as assumed in (15).The time derivative term ∂ t v i r in (15) is discretized by the standard central finite difference scheme with v a r given by the gridded observations (obtained from step 1) at t 0 − τ and t 0 (or t 0 and t 0 + τ).All the remaining terms in (11) are at the middle time level t 0 − τ/2 (or t 0 + τ/2), so v a r is interpolated to the middle time level in these terms.A further simplified version of (15) was used by Qiu and Xu [18] as a weak constraint with zero background velocity to retrieve lowaltitude winds from single-Doppler radar scans.The radial momentum equation in (15) and cost function in (11) are extensions and improvements relative to their counterparts used in ( 1) and ( 2) of Qiu and Xu [18].
The radial momentum equation in ( 15) is formulated at the middle time level and the equation errors are assumed to be spatially uncorrelated, so the associated equation term in (4) reduces to J rm and Q k reduces to Iσ 2 rm , as shown in (14).The second and third terms on the right-hand side of ( 14) are the horizontal advection increment, and the associated RMS error can be roughly estimated by σ vr σ rot /L rot +σ vr σ div /L div = 2 √ 2σ 2  vr /L χ (see (A.15)-(A.16) and (A.18)).The remaining terms in (4) are relatively small and so are their RMS errors.Based on this consideration, vr /L χ is estimated.As in step 1, by using ( 7)-( 8) and the B-spline representation, v i and v i r in (15) are all expressed as linear functions of the B-spline coefficients.With the use of the standard conjugate gradient algorithm, the cost function J 2 in ( 11) is minimized efficiently in the reduced space spanned by the Bspline coefficients-the final control variables.The estimated B-spline coefficients are then transformed back to (ψ, χ, w) and finally to v i .The total vector velocity analysis is given by v a = v f + v i on the model grid.

Thermodynamic Analysis
5.1.Perturbation Pressure Analysis.The estimated fields of v i at t 0 − τ/2 and t 0 + τ/2 from step 2 are used as input data to estimate the thermodynamic fields at t 0 in step 3. Note that the Exner function π = (p / p 00 ) R/Cp is used in place of the perturbation pressure p in COAMPS, where p 00 is a constant reference pressure, R is the gas constant for dry air, and C p is the specific heat at constant pressure for the atmosphere.In this paper, we simply call π "perturbation pressure".The perturbation pressure increment π i is estimated by minimizing the following cost function: Here, |π| 2 /σ 2 π = π iT C −1 π π i /σ 2 π is the background term and π(= C −1/2 π π i ) is the control variable.In this term, π i is the state vector of the grid field of π i , C π is the error correlation matrix for the background perturbation pressure π f , σ 2 π is the error variance for π f , and thus σ 2 π C π is the error covariance matrix for π f .The second term on the right-hand side of ( 16) is the equation term.In this term, E u and E v are the two horizontal momentum equations retained from COAMPS (see Hodur [22, (2)-(3)]) in the following incremental forms: and σ 2 uv is the equation error variance for E u and E v .Here, ∂ t u i and ∂ t v i are computed by the standard central finite difference scheme from the grid fields of (u i , v i ) at t 0 − τ/2 and t 0 + τ/2 produced in step 2. All the remaining terms in (17) are at the middle time level t 0 .
As explained in Section 2, the cost function in ( 16) is derived from ( 4) by setting M = 0 ( m = 0) and K = 1 with λ = 0, retaining only E u and E v in q k , and reducing Q k to Iσ 2 uv .The two error variances σ 2 π and σ 2 uv in ( 16) are difficult to estimate.Their ratio, however, can be tuned adaptively to make the cost function J π approximately equally partitioned by the two terms on the right-hand side of ( 16) when J π reaches the minimum.This approach is used in this paper.In (16), the grid field of π(= C −1/2 π π i ) is further expressed by a B-spline expansion on coarse finite element meshes (Xu et al. [19]), so the final control variables are the Bspline coefficients of π and J π is minimized efficiently in the reduced B-spline space.The estimated B-spline coefficients are then transformed back to π and finally to π i (= C 1/2 π π).The transformation represented by C 1/2 π is simulated by a recursive filter (Purser et al. [31]), in which the error covariance for π f is modeled by a horizontally isotropic Gaussian function.In this filter, the horizontal decorrelation length is set to L ψ as for (ψ f , χ f ) in (A.3) and (A.10), while the vertical decorrelation depth is set to D w as for w f in (A.4).

Perturbation Temperature Analysis.
After π i is computed, the perturbation potential temperature increment θ i is estimated by minimizing the following cost function: Here, is the background term and θ(= C −1/2 θ θ i ) is the control variable.In this term, θ i is the state vector of the grid field of θ i , C θ is the error correlation matrix for the background perturbation potential temperature θ f , σ 2 θ is the error variance for θ f , and thus σ 2 θ C θ is the error covariance matrix for θ f .The second term on the right-hand side of ( 18) is the equation term.In this term, Ew is the vertical momentum equation retained from COAMPS in the following incremental form: Advances in Meteorology 7 and σ 2 Ew is the equation error variance for E w .In (19), Θ is the basic state potential temperature and gθ i /Θ is the buoyancy increment caused by the perturbation potential temperature increment.Note that the mixing ratios of water vapor and hydrometeors are not updated yet, so they have zero increment at this step and thus no contribution to the buoyancy increment, as assumed in (16) (see Hodur [22, (4)]).Here, ∂ t w i are computed by the standard central finite difference scheme from the grid fields of w i at t 0 − τ/2 and t 0 + τ/2 produced in step 2. All the remaining terms in (19) are at the middle time level t 0 .
The cost function in ( 18) is derived from ( 4) by setting M = 0 ( m = 0) and K = 1 with λ = 0, retaining only E w in q k , and reducing Q k to Iσ 2  Ew .Again, as those in ( 16), the two error variances σ 2 θ and σ 2 Ew in ( 18) are hard to estimate, but their ratio can be tuned adaptively to make the cost function J θ approximately equally partitioned by the two terms on the right-hand side of (18) when J θ reaches the minimum.In (18), the grid field of θ is further expressed by a B-spline expansion on coarse finite element meshes similarly to that for π in (16).The final control variables are the B-spline coefficients of θ, so J θ can be minimized efficiently in the reduced B-spline space.The estimated B-spline coefficients are then transformed back to θ and finally to θ i (= C 1/2 θ θ).The transformation represented by C 1/2 θ is simulated by a recursive filter with the decorrelation length and depth set to be the same as those for (ψ f , χ f ).
The analyzed total perturbation pressure and potential temperature are given by π f + π i and θ f + θ i , respectively, on model grid.The basic idea for the above thermodynamic analysis is the same as Gal-Chen [20] and Hane and Scott [21], but the formulation is derived by considering not only the equation constraints but also the background constraint in connection with the general variational formulation in (4).

Moisture and Hydrometeor Adjustments.
After the perturbation pressure and potential temperature fields are updated by the above incremental analyses, the relative humidity is altered and thus needs to be recovered by adjusting the water vapor mixing ratio q v .After this, q v is further adjusted based on the radar observed reflectivity through the following three substeps: (i) Interpolate the radar observed reflectivity η o onto the model grid by minimizing the following cost function: where η = B −1/2 η η a , B η is the error covariance matrix for zero reflectivity background, η a is the state vector of the analyzed reflectivity η a , and σ ηo (≈1 dBZ) is the reflectivity observation error standard deviation.Again, the grid field of η is further expressed by a Bspline expansion on coarse finite element meshes.As the final control variables are the B-spline coefficients of η, J π can be minimized efficiently in the reduced B-spline space.The estimated B-spline coefficients are then transformed back to η and finally to η a (= B 1/2 η η).The transformation represented by B 1/2  η is simulated by a recursive filter, in which the decorrelation length and depth are set to be as small as Δx and Δz, respectively, and the error variance for zero reflectivity background is estimated by the spatial variance of the model predicted reflectivity η f .With the above parameter settings, the analyzed reflectivity η a obtained by minimization J η in ( 20) is an optimal interpolation of the observed reflectivity η o onto the model grid.After this interpolation, η a is further extrapolated downward from the lower radar beam height to the surface level (to fill the reflectivity data void area below the lowest beam in the boundary layer).
(ii) Check each grid point for the conditions of η a > 10 dBZ and η f < 10 dBZ.If these two conditions are both satisfied, then adjust q v to the saturated value if the analyzed vertical velocity is nonnegative (w ≥ 0) or to 80% of the saturated value if the analyzed vertical velocity is negative (w < 0) and the relative humidity is below 80% at this grid point.
(iii) Check each grid point for the conditions of η a < 10 dBZ and η f >10 dBZ.If these two conditions are both satisfied, then adjust q v to the value interpolated in xand y-directions from the nearest grid points where the two conditions are not both satisfied.This adjustment is designed to reduce q v and thus correct the model's tendency of overpredicting precipitation locally in areas detected by the above two conditions.
In companion with the q v adjustment, the rain mixing ratio q r (or snow mixing ratio q s ) is adjusted below (or above) the freezing level to make the computed reflectivity match η a at each grid point.The above moisture and hydrometeor adjustments have not been extensively tested and currently are used only as an option for case studies (including the example presented in the next section).For operational applications, the radar reflectivity observations are used together with the GOES satellite observations to adjust the moisture field together with the cloud field in a sophisticated cloud analysis package (Wei et al. [36], Zhao et al. [37,38]).

Application to Phased Array Radar Data
The 3.5dVar is applied to the phased array radar radial velocity and reflectivity observations collected during the period of 2100-2200 UTC 2 June 2004 when a four-quadrant electronic scan strategy was tested.During this period, a squall line moved southeastward through the central Oklahoma area in the radial range (140 km) of the phased array radar scans (Figure 2).The radar scanned roughly every two minutes per volume.Total 26 volume scans were collected.Among these 26 volume scans, there is one volume scan that covers only a single quadrant and this volume is not used.The remaining 25 volume scans cover all the four quadrants and are used in this study.Each volume scan has Advances in Meteorology  7 tilts with elevation angles of 0.75, 2.27, 3.78, 5.28, 6.78, 8.28, and 9.28 degree.On each tilt, the spatial resolutions are 240 m in the radial direction and approximately 1.5 • in the azimuthal direction.
The COAMPS is used to produce the background fields.The model is configured with three nested domains centered over the state of Oklahoma with resolutions of 54, 18, and 6 km for the coarse, medium, and fine grids, respectively, and 30 levels in the vertical.The three nested domains are shown in Figure 3.All other parameters are set to be the same as in Zhao et al. [38].For the control run, the model is initialized (cold start) at 0000 UTC 2 June 2004.After the first 12-hr model run, the conventional observations are assimilated, and then another 12-hr run is launched (warm start).The predicted wind fields on the 6 km grid are used as the reference by the dealiasing technique described in Section 3.1 to detect and correct alias errors in the radar radial velocity observations.The technique is found to be effective (as shown Figure 2) and better than the operationally used technique.The radial velocity innovations produced from the 25 volume scans of dealiased radial velocities and COAMPS predicted background fields were used in Xu et al. [34] to estimate the background and observation error variances.The estimated observation error variance is used in this paper to estimate the background error variance adaptively in each assimilation cycle, as explained in Section 3.2.
By using the 3.5dVar, three consecutive volume scans of the dealiased radial velocity and reflectivity data are assimilated through a single cycle around 2108 UTC, and then a test run, called Test 1, is launched for 58-minute forecast to 2200 UTC.Another test run, called Test 3, is also performed by assimilating the first nine volume scans of the dealiased radial velocity and reflectivity data in three cycles from 2108 to 2118 UTC.The velocity and reflectivity fields at z = 3.1 km produced by the prediction valid at 2108 UTC in the control run and by the analysis at 2108 UTC in Test 1 are plotted against the observed reflectivity (analyzed onto the grid) at 2108 UTC in Figures 4(a) and 4(b), respectively.The predicted velocity and reflectivity fields in Figure 4(a) are the background fields for the analysis at 2108 UTC in Test 1.As shown in Figure 4(a), the background wind field from the control run is dynamically self-consistent and is consistent with the reflectivity field of its predicted squall line, but the predicted reflectivity reveals a significant location error for the predicted squall line in comparison with the observed reflectivity field.In particular, the leading edge of the predicted squall line is lagged by about 30 km behind the observed.This location error is largely corrected by the analysis in Test 1 as shown in Figure 4    Figure 6(a) plots the incremental horizontal velocity and potential temperature (at z = 0.5 km) produced by the 10minute prediction from the analysis in Test 1 with respect to those produced by the control run valid at 2118 UTC.As shown, the potential temperature increment is negative and the horizontal wind increment is divergent over an elongated region along the leading edge of the squall line produced in Test 1.Note that this elongated region is to the south immediately ahead of the background squall line (produced by the control run valid at 2118 UTC), which implies that the under-predicted southward movement of the squall line in the control run is largely corrected at 2118 UTC in Test 1.The above features become slightly more distinct in the incremental fields produced at 2118 UTC by the three assimilation cycles in Test 3 as shown in Figure 6(b), which implies that the under-predicted southward movement of the squall line in the control run is further corrected at 2118 UTC after the second cycle and the third analysis in Test 3.
It is necessary to mention that the incremental potential temperature produced by the thermodynamic analysis is relatively small (within ±0.5 • K in the boundary layer) and nearly hydrostatically balanced with the perturbation pressure increment, while the latter is related to the horizontal wind increment through the horizontal momentum equations (see (17)).Thus, the potential temperature increment in either Figure 6(a) or Figure 6(b) is caused mainly by the moisture adjustment and subsequent model integration.Due to the very limited coverage of radial velocities from only one Doppler radar, the thermodynamic analysis is not well constrained and thus does not have the desired accuracy and significance to directly and adequately correct the background thermodynamic field.Instead, the main benefit of the thermodynamic analysis performed in this paper is to improve the overall dynamic and thermodynamic balance of the analyzed total field, and this explains why the thermodynamic analysis improves the subsequent forecast slightly in comparison with the forecast without thermodynamic analysis (not shown in this paper).
The velocity and reflectivity forecasts from the control run and Test 1 valid at 2200 UTC are plotted against the observed reflectivity at z = 3.1 km in Figures 7(a) and 7(b), respectively.As shown, the reflectivity field predicted by Test 1 is much closer to the observations than that from the control run, but slightly less close to the observations than that predicted by Test 3 (not shown).The spatial correlation coefficients between the predicted and observed reflectivity fields are plotted as functions of time for the control run and the two test runs in Figure 8(a).As shown by the plotted correlation coefficients, the predicted reflectivity and precipitation are significantly improved in Test 1 and further improved in Test 3.
To verify the three-dimensional winds, the modelproduced (analyzed and predicted) radial velocity fields are interpolated (using the observational operator) back to the phased-array radar observation points and compared with their respective observed values.The RMS differences between the predicted and observed radial velocities are plotted as functions of time for the control run and the two test runs in Figure 8(b).As one can see, the RMS difference is reduced by the analysis in each assimilation cycle, although the reductions by the analyses in the second and third assimilation cycles are relatively small and last only for about 0.5 hour.Clearly, the predicted reflectivity (associated with predicted precipitation) and velocity are both improved in Test 1 and both further improved in Test 3.

Conclusions
This paper presents a variational approach for Doppler radar data assimilation.This method analyzes four-dimensional radar observations (three consecutive volume scans) but updates the model state only in three-dimensional space at the central time level in each assimilation cycle, and therefore we call it three-and-half-dimensional variational (3.5dVar) method.The method can be considered as an upgraded combination of the previous retrieval methods (Qiu and Xu [18]; Gal-Chen [20]; Hane and Scott [21]) with the background information from NWP model forecasts.The finite element B-spline representations (Xu et al. [19]) and recursive filter (Purser et al. [31]) are used to reduce the dimension of the analysis space and enhance the computational efficiency.The method is tested with real radar data collected by the phasedarray radar at Norman, Oklahoma for the squall line event on June 2, 2004.The effectiveness of the 3.5dVar is demonstrated by the improved short-term predictions of the squall line (see Section 6).
The 3.5dVar method has also been tested and applied to many other cases, either as a stand-alone package (Gu et al. [15]; Xu et al. [34,39]; Zhao et al. [40]) or in combination with a cloud analysis package (Zhao et al. [37,38]).The 3.5dVar method is currently being further refined in combination with a cloud analysis package for the Navy's nowcasting applications.According to recent realtime tests, the current 3.5dVar algorithm code can be easily setup and run on a PC workstation.It takes only about 5 minutes to assimilate 9 volume scans from three radars (each radar provides three consecutive volume scans) on a high-resolution (100 × 100 × 30) COAMPS grid that covers a mesoscale area (600 × 600 km 2 ).The 3.5dVar is thus sufficiently efficient for real-time applications.
Because simplifications are made to the general variational formulation in each of the three steps in the method to achieve the desired computational efficiency for real-time applications, the resulting analyses are suboptimal (relative to the desired maximum likelihood estimates).This implies that the method can be improved by reducing the involved simplifications and refining the error covariance estimation and representation.Such an improvement is made in the current 3.5dVar over the early version (Gu et al. [15]) by recovering the matrix forms of background error covariances.Room still exists for further improvements.For example, including the vorticity equation as an additional constraint could substantially improve the thermodynamic analysis (at least in areas covered by two or more Doppler radars, as shown by Protat and Zawadzki [41]).The analysis may also be improved by considering nonisotropic flowdependent forms of background error covariances, which is under our current investigation in connection with the ensemble Kalman filter.

Appendix Background Error Variances and Decorrelation Lengths for Derived Variables
The mass continuity equation in (13) can be applied approximately to the background velocity v f .By neglecting the effect of terrain slope, this equation can be written into the following form: where Δ h = ∂ 2 x + ∂ 2 y is the Laplacian.Assume that the random error fields of χ f and w f are also constrained by (A.1), so their covariance functions (see Xu and Wei [35, (2.3) and (A.1)]) satisfy the following relationship: where r and η are the horizontal and vertical distances, respectively, between the two correlated points, and the two covariance functions are assumed to be horizontally isotropic and separable between the horizontal and vertical directions.Assume that C χχ (r, η) can be modeled by a Gaussian function in the horizontal direction, so Note that R χχ (η) and R ww (r) are obtained explicitly in (A.5) and (A.6), respectively.Substituting (A.5) into (A.3)gives the explicit form of C χχ (r, η), and substituting (A.6) into (A.4) results in another explicit form for C ww (r, η).Using these explicit forms, we can verify that the vertical decorrelation depth for the random error field of χ f is given by , (A.8) while the horizontal decorrelation length for the random error field of w f is given by Denote by v h the horizontal component of v f and represent by v rot and v div the rotational and divergent parts of v h , respectively.Then we have v rot = k × ∇ h ψ f , v div = ∇ h χ f , and v h = v rot + v div .Also assume that the error covariance function for ψ f can be modeled by a Gaussian function in the horizontal, so where σ 2 ψ and R ψψ (η) denote, respectively, the error variance and vertical error correlation for ψ f , and L ψ is the horizontal decorrelation length.C rot (r, η) and C div (r, η) represent the traces of the error covariance tensors for v rot and v div , respectively.By using (2.14) of Xu and Wei [35] with (A.4) and (A.10), one can get

Figure 1 :
Figure 1: Flowchart for the three steps of the method.Step 1 controls the quality of radial velocity observations v o r and produces radial velocity analyses v a r (on model grid) by minimizing J 1 in (6) at t 0 − τ, t 0 and t 0 + τ, respectively, for three consecutive radar volume scans.Step 2 produces the vector velocity increments v i by minimizing J 2 in (11) at t 0 − τ/2 and t 0 + τ/2, respectively.Step 3 produces the perturbation pressure increment π i by minimizing J π in (16) first, and then produces the perturbation potential temperature increment θ i by minimizing J θ in (18) at t 0 , and finally adjusts the water vapor mixing ratio q v and the rain (or snow) mixing ratio q r (or snow mixing ratio q s ) below (or above) the freezing level.
2 vr + σ 2 ob = σ 2 d , where σ 2 d denotes the radial velocity innovation variance, σ 2 vr and σ 2 ob denote the radial velocity background and observation error variances, respectively.Since zero background is used in step 1, the spatial variance of the innovation (v o r minus v f r in the observation space) is used to estimate σ 2 d in each assimilation cycle.The background error variance is estimated by σ 2 vr = σ 2
(b).The 5-minute forecasts of velocity and reflectivity fields valid at 2118 UTC from the second cycle (at 2113 UTC) in Test 3 are given Figure 5(a) while the analyses from the

Figure 4 :Figure 5 :
Figure4: Velocity (arrows for the horizontal wind and green contours for the vertical velocity) and reflectivity (color field) at z = 3.1 km produced by (a) the prediction in the control run and (b) the analysis in Test 1 at 2108 UTC.The black contours are plotted every 10 dBZ for the observed reflectivity (analyzed onto the grid).The vertical velocity contours are plotted in green for ±0.5, ±2, and ±5 ms −1 .The arrow (=30 m s −1 ) at the lower-left corner is the vector scale for the horizontal velocity.

Figure 6 :
Figure 6: As in Figure 4 but for the incremental horizontal velocity (arrows) and incremental potential temperature (contours every 1 degree, dashed for negative values) in the boundary layer (at z = 0.5 km) from (a) the 10-minute prediction (integrated from the analysis at 2108 UTC) in Test 1 and (b) the third analysis in Test 3 with respect to those produced by the control run valid at 2118 UTC.

Figure 7 :
Figure 7: As in Figure 4 but for velocity and reflectivity produced by the predictions valid at 2200 UTC in (a) the control run and (b) Test 1.

Figure 8 :
Figure 8: (a)  Correlation coefficient between the predicted and observed reflectivity fields.(b) RMS differences between the predicted and observed radial velocities.The red, blue, and green curves are for the results obtained from the control run, Test 1 and Test 3, respectively.