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.
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: 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 where denotes the probabilistic mean (expectation) of ( ), and b the random part of x0. A prior estimate of 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 ym is the vector of observations, h the observation operator, and rm the observation error, and the subscript 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 qk 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 ]), given observations , the conditional maximum likelihood estimate of xk is also the conditional mean of , and is given by the minimizer of the following cost function: where B = , Rm= and Qk = are the error covariance matrices for the background, observations, and model equations, respectively, ∑m denotes the summation over (from 1 to ), ∑k denotes the summation over (from 1 to ), denotes the transpose of ( ), and 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 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 only qk = 0 but also M = 1, K = 0 (∑k = 0) and 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 grid for a 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 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 ) 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: , and (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 − and (or and + ) in four-dimensional space, the vector velocity increment can be estimated at − /2 (or + /2). Once the vector velocity fields are updated at − /2 and + /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 . 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.
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 velocity . Since the model predicted velocity = (, , ) 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)Interpolate 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 ]). Project the interpolated onto the radial direction (along the radar beam) to obtain the background radial velocity , where 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, () and ()]). Calculate the innovation at each observation point, where denotes the observation.(ii)Take the integer n nearest to (where is the Nyquist velocity) and adjust the observational radial velocity to so that the adjusted innovation (i.e., the adjusted observation minus the background value) is between ±.(iii)If the absolute value of the adjusted innovation is less than 0.5, then the adjusted observation is accepted as a “good” one and used to replace the original . Otherwise, perform the continuity check (buddy check) which is similar to the third step described in Gong et al. [28, Section ] but the threshold value is set to 0.25 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, , 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 where α is the local tilt angle of the radar beam. Thus, with the above two types of correction, the originally observed is replaced by at each observation point where n 0 or . For simplicity, the corrected radial velocity observation will still be denoted by in the remaining part of this paper.
3.2. Radial Velocity Analysis
The radial velocity analysis is performed by minimizing the following cost function: where is the background term and 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 is the analysis increment defined by where and 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 () but the background error covariance is nonzero in the background term , and the radio between the background and observation error variance is set to , where and denote the radial velocity background and observation error variances, respectively.
To facilitate the formulation and computation of the background term, or, equivalently, (since in step 1) is treated as an intermediate control variable with its two horizontal components expressed in terms of stream function and velocity potential by The grid field of () can be computed from the grid field of () by using the standard central finite difference scheme (in the horizontal). Denote by () the state vectors of the grid fields of (). These three state vectors can be related to three new intermediate vectors (, , w) by where B1, B2, and B3 are the univariate error covariance matrices for the background fields , and , respectively. The background term is then given by Here, it is assumed that the background errors are not cross-correlated between , and . With this assumption, , and are decoupled in the background term as shown by (9), so the computation is simplified.
The background error covariances for () 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 2x waves. Based on these considerations, we set the horizontal decorrelation length to (which is 1/2 of that used in step 2 as shown Section 4) and the vertical decorrelation depth to = = for () in step 1. This gives = and = for according to (A.8) and (A.9), where (=6 km in this study) is the horizontal grid spacing and 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 () and is then related to by and , 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 innovation, is the incremental radial velocity at each observational point, and 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, () and ()]). As explained earlier, since , we have and thus and 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 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 as explained earlier, there is no need to estimate and separately for the analysis in step 1.
As in (9), and 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 is linearly related to the B-spline coefficients, is readily estimated once the cost function is minimized. Since only one time level of radial velocity observations is used, the estimated is not accurate, but its radial component (= 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: This cost function is derived from (4) by setting and to retain only the mass continuity equation and the radial momentum equation in , so the equation term in (4) reduces to the sum of the mass continuity equation term and radial momentum equation term 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 , where 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 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 (). Thus, the decorrelation lengths for () and can be estimated by and (see (A.15)–(A.17)). The vertical decorrelation depth is set empirically to , so (see (A.8) and (A.10)). The analysis time level for in (11) is at —the middle time level between and (or and ), so the incremental field () and background velocity field are also all at .
The observation term for in (11) has the same form as that in (10) for , but the background velocity is given by (instead of as in step 1) and the observations are binned from the two time levels at and (or and ) to the middle-time level at . Since the background and observation errors are uncorrelated (Xu and Wei [35]), the sum of their error variances is given by the innovation variance , where denotes the radial velocity innovation variance, and denote the radial velocity background and observation error variances, respectively. Since zero background is used in step 1, the spatial variance of the innovation ( minus in the observation space) is used to estimate in each assimilation cycle. The background error variance is estimated by with given by a pre-estimated value. For the phased-array radar data used in this paper, 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, is the mass continuity equation constraint for , and is the equation error variance or, specifically, the variance of the residual error of . 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, () and ()]). From the COAMPS pressure tendency equation, the extracted mass continuity equation has the following incremental form: where , and are the heights of the model top and bottom terrain boundaries, respectively, the basic state virtual potential temperature, and the basic state density. Here, the equation errors are assumed to be not spatially correlated, so (with ) reduces to and the equation term in (4) reduces to . The first two terms on the right-hand side of (13) are the horizontal divergence and the associated RMS error can be estimated by (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, is estimated.
The last term in (11) is given by where is the radial momentum equation in an incremental form, and is the variance of the residual error of . The radial momentum equation is obtained by projecting the model vector momentum equation (Hodur [22, ()–()]) onto the radial direction along the radar beam, and then converted into the following incremental form: where , f is the Coriolis parameter, is the radial distance from the radar, is the gradient operator in the -coordinates, is the diffusion operator for subgrid-scale turbulent mixing as in COAMPS, and 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 in (15) is discretized by the standard central finite difference scheme with given by the gridded observations (obtained from step 1) at and (or and ). All the remaining terms in (11) are at the middle time level (or ), so 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 and reduces to , 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 (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, is estimated.
As in step 1, by using (7)-(8) and the B-spline representation, and 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 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 () and finally to . The total vector velocity analysis is given by on the model grid.
5. Thermodynamic Analysis
5.1. Perturbation Pressure Analysis
The estimated fields of at and from step 2 are used as input data to estimate the thermodynamic fields at in step 3. Note that the Exner function is used in place of the perturbation pressure in COAMPS, where is a constant reference pressure, is the gas constant for dry air, and is the specific heat at constant pressure for the atmosphere. In this paper, we simply call “perturbation pressure”. The perturbation pressure increment is estimated by minimizing the following cost function: Here, is the background term and is the control variable. In this term, is the state vector of the grid field of is the error correlation matrix for the background perturbation pressure is the error variance for and thus is the error covariance matrix for . The second term on the right-hand side of (16) is the equation term. In this term, and are the two horizontal momentum equations retained from COAMPS (see Hodur [22, ()-()]) in the following incremental forms: and is the equation error variance for and . Here, and are computed by the standard central finite difference scheme from the grid fields of () at and produced in step 2. All the remaining terms in (17) are at the middle time level .
As explained in Section 2, the cost function in (16) is derived from (4) by setting and with , retaining only and in , and reducing to . The two error variances and in (16) are difficult to estimate. Their ratio, however, can be tuned adaptively to make the cost function approximately equally partitioned by the two terms on the right-hand side of (16) when reaches the minimum. This approach is used in this paper. In (16), the grid field of 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 is minimized efficiently in the reduced B-spline space. The estimated B-spline coefficients are then transformed back to and finally to . The transformation represented by is simulated by a recursive filter (Purser et al. [31]), in which the error covariance for is modeled by a horizontally isotropic Gaussian function. In this filter, the horizontal decorrelation length is set to as for () in (A.3) and (A.10), while the vertical decorrelation depth is set to as for in (A.4).
5.2. Perturbation Temperature Analysis
After is computed, the perturbation potential temperature increment is estimated by minimizing the following cost function: Here, is the background term and ) is the control variable. In this term, is the state vector of the grid field of , is the error correlation matrix for the background perturbation potential temperature , is the error variance for , and thus is the error covariance matrix for . The second term on the right-hand side of (18) is the equation term. In this term, is the vertical momentum equation retained from COAMPS in the following incremental form: and is the equation error variance for . In (19), is the basic state potential temperature and 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, ()]). Here, are computed by the standard central finite difference scheme from the grid fields of at and produced in step 2. All the remaining terms in (19) are at the middle time level .
The cost function in (18) is derived from (4) by setting and with , retaining only in , and reducing to . Again, as those in (16), the two error variances and in (18) are hard to estimate, but their ratio can be tuned adaptively to make the cost function approximately equally partitioned by the two terms on the right-hand side of (18) when 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 can be minimized efficiently in the reduced B-spline space. The estimated B-spline coefficients are then transformed back to and finally to . The transformation represented by is simulated by a recursive filter with the decorrelation length and depth set to be the same as those for .
The analyzed total perturbation pressure and potential temperature are given by and , 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 ratio . After this, is further adjusted based on the radar observed reflectivity through the following three substeps:Interpolate the radar observed reflectivity onto the model grid by minimizing the following cost function: where is the error covariance matrix for zero reflectivity background, is the state vector of the analyzed reflectivity , and (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 , can be minimized efficiently in the reduced B-spline space. The estimated B-spline coefficients are then transformed back to and finally to . The transformation represented by is simulated by a recursive filter, in which the decorrelation length and depth are set to be as small as and , respectively, and the error variance for zero reflectivity background is estimated by the spatial variance of the model predicted reflectivity . With the above parameter settings, the analyzed reflectivity obtained by minimization in (20) is an optimal interpolation of the observed reflectivity onto the model grid. After this interpolation, 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). Check each grid point for the conditions of dBZ and dBZ. If these two conditions are both satisfied, then adjust to the saturated value if the analyzed vertical velocity is nonnegative () or to 80% of the saturated value if the analyzed vertical velocity is negative () and the relative humidity is below 80% at this grid point.Check each grid point for the conditions of 10 dBZ and 10 dBZ. If these two conditions are both satisfied, then adjust to the value interpolated in - and -directions from the nearest grid points where the two conditions are not both satisfied. This adjustment is designed to reduce and thus correct the model’s tendency of overpredicting precipitation locally in areas detected by the above two conditions.
In companion with the adjustment, the rain mixing ratio (or snow mixing ratio ) is adjusted below (or above) the freezing level to make the computed reflectivity match 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 (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 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 in the azimuthal direction.
(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 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 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).
(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 Figure 5(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).
(a)
(b)
Figure 6(a) plots the incremental horizontal velocity and potential temperature (at 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.
(a)
(b)
It is necessary to mention that the incremental potential temperature produced by the thermodynamic analysis is relatively small (within ± 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 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.
(a)
(b)
(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 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.
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 () COAMPS grid that covers a mesoscale area (). 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.
Appendix
Background Error Variances and Decorrelation Lengths for Derived Variables
The mass continuity equation in (13) can be applied approximately to the background velocity . By neglecting the effect of terrain slope, this equation can be written into the following form: where is the Laplacian. Assume that the random error fields of and are also constrained by (A.1), so their covariance functions (see Xu and Wei [35, () and ()]) satisfy the following relationship: where 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 can be modeled by a Gaussian function in the horizontal direction, so where and denote, respectively, the error variance and vertical error correlation for denotes the Gaussian function form, and is the horizontal decorrelation length defined by as in Daley [42, ()]. By assuming that can be modeled by a Gaussian function in the vertical, we get where and denote, respectively, the error variance and vertical error correlation for , and is the vertical decorrelation depth defined by . Substituting (A.3) and (A.4) into (A.2) and setting with , we have Substituting (A.3) and (A.4) into (A.2) again and setting with , we have another equation Setting in (A.5) or in (A.6) gives the ratio of the two variances
Note that and are obtained explicitly in (A.5) and (A.6), respectively. Substituting (A.5) into (A.3) gives the explicit form of , and substituting (A.6) into (A.4) results in another explicit form for . Using these explicit forms, we can verify that the vertical decorrelation depth for the random error field of is given by while the horizontal decorrelation length for the random error field of is given by
Denote by the horizontal component of and represent by and the rotational and divergent parts of , respectively. Then we have , and . Also assume that the error covariance function for can be modeled by a Gaussian function in the horizontal, so where and denote, respectively, the error variance and vertical error correlation for , and is the horizontal decorrelation length. and represent the traces of the error covariance tensors for and , respectively. By using () of Xu and Wei [35] with (A.4) and (A.10), one can get The error variances for and can be then obtained by setting and in (A.11) and (A.12), that is, Using (A.11) and (A.12), one can verify that the horizontal decorrelation lengths for the random error fields of and are given, respectively, by
For mesoscale flows, we may assume that and have roughly the same error variances and decorrelation lengths. In this case, we have where is the error variance for (the horizontal component of ), and is the error variance for . Note that , and therefore . Also note that is the radial component of and approximately the radial component of , so is a reasonable approximation (Xu and Gong [43]). These relationships are used in the derivation of (A.18). Substituting (A.18) into (A.13) and (A.7) gives
Acknowledgments
The authors are thankful to Carl Hane, Jidong Gao, and the anonymous reviewer for their comments and suggestions that improved the presentation of the paper and to Douglas Forsyth, Kurt Hondl, Richard Adams, Pengfei Zhang and Kang Nai for their help in obtaining and processing the phased-array radar data. The research work was supported by the NOAA HPCC program, the FAA contract IA# DTFA03-01-X-9007 to NSSL, the ONR Grants N000140410312 to the University of Oklahoma, and the NOAA-University of Oklahoma Cooperative Agreement no. NA17RJ1227.