# WRF-ARW Variational Storm-Scale Data Assimilation: Current Capabilities and Future Developments

Authors: Juanzhen Sun; Hongli Wang
Journal: Advances in Meteorology
 (2013)

---

## Abstract

The variational radar data assimilation system has been developed and tested for the Advanced Research Weather Research and Forecasting (WRF-ARW) model since 2005. Initial efforts focused on the assimilation of the radar observations in the 3-dimensional variational framework, and recently the efforts have been extended to the 4-dimensional system. This article provides a review of the basics of the system and various studies that have been conducted to evaluate and improve the performance of the system. Future activities that are required to further improve the system and to make it operational are also discussed.


---

## Body

## 1. Introduction

In the past two decades active research was conducted on the development of techniques to initialize storm-scale numerical prediction models. It has been recognized that the success will critically depend on the optimal use of the national operational WSR-88D radar network that covers the United States with single Doppler coverage in most areas. Although the network provides observations at a resolution that is able to resolve atmospheric convection, they are only limited to radial wind and reflectivity. Therefore several early studies focused on the feasibility of retrieving meteorological fields from these single Doppler observations. Techniques with different complexities have been developed which aim at obtaining the unobserved meteorological variables such as 3-dimensional (3D) wind, temperature, and microphysical fields from the radar observations of radial velocity and reflectivity (e.g., [1–5]).The techniques that make use of a numerical model in a data assimilation (DA) context received particular attention because they combine the retrieval, initialization, and forecast in one system. The first radar DA system for the storm-scale was developed based on the 4-dimensional variational data assimilation (4D-Var) technique and a boundary layer fluid dynamics model for the retrieval of the 3D wind and temperature [1]. This system, known as VDRAS (Variational Doppler Radar Analysis System), was later expanded to include microphysical retrieval, as well as short-term forecasts initialized by these retrieved fields [6–9]. Another variational-based radar DA system was developed by Gao et al. [4] using a 3-dimensional variational data assimilation (3D-Var) technique in the framework of the ARPS (Advanced Research and Prediction System [10]) model. A so-called 3.5-dimensional variational radar data assimilation based on Navy’s COAMPS (The Coupled Ocean/Atmosphere Mesoscale Prediction System) was developed and demonstrated through a number of studies [11, 12]. These variational systems showed great potentials in the use of radar observations for initializing high-resolution numerical models through several case studies and real-time demonstrations [13, 14].The demonstrated potential motivated the development of a radar DA scheme in the variational DA system WRFVAR of the community model ARW-WRF (Advanced Research Weather Research and Forecasting; hereafter refer to WRF). WRFVAR includes both 3D-Var [15, 16] and 4D-Var [17] components. The radar DA scheme was first developed for WRF 3D-Var [18, 19] and recently expanded to 4D-Var [20, 21]. Unlike the ARPS 3D-Var system and VDRAS whose developments were motivated by the assimilation of radar observations for convective-scale analysis, the WRFVAR system was started with conventional data assimilation without the convective-scale analysis as a priority. It thus emphasized the synoptic-scale balance through the use of model forecast background and its error covariance. It also used the stream function and velocity potential as the momentum control variables instead of the direct velocities as in VDRAS and ARPS 3D-Var. In addition, the incremental formulation [22] that is commonly used in large-scale data assimilation systems was adopted by WRFVAR, which might not be suitable for the storm-scale applications because of the linearization required in the formulation. Therefore, the addition of a radar data assimilation capability to WRF, which is a forecast system aiming at broad research and operational applications, presents new challenges that do not exist in the systems that focus on the storm-scale analysis and forecasting.In recent years several operational centers throughout the world have also been developing and testing capabilities of assimilating radar observations to initialize high-resolution convection-permitting NWP models. One of the notable efforts was taken by the Met Office of the United Kingdom using the unified model and its variational data assimilation system [23]. The lessons learnt in understanding and finding solutions to the aforementioned challenges in WRFVAR have a wide application to the operational community because some operational data assimilation systems use the similar framework as WRFVAR and will hence face the same challenges.The purpose of this paper is to provide a review of the progress of the WRFVAR system and its planned future developments. Since the initial development of the system, a good number of studies have been conducted which either aim at the further development of the system or at answering some challenging scientific questions. A review of these studies will benefit the future users of the community system. In Section2, we provide a review on the overall system design, including the fundamental framework of WRFVAR, radar observation operators, control variables and their background error (BE) estimates, and automated radar data quality control. The performance of the 3D-Var system is reviewed in Section 3 and of the 4D-Var system in Section 4. At the last section, we discuss the planned future efforts for the DA system.


## 2. System Design of WRFVAR Radar DA

### 2.1. Fundamental Framework of WRFVAR

The WRFVAR follows the incremental variational formulation [22] that is commonly used in operational systems. In the standard variational formulation, the optimal analysis is obtained by minimizing a cost function measuring the misfit between observation and prediction by a nonlinear model. In the incremental variational formulation, however, a linearized model is used to generate the prediction needed in the cost function. The advantage of the incremental approach is that it not only reduces the computational cost of 4D-Var but also improves the mathematical conditioning of the cost function (i.e., a quadratic cost function) because of the linearization of the forward operator. It facilitates the use of a coarser resolution tangent linear model (TLM) and adjoint model (ADM) to reduce the computational cost and a simplified representation of physical processes in TLM and ADM in the data assimilation cycle.The cost function of the incremental formulation can be derived from the standard cost function with some assumptions that will be described later. First we express the standard cost function as follows:(1)J(x0)=12(x0-x0b)TB-1(x0-x0b)+12∑k=0K{yko-Hk[Mk(x0)]}T×R-1{yko-Hk[Mk(x0)]}.
This cost function assumes that the assimilation window covers K observation windows with each represented by the subscript k. The variables of x0,  x0b, and yko represent the initial atmospheric state, the background state which can be previous model forecast, and the observed state, respectively. Mk is the nonlinear prediction model to propagate the initial atmospheric state to that at the kth observation time in the case of 4D-Var. It should be noted that 3D-Var is regarded as a special case in which K can simply set to 0 and Mk set to the identity matrix. Hk is the nonlinear observation operator.  B  and  R  are the background and observation error covariance matrices, respectively.To derive the cost function of the incremental formulation, we first introduce the innovation variabled at the kth observation time as
(2)dk=yko-Hk[Mk(x0)]
and then linearize the operators Hk and Mk in (1) as follows:
(3)Hk[Mk(x0n)]≈Hk[Mk(x0n-1)]+HkMk[x0n-1]δx0n,
where the superscripts n and n-1 represent the current and previous outer loop iterations, respectively. Hk and Mk are the tangent linear operator of Hk and Mk. Using (2) and (3), the cost function (1) can be changed to
(4)Jn[δx0n]=12{δx0n-[x0b-x0n-1]}TB-1{δx0n-[x0b-x0n-1]}+12∑k=0K{HkMk[x0n-1]δx0n-d0n-1}T×R-1{HkMk[x0n-1]δx0n-dkn-1}.
When the increment δx0n at the nth outer loop iteration is obtained, the estimate of the atmospheric state is updated by x0n=x0n-1+δx0n and used to produce the first guess trajectory Mk[x0n] for the next outer loop n+1. The basic assumption of the incremental approach is that the solution of the cost function (4) gradually approaches that of the nonlinear cost function (1) given enough number of outer loop iterations.Assuming the background error covariance matrix is given byB=UUT and a control variable transform by δx0=Uv, the background term of the cost function (4) is simplified to (1/2)(∑i=1Nvi)T(∑i=1Nvi), where N is the total number of outer loops, whereby avoiding the difficulty in computing the inverse of B. The control variable transform δx0=Uv is implemented through a series of operations U=UpUvUh [15]. The horizontal transform Uh is performed using a recursive filter [24]. The vertical transform Uv applies an EOF (empirical orthogonal function) decomposition on the vertical component of the background error covariance. The physical transform Up converts the increment in control variable space to analysis variable space. The control variables of the WRFVAR system will be described later in this section.The optimal analysis is obtained by minimizing the cost function (4) by iterating in both the outer and inner loops. The number of inner loop iterations is determined either by a convergence criterion or a prespecified fixed number. The number of the outer loop iterations is typically 2 to 6. The gradient of the cost function required in the minimization process is computed by the adjoint method. For the 4D-Var system with the incremental formulation, a forward model and its TLM and ADM are needed. The forward model is the same as the WRF forecast model but currently only with the physical schemes of the Kessler microphysics and the diffusion. More schemes will be added to the forward model in the future as the adjoint of these schemes has been developed.


### 2.2. Radar Observation Operators

To assimilate radar observations, the following observation terms are added to the existing cost function:(5)J=Jold+12∑k=0k[(Vrk-Vrkob)TRv-1(Vrk-Vrkob)nnnnnnnn+12(Zk-Zkob)TRz-1(Zk-Zkob)(Vrk-Vrrkob)T],
where Jold is used to represent the existing cost function before radar data assimilation is developed. The variables Vr stand for the radial velocity and Z for the reflectivity factor. The superscript “ob” indicates the observations. The symbols Rv-1 and Rz-1 stand for observation error covariance matrices for radial velocity and reflectivity, respectively. Note that the summation over the observation time levels k is not needed in the case of 3D-Var. The observation operator Hk in the cost function (1) links the model variables in a model coordinate to the observation variables in an observation space. For the radar radial velocity, this linkage is formulated with the 3D wind field (u, v, and w), the hydrometeor fall speed Vt, and the distance D between the location of a data point and the radar antenna:
(6)Vr=1D[(xd-xr)u+(yd-yr)v+(zd-zr)(w-Vt)],(7)D=(xd-xr)2+(yd-yr)2+(zd-zr)2,
where (xd,yd,zd) represents the location of the observation point and (xr,yr,zr) represents the location of the radar station. Vt is calculated from the rainwater mixing ratio with a height correction following Sun and Crook [6]. Relation (6) is linear except that Vt is nonlinearly dependent on qr, which needs to be linearized to obtain the linear observation operator for Vr. Note that the earth curvature effect [25] must be considered when mapping the radar observations to the model grid, which is done in the radar data preprocessing package that is described in the next subsection.The formulation of the reflectivity operator is not as straightforward because it depends on the assumption of drop size distribution in a microphysical parameterization scheme and the classification of hydrometeors. Following Sun and Crook [6], Xiao et al. [19] used the following relation for WRF 3D-Var by assuming the Marshall-Palmer drop size distribution for rain:
(8)Z=43.1+17.5log10(ρqr),
where ρ is the air density (kg m−3) and qr is the rainwater mixing ratio (g kg−1). Zhang et al. [26] examined a different relation by assuming a constrained gamma drop size distribution and found some improvement in the analysis of low precipitation.Equation (8) was used in the cost function (5) to assimilate reflectivity in the WRF 3D-Var developed by Xiao et al. [19]. However, Wang et al. [27] showed that the linearization of (8) as required by the incremental formulation of the cost function can result in a dry bias in rainwater analysis. Thus they proposed to indirectly assimilate the derived rainwater mixing ratio qr from (8) by replacing Z with qr in (5).


### 2.3. Control Variables and Background Error Estimate

The standard control variables (CV) [15] for WRF 3D-Var and 4D-Var are the stream function ψ, the unbalanced component of velocity potential χu, the unbalanced component of temperature Tu, the unbalanced component of surface pressure Psu, and “pseudo” relative humidity Rhs. The unbalanced component for each of the three variables is the residual after subtracting the respective correlated component with ψ which is obtained by a statistical fitting. The “pseudo” relative humidity is the water vapor mixing ratio divided by its saturated value in the background state. To assimilate radar radial velocity and reflectivity associated with warm rain, vertical velocity and the microphysical parameters of cloud water and rainwater must be added as control variables. In the current WRFVAR system, no cross-correlation with other variables is considered for relative humidity, cloud water mixing ratio, rain water mixing ratio, and vertical velocity. The 4D-Var technique through the use of the WRF model can implicitly produce multivariate correlations for the analysis.One issue that raises concerns for the convective-scale DA is the selection of momentum control variables. The use of stream function and velocity potential (ψ and χ) was widely accepted by global operational models (e.g., National Center of Environmental Prediction’s Gridpoint Statistical Interpolation and United Kingdom Met Office Unified Model 3D-Var system). However, a recent study by Xie and MacDonals [28] suggested that the use of stream function and velocity potential as momentum control variables might not be suitable for analysis over a regional domain for small-scale problems. Other mesoscale DA systems that involve radar observations, such as VDRAS and ARPS 3D-Var, used u-wind along x-grid direction and v-wind along y-grid direction as momentum control variables. An effort is being undertaken to compare the impact of the choice of momentum control variables on high-resolution analysis. Some preliminary results will be shown in Section 5.For control variables with option 5 (CV5) in WRFVAR, a recursive filter and EOFs (empirical orthogonal functions) are used to model background error covariance. The horizontal autocovariance, which is assumed to be spatially homogeneous and isotropic, is modeled by the recursive filter. The vertical error correlations are modeled using the diagnosed EOFs (empirical orthogonal functions) from estimated background errors. Since convective-scale data assimilation using radar observations is commonly done on limited area domains, it is preferable to derive background error statistics over the domain of interest as well as the season of interest. There are two approaches to estimate the background error statistics of the control variables in WRFVAR: the ensemble method [20, 29] and the NMC method [30]. Both methods are included in the utility GEN_BE [15] of the standard released WRF data assimilation system. Both algorithms were used in previous WRFVAR studies.An example of the estimated background error statistics is given in Figure1, which shows the horizontal length scales of the control variables with respect to the vertical EOF modes. The error statistics is computed on a domain covering a large part of the USA (N24°–N49° and W116°–W74°) using the NMC method and the June 2012 data from the real-time WRF 3 km runs. The plot shows that the length scales are a few hundred kilometers for the dynamical variables and temperature in the first few modes and reduced to below 100 kilometers in high frequency modes. The length scale of relative humidity is below 50 kilometers for all EOF modes.Figure 1
Length scales of forecast errors of the WRFVAR control variables of stream function (ψ), unbalanced velocity potential (χu), unblanced temperature (Tu), and relative humidity (rh), with respect to the vertical EOF mode, computed using 3 km WRF forecasts over a domain that covers a region of N24°–N49° and W116°–W74° in the USA. Note that the surface pressure has only one mode (indicated by the purple *) because it is a 2D field.


### 2.4. Automated Radar Data Quality Control

An automated radar data quality control and error estimate system that can be executed in real-time is an essential component of an operational radar DA system. The automated radar data preprocessing system embedded in NCAR’s real-time 4D-Var system VDRAS that is based on a cloud-scale model (different from WRF) [6] produces quality controlled radar data for use by WRFVAR. The system along with its radar data quality control package has been used in real-time since 2001. The VDRAS data preprocessing system includes several quality control algorithms as well as filtering, superobbing, and error computation. In Sun [31], some of these preprocessing methods were described. One of the challenging issues in the use of radar radial velocity for various applications is the data quality control of aliased velocities, known as velocity dealiasing. Radial velocity is aliased when the true radial velocity is larger than the unambiguous velocity of a radar or Nyquist velocity [25]. A unique feature of the VDRAS data preprocessing system is a velocity dealiasing algorithm coupled with the 4D-Var radar data assimilation in which the frequently updated wind analysis through a continuous cycling is used as the reference wind for the dealiasing. Lim and Sun [32] demonstrated that the coupling of the velocity dealiasing algorithm with the radar DA system improved the accuracy of the reference wind and hence resulted in an improved dealiased velocity field by significantly reducing the percentage of the improperly dealiased data from using a reference wind field that had no radar data assimilation.The velocity dealiasing can be particularly challenging during hurricane events because of the strong wind and the large error in the reference wind, especially near the center where a slight offset of center location can cause a large difference in wind direction. Recent applications of the coupled dealiasing scheme in VDRAS suggest that the scheme is effective to automatically unfold severely aliased velocity observations. Figure2 shows an example from a large number of radar volumes that were successfully dealiased for Typhoon Morakot that affected Taiwan during 7–9 August 2008.Radial velocity observations at 0.5° elevation from RCHL radar of the Central Weather Bureau (CWB) of Taiwan before (a) and after (b) the dealiasing scheme in VDRAS was applied. The thin white contours are costal lines.
(a)
(b)


## 2.1. Fundamental Framework of WRFVAR

The WRFVAR follows the incremental variational formulation [22] that is commonly used in operational systems. In the standard variational formulation, the optimal analysis is obtained by minimizing a cost function measuring the misfit between observation and prediction by a nonlinear model. In the incremental variational formulation, however, a linearized model is used to generate the prediction needed in the cost function. The advantage of the incremental approach is that it not only reduces the computational cost of 4D-Var but also improves the mathematical conditioning of the cost function (i.e., a quadratic cost function) because of the linearization of the forward operator. It facilitates the use of a coarser resolution tangent linear model (TLM) and adjoint model (ADM) to reduce the computational cost and a simplified representation of physical processes in TLM and ADM in the data assimilation cycle.The cost function of the incremental formulation can be derived from the standard cost function with some assumptions that will be described later. First we express the standard cost function as follows:(1)J(x0)=12(x0-x0b)TB-1(x0-x0b)+12∑k=0K{yko-Hk[Mk(x0)]}T×R-1{yko-Hk[Mk(x0)]}.
This cost function assumes that the assimilation window covers K observation windows with each represented by the subscript k. The variables of x0,  x0b, and yko represent the initial atmospheric state, the background state which can be previous model forecast, and the observed state, respectively. Mk is the nonlinear prediction model to propagate the initial atmospheric state to that at the kth observation time in the case of 4D-Var. It should be noted that 3D-Var is regarded as a special case in which K can simply set to 0 and Mk set to the identity matrix. Hk is the nonlinear observation operator.  B  and  R  are the background and observation error covariance matrices, respectively.To derive the cost function of the incremental formulation, we first introduce the innovation variabled at the kth observation time as
(2)dk=yko-Hk[Mk(x0)]
and then linearize the operators Hk and Mk in (1) as follows:
(3)Hk[Mk(x0n)]≈Hk[Mk(x0n-1)]+HkMk[x0n-1]δx0n,
where the superscripts n and n-1 represent the current and previous outer loop iterations, respectively. Hk and Mk are the tangent linear operator of Hk and Mk. Using (2) and (3), the cost function (1) can be changed to
(4)Jn[δx0n]=12{δx0n-[x0b-x0n-1]}TB-1{δx0n-[x0b-x0n-1]}+12∑k=0K{HkMk[x0n-1]δx0n-d0n-1}T×R-1{HkMk[x0n-1]δx0n-dkn-1}.
When the increment δx0n at the nth outer loop iteration is obtained, the estimate of the atmospheric state is updated by x0n=x0n-1+δx0n and used to produce the first guess trajectory Mk[x0n] for the next outer loop n+1. The basic assumption of the incremental approach is that the solution of the cost function (4) gradually approaches that of the nonlinear cost function (1) given enough number of outer loop iterations.Assuming the background error covariance matrix is given byB=UUT and a control variable transform by δx0=Uv, the background term of the cost function (4) is simplified to (1/2)(∑i=1Nvi)T(∑i=1Nvi), where N is the total number of outer loops, whereby avoiding the difficulty in computing the inverse of B. The control variable transform δx0=Uv is implemented through a series of operations U=UpUvUh [15]. The horizontal transform Uh is performed using a recursive filter [24]. The vertical transform Uv applies an EOF (empirical orthogonal function) decomposition on the vertical component of the background error covariance. The physical transform Up converts the increment in control variable space to analysis variable space. The control variables of the WRFVAR system will be described later in this section.The optimal analysis is obtained by minimizing the cost function (4) by iterating in both the outer and inner loops. The number of inner loop iterations is determined either by a convergence criterion or a prespecified fixed number. The number of the outer loop iterations is typically 2 to 6. The gradient of the cost function required in the minimization process is computed by the adjoint method. For the 4D-Var system with the incremental formulation, a forward model and its TLM and ADM are needed. The forward model is the same as the WRF forecast model but currently only with the physical schemes of the Kessler microphysics and the diffusion. More schemes will be added to the forward model in the future as the adjoint of these schemes has been developed.


## 2.2. Radar Observation Operators

To assimilate radar observations, the following observation terms are added to the existing cost function:(5)J=Jold+12∑k=0k[(Vrk-Vrkob)TRv-1(Vrk-Vrkob)nnnnnnnn+12(Zk-Zkob)TRz-1(Zk-Zkob)(Vrk-Vrrkob)T],
where Jold is used to represent the existing cost function before radar data assimilation is developed. The variables Vr stand for the radial velocity and Z for the reflectivity factor. The superscript “ob” indicates the observations. The symbols Rv-1 and Rz-1 stand for observation error covariance matrices for radial velocity and reflectivity, respectively. Note that the summation over the observation time levels k is not needed in the case of 3D-Var. The observation operator Hk in the cost function (1) links the model variables in a model coordinate to the observation variables in an observation space. For the radar radial velocity, this linkage is formulated with the 3D wind field (u, v, and w), the hydrometeor fall speed Vt, and the distance D between the location of a data point and the radar antenna:
(6)Vr=1D[(xd-xr)u+(yd-yr)v+(zd-zr)(w-Vt)],(7)D=(xd-xr)2+(yd-yr)2+(zd-zr)2,
where (xd,yd,zd) represents the location of the observation point and (xr,yr,zr) represents the location of the radar station. Vt is calculated from the rainwater mixing ratio with a height correction following Sun and Crook [6]. Relation (6) is linear except that Vt is nonlinearly dependent on qr, which needs to be linearized to obtain the linear observation operator for Vr. Note that the earth curvature effect [25] must be considered when mapping the radar observations to the model grid, which is done in the radar data preprocessing package that is described in the next subsection.The formulation of the reflectivity operator is not as straightforward because it depends on the assumption of drop size distribution in a microphysical parameterization scheme and the classification of hydrometeors. Following Sun and Crook [6], Xiao et al. [19] used the following relation for WRF 3D-Var by assuming the Marshall-Palmer drop size distribution for rain:
(8)Z=43.1+17.5log10(ρqr),
where ρ is the air density (kg m−3) and qr is the rainwater mixing ratio (g kg−1). Zhang et al. [26] examined a different relation by assuming a constrained gamma drop size distribution and found some improvement in the analysis of low precipitation.Equation (8) was used in the cost function (5) to assimilate reflectivity in the WRF 3D-Var developed by Xiao et al. [19]. However, Wang et al. [27] showed that the linearization of (8) as required by the incremental formulation of the cost function can result in a dry bias in rainwater analysis. Thus they proposed to indirectly assimilate the derived rainwater mixing ratio qr from (8) by replacing Z with qr in (5).


## 2.3. Control Variables and Background Error Estimate

The standard control variables (CV) [15] for WRF 3D-Var and 4D-Var are the stream function ψ, the unbalanced component of velocity potential χu, the unbalanced component of temperature Tu, the unbalanced component of surface pressure Psu, and “pseudo” relative humidity Rhs. The unbalanced component for each of the three variables is the residual after subtracting the respective correlated component with ψ which is obtained by a statistical fitting. The “pseudo” relative humidity is the water vapor mixing ratio divided by its saturated value in the background state. To assimilate radar radial velocity and reflectivity associated with warm rain, vertical velocity and the microphysical parameters of cloud water and rainwater must be added as control variables. In the current WRFVAR system, no cross-correlation with other variables is considered for relative humidity, cloud water mixing ratio, rain water mixing ratio, and vertical velocity. The 4D-Var technique through the use of the WRF model can implicitly produce multivariate correlations for the analysis.One issue that raises concerns for the convective-scale DA is the selection of momentum control variables. The use of stream function and velocity potential (ψ and χ) was widely accepted by global operational models (e.g., National Center of Environmental Prediction’s Gridpoint Statistical Interpolation and United Kingdom Met Office Unified Model 3D-Var system). However, a recent study by Xie and MacDonals [28] suggested that the use of stream function and velocity potential as momentum control variables might not be suitable for analysis over a regional domain for small-scale problems. Other mesoscale DA systems that involve radar observations, such as VDRAS and ARPS 3D-Var, used u-wind along x-grid direction and v-wind along y-grid direction as momentum control variables. An effort is being undertaken to compare the impact of the choice of momentum control variables on high-resolution analysis. Some preliminary results will be shown in Section 5.For control variables with option 5 (CV5) in WRFVAR, a recursive filter and EOFs (empirical orthogonal functions) are used to model background error covariance. The horizontal autocovariance, which is assumed to be spatially homogeneous and isotropic, is modeled by the recursive filter. The vertical error correlations are modeled using the diagnosed EOFs (empirical orthogonal functions) from estimated background errors. Since convective-scale data assimilation using radar observations is commonly done on limited area domains, it is preferable to derive background error statistics over the domain of interest as well as the season of interest. There are two approaches to estimate the background error statistics of the control variables in WRFVAR: the ensemble method [20, 29] and the NMC method [30]. Both methods are included in the utility GEN_BE [15] of the standard released WRF data assimilation system. Both algorithms were used in previous WRFVAR studies.An example of the estimated background error statistics is given in Figure1, which shows the horizontal length scales of the control variables with respect to the vertical EOF modes. The error statistics is computed on a domain covering a large part of the USA (N24°–N49° and W116°–W74°) using the NMC method and the June 2012 data from the real-time WRF 3 km runs. The plot shows that the length scales are a few hundred kilometers for the dynamical variables and temperature in the first few modes and reduced to below 100 kilometers in high frequency modes. The length scale of relative humidity is below 50 kilometers for all EOF modes.Figure 1
Length scales of forecast errors of the WRFVAR control variables of stream function (ψ), unbalanced velocity potential (χu), unblanced temperature (Tu), and relative humidity (rh), with respect to the vertical EOF mode, computed using 3 km WRF forecasts over a domain that covers a region of N24°–N49° and W116°–W74° in the USA. Note that the surface pressure has only one mode (indicated by the purple *) because it is a 2D field.


## 2.4. Automated Radar Data Quality Control

An automated radar data quality control and error estimate system that can be executed in real-time is an essential component of an operational radar DA system. The automated radar data preprocessing system embedded in NCAR’s real-time 4D-Var system VDRAS that is based on a cloud-scale model (different from WRF) [6] produces quality controlled radar data for use by WRFVAR. The system along with its radar data quality control package has been used in real-time since 2001. The VDRAS data preprocessing system includes several quality control algorithms as well as filtering, superobbing, and error computation. In Sun [31], some of these preprocessing methods were described. One of the challenging issues in the use of radar radial velocity for various applications is the data quality control of aliased velocities, known as velocity dealiasing. Radial velocity is aliased when the true radial velocity is larger than the unambiguous velocity of a radar or Nyquist velocity [25]. A unique feature of the VDRAS data preprocessing system is a velocity dealiasing algorithm coupled with the 4D-Var radar data assimilation in which the frequently updated wind analysis through a continuous cycling is used as the reference wind for the dealiasing. Lim and Sun [32] demonstrated that the coupling of the velocity dealiasing algorithm with the radar DA system improved the accuracy of the reference wind and hence resulted in an improved dealiased velocity field by significantly reducing the percentage of the improperly dealiased data from using a reference wind field that had no radar data assimilation.The velocity dealiasing can be particularly challenging during hurricane events because of the strong wind and the large error in the reference wind, especially near the center where a slight offset of center location can cause a large difference in wind direction. Recent applications of the coupled dealiasing scheme in VDRAS suggest that the scheme is effective to automatically unfold severely aliased velocity observations. Figure2 shows an example from a large number of radar volumes that were successfully dealiased for Typhoon Morakot that affected Taiwan during 7–9 August 2008.Radial velocity observations at 0.5° elevation from RCHL radar of the Central Weather Bureau (CWB) of Taiwan before (a) and after (b) the dealiasing scheme in VDRAS was applied. The thin white contours are costal lines.
(a)
(b)


## 3. Review of the 3D-Var Performances

The performance and impact of radar DA using WRF 3D-Var were examined through studies of individual cases as well as multiple cases/days under operational settings. The first case study with radar radial velocity data assimilation in WRF 3D-Var was reported by Xiao et al. [18] using a heavy rainfall case that occurred over the Korean Peninsula. They found that the assimilation of the radial velocity data from Jindo radar operated by the Korean Meteorological Administration improved the wind analysis (especially the northward component) resulting in substantial improvement for 3 h rainfall forecasts and noticeable improvement for 6 h forecasts. In a follow-up paper, Xiao et al. [19] described a methodology to assimilate radar reflectivity data into the WRF 3D-Var. Although the reflectivity is linked to the rain water mixing ratio, they chose the total liquid water qt as the control variable instead of qr not only because the distribution of total water is more Gaussian-like but also the multivariate analyses can be achieved by use of a partition scheme. In this scheme, the analysis increment of qt is partitioned into the water vapor mixing ratio qv, the cloud water qc, and the rain water qr using a warm rain microphysical scheme. The method was first tested with the landfalling typhoon Rusa (2002) over South Korea. Their results indicated that the positive impact of radar reflectivity lasted only for 3 hours. However, when both radial velocity and reflectivity radar observations were assimilated, noticeable improvement lasted up to 12 hours.The above 3D-Var radar data assimilation system was further evaluated by Pu et al. [33] in a case study of the Hurricane Dennis near landfall that occurred during 2005 in southeastern Cuba. In this study, the eastward and northward velocities synthesized from an airborne Doppler radar as well as reflectivity were assimilated into the WRF 3D-Var. They demonstrated that the radar data significantly improved the short-term (18 hours) forecasts of the intensity, track, and precipitation. However, similarly as in Xiao et al. [19], the impact was mainly attributed to the assimilation of the Doppler wind observations.The relatively small impact of the reflectivity data in the scheme developed by Xiao et al. [19] motivated a further examination of the reflectivity assimilation. Wang et al. [27] pointed out some issues with this scheme, namely, the use of a linearized form of the Z-qr equation (8) as the observation operator and the warm-rain partition scheme. As aforementioned, the WRF 3D-Var applies the incremental approach that requires the linearization of the forward model. The logarithm Z-qr equation (8) has a high degree of nonlinearity especially when qr is small. They found a dry bias was introduced because of the errors in the linearization. In addition, the performance of the partition scheme largely depends on quality of relative humidity in the first guess state because the switches that initiate the warm rain processes may never be turned on during the minimization of the cost function. Moreover, as noted by Xiao et al. [19], the warm-rain partition scheme is not suitable for winter storms in which cold rain process may play an important role. Wang et al. [27] described a new scheme for the assimilation of reflectivity in which the derived rain water mixing ratio from (8) was assimilated instead of the direct assimilation of the reflectivity observations. In addition, an extra observation term was added to the cost function (5) that measures the misfit between the water vapor mixing ratios from the model and an estimated field from radar reflectivity observations. The estimated water vapor was obtained by assuming saturation where radar reflectivity is higher than a specified threshold above the cloud base. In this new scheme, the control variables for the microphysics are pseudorelative humidity and rain water mixing ratios. The partition scheme used in Xiao et al. [19] is not needed.The performance of the new reflectivity assimilation scheme was demonstrated using an operational WRF 3D-Var DA and forecasting system with rapid update cycles by Wang et al. [27]. Four heavy-rain-producing convective cases that occurred during summer 2009 in Beijing, China, were studied using the new system. Results showed that the indirect assimilation of reflectivity significantly improved the short-term precipitation forecast skill up to 6 hours, extended from the 3 hours shown in Xiao et al. [19]. Figure 3 compares the 4-hour forecasts from two experiments that assimilate radial velocity alone and both radial velocity and reflectivity, respectively, for one of the four cases studied in Wang et al. [27]. It is evident from this comparison that the addition of the reflectivity data assimilation improves the forecast by making the area of the precipitation more compact. The improved positive impact of the new reflectivity assimilation scheme was also confirmed by Sun et al. [34] by running the 3D-Var over a consecutive one-week period over the US Great Plains region.Hourly accumulated precipitation (unit: mm) for the convective case that occurred on July 23, 2009, in Beijing, China. The 3-hour (left column) and 4-hour (right column) forecasts are initialized at 0600 UTC by a 3 hourly cycled WRF 3D-Var. The radar QPE ((a) and (b)) is used to verify the experiments that assimilate radar radial velocity ((c) and (d)) only and that assimilate both radial velocity and rain water mixing ratio derived from reflectivity ((e) and (f)) (from [27]).
(a)
(b)
(c)
(d)
(e)
(f)All of the studies reviewed previously assimilated radar observations with one cycle or multiple 3 hourly cycles. Sun et al. [34] demonstrated that the improvement of the short-term precipitation forecast could be achieved by running the WRF 3D-Var with 3 hourly rapid update cycles even without radar observations, although adding radar data further improved the forecasts. It is desirable to run 3D-Var with more frequent cycles (less than one hour) for nowcasting applications, given the fact that radars provide frequent observations. However, initial experiments yielded mixed results with hourly cycles. Digital filter initialization [35] was used in an attempt to reduce noises that may have been the reason for the difficulty with the hourly cycles, but no significant impact was found. Some technical issues need to be resolved and will be discussed in the last section.


## 4. Some Early Results of the 4D-Var System

The ultimate goal of the WRFVAR radar data assimilation system is the development of a 4D-Var scheme that is capable of using the frequent radar observations at multiple time levels. In their OSSE (Observation System Simulation Experiments) study, Sugimoto et al. [36] showed that the 3D-Var had only limited ability in retrieving the tangential wind component that was not observed by radar. Using the 4D-Var of a cloud-scale model, Sun and Crook [37], Sun and Crook [38], and Crook and Sun [39] all showed that the 4D-Var technique was able to retrieve the tangential wind with reasonable accuracy. We therefore anticipate that the WRF 4D-Var can improve the 3D wind retrieval and hence the convective forecasting over its 3D-Var.The 4D-Var technique is used in the global operational systems by several major operational weather service centers throughout the world. The European Center for Medium-Range Weather Forecasts (ECMWF) is the first center that implemented a 4D-Var operational system [40–42]. Following the success of ECMWF, other 4D-Var systems based on operational NWP models were developed and implemented at several operational centers throughout the world as their data assimilation systems for global or regional NWP without the convective-scale radar data assimilation. Recently the Met Office has run a real-time demonstration with an hourly cycling 4D-Var system and a 1.5 km resolution over a domain covering Southern England and Wales. This system currently assimilates radar radial velocity in the 4D-Var, but the reflectivity is assimilated with a diabatic initialization following Jones and Macpherson [43].The basic framework of the WRF 4D-Var was described by Huang et al. [17] and tested using a Typhoon case. Recently, the TLM and its ADM of the system were upgraded to match the updated WRF version 3.3. These upgrades along with some other software enhancements were described in Zhang et al. [44]. Wang et al. [20] described the latest development of the system in regard to radar data assimilation and showed the results of a case study from a squall line precipitation event. Using the same case, Sun and Wang [21] compared the performance of the 4D-Var with that of 3D-Var and found that the precipitation forecast initialized by the 4D-Var noticeably improved over that of 3D-Var.One question sought to answer by Wang et al. [20] was whether the incremental formulation would be feasible for the application at the convective-scale. They concluded that the incremental formulation worked well in the context of radar data assimilation with a high-resolution configuration when the typical assimilation window lengths of less than 30 minutes were used. Figure 4 shows the reduction of the cost function (5) when the window lengths of 5, 10, and 30 minutes are used. Note that the frequency of the data is the same (5 min) in all the three experiments. In these experiments, six outer loop iterations for the nonlinear base state update were used, and 20 inner loop iterations for the minimization of the cost function were performed. It is seen that the cost function is smoothly reduced when the 5-minute window was used. Some discontinuities are shown as the base state is updated in the case of 15- and 30-minute window lengths, but the cost function is steadily reduced at the end of each outer loop from the previous loop. Figure 4 also suggests that the 30-minute assimilation window results in a largest reduction of the cost function despite of the discontinuity caused by the mismatch between the linearized forward model and the nonlinear base state. Wang et al. [20] showed that the forecast initialized by the 4D-Var with the 30-minute window produced precipitation with smallest bias.Figure 4
Cost functions with respect to the number of total iterations accumulated by the 6 outer loops indicated on top of the solid curve. The three cost functions are from three experiments that use the assimilation window length of 5 minutes (4DV_T05), 15 minutes (4DV_T15), and 30 minutes (4DV_T30) (from [20]).The performance of WRF 4D-Var was compared with 3D-Var by Sun and Wang [21] in a case study. They showed that the 4D-Var substantially improved the precipitation skill over the 3D-Var for the 0–6 h forecasts of the squall line studied. Figure 5 compares the 6-hour precipitation forecasts (hourly accumulation) of a 4D-Var experiment with two 3D-Var experiments from their study. The difference between the two 3D-Var experiments 3DV and 3DVQV is that the latter accounts for the effect of latent heat by assimilating saturated in-cloud humidity as described in Wang et al. [20, 27]. The improvement of the precipitation forecast from the 4D-Var over the basic 3D-Var (Figure 5(b)) is clearly shown when compared with the Stage IV data (Figure 5(a)). Although the line structure is forecasted in the experiment 3DVQV (Figure 5(c)), the squall line has slower propagation speed and the precipitation is overforecasted. It is clearly seen that the 4D-Var experiment (Figure 5(d)) has the closest agreement with the observations. The qualitative verifications of the forecasts are compared in Figure 6 using the Fractions Skill Score (FSS) [45] with the radii of influence of 8 km, 24 km, and 48 km. The superior performance of the 4D-Var is shown by the improved FSS and the elimination of the initial skill drop that is present in the 3D-Var techniques. Sun and Wang [21] also found that the radial velocity observations had greater impact on the 4D-Var analyses and forecasts compared to the 3D-Var, suggesting that the 4D-Var has a better ability in retrieving the 3D wind as expected. Diagnostic analyses in their study suggested that the 4D-Var produced more realistic low-level cold pool, its leading edge convergence, and mid-level latent heating when compared with the radar analyses produced by VDRAS.6-hour forecasts of hourly accumulated precipitation (mm) from (b) the basic WRF 3D-Var (3DV), (c) the enhanced WRF 3D-Var (3DVQV), and (d) the 4D-Var (4DV). The Stage IV precipitation analysis is shown by (a) for verification.
(a)
Stage IV
(b)
3DV
(c)
3DVQV
(d)
4DVFractions scale score (FSS) versus forecast hour of hourly accumulated precipitation (mm h−1) for the thresholds 5 mm with the radius of influence (ROI) of 8 km (a), 24 km (b), and (c) 48 km, respectively.
(a)
5 mm ROI = 8 km
(b)
5 mm ROI = 24 km
(c)
5 mm ROI = 48 km


## 5. Future Developments for WRFVAR Radar DA

The encouraging results produced in the past studies warrant further development and improvement of the WRFVAR radar data assimilation system. An immediate next step is to investigate the impact of the analysis and forecast with regard to the selection of momentum control variables. One of the motivations for investigating the impact of the momentum control variables is that we have found that WRF 3D-Var tended to produce degraded results when an hourly cycle frequency was used. A new option that uses theu-wind and v-wind is being added to the WRFVAR system and compared with the existing scheme using ψ and χ momentum variables. Note that the unbalance component of the velocity potential is used in WRF 3D-Var although, for simplicity, we use the notation χ in the following description. Preliminary results indicate that the two options can yield significant differences in the characteristics of background error covariance and hence affect the final analysis. Figure 7 compares the analysis increments of u-wind and v-wind from single observation tests using the two momentum control variable options. These single observation tests were conducted by assuming an observation of u (left column) and v (right column) at the point marked by “+” in Figure 7 with an innovation and observation error of 1 ms−1. The background error covariance values for the two tests were calculated using an ensemble of short-term WRF forecasts. Apparently, the test using the ψ-χ control variables (Figures 7(a) and 7(d)) gives larger horizontal spread than the test using the u-v control variables (Figures 7(b) and 7(e)). By plotting the increments along the x-direction on the observation level (Figures 7(c) and 7(f)) it is clearly shown that the magnitudes of the wind increments are substantially different and much larger values have resulted from the use of the u-v option than the ψ-χ option. In addition, the ψ-χ control variable scheme produces negative increments in the distance far away from the observation point, which is questionable for the convective-scale. The impact of the different momentum control variables in the context of radar data assimilation for convective forecasting is being currently studied in more depth to investigate what are the causes of the different error characteristics and is also being studied using real observations.Figure 7
Vertical sections of analysis increments from experiments that assume a single observation ofu-wind (left columns) and v-wind (right columns) with 1 ms−1 innovation located at the point marked by “+” in (a) and (d). The u and v increments (unit: ms−1) resulting from the analyses that use ψ and χ as momentum control variables are shown in (a) and (d), respectively, and from the direct use of u and v as control variables that are shown in (b) and (e), respectively. The magnitudes of u and v increments on the observation level (6 km) are plotted as a function of x- and y-grid points. The horizontal grid spacing is 4 km.Although the initial results from Sun and Wang [21] are encouraging, it is obvious that the WRF 4D-Var should be examined with more cases to confirm its superiority over the 3D-Var. One important issue that was not yet examined is the continuous cycle in WRF 4D-Var. In both Wang et al. [20] and Sun and Wang [21], only one assimilation cycle was performed. Continuous cycles should be designed and tested before the system can be implemented operationally. One issue that requires attention is how to maintain the continuity between the assimilation cycle and the short-forecast cycle in between two assimilation cycles. Difficulty may arise from the use of different microphysics schemes during the assimilation period and the short-forecast period. Another obstacle for operational implementation of 4D-Var is its high computational cost. The incremental formulation of WRF 4D-Var can be implemented with a coarser resolution forward model for the data assimilation while the high-resolution nonlinear model is used for the outer loop background run as well as the forecast run. This strategy will be tested in the near future, and its performance both in terms of computational efficiency and analysis/forecast quality will be evaluated.Another planned future development is the use of dual-polarization observations to improve the initialization of microphysics variables. Li and Mecicalski [46] examined the impact of assimilating an improved estimate of the rain water mixing ratio through the use of differential reflectivity and specific differential phase. Some impacts were found both in the analyses and subsequent forecasts. A microphysics scheme that includes the ice phase will be developed to facilitate the assimilation and initialization of ice microphysics. Although the preliminary study by Li and Mecicalski [46] was encouraging, they raised the question how much additional information can be retrieved from the dual-polarization observations and to what extent it can improve the NWP model forecast. It is certainly an area requiring more studies.In recent years, a great deal of research has been devoted to the development and testing of methodologies of hybrid data assimilation that make use of both variational and ensemble data assimilation techniques. Wang et al. [47] developed a hybrid ensemble 3D-Var system for the WRF model. This system was applied to a hurricane prediction, and promising results were obtained [48]. Although the hybrid concept is applicable to WRF 4D-Var, its effectiveness is yet to be tested with radar data assimilation for convective forecasting.



---

*Source: 815910-2013-07-28.xml*