# A 3.5-Dimensional Variational Method for Doppler Radar Data Assimilation and Its Application to Phased-Array Radar Observations

Authors: Qin Xu; Li Wei; Wei Gu; Jiandong Gong; Qingyun Zhao
Journal: Advances in Meteorology
 (2010)

---

## Abstract

A 3.5-dimensional variational method is developed for Doppler radar data assimilation. In this method, incremental analyses are performed in three steps to update the model state upon the background state provided by the model prediction. First, radar radial-velocity observations from three consecutive volume scans are analyzed on the model grid. The analyzed radial-velocity fields are then used in step 2 to produce incremental analyses for the vector velocity fields at two time levels between the three volume scans. The analyzed vector velocity fields are used in step 3 to produce incremental analyses for the thermodynamic fields at the central time level accompanied by the adjustments in water vapor and hydrometeor mixing ratios based on radar reflectivity observations. The finite element B-spline representations and recursive filter are used to reduce the dimension of the analysis space and enhance the computational efficiency. The method is applied to a squall line case observed by the phased-array radar with rapid volume scans at the National Weather Radar Testbed and is shown to be effective in assimilating the phased-array radar observations and improve the prediction of the subsequent evolution of the squall line.


---

## Body

## 1. 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 in Sections 3–5. The method is applied to phased-array radar observations for a squall line case in Section 6. Conclusions follow in Section 7.


## 2. Estimate Theory and Variational Methods

Equations in an NWP model can be expressed symbolically in the following vector form:(1)xk-f(xk-1)=qk,
where xk is the state vector representing the discrete fields of the prognostic variables, f the forward model operator, qk the model error, and the subscript k denotes the kth time level of the model integration. A prior estimate of the state vector x0 at the initial time (k = 0) is given by the model prediction from the previous data assimilation cycle. The state vector x0 at the initial time (k = 0) can be expressed by(2)x0=〈x0〉+b,
where 〈()〉 denotes the probabilistic mean (expectation) of ( ), and b the random part of x0. A prior estimate of 〈x0〉 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(3)ym-h(xm)=rm,
where ym is the vector of observations, h the observation operator, and rm 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 thatqk is a white Gaussian sequence and xk generated by (1) is Markov, which would be the case with linear dynamics. Assume also that b and rm are Gaussian random, unbiased (with zero mean), not cross-correlated and not correlated with qk. Then, according to Sherman and Bayes theorems (see Jazwinski [23, Chapter 5]), given observations ym(m=1,2,…M), the conditional maximum likelihood estimate of xk is also the conditional mean of xk(k=0,1,…K), and is given by the minimizer of the following cost function:(4)J=(x0-〈x0〉)TB-1(x0-〈x0〉)+∑m[ym-h(xm)]TRm-1[ym-h(xm)]+∑kqkTQk-1qk+∑kλkT[xk-f(xk-1)-qk],
where B = 〈bbT〉, Rm= 〈rmrmT〉  and Qk = 〈qkqkT〉  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 qk = 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 Rm-1 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 106 observations. To assimilate high-resolution radar observations over a mesoscale area, the model state vector xk often needs to contain at least as many as 106 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 onlyqk = 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 real-time 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 km3 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 Figure1, 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 their respective central time levels: t0-τ, t0 and t0+τ  (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 qk 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 t0 − τ and t0 (or t0 and t0 + τ) in four-dimensional space, the vector velocity increment can be estimated at t0 − τ/2 (or t0 + τ/2). Once the vector velocity fields are updated at t0 − τ/2 and t0 + τ/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 t0. 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.Figure 1
Flowchart for the three steps of the method. Step1 controls the quality of radial velocity observations vro and produces radial velocity analyses vra (on model grid) by minimizing J1 in (6) at t0 − τ, t0 and t0 + τ, respectively, for three consecutive radar volume scans. Step 2 produces the vector velocity increments vi by minimizing J2 in (11) at t0 − τ/2 and t0 + τ/2, respectively. Step 3 produces the perturbation pressure increment πi by minimizing Jπ in (16) first, and then produces the perturbation potential temperature increment θiby minimizing Jθ in (18) at t0, and finally adjusts the water vapor mixing ratio qv and the rain (or snow) mixing ratio qr (or snow mixing ratio qs) below (or above) the freezing level.


## 3. Data Quality Control and Radial Velocity Analysis

### 3.1. Data Quality Control

The quality control procedure checks for and corrects two types of errors in raw level II radial velocity data. First, it corrects the possible velocity alias error (if any) caused by the finite range of radar velocity measurements limited by the Nyquest velocityvN. Since the model predicted velocity vf = (uf, vf, wf) is used as the reference field, the dealiasing technique used in this paper is a simplified version of the three-step dealiasing of Gong et al. [28]. The dealiasing procedure is described as follows: (i)
Interpolatevf from model grid to each observation point with the radar beam height computed by using the equivalent Earth model (see Doviak and Zrnić [29, Section 2.2.3.1]). Project the interpolated vf onto the radial direction (along the radar beam) to obtain the background radial velocity vrf=ro·vf, where ro 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 vrd=vro-vrf at each observation point, where vro denotes the observation.(ii)
Take the integern nearest to vrd/(2vN) (where vN is the Nyquist velocity) and adjust the observational radial velocity to vro-2nvN so that the adjusted innovation (i.e., the adjusted observation minus the background value) is between ±vN.(iii)
If the absolute value of the adjusted innovation is less than 0.5vN, then the adjusted observation is accepted as a “good” one and used to replace the original vro. 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.25vN 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,wT,  is estimated empirically (Kessler [30]) by(5)wT=-3.088×100.00957Z,
where Z is the observed reflectivity in dBZ. The projection of the terminal velocity on the radial direction is wTsinα where α is the local tilt angle of the radar beam. Thus, with the above two types of correction, the originally observed vro is replaced by vro-2nvN-wTsinα at each observation point where n≠ 0 or wT≠0. For simplicity, the corrected radial velocity observation will still be denoted by vro in the remaining part of this paper.


### 3.2. Radial Velocity Analysis

The radial velocity analysis is performed by minimizing the following cost function:(6)J1=Jbk+Job,
where Jbk is the background term and Job 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 J1 is the analysis increment defined by vi=(ui,vi,wi)=va-vb where va and vb 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 (vb=0) but the background error covariance is nonzero in the background term Jbk, and the radio between the background and observation error variance is set to σvr2/σob2=10, where σvr2 and σob2 denote the radial velocity background and observation error variances, respectively.To facilitate the formulation and computation of the background term,va or, equivalently, vi (since vb=0 in step 1) is treated as an intermediate control variable with its two horizontal components (ui,vi) expressed in terms of stream function and velocity potential by (7)ui=-∂ψi∂y+∂χi∂x,vi=∂ψi∂x+∂χi∂y.
The grid field of (ui,vi) 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,wi) the state vectors of the grid fields of (ψi,χi,wi). These three state vectors can be related to three new intermediate vectors (ψ, χ, w) by (8)(ψi,χi,wi)=(B11/2ψ,B21/2χ,B31/2w),
where B1, B2, and B3 are the univariate error covariance matrices for the background fields ψf,χf, and wf, respectively. The background term is then given by(9)Jbk=|ψ|2+|χ|2+|w|2.
Here, it is assumed that the background errors are not cross-correlated between ψf,χf, and wf. With this assumption, ψi,χi, and wi are decoupled in the background term as shown by (9), so the computation is simplified.The background error covariances for (ψb,χb,wb) 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 Dψ = Dχ = Δz/√3 for (ψb,χb) in step 1. This gives Lw=Lχ/√3 = √6Δx and Dw = Δz for wb 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 wb is then related to σvr2 by σψ2=σχ2=0.5σvr2Lχ2 and σw2=4σvr2Dw2/Lχ2, 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)(10)Job={{(vri-vrd)2}}σob2,
where {{}} denotes the summation over all observation points, vrd=vro-vrb is the innovation, vri=ro·vi is the incremental radial velocity at each observational point, and ro 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 vb=0, we have vrb=ro·vb=0 and thus vrd=vro and vri=vra=ro·va in step 1. In (10), the observation errors are assumed to be not correlated between different observation points, so Rm (with m = M = 1) reduces to Iσob2 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 σvr2/σob2=10 as explained earlier, there is no need to estimate σob2 and σvr2 separately for the analysis in step 1.As in (9), vi and vri 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 vi is linearly related to the B-spline coefficients, vi is readily estimated once the cost function is minimized. Since only one time level of radial velocity observations is used, the estimated vi is not accurate, but its radial component vri=ro·vi (=vra 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.


## 3.1. Data Quality Control

The quality control procedure checks for and corrects two types of errors in raw level II radial velocity data. First, it corrects the possible velocity alias error (if any) caused by the finite range of radar velocity measurements limited by the Nyquest velocityvN. Since the model predicted velocity vf = (uf, vf, wf) is used as the reference field, the dealiasing technique used in this paper is a simplified version of the three-step dealiasing of Gong et al. [28]. The dealiasing procedure is described as follows: (i)
Interpolatevf from model grid to each observation point with the radar beam height computed by using the equivalent Earth model (see Doviak and Zrnić [29, Section 2.2.3.1]). Project the interpolated vf onto the radial direction (along the radar beam) to obtain the background radial velocity vrf=ro·vf, where ro 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 vrd=vro-vrf at each observation point, where vro denotes the observation.(ii)
Take the integern nearest to vrd/(2vN) (where vN is the Nyquist velocity) and adjust the observational radial velocity to vro-2nvN so that the adjusted innovation (i.e., the adjusted observation minus the background value) is between ±vN.(iii)
If the absolute value of the adjusted innovation is less than 0.5vN, then the adjusted observation is accepted as a “good” one and used to replace the original vro. 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.25vN 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,wT,  is estimated empirically (Kessler [30]) by(5)wT=-3.088×100.00957Z,
where Z is the observed reflectivity in dBZ. The projection of the terminal velocity on the radial direction is wTsinα where α is the local tilt angle of the radar beam. Thus, with the above two types of correction, the originally observed vro is replaced by vro-2nvN-wTsinα at each observation point where n≠ 0 or wT≠0. For simplicity, the corrected radial velocity observation will still be denoted by vro in the remaining part of this paper.


## 3.2. Radial Velocity Analysis

The radial velocity analysis is performed by minimizing the following cost function:(6)J1=Jbk+Job,
where Jbk is the background term and Job 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 J1 is the analysis increment defined by vi=(ui,vi,wi)=va-vb where va and vb 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 (vb=0) but the background error covariance is nonzero in the background term Jbk, and the radio between the background and observation error variance is set to σvr2/σob2=10, where σvr2 and σob2 denote the radial velocity background and observation error variances, respectively.To facilitate the formulation and computation of the background term,va or, equivalently, vi (since vb=0 in step 1) is treated as an intermediate control variable with its two horizontal components (ui,vi) expressed in terms of stream function and velocity potential by (7)ui=-∂ψi∂y+∂χi∂x,vi=∂ψi∂x+∂χi∂y.
The grid field of (ui,vi) 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,wi) the state vectors of the grid fields of (ψi,χi,wi). These three state vectors can be related to three new intermediate vectors (ψ, χ, w) by (8)(ψi,χi,wi)=(B11/2ψ,B21/2χ,B31/2w),
where B1, B2, and B3 are the univariate error covariance matrices for the background fields ψf,χf, and wf, respectively. The background term is then given by(9)Jbk=|ψ|2+|χ|2+|w|2.
Here, it is assumed that the background errors are not cross-correlated between ψf,χf, and wf. With this assumption, ψi,χi, and wi are decoupled in the background term as shown by (9), so the computation is simplified.The background error covariances for (ψb,χb,wb) 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 Dψ = Dχ = Δz/√3 for (ψb,χb) in step 1. This gives Lw=Lχ/√3 = √6Δx and Dw = Δz for wb 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 wb is then related to σvr2 by σψ2=σχ2=0.5σvr2Lχ2 and σw2=4σvr2Dw2/Lχ2, 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)(10)Job={{(vri-vrd)2}}σob2,
where {{}} denotes the summation over all observation points, vrd=vro-vrb is the innovation, vri=ro·vi is the incremental radial velocity at each observational point, and ro 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 vb=0, we have vrb=ro·vb=0 and thus vrd=vro and vri=vra=ro·va in step 1. In (10), the observation errors are assumed to be not correlated between different observation points, so Rm (with m = M = 1) reduces to Iσob2 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 σvr2/σob2=10 as explained earlier, there is no need to estimate σob2 and σvr2 separately for the analysis in step 1.As in (9), vi and vri 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 vi is linearly related to the B-spline coefficients, vi is readily estimated once the cost function is minimized. Since only one time level of radial velocity observations is used, the estimated vi is not accurate, but its radial component vri=ro·vi (=vra 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.


## 4. Vector Velocity Analysis

The vector velocity analysis is performed by minimizing the following cost function:(11)J2=Jbk+Job+Jrm+Jms.
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 qk, so the equation term in (4) reduces to the sum of the mass continuity equation term Jms and radial momentum equation term Jrm 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 vb=vf, where vf 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 wf can be estimated by Lψ=Lχ=6√2Δx and Lw=Lχ/√3 (see (A.15)–(A.17)). The vertical decorrelation depth is set empirically to Dw=Δzforwf, so Dψ=Dχ=Δz/√3for(ψf,χf) (see (A.8) and (A.10)). The analysis time level for Jbk in (11) is at t0-τ/2(ort0+τ/2)—the middle time level between t0-τ and t0 (or t0 and t0+τ), so the incremental field vi (=va-vf) and background velocity field vf are also all at t0-τ/2(ort0+τ/2).The observation termJob for J2 in (11) has the same form as that in (10) for J1, but the background velocity is given by vb=vf (instead of vb=0 as in step 1) and the observations are binned from the two time levels at t0-τ  and t0 (or t0 and t0+τ) to the middle-time level at t0-τ/2(ort0+τ/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 σvr2+σob2=σd2, where σd2 denotes the radial velocity innovation variance, σvr2 and σob2 denote the radial velocity background and observation error variances, respectively. Since zero background is used in step 1, the spatial variance of the innovation (vro minus vrf in the observation space) is used to estimate σd2 in each assimilation cycle. The background error variance is estimated by σvr2=σd2-σob2 with σob given by a pre-estimated value. For the phased-array radar data used in this paper, σob=2ms-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(12)Jms=[[Ems2]]σms2,
where [[]] denotes summation over all grid points, Ems is the mass continuity equation constraint for vi, and σms2 is the equation error variance or, specifically, the variance of the residual error of Ems. 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:
(13)Ems=(∂x+Gx∂σ)(ρΘvui)+(∂y+Gy∂σ)(ρΘvvi)+Gz∂σ(ρΘvwi),where Gx=a∂xzs,Gy=a∂yzs,Gz=∂zσ, a=(s-H)/(H-zs),σ=H(z-zs)/(H-zs),H and zs 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 Qk (with k=K=1) reduces to Iσms2 and the equation term in (4) reduces to Jms. 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/Ldiv=√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(14)Jrm=[[Erm2]]σrm2,
where Erm is the radial momentum equation in an incremental form, and σrm2 is the variance of the residual error of Erm. 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:(15)Erm=∂tvri+va·∇vri+vi·∇vrf-[vi·(vi+2vf)-(vra)2+(vrf)2]r+ro·(fk×vi)-Dvri,
where va=vf+vi,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 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∂tvri in (15) is discretized by the standard central finite difference scheme with vra given by the gridded observations (obtained from step 1) at t0-τ  and t0 (or t0 and t0+τ). All the remaining terms in (11) are at the middle time level t0-τ/2 (or t0+τ/2), so vra 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 low-altitude 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 Jrm and Qk reduces to Iσrm2, 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/Lrot+σvrσdiv/Ldiv=2√2σvr2/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, σrm=3√2σvr2/Lχ is estimated.As in step 1, by using (7)-(8) and the B-spline representation, vi and vri 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 J2 in (11) is minimized efficiently in the reduced space spanned by the B-spline coefficients—the final control variables. The estimated B-spline coefficients are then transformed back to (ψ,χ,w) and finally to vi. The total vector velocity analysis is given by va=vf+vi on the model grid.


## 5. Thermodynamic Analysis

### 5.1. Perturbation Pressure Analysis

The estimated fields ofvi at t0-τ/2 and t0+τ/2 from step 2 are used as input data to estimate the thermodynamic fields at t0 in step 3. Note that the Exner function π=(p′/p00)R/Cp is used in place of the perturbation pressure p′ in COAMPS, where p00 is a constant reference pressure, R is the gas constant for dry air, and Cp 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: (16)Jπ=|π|2σπ2+[[Eu2+Ev2]]σuv2.
Here,|π|2/σπ2=πiTCπ-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 σπ2Cπ is the error covariance matrix for πf. The second term on the right-hand side of (16) is the equation term. In this term, Eu and Ev are the two horizontal momentum equations retained from COAMPS (see Hodur [22, (2)-(3)]) in the following incremental forms:(17)Eu=∂tui+va·∇ui+vi·∇uf-fvi-Dui+CpΘv(∂xπi+Gx∂σπi),Ev=∂tvi+va·∇vi+vi·∇vf+fui-Dvi+CpΘv(∂yπi+Gy∂σπi),
and σuv2 is the equation error variance for Eu and Ev. Here, ∂tui and ∂tvi are computed by the standard central finite difference scheme from the grid fields of (ui,vi) at t0-τ/2 and t0+τ/2 produced in step 2. All the remaining terms in (17) are at the middle time level t0.As explained in Section2, the cost function in (16) is derived from (4) by setting M=0(∑m=0) and K=1 with λ=0, retaining only Eu and Ev in qk, and reducing Qk to Iσuv2. The two error variances σπ2 and σuv2 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 B-spline 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 Dw as for wf in (A.4).


### 5.2. Perturbation Temperature Analysis

Afterπi is computed, the perturbation potential temperature increment θi is estimated by minimizing the following cost function:(18)Jθ=|θ|2σθ2+[[Ew2]]σEw2.
Here, |θ|2/σθ2=θiTCθ-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 potential temperature θf, σθ2 is the error variance for θf, and thus σθ2Cθ 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:(19)Ew=∂twi+CpΘvGz∂σπi-gθiΘ+(vi+vf)·∇wi+vi·∇wf-Dwi,
and σEw2 is the equation error variance for Ew. 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, ∂twi are computed by the standard central finite difference scheme from the grid fields of wi at t0-τ/2 and t0+τ/2 produced in step 2. All the remaining terms in (19) are at the middle time level t0.The cost function in (18) is derived from (4) by setting M=0(∑m=0) and K=1 with λ=0, retaining only Ew in qk, and reducing Qk to IσEw2. Again, as those in (16), the two error variances σθ2 and σEw2 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).


### 5.3. 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 ratioqv. After this, qv 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:
(20)Jη=|η|2+{{(ηa-ηo)2}}σηo2,
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 B-spline 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 qv 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 qv to the value interpolated in x- and y-directions from the nearest grid points where the two conditions are not both satisfied. This adjustment is designed to reduce qv and thus correct the model’s tendency of overpredicting precipitation locally in areas detected by the above two conditions.In companion with theqv adjustment, the rain mixing ratio qr (or snow mixing ratio qs) 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]).


## 5.1. Perturbation Pressure Analysis

The estimated fields ofvi at t0-τ/2 and t0+τ/2 from step 2 are used as input data to estimate the thermodynamic fields at t0 in step 3. Note that the Exner function π=(p′/p00)R/Cp is used in place of the perturbation pressure p′ in COAMPS, where p00 is a constant reference pressure, R is the gas constant for dry air, and Cp 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: (16)Jπ=|π|2σπ2+[[Eu2+Ev2]]σuv2.
Here,|π|2/σπ2=πiTCπ-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 σπ2Cπ is the error covariance matrix for πf. The second term on the right-hand side of (16) is the equation term. In this term, Eu and Ev are the two horizontal momentum equations retained from COAMPS (see Hodur [22, (2)-(3)]) in the following incremental forms:(17)Eu=∂tui+va·∇ui+vi·∇uf-fvi-Dui+CpΘv(∂xπi+Gx∂σπi),Ev=∂tvi+va·∇vi+vi·∇vf+fui-Dvi+CpΘv(∂yπi+Gy∂σπi),
and σuv2 is the equation error variance for Eu and Ev. Here, ∂tui and ∂tvi are computed by the standard central finite difference scheme from the grid fields of (ui,vi) at t0-τ/2 and t0+τ/2 produced in step 2. All the remaining terms in (17) are at the middle time level t0.As explained in Section2, the cost function in (16) is derived from (4) by setting M=0(∑m=0) and K=1 with λ=0, retaining only Eu and Ev in qk, and reducing Qk to Iσuv2. The two error variances σπ2 and σuv2 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 B-spline 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 Dw as for wf in (A.4).


## 5.2. Perturbation Temperature Analysis

Afterπi is computed, the perturbation potential temperature increment θi is estimated by minimizing the following cost function:(18)Jθ=|θ|2σθ2+[[Ew2]]σEw2.
Here, |θ|2/σθ2=θiTCθ-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 potential temperature θf, σθ2 is the error variance for θf, and thus σθ2Cθ 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:(19)Ew=∂twi+CpΘvGz∂σπi-gθiΘ+(vi+vf)·∇wi+vi·∇wf-Dwi,
and σEw2 is the equation error variance for Ew. 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, ∂twi are computed by the standard central finite difference scheme from the grid fields of wi at t0-τ/2 and t0+τ/2 produced in step 2. All the remaining terms in (19) are at the middle time level t0.The cost function in (18) is derived from (4) by setting M=0(∑m=0) and K=1 with λ=0, retaining only Ew in qk, and reducing Qk to IσEw2. Again, as those in (16), the two error variances σθ2 and σEw2 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).


## 5.3. 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 ratioqv. After this, qv 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:
(20)Jη=|η|2+{{(ηa-ηo)2}}σηo2,
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 B-spline 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 qv 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 qv to the value interpolated in x- and y-directions from the nearest grid points where the two conditions are not both satisfied. This adjustment is designed to reduce qv and thus correct the model’s tendency of overpredicting precipitation locally in areas detected by the above two conditions.In companion with theqv adjustment, the rain mixing ratio qr (or snow mixing ratio qs) 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]).


## 6. 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 (Figure2). 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 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.Phased-array radar observed reflectivity (a), radial velocity (b) at0.75°, and dealiased radial velocity (c) at 2107 UTC on 2 June 2004. Aliased velocity areas are highlighted in (b).
(a)(b)(c)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 Figure3. 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.Figure 3
COAMPS domains nested with 54, 18, and 6 km grids.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 atz=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(b).Velocity (arrows for the horizontal wind and green contours for the vertical velocity) and reflectivity (color field) atz = 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.
(a)(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 Figure5(a) while the analyses from the third cycle at 2118 UTC are plotted in Figure 5(b). For comparison, radar observations of reflectivity at 2118 UTC are also shown (by solid contours) in Figures 5(a) and 5(b). It should be mentioned again that the analyzed fields in Figure 5(b) use the predicted velocity and reflectivity fields in Figure 5(a) as background fields. It is apparent that the predicted background reflectivity in Figure 5(a) already becomes quite close to the observations, and the minor miss-matches are further corrected by the analysis as shown by the analyzed reflectivity in Figure 5(b).As in Figure4 but for velocity and reflectivity produced by (a) the 5-minute prediction (from the second cycle at 2113 UTC) and (b) analysis of the third cycle at 2118 UTC in Test 3.
(a)(b)Figure6(a) plots the incremental horizontal velocity and potential temperature (at z=0.5 km) produced by the 10-minute 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.As in Figure4 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.
(a)(b)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 atz=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.As in Figure4 but for velocity and reflectivity produced by the predictions valid at 2200 UTC in (a) the control run and (b) Test 1.
(a)(b)(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.
(a)(b)To verify the three-dimensional winds, the model-produced (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 Figure8(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.


## 7. 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 phased-array 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 real-time 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×600km2). 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 flow-dependent forms of background error covariances, which is under our current investigation in connection with the ensemble Kalman filter.



---

*Source: 797265-2010-06-27.xml*