Articles | Volume 12, issue 9
Geosci. Model Dev., 12, 4031–4051, 2019
Geosci. Model Dev., 12, 4031–4051, 2019

Model description paper 13 Sep 2019

Model description paper | 13 Sep 2019

A radar reflectivity operator with ice-phase hydrometeors for variational data assimilation (version 1.0) and its evaluation with real radar data

A radar reflectivity operator with ice-phase hydrometeors for variational data assimilation (version 1.0) and its evaluation with real radar data
Shizhang Wang1,2 and Zhiquan Liu2 Shizhang Wang and Zhiquan Liu
  • 1Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters, Key Laboratory of Meteorological Disaster of Ministry of Education, Nanjing University of Information Science and Technology, Nanjing, 210044, China
  • 2National Center for Atmospheric Research, Boulder, CO 80301, USA

Correspondence: Zhiquan Liu (


A reflectivity forward operator and its associated tangent linear and adjoint operators (together named RadarVar) were developed for variational data assimilation (DA). RadarVar can analyze both rainwater and ice-phase species (snow and graupel) by directly assimilating radar reflectivity observations. The results of three-dimensional variational (3D-Var) DA experiments with a 3 km grid mesh setting of the Weather Research and Forecasting (WRF) model showed that RadarVar was effective at producing an analysis of reflectivity pattern and intensity similar to the observed data. Two to three outer loops with 50–100 iterations in each loop were needed to obtain a converged 3-D analysis of reflectivity, rainwater, snow, and graupel, including the melting layers with mixed-phase hydrometeors. It is shown that the deficiencies in the analysis using this operator, caused by the poor quality of the background fields and the use of the static background error covariance, can be partially resolved by using radar-retrieved hydrometeors in a preprocessing step and tuning the spatial correlation length scales of the background errors. The direct radar reflectivity assimilation using RadarVar also improved the short-term (2–5 h) precipitation forecasts compared to those of the experiment without DA.

1 Introduction

Over the past several decades, radar reflectivity observations have been used in many data assimilation (DA) studies (Borderies et al., 2018; Caumont et al., 2010; Gao and Stensrud, 2012; Hu et al., 2006; Jung et al., 2010, 2008a; Liu et al., 2019; Putnam et al., 2014; Snook et al., 2012, 2015; Sun and Crook, 1997; Sun and Wang, 2013; Tong and Xue, 2005; Wang et al., 2013b; Wang and Wang, 2017; Wattrelot et al., 2014; Xiao et al., 2007; Xue et al., 2006) and they have demonstrated that assimilating this radar reflectivity improves the initial conditions of the convective scale and benefits the subsequent forecasts. To assimilate the reflectivity, it is necessary to transform the model's prognostic variables (e.g., rainwater, snow, and graupel) to the observed radar reflectivity. To perform this transformation, early studies (e.g., Sun and Crook, 1997; Xiao et al., 2007) used the Marshall–Palmer distribution of raindrop size (ZR relationship). However, this relationship is only valid in precipitation areas without ice-phase species; thus, its applications (e.g., Schwitalla and Wulfmeyer, 2018) are often limited to layers lower than 4 km or 8 km above ground level (a.g.l.). To overcome this deficiency, more comprehensive operators that involve snow and graupel have been developed (Gao and Stensrud, 2012; Tong and Xue, 2005). Several studies (e.g., Gao and Stensrud, 2012; Wang and Wang, 2017) have demonstrated that involving these ice species in the reflectivity operator improves the analysis of hydrometeors in terms of their spatial distribution, especially in the vertical direction.

Although these operators have been successfully applied in many convective-scale DA studies, they were developed for a specific band of radar (e.g., S-band radar; Gao and Stensrud, 2012; Sun and Crook, 1997) and specific microphysics characteristics (e.g., with a fixed intercept parameter; Sun and Crook, 1997). However, mixed-phase species such as wet snow and wet graupel have not been considered in these operators. Recently, the contributions from mixed-phase species have been studied (Jung et al., 2008a, hereafter, J08; Posselt et al., 2015). To compute the mixed-phase species' contributions, J08 proposed an operator that was based on the expressions given by Zhang et al. (2001). Their expressions were derived according to the scattering amplitudes that were estimated through the T-matrix method and the Rayleigh scattering approximation (J08). For computational efficiency, these expressions were rewritten in the polynomial form that was only valid for S-band radar in which the Rayleigh assumption was satisfied. Later, a more general and exact operator that fully used the T-matrix scattering method was proposed (Jung et al., 2010). This operator was given as the integral of the complex backscattering amplitudes over the size distribution of the precipitation particles (i.e., rainwater, snow, and graupel). In addition to these operators, several complex reflectivity operators in the integral form have also been proposed (Borderies et al., 2018; Caumont et al., 2006; Pfeifer et al., 2008; Ryzhkov et al., 2011; Wattrelot et al., 2014). Some were designed for a specific band of radar (e.g., W band; Borderies et al., 2018), whereas some were designed for the bin microphysics scheme (e.g., Ryzhkov et al., 2011).

Reflectivity operators have been developed both for the variational method (Caumont et al., 2010; Gao and Stensrud, 2012; Hu et al., 2006; Sun and Crook, 1997; Sun and Wang, 2013; Wang et al., 2013b; Wattrelot et al., 2014; Xiao et al., 2007) and for the ensemble Kalman filter method (EnKF; Dawson et al., 2010; Jung et al., 2010, 2008a, b; Putnam et al., 2014; Snook et al., 2011, 2015; Tong and Xue, 2005; Xue et al., 2006). The variational method requires the tangent linear (TL) and adjoint (AD) operators, which are not required by the EnKF (Evensen, 2003). Therefore, complex operators such as those proposed by J08 are often employed in EnKF DA applications. For the variational method, a common approach to avoid using the TL/AD operators is to assimilate the reflectivity-retrieved hydrometeor profiles (Caumont et al., 2010; Wang et al., 2013a; Wattrelot et al., 2014); an alternative is to use the reflectivity as an additional control variable with the ensemble–variational DA approach (Wang and Wang, 2017).

Despite the difficulty, some efforts have been undertaken for reflectivity assimilation with the TL/AD operators (Gao and Stensrud, 2012; Kawabata et al., 2018; Liu et al., 2019; Xiao et al., 2007), and reasonable results have been obtained in terms of hydrometeor analysis and precipitation forecasts. However, none of these studies employed operators as complex as those proposed by J08. Kawabata et al. (2018) adopted the expressions of Zhang et al. (2001) and developed the TL/AD operator for C-band radar but without taking into account the contributions from ice-phase species.

The main purpose of this study is to develop a TL/AD operator based on Jung et al. (2008a) with the contributions of ice-phase precipitation and apply it in a variational DA framework. For convenience, the operator implemented in this study is called RadarVar to represent that it was developed for variational DA and contains ice-phase species. The original J08 operator is called J08orig. The reminder of this paper is organized as follows. In Sect. 2, the J08 operator is reviewed, and its TL and AD operators are derived. The experimental design is given in Sect. 3, and the new operators are verified in Sect. 4. The performance of RadarVar is discussed in Sect. 5, and the conclusions are presented in Sect. 6.

2 Reflectivity operator

2.1 Review of the J08 operator

The radar-observed reflectivity, Z, is given in logarithmic form as

(1) Z = 10 log 10 Z e ,

where Ze is the equivalent reflectivity factor, which is the sum of the contributions from pure rainwater (Zr), dry snow (Zds), dry graupel (Zdg), wet snow (Zws), and wet graupel (Zwg) as follows:

(2) Z e = Z r + Z ds + Z dg + Z ws + Z wg .

To compute Eq. (2), the mixing ratios of mixed-phase species (wet snow and wet graupel) are required. However, many widely used microphysics schemes, such as the Lin, WSM6, and Morrison schemes, do not predict or diagnose the mixed-phase species; thus, the amount of rainwater in wet snow or graupel cannot be directly extracted from the model output. To solve this issue, J08 modeled the rain–snow (rain–graupel) mixture using a fraction that is given by

(3) F = [ min ( q r / q x , q x / q r ) ] 0.3 F max ,

where Fmax is the maximum fraction, which is 0.5 (0.3) for rain–graupel (rain–snow) mixtures; qr is the mixing ratio of rainwater; and qx is the general form of the mixing ratio of ice-phase species. The subscript “x” can be either “s” for snow or “g” for graupel. With this fraction, the mixing ratios of pure rainwater, dry snow, dry graupel, and mixed-phase species are given by

(4) q pr = ( 1 - F ws - F wg ) q r q ds = ( 1 - F ws ) q s q dg = ( 1 - F wg ) q g q ws = F ws ( q s + q r ) q wg = F wg ( q g + q r ) ,

where the subscripts “ws” and “wg” are added to F to represent the fractions of wet snow and wet graupel, respectively, and the subscripts “pr”, “ds”, and “dg” represent pure water, dry snow, and dry graupel, respectively. The mixed-phase density, ρwx, is not a constant and is parameterized by

(5) ρ w x = ( 1 - f w x 2 ) ρ x + f w x 2 ρ r ,


(6) f w x = q r q r + q x .

The subscript “x” in ρwx, ρx, and fwx represents either snow (s) or graupel (g), and fwx is called the water fraction.

2.1.1 Contribution from rainwater

In accordance with J08, all of the contributions are computed by integrations over the drop size distribution (DSD) weighted by the scatter cross section determined by the density, shape, and DSD. The DSD is modeled by an exponential distribution. After performing the integration, the contribution from pure rainwater, Zr, is written in a simple form (Zhang et al., 2001; Posselt et al., 2015; Kawabata et al., 2018) as follows:

(7) Z r = 4 λ 4 α ra 2 N 0 r π 4 K w 2 Λ r - ( 2 β ra + 1 ) Γ ( 2 β ra + 1 ) ,

where λ is the wavelength of the radar, which is 107 mm for S-band radar, and N0r is the intercept parameter of rainwater, which is 8×106 m−4 in this study. N0r values are typically fixed (or constant) in single-moment microphysics schemes; for a two-moment scheme, this value should be determined using the predicted number concentration. |Kw|2 is the dielectric factor for rainwater and is equal to 0.93, and αra and βra are 4.28×10-4 and 3.04, respectively. The complete gamma function is written as Γ(…). The slope parameter of rain, Λr, is

(8) Λ r = ( π ρ r N 0 r ρ a q pr ) 1 4 ,

where ρr=1000 kg m−3 is the rain density, qpr is given by Eq. (4), and ρa is the density of air. By substituting Eq. (8) and the constant parameters into Eq. (7), we can rewrite Eq. (7) as a function of qpr as follows:

(9) Z r ( q pr ) = P r ( q pr ) 1.77 ,


(10) P r = 4 λ 4 α ra 2 π 4 K w 2 ( π ρ r ρ a ) - 2 β ra + 1 4 ( N 0 r ) 1 - 2 β ra + 1 4 Γ ( 2 β ra + 1 ) .

The value of Pr is approximately 4.8×109 with an air density of 1.0 kg m−3. This value has the same magnitude as those proposed by Sun and Crook (1997) and Gao and Stensrud (2012).

2.1.2 Contribution from dry snow/graupel

The contributions from both dry and mixed-phase ice species after integration have the same form but differ in their coefficients. For dry ice species, the contribution is given by

(11) Z d x = 4 Γ ( 7 ) λ 4 N 0 x π 4 K w 2 Λ d x - 7 ( A α d x a 2 + B α d x b 2 + 2 C α d x a α d x b ) ,

where the subscript x represents either snow (s) or graupel (g). The intercept parameters of these species are denoted by N0x, the values of which are 3×106 and 4×105 m−4 for snow and graupel, respectively. Both values are consistent with the default values of Advanced Regional Prediction System (ARPS) EnKF where the J08 operator was implemented. The slope parameter in Eq. (11) for either dry snow or dry graupel is written as

(12) Λ d x = ( π ρ x N 0 x ρ a q d x ) 1 4 ,

where ρx is either the density of snow (100 kg m−3) or graupel (400 kg m−3). The parameters A, B, and C in Eq. (11) are functions of the mean (ϕ¯) and the standard deviation (σ) of the canting angle and are given by

(13) A = 1 8 ( 3 + 4 cos 2 ϕ ¯ e - 2 σ 2 + cos 4 ϕ ¯ e - 8 σ 2 ) B = 1 8 ( 3 - 4 cos 2 ϕ ¯ e - 2 σ 2 + cos 4 ϕ ¯ e - 8 σ 2 ) C = 1 8 ( 1 - cos 4 ϕ ¯ e - 8 σ 2 ) .

According to J08, ϕ¯ is zero for all hydrometeors, and σ differs for snow (20) and hail (60). Here, we assume that σfor graupel is also 60. The horizontal reflectivity that is regarded in this study is not sensitive to the canting angle (will be demonstrated in Sect. 2.4), although the differential reflectivity is sensitive to canting angle (Aydin and Seliga, 1984). In Eq. (11), A, B, and C are multiplied by coefficients (αdxa and αdxb) that describe the backscattering amplitudes. For dry ice-phase species, these coefficients are precalculated constants and are listed in Table 1. The coefficients of graupel are calculated using the backscattering amplitudes (for particle size <10 mm) in the pyCAPS-PRS v1.1 software (Dawson et al., 2014; Johnson et al., 2016; Jung et al., 2010; Jung et al., 2008a) provided by the Center for Analysis and Prediction of Storms (CAPS), and fitted to the polynomial function of fwg. αdga and αdgb are the coefficients when fwg is zero, which indicates no rainwater.

Table 1Values of αdsa, αdsb, αdga, and αdgb.

Download Print Version | Download XLSX

Table 2Values of Pwxak and Pwxbk in Eq. (16).

Download Print Version | Download XLSX

For brevity, Eq. (11) is rewritten as a function of qdx for dry snow and dry graupel and is given by

(14) Z d x ( q d x ) = P d x q d x 1.75 ,


(15) P d x = 4 Γ ( 7 ) λ 4 N 0 x - 0.75 π 4 K w 2 ( π ρ x ρ a ) - 1.75 ( A α d x a 2 + B α d x b 2 + 2 C α d x a α d x b ) .

2.1.3 Contribution from wet snow/graupel

The equation for the contribution from mixed-phase species has the same form as Eq. (11), except that the subscript “d” is replaced by “w” to represent wet species. The slope parameter for mixed-phase species also has the same form as Eq. (12), except that the subscript “d” is replaced by “w”, and ρwx is substituted for ρx. For wet species, σ in A, B, and C is a function of fw and qwg. Additional details are given in Sect. 3c of J08. The coefficients that are multiplied by A, B, and C are functions of fw and are written as

(16) α w x a = ε x k = 0 n P w x a k f w x k α w x b = ε x k = 0 n P w x b k f w x k ,

where “x” is “s” (g) for snow (graupel), εx is 10−4 (10−3) for snow (graupel), Pwxak and Pwxbk are precalculated constants for S-band radar, the value of n is 6, and the superscript k denotes the index of these constants. All these values are based on J08, except for those of graupel, which are computed using the same method mentioned in Sect. 2.1.2. The values of Pwxak and Pwxbk are listed in Table 2.

To simplify the derivation of the TL/AD operators, we rewrite Eq. (11) as a function of the mixed-phase mixing ratio (qwx) and the water fraction (fw) as follows:

(17) Z w x ( q w x , f w x ) = P w x q w x 1.75 ε x 2 k = 0 2 n P x k f w x k ,

where the coefficients Pwx and Pxk are given by

(18) P w x = 4 Γ ( 7 ) λ 4 N 0 x - 0.75 π 4 K w 2 ( π ρ w x ρ a ) - 1.75 ,


(19) P x k = A P A x k + B P B x k + 2 C P C x k ,

respectively. PAxk, PBxk, and PCxk in Eq. (19) are given by

(20) P A x k = i = 0 k P w x a i P w x a ( k - i ) ( k n ) i = k - n n P w x a i P w x a ( k - i ) , ( k > n ) P B x k = i = 0 k P w x b i P w x b ( k - i ) , ( k n ) i = k - n n P w x b i P w x b ( k - i ) ( k > n ) P C x k = i = 0 k P w x a i P w x b ( k - i ) , ( k n ) i = k - n n P w x a i P w x b ( k - i ) , ( k > n ) ,

where Pwxai and Pwxbi are precalculated constants in Eq. (16) and are listed in Table 2. The subscript “x” in these constant coefficients represents either snow (s) or graupel (g). The derivations of Eqs. (17) to (20) are given in the Appendix.

2.2 Tangent linear operator

Because RadarVar is highly nonlinear and complex, performing the derivation for every nonconstant variable in this operator is difficult. In addition, some components of RadarVar (e.g., the fraction F in Eq. 3) are discontinuous and may cause serious convergence problems in the minimization (Janisková and Lopez, 2013; Janisková et al., 1999). Although this issue can be addressed by performing regularization for the discontinuous components in RadarVar, it is beyond the scope of this study. Here, we assumed that five variables in RadarVar are not changed in the minimization: (i) the air density, (ii) the fraction F, (iii) the intercept parameter, (iv) the standard deviation of the canting angle σ, and (v) the density of the mixed-phase species. The air density is held constant for simplification and to focus on the impact of the changes in the hydrometeors. The fraction F is assumed to remain unchanged because of its second-order discontinuity at the point that qx is equal to qr. The intercept parameter is constant because RadarVar was currently designed for single-moment microphysics schemes. Although multi-moment microphysics schemes are more realistic because they predict the number density of hydrometeors, single-moment microphysics schemes are still widely used to compute reflectivity and the polarization variables (e.g., Posselt et al., 2015). For simplification, the standard deviation of the canting angle and the density of the mixed-phase species remain unchanged in the minimization. These assumptions are discussed in Sect. 4.

2.2.1 Linearization for rain

Considering the assumptions presented above, Pr in Eq. (9) becomes a constant in the minimization. Therefore, the linearized form of Eq. (9) for qpr is given by

(21) δ Z r ( q pr ) = Z r ( q pr ) q pr q pr q r δ q r = 1.77 P r ( q pr ) 0.77 ( 1 - F ws - F wg ) δ q r .

2.2.2 Linearization for dry and wet snow/graupel

For dry snow and graupel, the linearized form of Eq. (14) is given by

(22) δ Z d x ( q d x ) = Z d x q d x q d x q x δ q x = 1.75 P d x q d x 0.75 ( 1 - F w x ) δ q x .

The linearization of Eq. (17) can be categorized into two parts, which represent the variations in Zx caused by changes in qwx and fwx, respectively. The linear equation of Eq. (17) is written as

(23) δ Z w x ( q w x , f w x ) = Z w x q w x ( q w x q x δ q x + q w x q r δ q r ) + Z w x f w x ( f w x q x δ q x + f w x q r δ q r ) = ( Z w x q w x q w x q r + Z w x f w x f w x q r ) δ q r + ( Z w x q w x q w x q x + Z w x f w x f w x q r ) δ q x = { 1.75 P w x q w x 0.75 F w x ( ε x 2 k = 0 2 n P x k f w x k ) + P w x q w x 1.75 [ ε x 2 k = 0 2 n P x k k f w x k ( 1 q r - 1 q r + q x ) ] } δ q r + { 1.75 P w x q w x 0.75 F w x ( ε x 2 k = 0 2 n P x k f w x k ) + P w x q w x 1.75 [ ε x 2 k = 0 2 n P x k k f w x k ( - 1 q r + q x ) ] } δ q x ,

where the subscript “x” represents either snow (s) or graupel (g). Using s (g) to replace the subscript “x” in Eq. (23), we can obtain the tangent linear operator of the contribution from wet snow (wet graupel).

The linearization of Eq. (1) is given by

(24) δ Z = 10 1 Z e ln 10 δ Z e ,

where δZe is the sum of δZr, δZds, δZdg, δZws, and δZwg. Note that Ze cannot be zero. Moreover, a too-small Ze may result in an extremely large gradient during the minimization, which is undesirable. Therefore, the accepted minimum Ze is set to 1.0 in the current RadarVar. When the background value is smaller than this minimum value, DA system will discard the corresponding observation to prevent a large gradient. This minimum value can be smaller to ingest more observation and will be tuned in future work.

Figure 1(a) The reflectivity as (a) a function of qr (red), qs (green), and qg (blue) for pure water and dry snow/graupel and (b) a function of qrqs (colors) and qrqg (contours) for wet snow/graupel.


2.3 Adjoint operator

The adjoint operator is the transpose of the tangent linear operator. Because the tangent linear operator is applied to the model variables qr, qs, and qg, the adjoint operator is written for these variables. First, the adjoint operator is written for Eq. (24). This operator has the following form:

(25) δ Z e A = 10 1 Z e ln 10 δ Z ,

where the superscript “A” means adjoint.

For rainwater in Eq. (21), the adjoint operator is given by

(26) δ q r A = δ q r A + 1.77 P r ( q pr ) 0.77 ( 1 - F ws - F wg ) δ Z e A .

The parameter qrA on the right-hand side of Eq. (26) is the accumulated qrA before computing Eq. (26). This rule is also valid for qx.

The adjoint operator of Eq. (22) is given by

(27) δ q x A = δ q x A + 1.75 P d x q d x 0.75 ( 1 - F w x ) δ Z e A .

Because Eq. (23) is the derivation with respect to both qr and qx, the adjoint operator of Eq. (23) contains two parts: one for rainwater and the other for ice species. For rainwater involved in Eq. (23), the adjoint operator is given by

(28) δ q r A = δ q r A + 1.75 P w x q w x 0.75 F w x ( ε x 2 k = 0 2 n P x k f w x k ) δ Z e A + P w x q w x 1.75 [ ε x 2 k = 0 2 n P x k k f w x k ( 1 q r - 1 q r + q x ) ] δ Z e A .

For ice species in Eq. (23), the adjoint operator is given by

(29) δ q x A = δ q x A + 1.75 P w x q w x 0.75 F w x ( ε x 2 k = 0 2 n P x k f w x k ) δ Z e A + P w x q w x 1.75 [ ε x 2 k = 0 2 n P x k k f w x k ( - 1 q r + q x ) ] δ Z e A .

2.4 Sensitivity of RadarVar to changes in hydrometeors

Since RadarVar is highly nonlinear, understanding its response to changes in qr, qs, and qg assists in the analysis of DA results using this operator. The response functions for Eqs. (9), (14), and (17) are plotted in Fig. 1. Figure 1b is a two-dimensional plot because Eq. (17) involves two kinds of hydrometeors. Figure 1a shows that the reflectivity changes more rapidly for all three hydrometeors when the mixing ratios are less than 0.5 g kg−1. As the mixing ratios increase to 2.0 g kg−1 or greater, the relationship between the reflectivity and the mixing ratio is approximately linear. This feature indicates that the reflectivity is more sensitive at small mixing ratios, and it also implies that the tangent linear approximation may give larger errors when the background reflectivity is small.

Figure 2The simulation domain with radar sites marked by radar icons and names. The areas of precipitation greater than 5 mm h−1 are plotted for 00:00 Z (red), 02:00 Z (green), and 04:00 Z (blue) on 2 June 2018.

In addition, the reflectivity contribution from wet snow increases more substantially than that from dry snow when qs or qr are in the range of 0–0.5 g kg−1. Figure 1b shows that the reflectivity reaches 35 dBZ when both qs and qr are approximately 0.2 g kg−1, while this reflectivity value requires qs of ∼1.2 g kg−1 for dry snow. This result is expected because many observation studies (e.g., Zhang et al., 2008) have shown that wet snow causes a bright band (large reflectivity) in the melting layer. This result also implies that the approximation error in the melting layer could be large.

Figure 3(a–c) The standard deviations of the background errors at different vertical levels, (d–f) the horizontal correlation length scale as a function of EOF mode, and (g–i) the vertical correlation coefficients for qr, qs, and qg.


3 Data and experimental design

3.1 Case review

This study focuses on a precipitation case that occurred in the northern US. The precipitation initiated at approximately 21:00 Z on 1 June 2018, when convective cells formed near the border between South Dakota and Nebraska. By 00:00 Z on 2 June 2018, these cells had developed into a linear convective system that was approximately 300 km long in the northeast–southwest direction and stretched from south of South Dakota to northern Nebraska, as shown in Fig. 2. Note that there is also a weaker precipitation system near the north boundary of the domain. The top of the convective system at this time, identified by reflectivity greater than 5 dBZ, reached 16 km a.g.l. This convective line developed further and moved to southeast Nebraska from 00:00 to 04:00 Z. During this period, a bow echo was observed on the radar mosaic provided by the National Centers for Environmental Information (NCEI). By 08:00 Z, the convective line had moved to the northern border of Kansas, followed by a large area of stratiform clouds that covered eastern Nebraska. The precipitation caused by this convective system lasted for nearly 20 h and ended at approximately 18:00 Z on 2 June 2018.

3.2 Settings of the forecast model

Version 3.9.1 of the Weather Research and Forecasting (WRF; Skamarock et al., 2018) model was employed. All of the experiments were performed on a 450×450×42 domain centered at 41 N, 96 W, in eastern Nebraska (Fig. 2). The horizontal grid spacing was 3 km. The terrain-following vertical grid was employed with the model top at 50 hPa. All of the experiments used the same physical parameterizations: no cumulus parameterization, the Thompson microphysics scheme (Thompson et al., 2004), the Rapid Radiative Transfer Model for general circulation models (RRTMG) longwave and shortwave radiation scheme (Iacono et al., 2008; Mlawer et al., 1997), and the Unified Noah land-surface model (Chen and Dudhia, 2001). Note that in the current RadarVar implementation, the intercept parameters are fixed, while they spatiotemporally vary in the Thompson scheme. This inconsistency may increase the adjustment time for model initialization. However, this issue is secondary in the present work because no improvement in terms of forecast skill which will be introduced was found in our early tests using a single-moment microphysics scheme with the density and intercept parameters being identical to RadarVar. The primary DA issue in this study is the poor background quality (due to no hydrometeor in Global Forecast System (GFS) analysis or the precipitation displacement). The inconsistency between the operator and microphysics scheme will be considered in the future. The initial conditions and the lateral boundary conditions were generated with the GFS data at 00:00 Z on 2 June 2018.

3.3 Generation of the background error covariance

RadarVar was implemented in the WRF data assimilation (WRFDA) (Barker et al., 2012) system. To perform the variational DA with this newly developed radar operator, it is necessary to generate the background error covariance matrix with qr, qs, and qg as part of the control variables. The generalized software package for the background error covariance statistics (GEN_BE) developed by Descombes et al. (2015) was used. The GEN_BE package can generate the univariate background error statistics for 11 variables, including these three hydrometeors. The background error statistics were computed using the National Meteorological Center (NMC) method (Parrish and Derber, 1992), which uses pairs of the differences between the 12 and 24 h forecasts. In total, 27 d forecasts from 20 May to 15 June 2018 were employed to generate the background error covariance.

The background error statistics for qr, qs, and qg are shown in Fig. 3. The vertical distributions of the background error's standard deviation (Fig. 3a, b, c) are consistent with those of the corresponding hydrometeor profiles: the error of the rainwater mainly appears in the lower levels, while those of snow and graupel mainly exist in the upper levels. The graupel may fall from the upper levels into the lower levels, so the graupel error has a broad vertical distribution. The horizontal correlation length scales of the background errors are often less than four grids (<12 km), and the vertical correlation of each hydrometeor can be large at the associated precipitation levels. These spatial correlations of the hydrometeor errors determine the remote horizontal and vertical influences of the observed reflectivity.

3.4 Observation data and verification data

Radar data at 00:00 and 01:00 Z on 2 June 2018 were selected to evaluate the radar operator in a variational analysis framework because the convection was sufficiently deep to contain ice-phase species. To fully cover the convective system at 00:00 Z, data from KABR, KFSD, KLNX, KOAX, KUDX, and KUEX were used, and KLNX was the closest radar from the convective line (Fig. 2). These radars were located in Nebraska and South Dakota. The radar data were stored in level-II format and converted to WRFDA format using a modified 88D2ARPS package, which is widely used in radar DA studies (Putnam et al., 2014; Snook et al., 2011, 2012). During this conversion, the radar data were also horizontally remapped to the model grids but remained at the radar elevations in the vertical direction; in other words, the horizontal resolution of the radar data after the conversion was consistent with that of the model. To be consistent with the work using the J08 operator (Jung et al., 2008b), the observation error of the reflectivity was set to 2.0 dBZ. Our early tests using different observation errors indicated that the errors of the analysis reflectivity were comparable when using reflectivity errors ranging from 0.5 to 2.0 dBZ. For computational efficiency, we selected the remapped data every two grids in both the x and y directions for the DA.

The NCEP gridded stage IV (ST4) dataset (Lin and Mitchell, 2005) was used for the precipitation forecast verification. ST4 data with a horizontal resolution of 4 km were interpolated into the 3 km model grid mesh to evaluate the precipitation prediction. At each model grid, the interpolated value was the weighted average of the ST4 data within 10 km of the grid; these data were weighted by the square of the inverse of the distance between the model grid and the ST4 data location.

Figure 4Scatter plot of the reflectivity for J08orig (x axis) versus RadarVar (y axis). The bias, standard deviation (SD), and number of samples are listed in the plot.


Figure 5(a) Vertical distributions of the maximum absolute values of the perturbed qr (red), qs (green), and qg (blue). Reflectivity scatter plots of all_prt (x axis) versus (b) F_unprt (y axis), (c) rhom_unprt (y axis), and (d) SD_unprt (y axis) at model levels 1–12 (black, lower than 3 km a.g.l.), 12–20 (red, between 3 and 7 km a.g.l.), and 21–30 (green, above 7 km a.g.l.).


3.5 Experimental design

As the first attempt to implement and apply RadarVar in WRFDA, this study focused on the quality of the analysis using the univariate three-dimensional variational (3D-Var) DA method in terms of the root mean square difference (RMSD) against the observed reflectivity and the similarity between the observed reflectivity distributions and the analysis. The forecast performance is the secondary concern and will be explored more thoroughly in a future study with multivariate analysis using more advanced DA techniques.

All of the DA experiments analyzed only rainwater, snow, and graupel by assimilating only the radar reflectivity observations. The first DA experiment, called Exp_ref, is considered the benchmark experiment; it mostly used the default configurations of WRFDA-3D-Var, except the number of outer loops was set to six, and the number of maximum iterations in each outer loop was set to 150. More outer loops were necessary due to inaccurate background hydrometeors and the high nonlinearity of the radar operator. To determine the tradeoff between the analysis quality and computational cost, two variants of Exp_ref were conducted with 50 and 100 inner iterations. In each experiment, the radar DA analyses were performed at 00:00 and 01:00 Z. The background for the 00:00 Z analysis was interpolated from the GFS analysis with zeros for the hydrometeor fields, and the background for the 01:00 Z analysis was the 1 h WRF forecast from the 00:00 Z analysis with more realistic hydrometeor fields.

Note that TL/AD of RadarVar will not be able to create reflectivity increments with the zero-hydrometeor background (serves as the base state of TL/AD in the first outer loop), which is the case for the 00:00 Z analysis. An approach to address this issue is to reset the zero background values of qr, qs, and qg to small values that can range from 10−9 to 10−6 kg kg−1 (e.g., Wang and Wang, 2017). However, this approach will result in the fraction F being a constant, while it is expected that F will peak near the middle of the melting layer. Therefore, we introduced a hydrometeor preprocessing step before performing the analysis, which constructs a new background with the weighted sum of the radar-retrieved hydrometeor mixing ratios and their background counterparts. The hydrometeor retrieval followed the procedure that is available in WRFDA, which is based on Gao and Stensrud (2012). The weight coefficients are arbitrarily set to 0.1 for the retrieval part and to 0.9 for the background. The small weight for the retrieval part was used to minimize the impact of the retrieval, mainly to ensure a nonzero background. In addition, current RadarVar cannot work with the weight coefficient of the retrieval part being smaller than 6×10-4 if the background contains no hydrometeor. A larger weight coefficient of the retrieval part (e.g., 0.5) reduces the difference between the background and the observation, which weakens the impact of direct DA using RadarVar and is contradictory to the purpose of this study. The hydrometeor preprocessing was performed at both 00:00 and 01:00 Z in Exp_ref as a reference.

Table 3Perturbation samples for the verification of the tangent linear operator.

Download Print Version | Download XLSX

Figure 6Scatter plots of the observed (x axis) versus (a, c) background and (b, d) analysis reflectivity at (a, b) 00:00 Z and (c, d) 01:00 Z in Exp_ref. The mean bias and SD between the observations and the background (or analysis) are listed in each plot.


Figure 7The (a, e, i, c, g, k) observed (KLNX) and (b, f, j, d, h, l) analysis reflectivities at (a–d) 2 km, (e–h) 4 km, and (i–l) 6 km a.g.l. at (a, e, i, b, f, j) 00:00 Z and (c, g, k, d, h, l) 01:00 Z in Exp_ref.


To examine the analysis performance with a bad background, the Exp_ref analysis at 00:00 Z was also run with a very small retrieval weight of 6×10-4. In addition, the impact of the hydrometeor preprocessing on a relatively “good” background (01:00 Z) was examined by comparing it with an experiment without hydrometeor preprocessing.

In several previous studies (Ban et al., 2017; Choi et al., 2017; e.g., Shen and Min, 2015), the horizontal correlation length scale factor of the background error had a large impact on the analysis. Therefore, two additional experiments, Exp_ls0.5 and Exp_ls0.125, were performed at 00:00 Z with the length scale factors of 0.5 and 0.125, respectively.

Short-term forecasts from some of these analyses and from the GFS analysis (referred to as the “noDA” experiment) were also carried out and evaluated.

4 Radar operator validation

Before evaluating the performance of RadarVar in WRFDA-3D-Var, we first examined the consistency between J08orig and RadarVar. Because RadarVar follows the J08 operator, the operator implemented in ARPS EnKF (Jung et al., 2008a) was used to serve as the J08orig operator. To be comparable to J08orig operator, the coefficients of graupel listed in Tables 1 and 2, as well as the intercept parameter and hydrometeor density, are replaced by those listed in J08, which were designed for hail. These hail-associated coefficients are used only in Sect. 4 for comparison. The noDA 4 h forecast initialized at 00:00 Z was used as the input of the radar operators. The results, which are shown in Fig. 4, indicate that the difference between J08orig and RadarVar is small and acceptable. The small difference is likely caused by the subtle programming differences between the two software packages.

To verify whether RadarVar is sensitive to the unchanged variables during the minimization (see Sect. 2b), four tests were conducted. The same noDA 4 h forecast was used. In these tests, qr, qs, and qg were randomly perturbed with standard deviations proportional to their input values; in other words, a large background mixing ratio caused a large perturbation. The tests called F_unprt, rhom_unprt, and SD_unprt represent that the fraction F, the mixed form density, and the standard deviation of canting angle, respectively, were unperturbed (i.e., fixed input calculated from the forecast input). The test called all_prt denotes that all of the variables in RadarVar were calculated using perturbed mixing ratios. The standard deviation of the perturbation was set to 10 % of the background value; thus, the maximum perturbation could be large, as shown in Fig. 5a. Figure 5 shows that the reflectivities computed in F_unprt, rhom_unprt, and SD_unprt do not differ significantly from those computed in all_prt. This result indicates that keeping these variables unchanged during the minimization is an acceptable approximation. The most noticeable difference occurs in the middle layer, which is marked in red circles that plot off the diagonal line by several dBZ (Fig. 5c), but there are few of these samples.

Figure 8The (a, b, c) 00:00 Z and (d, e, f) 01:00 Z analyses of qr (red), qs (green), and qg (blue) at model levels 5 (left column), 15 (middle column), and 20 (right column) for outer loop 6 in Exp_ref. Model levels 5, 15, and 20 approximately correspond to 0.7, 4.0, and 8.0 km a.g.l., respectively.


The tangent linear operator of RadarVar was verified through a ratio, which is given by

(30) H ( x + ε δ x ) - H ( x ) ε H ( δ x ) = 1 + O ( ε ) ,

where H and H represent the nonlinear operator and the corresponding TL version, respectively, x is the vector of the model state variables (qr, qs, qg), whose perturbation is denoted by δx, and ε is the perturbation magnitude and is greater than zero; δx used the perturbations generated for all_prt. Table 3 shows that the tangent linear operator is sufficiently accurate with a ratio close to 1.0.

A correct adjoint operator must satisfy the relationship that is given by

(31) ( H δ x ) T H δ x = δ x T H T ( H δ x ) .

The meanings of all of the symbols in this equation are consistent with those in Eq. (30), and HT is the adjoint operator. All of the perturbations used for the tangent linear test were adopted in the adjoint test. In the double precision test, there were 14 identical digits on both sides of Eq. (31); in the single precision test, there were seven identical digits.

5 Results

5.1 Analyses in the observation space and the model space

The Exp_ref analyses (six outer loops, 150 inner iterations, and 0.1 weight for the retrieval part) are first examined in terms of the radar reflectivity and the mixing ratios of rain, snow, and graupel. As expected, the analyzed reflectivity agrees much more closely with the observed reflectivity than the background reflectivity for both analysis times (00:00 and 01:00 Z) (Fig. 6). In addition, considering that the 00:00 Z background bias is much larger than that at 01:00 Z, the comparable analysis error for both times indicates that the analysis is generally relatively insensitive to the initial background. However, a small number of points in the analysis have zero reflectivity, while the corresponding observations can be greater than 10 dBZ (Fig. 6b, d). Further examination indicates that these failed points are related to the locations with very small background values of qr, qs, and qg, where the nonlinearity of the radar operator and the deficiency of the TL/AD operator are more pronounced.

Figure 9The cost function and normalized gradient norm as functions of the iteration step during the minimization at (a, c) 00:00 Z and (b, d) 01:00 Z in Exp_ref.


Figure 7 shows the horizontal distributions of the observed and analyzed reflectivities at 2 km, 4 km (melting layer), and 6 km a.g.l. In general, they match well in observed areas with a weaker analysis at 00:00 Z, which is likely caused by the bad background. Spurious echoes appear over unobserved areas in the analyses at both 00:00 and 01:00 Z and most likely resulted from the spatial correlations in the background error covariance; these correlations allow the propagation of information from observed areas to unobserved areas both horizontally and vertically. This result implies that the statistical correlation length scales obtained from a 1-month forecast sample may be too large for this squall line case. The results of tuning the horizontal correlation length scales will be given in Sect. 5.4. Another cause of these spurious echoes is that non-precipitation echoes were not assigned in the observation data in this study such that DA has no impact outside the observed convective area. An approach to suppress the spurious echoes is to determine the non-precipitation points, assign a specific value like 0 dBZ to these points, and assimilate these non-precipitation echoes. The non-precipitation echoes will be considered in the future.

The 00:00 and 01:00 Z analyses of qr, qs, and qg at three model levels are shown in Fig. 8. The direct assimilation of the reflectivity data using RadarVar successfully retrieved the lower level rain, the upper level snow or graupel, and mixed rain/snow and/or rain/graupel in the melting layer (model level 15). Note that the analysis increment can be created only at levels where the SD of the background error is larger than zero (see Fig. 3).

Figure 10The correlation coefficients of qr (red), qs (green), qg (blue), and reflectivity (gray) between the Exp_ref analysis (with six outer loops and 150 inner iterations each loop) and that with different numbers of outer loops and inner iterations at (a) 00:00 Z and (b) 01:00 Z.


5.2 Convergence of the minimization

Figure 9 shows the cost function and the norm of its gradient as a function of the number of inner iterations for the first, third, and sixth outer loops. For the 00:00 Z Exp_ref analysis, both the cost function and the gradient norm decrease rapidly in the first outer loop due to the large adjustment from the bad initial background. Using the improved guess after two outer loops, the third outer loop starts from an ∼85 % smaller cost function, which is then reduced more gradually with increasing iterations. Performing more outer loops does not further reduce the cost function substantially, and outer loop 6 even has a slightly larger cost function than that of outer loop 3. Similar behavior can also be observed in the 01:00 Z analysis, but performing more than three outer loops can further reduce the cost function to a small extent. This indicates that performing six outer loops may not be necessary and that two or three outer loops may be sufficient.

Figure 11Same as Fig. 6 but for the experiments without hydrometeor preprocessing at (a, b) 00:00 and (c, d) 01:00 Z.


To more quantitatively measure the impacts of the number of outer loops and inner iterations, the correlation coefficients between the Exp_ref analysis (i.e., six outer loops and 150 iterations in each loop) and the analyses with one to five outer loops and 50/100/150 inner iterations in each loop are calculated and shown in Fig. 10. In the 00:00 Z analysis, the correlation coefficients of qr and qg are already greater than 0.9 after the first outer loop, and using 50 or 100 inner iterations also results in good agreement with the Exp_ref analysis for qr and qg. However, this is not the case for qs, which is likely an indication of a bad qs background (0.1 of the retrieved values). Given that the correlation coefficients of reflectivity reach 0.9 after the second outer loop, the low correlation coefficients of qs also imply that it is difficult to distinguish the contributions from ice-phase species. Even with a more realistic background at 01:00 Z, the issue of the qs analysis with less iterations and outer loops (Fig. 10b) still exists, but the higher correlation coefficient of the run with 150 iterations after four outer loops, compared to the counterparts in the 00:00 Z analysis, indicates the importance of the background quality. Note that the hydrometeor's phase information could be better determined by assimilating polarization measurements from dual-polarization radar, which will be explored in a future study.

Figure 12(a, c) The analysis reflectivity at 4 km a.g.l. and (b, d) the analyses of qr (red), qs (green), and qg (blue) at model level 15 at 00:00 and 01:00 Z for the experiments without hydrometeor preprocessing.


5.3 Importance of hydrometeor preprocessing

To further investigate the sensitivity of the analysis to the background quality, additional analyses without the hydrometeor preprocessing step were performed at 00:00 and 01:00 Z. At 00:00 Z, the nonzero background was constructed by using a very small weight (7×10-4) for the retrievals. At 01:00 Z, the background is simply taken from the 1 h WRF forecast. Compared to Fig. 6b, the analysis at 00:00 Z (Fig. 11b) has many more failed points (i.e., those with zero analyzed reflectivity) due to the excessively small background hydrometeor values. In contrast, the 01:00 Z analysis (Fig. 11d) without preprocessing has a comparable error bias and SD to those with the preprocessing (Fig. 6d). These results indicate that the preprocessing step is mostly needed to address the bad background.

Figure 13Same as Fig. 7 but for (a) the observations and the analyses from (b) Exp_ref, (c) Exp_ls0.5, and (d) Exp_ls0.125.


Figure 12a shows that the analyzed reflectivity in the melting layer has a reasonable fit to the observations even though it started from very small values of the background reflectivity. The analysis in the model space (Fig. 12b) is also similar to that with the preprocessing (Fig. 8b). However, the reflectivity coverage (>35 dBZ) in Fig. 12a is smaller than those in the observation and Exp_ref, which corresponds to the horizontally distributed points on the bottom of Fig. 11b. In contrast, the 01Z analysis without the preprocessing step (Fig. 12c, d) closely resembles that with the preprocessing (Figs. 7h and 8e) in terms of the radar reflectivity and hydrometeor mixing ratios, especially over the strong convective line (dashed black line in Fig. 12d). However, removing this preprocessing step results in a large spurious echo (marked by “A” in Fig. 12c) that is smaller in Fig. 7h. This implies that the preprocessing step will still be helpful for some “bad” locations (e.g., mismatched convective cells between the model background and observation) for a generally “good” background.

5.4 Impact of the spatial correlation scale

Figure 13 shows the 00:00 Z analysis reflectivity at 4 km a.g.l. with the horizontal correlation length scales of the background errors reduced by factors of 2 and 4. The spurious echoes weaken as the length scale decreases, especially for the spurious echoes marked by “A”. In addition, reducing the length scale improves the intensity analysis at certain locations (e.g., the convective cell marked by “B”; also see Fig. 7e, f). Note that the background error statistics used for these 3D-Var analyses were obtained from forecast samples over a month-long period and are likely not optimal for this particular squall line case. A better solution would be the use of the flow-dependent background error covariance, which will be investigated in the future using more advanced DA methods, such as 4D-Var and hybrid-3D/4DEnVar that are available in WRFDA.

5.5 Impact on the precipitation forecast

The performance of the precipitation forecasts was quantitatively evaluated using the stage IV dataset with the fractions skill score (FSS; Roberts and Lean, 2008). Hourly precipitation forecasts between 02:00 and 05:00 Z were evaluated because assimilating only the reflectivity is expected to have an impact mostly on the short-term forecast. Similar to Schwartz et al. (2014), the aggregated FSS over the period from 02:00 to 05:00 Z is shown in Fig. 14 for the forecasts initialized at 00:00 and 01:00 Z and for different rainfall thresholds. For the forecast from 00:00 Z (Fig. 15a), Exp_ref obtained greater FSSs than those of noDA for thresholds greater than 5 mm h−1 with a radius of influence smaller than 40 km. With a larger radius, this difference remains for the lighter rain (<20 mm h−1) but disappears for the heavy rain (≥20 mm h−1). Similar behavior can be observed for the forecast at 01:00 Z (Fig. 15b) but with a more positive impact from the radar DA for heavier rainfall and for a radius of influence smaller than 50 km. Further examination (not shown) indicated that the improvement in the moderate rain (>5 mm h−1) prediction was associated with the larger snow area in the analysis, so the light rain missed by noDA was better captured in Exp_ref. This examination also showed that the higher FSSs obtained by Exp_ref for the heavier rainfall were mostly associated with the smaller displacement error.

Figure 14The FSSs of the hourly precipitation forecasts aggregated over the period from 02:00 to 05:00 Z as functions of the radius of influence for forecasts launched from the (a) 00:00 and (b) 01:00 Z analyses for Exp_ref (solid lines) and noDA (dashed lines). The hourly precipitation thresholds are denoted by green (1 mm h−1), blue (5 mm h−1), orange (10 mm h−1), and red (20 mm h−1) lines. The straight dashed line represents the skillful FSS at 1 mm h−1.


6 Conclusions

This study developed tangent linear and adjoint operators based on the J08 reflectivity operator and implemented them in WRFDA. This new operator can compute the reflectivity contributed by ice-phase species and is called RadarVar. RadarVar is effective at analyzing rainwater, snow, and graupel by directly assimilating US WSR-88D S-band radar reflectivity with WRFDA's 3D-Var. The analysis accuracy is somewhat sensitive to the numbers of outer loops and inner iterations. The results indicated that two to three outer loops with 50–100 iterations in each loop are needed to obtain a sufficiently accurate analysis. Two problems of RadarVar were found in our test. One issue is the analysis failures at locations with observed radar echoes but with zero or excessively small model background values of the hydrometeors. This issue can be partially resolved using a preprocessing step with radar-retrieved hydrometeors to improve the “bad” background before the analysis. Another issue is the spurious radar echoes (i.e., precipitation) in the analysis caused by the spatial correlations in the background error covariance. These can be reduced by decreasing the correlation length scales. In addition, the short-term (2–5 h) precipitation forecast is improved by the direct reflectivity DA even though the inexpensive univariate 3D-Var technique is used in this first attempt of applying RadarVar.

A more thorough evaluation of reflectivity DA with RadarVar will be examined in a future study using a more advanced hybrid ensemble–variational DA technique, which allows flow-dependent background errors with multivariate correlations and is expected to further reduce the aforementioned deficiencies. Moreover, RadarVar will be extended to include the computation of polarimetric quantities to better determine the phases of the hydrometeors, especially in the melting layers.

Code and data availability

The RadarVar v1.0 operator is integrated into the community WRFDA software and will be publicly available in a future release. The code of RadarVar v1.0 and the scripts for running experiments in this study can be obtained at (Wang, 2019). The GFS data are available at (National centers for environmental information, 2019a) and the radar data can be downloaded at (National centers for environmental information, 2019b).

Appendix A

This section provides details about the reorganization of all of the terms in the brackets in Eq. (11) in terms of the fwx power. For convenience, the expressions in these brackets are represented by G(fwx).

Applying Eq. (16) to the brackets in Eq. (11) results in

(A1) G ( f w x ) = ε x 2 [ A ( k = 0 n P w x a k f w x k ) 2 + B ( k = 0 n P w x b k f w x k ) 2 + 2 C ( k = 0 n P w x a k f w x k ) ( k = 0 n P w x b k f w x k ) ] .

Expanding all of the terms of Eq. (A1) results in

(A2) G ( f w x ) = ε x 2 { A i = 0 n [ P w x a i ( j = 0 n P w x a j f w x i + j ) ] + B i = 0 n [ P w x b i ( j = 0 n P w x b j f w x i + j ) ] + 2 C i = 0 n [ P w x a i ( j = 0 n P w x b j f w x i + j ) ] } .

Reorganizing the third term (omitting 2C) of Eq. (A2) in terms of the fwx power (i+j) results in

(A3) i = 0 n [ P w x a i ( j = 0 n P w x b j f w x i + j ) ] = P wxa 0 P wxb 0 f w x 0 + P wxa 0 P wxb 1 f w x 1 + P wxa 1 P wxb 0 f w x 1 + P wxa 0 P w x b n f w x n + P wxa 1 P w x b ( n - 1 ) f w x n + + P w x a n P wxb 0 f w x n + P wxa 1 P w x b n f w x n + 1 + P wxa 2 P w x b ( n - 1 ) f w x n + 1 + + P w x a n P wxb 1 f w x n + 1 + P w x a ( n - 1 ) P w x b n f w x 2 n - 1 + P w x a n P w x b ( n - 1 ) f w x 2 n - 1 + P w x a n P w x b n f w x 2 n = k = 0 n [ f w x k i = 0 k P w x a i P w x b ( k - i ) ] + k = n + 1 2 n [ f w x k i = k - n n P w x a i P w x b ( k - i ) ] .

The sum functions in the square brackets on the right-hand side of Eq. (A3) correspond to the third expression of Eq. (20). Using the third expression of Eq. (20), the third term of Eq. (A2) is rewritten as follows:

(A4) 2 C i = 0 n [ P w x a i ( j = 0 n P w x b j f w x i + j ) ] = 2 C { k = 0 n [ f w x k i = 0 k P w x a i P w x b ( k - i ) ] + k = n + 1 2 n [ f w x k i = k - n n P w x a i P w x b ( k - i ) ] } = 2 C k = 0 2 n P C x k f w x k ,

where PCxk has the same meaning as in Eq. (20). Similarly, we can rewrite the other two terms in Eq. (A2) as follows:

(A5) A i = 0 n [ P w x a i ( j = 0 n P w x a j f w x i + j ) ] = A { k = 0 n [ f w x k i = 0 k P w x a i P w x a ( k - i ) ] + k = n + 1 2 n [ f w x k i = k - n n P w x a i P w x a ( k - i ) ] } = A k = 0 2 n P A x k f w x k B i = 0 n [ P w x b i ( j = 0 n P w x b j f w x i + j ) ] = B { k = 0 n [ f w x k i = 0 k P w x b i P w x b ( k - i ) ] + k = n + 1 2 n [ f w x k i = k - n n P w x b i P w x b ( k - i ) ] } = B k = 0 2 n P B x k f w x k .

Because these three expressions (Eqs. A4 and A5) contain the same sum function with respect to k from 0 to 2n, Eq. (A2) can be rewritten as follows:

(A6) G ( f w x ) = ε x 2 ( A k = 0 2 n P A x k f w x k + B k = 0 2 n P B x k f w x k + 2 C k = 0 2 n P C x k f w x k ) = ε x 2 k = 0 2 n ( A P A x k f w x k + B P B x k f w x k + 2 C P C x k f w x k ) = ε x 2 k = 0 2 n ( A P A x k + B P B x k + 2 C P C x k ) f w x k = ε x 2 k = 0 2 n P x k f w x k .

Author contributions

SW performed the coding and designed the data assimilation experiments. ZL supervised this study. All authors contributed to the writing of the paper.


This work was jointly sponsored by the National Key Research and Development Program of China (2017YFC1502103), the National Natural Science Foundation of China (41505089, 41875129, 41505090, 41430427, and 41805070), the Startup Foundation for Introducing Talent of NUIST (2014R007), and National Key Research and Development Program of China (2018YFC1506404). NCAR is sponsored by the US National Science Foundation.

Financial support

This work was jointly sponsored by the National Key Research and Development Program of China (grant no. 2017YFC1502103) and the National Natural Science Foundation of China (grant nos. 41875129, 41430427, and 41805070). NCAR is sponsored by the US National Science Foundation.

Review statement

This paper was edited by Richard Neale and reviewed by Juanzhen Sun and one anonymous referee.


Aydin, K. and Seliga, T. A.: Radar Polarimetric Backscattering Properties of Conical Graupel, J. Atmos. Sci., 41, 1887–1892, 1984. 

Ban, J., Liu, Z., Zhang, X., Huang, X.-Y., and Wang, H.: Precipitation data assimilation in WRFDA 4D-Var: implementation and application to convection-permitting forecasts over United States, Tellus A, 69, 1368310,, 2017. 

Barker, D., Huang, X.-Y., Liu, Z., Auligné, T., Zhang, X., Rugg, S., Ajjaji, R., Bourgeois, A., Bray, J., Chen, Y., Demirtas, M., Guo, Y.-R., Henderson, T., Huang, W., Lin, H.-C., Michalakes, J., Rizvi, S., and Zhang, X.: The Weather Research and Forecasting (WRF) Model's Community Variational/Ensemble Data Assimilation System: WRFDA, B. Am. Meteorol. Soc., 93, 831–843, 2012. 

Borderies, M., Caumont, O., Augros, C., Bresson, É., Delanoë, J., Ducrocq, V., Fourrié, N., Bastard, T. L., and Nuret, M.: Simulation of W©-band radar reflectivity for model validation and data assimilation, Q. J. Roy. Meteor. Soc., 144, 391–403, 2018. 

Caumont, O., Ducrocq, V., Delrieu, G., Gosset, M., Pinty, J.-P., Parent du Châtelet, J., Andrieu, H., Lemaître, Y., and Scialom, G.: A radar simulator for high-resolution nonhydrostatic models, J. Atmos. Ocean. Tech., 23, 1049–1067, 2006. 

Caumont, O., Ducrocq, V., Wattrelot, É., Jaubert, G., and Pradier-Vabre, S.: 1D+3DVar assimilation of radar reflectivity data: A proof of concept, Tellus A, 62, 173–187, 2010. 

Chen, F. and Dudhia, J.: Coupling an advanced land surface–hydrology model with the Penn State–NCAR MM5 modeling system. Part I: Model implementation and sensitivity, Mon. Weather Rev., 129, 569–585, 2001. 

Choi, Y., Cha, D. H., and Kim, J.: Tuning of lengthscale and observation error for radar data assimilation using four dimensional variational (4D‐Var) method, Atmos. Sci. Lett., 18, 441–448, 2017. 

Dawson, D. T., Ii, Xue, M., Milbrandt, J. A., and Yau, M. K.: Comparison of Evaporation and Cold Pool Development between Single-Moment and Multimoment Bulk Microphysics Schemes in Idealized Simulations of Tornadic Thunderstorms, Mon. Weather Rev., 138, 1152–1171, 2010. 

Descombes, G., Auligné, T., Vandenberghe, F., Barker, D. M., and Barré, J.: Generalized background error covariance matrix model (GEN_BE v2.0), Geosci. Model Dev., 8, 669–696,, 2015. 

Evensen, G.: The ensemble Kalman filter: Theoretical formulation and practical implementation, Ocean Dynam., 53, 343–367, 2003. 

Gao, J. and Stensrud, D. J.: Assimilation of reflectivity data in a convective-scale, cycled 3DVAR framework with hydrometeor classification, J. Atmos. Sci., 69, 1054–1065, 2012. 

Hu, M., Xue, M., Gao, J., and Brewster, K.: 3DVAR and cloud analysis with WSR-88D level-II data for the prediction of the Fort Worth, Texas, tornadic thunderstorms. Part II: Impact of radial velocity analysis via 3DVAR, Mon. Weather Rev., 134, 699–721, 2006. 

Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by long-lived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res.-Atmos., 113, D13103,, 2008. 

Janisková, M. and Lopez, P.: Linearized physics for data assimilation at ECMWF, in: Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications, Vol. II, Springer, 2013. 

Janisková, M., Thépaut, J.-N., and Geleyn, J.-F.: Simplified and regular physical parameterizations for incremental four-dimensional variational assimilation, Mon. Weather Rev., 127, 26–45, 1999. 

Johnson, M., Jung, Y., Dawson, D. T., and Xue, M.: Comparison of simulated polarimetric signatures in idealized supercell storms using two-moment bulk microphysics schemes in WRF, Mon. Weather Rev., 144, 971–996, 2016. 

Jung, Y., Zhang, G., and Xue, M.: Assimilation of simulated polarimetric radar data for a convective storm using the ensemble Kalman filter. Part I: Observation operators for reflectivity and polarimetric variables, Mon. Weather Rev., 136, 2228–2245, 2008a. 

Jung, Y., Xue, M., Zhang, G., and Straka, J. M.: Assimilation of simulated polarimetric radar data for a convective storm using the ensemble Kalman filter. Part II: Impact of polarimetric data on storm analysis, Mon. Weather Rev., 136, 2246–2260, 2008b. 

Jung, Y., Xue, M., and Zhang, G.: Simulations of polarimetric radar signatures of a supercell storm using a two-moment bulk microphysics scheme, J. Appl. Meteorol. Clim., 49, 146–163, 2010. 

Kawabata, T., Schwitalla, T., Adachi, A., Bauer, H.-S., Wulfmeyer, V., Nagumo, N., and Yamauchi, H.: Observational operators for dual polarimetric radars in variational data assimilation systems (PolRad VAR v1.0), Geosci. Model Dev., 11, 2493–2501,, 2018. 

Lin, Y. and Mitchell, K. E.: 1.2 the NCEP stage II/IV hourly precipitation analyses: Development and applications, 2005. 

Liu, C., Xue, M., and Kong, R.: Direct Assimilation of Radar Reflectivity Data using 3DVAR: Treatment of Hydrometeor Background Errors and OSSE Tests, Mon. Weather Rev., 147, 17–29, 2019. 

Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., and Clough, S. A.: Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated‐k model for the longwave, J. Geophys. Res.-Atmos., 102, 16663–16682, 1997. 

National centers for environmental information: NEXRAD Inventory, available at:, last access: 8 September 2019. 

National centers for environmental information: Global Forecast System (GFS), available at:, last access: 8 September 2019. 

Parrish, D. F. and Derber, J. C.: The National Meteorological Center's spectral statistical-interpolation analysis system, Mon. Weather Rev., 120, 1747–1763, 1992. 

Pfeifer, M., Craig, G., Hagen, M., and Keil, C.: A polarimetric radar forward operator for model evaluation, J. Appl. Meteorol. Clim., 47, 3202–3220, 2008. 

Posselt, D. J., Li, X., Tushaus, S. A., and Mecikalski, J. R.: Assimilation of dual-polarization radar observations in mixed-and ice-phase regions of convective storms: Information content and forward model errors, Mon. Weather Rev., 143, 2611–2636, 2015. 

Putnam, B. J., Xue, M., Jung, Y., Snook, N., and Zhang, G.: The Analysis and Prediction of Microphysical States and Polarimetric Radar Variables in a Mesoscale Convective System Using Double-Moment Microphysics, Multinetwork Radar Data, and the Ensemble Kalman Filter, Mon. Weather Rev., 142, 141–162, 2014. 

Ryzhkov, A., Pinsky, M., Pokrovsky, A., and Khain, A.: Polarimetric radar observation operator for a cloud model with spectral microphysics, J. Appl. Meteorol. Clim., 50, 873–894, 2011. 

Schwitalla, T. and Wulfmeyer, V.: Radar data assimilation experiments using the IPM WRF Rapid Update Cycle, Meteorol. Z., 23, 79–102, 2018. 

Shen, F. and Min, J.: Assimilating AMSU-A radiance data with the WRF Hybrid En3DVAR system for track predictions of Typhoon Megi (2010), Adv. Atmos. Sci., 32, 1231–1243, 2015. 

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Duda, M. G., Huang, X.-Y., Wang, W., and Powers, J. G.: A description of the advanced research WRF version 3, National Center For Atmospheric Research Boulder CoNCAR/TN-475+STR, 91 pp., 2018. 

Snook, N., Xue, M., and Jung, Y.: Analysis of a Tornadic Mesoscale Convective Vortex Based on Ensemble Kalman Filter Assimilation of CASA X-Band and WSR-88D Radar Data, Mon. Weather Rev., 139, 3446–3468, 2011. 

Snook, N., Xue, M., and Jung, Y.: Ensemble Probabilistic Forecasts of a Tornadic Mesoscale Convective System from Ensemble Kalman Filter Analyses Using WSR-88D and CASA Radar Data, Mon. Weather Rev., 140, 2126–2146, 2012.  

Snook, N., Xue, M., and Jung, Y.: Multiscale EnKF Assimilation of Radar and Conventional Observations and Ensemble Forecasting for a Tornadic Mesoscale Convective System, Mon. Weather Rev., 143, 1035–1057, 2015. 

Sun, J. and Crook, N. A.: Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part I: Model development and simulated data experiments, J. Atmos. Sci., 54, 1642–1661, 1997. 

Sun, J. and Wang, H.: Radar data assimilation with WRF 4D-Var. Part II: Comparison with 3D-Var for a squall line over the US Great Plains, Mon. Weather Rev., 141, 2245–2264, 2013. 

Thompson, G., Rasmussen, R. M., and Manning, K.: Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part I: Description and Sensitivity Analysis, Mon. Weather Rev., 132, 519–542, 2004. 

Tong, M. and Xue, M.: Ensemble Kalman Filter Assimilation of Doppler Radar Data with a Compressible Nonhydrostatic Model: OSS Experiments, Mon. Weather Rev., 133, 1789–1807, 2005. 

Wang, H., Sun, J., Fan, S., and Huang, X.-Y.: Indirect assimilation of radar reflectivity with WRF 3D-Var and its impact on prediction of four summertime convective events, J. Appl. Meteorol. Clim., 52, 889–902, 2013a. 

Wang, H., Sun, J., Zhang, X., Huang, X.-Y., and Auligné, T.: Radar data assimilation with WRF 4D-Var. Part I: System development and preliminary testing, Mon. Weather Rev., 141, 2224–2244, 2013b. 

Wang, S: WRFDA_gmdd, available at:, last access: 8 September 2019. 

Wang, Y. and Wang, X.: Direct assimilation of radar reflectivity without tangent linear and adjoint of the nonlinear observation operator in the GSI-based EnVar system: Methodology and experiment with the 8 May 2003 Oklahoma City tornadic supercell, Mon. Weather Rev., 145, 1447–1471, 2017. 

Wattrelot, E., Caumont, O., and Mahfouf, J.-F.: Operational implementation of the 1D+ 3D-Var assimilation method of radar reflectivity data in the AROME model, Mon. Weather Rev., 142, 1852–1873, 2014. 

Xiao, Q., Kuo, Y.-H., Sun, J., Lee, W.-C., Barker, D. M., and Lim, E.: An approach of radar reflectivity data assimilation and its assessment with the inland QPF of Typhoon Rusa (2002) at landfall, J. Appl. Meteorol. Clim., 46, 14–22, 2007. 

Xue, M., Tong, M., and Droegemeier, K. K.: An OSSE framework based on the ensemble square root Kalman filter for evaluating the impact of data from radar networks on thunderstorm analysis and forecasting, J. Atmos. Ocean. Tech., 23, 46–66, 2006. 

Zhang, G., Vivekanandan, J., and Brandes, E.: A method for estimating rain rate and drop size distribution from polarimetric radar measurements, IEEE T. Geosci. Remote S., 39, 830–841, 2001. 

Zhang, J., Langston, C., and Howard, K.: Brightband Identification Based on Vertical Profiles of Reflectivity from the WSR-88D, J. Atmos. Ocean. Tech., 25, 1859–1872, 2008. 

Short summary
A reflectivity operator was developed for directly assimilating radar reflectivity involving contributions from ice species with the variational data assimilation method. Its current version was implemented in WRFDA 3.9.1. This operator allows for not only the dry snow/graupel but also the wet species so that it can effectively obtain the rainwater, snow, and graupel analysis which improved the short-term precipitation forecasts compared to those of the experiment without DA.