A Radar Reflectivity Operator with Ice-Phase Hydrometeors for Variational Data Assimilation ( RadZIceVar v 1 . 0 ) and Its Evaluation with Real Radar Data

A reflectivity forward operator and its associated tangent linear and adjoint operators (together named RadZIceVar) were developed for variational data assimilation (DA). RadZIceVar can analyze both rainwater and ice-phase 10 species (snow and graupel) by directly assimilating radar reflectivity observations. The results of three-dimensional variational (3DVAR) DA experiments with a 3 km grid mesh setting of the Weather Research and Forecasting (WRF) model showed that RadZIceVar 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 3D analysis of rainwater, snow, and graupel, including the melting layers with mixed-phase hydrometeors. The deficiencies in the analysis 15 using this operator could be caused by the poor quality of the background fields and the use of the static background error covariance, and these issues 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 RadZIceVar also improved the short-term (2 h-5 h) precipitation forecasts compared to those of the experiment without DA.


3
The main purpose of this study is to develop a TL/AD operator based on Jung et al. (2008) with the contributions of icephase precipitation and apply it in a variational DA framework.For convenience, the operator implemented in this study is called RadZIceVar 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 Section 2, the J08 operator is reviewed, and its TL and AD operators are derived.The experimental design is given in Section 3, and the new operators are verified in Section 4. The performance of RadZIceVar is discussed in Section 5, and the conclusions are presented in Section 6.

Review of the J08 operator
The radar-observed reflectivity, Z, is given in logarithmic form as where Z e is the equivalent reflectivity factor, which is the sum of the contributions from pure rainwater (Z r ), dry snow (Z ds ), dry graupel (Z dg ), wet snow (Z ws ), and wet graupel (Z wg ) as follows: e r ds dg ws wg To compute Eq. ( 2), it is necessary to use the mixing ratios of mixed-phase species (wet snow and wet graupel).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 0.3 r x x r max [min( / , / )] F q q q q F  , where F max is the maximum fraction, which is 0.5 (0.3) for rain-graupel (rain-snow) mixtures; q r is the mixing ratio of rainwater; and q x 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 pr ws wg r ds ws s dg wg g ws w s r wg wg g r (1 ) (1 ) (1 ) ( ) ( ) s q F F q q F q q F q q F q q q F q q 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.(1 The subscript "x" in ρ wx , ρ x , and f wx represents either snow (s) or graupel (g), and f wx is called the water fraction.

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, Z r , is written in a simple form (Zhang et al., 2001;Posselt et al., 2015;Kawabata et al., 2018) as follows: where λ is the wavelength of the radar, which is 107 mm for S-band radar, and N 0r is the intercept parameter of rainwater, which is 8×10 6 m -4 in this study.A fixed N 0r value is only available for a single moment microphysics scheme; for a twomoment scheme, this value should be determined using the predicted number concentration.K w 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 where ρ r =1000 kg m -3 is the rain density, q pr 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 q pr as follows: 1.77 r pr r pr where The value of P r is approximately 4.8×10 9 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).

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 where the subscript x represents either snow (s) or graupel (g).The intercept parameters of these species are denoted by N 0x , the values of which are 3×10 6 m -4 and 4×10 5 m -4 for snow and graupel, respectively.The slope parameter in Eq. ( 11) for either dry snow or dry graupel is written as 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 According to J08,  is zero for all hydrometeors, and σ is different for snow (20°) and hail (60°).In Eq. ( 11), these functions 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.
For brevity, Eq. ( 11) is rewritten as a function of q dx for dry snow and dry graupel and is given by

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 substitutes for ρ x .For wet species, σ in A, B, and C is a function of f w and q wg .Additional details are given in Section 3c of J08.
where "x" is "s" (g) for snow (graupel), ε x is 10 -4 (10 -3 ) for snow (graupel), P wxak and P wxbk are precalculated constants for Sband radar, and the superscript k denotes the index of these constants.Based on J08, the value of n is 6.The values of P wxak and P wxbk 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 (q wx ) and the water fraction (f w ) as follows: 2 1.75 2 wx wx wx wx wx x x wx 0 ( , ) where where P wxai and P wxbi are precalculated constants listed in Table 2.The subscript "x" in these constant coefficients represents either snow (s) or graupel (g).The derivations of Eq. ( 17) to Eq. ( 20) are given in the appendix.

Tangent linear operator
Because RadZIceVar is highly nonlinear and complex, performing the derivation for every nonconstant variable in this operator is difficult.In addition, some components of RadZIceVar (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 RadZIceVar, it is beyond the scope of this study.Here, we assumed that five variables in RadZIceVar 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 mixedphase 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 q x is equal to q r .The intercept parameter is constant because RadZIceVar was currently designed for single-moment microphysics schemes.Although multimoment 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 Section 4.

Linearization for rain
Considering the assumptions presented above, P r in Eq. ( 9) becomes a constant in the minimization.Therefore, the linearized form of Eq. ( 9) for q pr is given by r pr pr r pr r pr r 0.77 r pr ws wg r () () =1.77P ( ) ( 1)

Linearization for dry and wet snow/graupel
For dry snow and graupel, the linearized form of Eq. ( 14) is given by The linearization of Eq. ( 17) can be categorized into two parts, which represent the variations in Z x caused by changes in q wx and f wx , respectively.The linear equation of Eq. ( 17) is written as P q P kf q q q q P q F P f 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 e e 1 10 ln10 where δZ e is the sum of δZ r , δZ ds , δZ dg , δZ ws , and δZ wg .

Adjoint operator
The adjoint operator is the transpose of the tangent linear operator.Because the tangent linear operator is applied to the model variables q r , q s , and q g , the adjoint operator is written for these variables.First, the adjoint operator is written for Eq. ( 24).This operator has the following form: A e e 1 10 ln10 where the superscript "A" means adjoint.
For rainwater in Eq. ( 21), the adjoint operator is given by The parameter A r q on the right-hand side of Eq. ( 26) is the accumulated A r q before computing Eq. ( 26).This rule is also valid for q x .
The adjoint operator of Eq. ( 22) is given by Because Eq. ( 23) is the derivation with respect to both q r and q x , 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 For ice species in Eq. ( 23), the adjoint operator is given by

Sensitivity of RadZIceVar to changes in hydrometeors
Since RadZIceVar is highly nonlinear, understanding its response to changes in q r , q s , and q g assists in the analysis of DA results using this operator.The response functions for Eq. ( 9), Eq. ( 14), and Eq. ( 17) are plotted in Fig. 1. Figure 1b is a twodimensional plot because Eq. ( 17) involves two kinds of hydrometeors.Fig. 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.
In addition, the reflectivity contribution from wet snow increases more substantially when q s or q r are in the range of 0 -0.5 g kg -1 .Fig. 1b shows that the reflectivity reaches 35 dBZ when both q s and q r are approximately 0.2 g kg -1 , while this reflectivity value requires q s 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.
3 Data and experimental design

Case review
This study focuses on a precipitation case that occurred in the northern U.S. The precipitation initiated at approximately 21Z on 1 June 2018, when convective cells formed near the border between South Dakota and Nebraska.By 00Z 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 the middle of South Dakota to the middle of Nebraska.The top of the convective system at this  (Thompson et al., 2004), the 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).The initial conditions and the lateral boundary conditions were generated with the Global Forecast System (GFS) data at 00Z on 2 June 2018.

Generation of the background error covariance
RadZIceVar 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 q r , q s , and q g 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 h and 24 h forecasts.A total of 27 days' forecasts from 20 May 2018 to 15 June 2018 were employed to generate the background error covariance.
The background error statistics for q r , q s , and q g 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 that of snow and graupel mainly exists in the upper levels.The graupel may fall from the upper levels into the lower levels, so the graupel error has a broad distribution.The horizontal correlation length scales of the background errors are often less than 4 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.

Observation data and verification data
Radar data at 00Z and 01Z 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 00Z, 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., 2011Snook et al., , 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.Based on other studies (e.g., Umemoto et al., 2009), the observation error of the reflectivity was set to 1 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 dBZ 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.

Experimental design
As the first attempt to implement and apply RadZIceVar in WRFDA, this study focused on the quality of the analysis using the univariate 3DVAR DA method in terms of the root-mean-square error (RMSE) 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-3DVAR, 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 is utilized to consider 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.Two Exp_ref analyses at 00Z and 01Z were performed.
The background for the 00Z analysis was interpolated from the GFS analysis with zeros for the hydrometeor fields, and the background for the 01Z analysis was the 1-h WRF forecast from the 00Z analysis with more realistic hydrometeor fields.
Note that TL/AD of RadZIceVar will not be able to create the change in reflectivity with hydrometeor perturbations with the zero-hydrometeor background (serves as the base state of TL/AD in the first outer loop), which is the case for the 00Z analysis.An approach to address this issue is to reset the zero background values of q r , q s , and q g 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 J08 expected F to 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 radarretrieved 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 nonvanishing background.The hydrometeor preprocessing was performed at both 00Z and 01Z in Exp_ref as a reference.
To observe the analysis performance with a bad background, the Exp_ref analysis at 00Z was also run with a very small retrieval weight of 5×10 -4 .In addition, the impact of the hydrometeor preprocessing on a relatively "good" background (01Z) 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 00Z with 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.

Radar operator validation
Before evaluating the performance of RadZIceVar in WRFDA-3DVAR, we first examined the consistency between J08orig and RadZIceVar.The pyCAPS-PRS v1.1 software (Dawson et al., 2014;Johnson et al., 2016;Jung et al., 2010;Jung et al., 2008) provided by the Center for Analysis and Prediction of Storms (CAPS) was used to serve as the J08orig operator.The noDA 4-hr forecast initialized at 00Z was used as the input of the radar operators.The results, which are shown in Fig. 4, indicate that the difference between J08orig and RadZIceVar is small and acceptable.The small difference is likely caused by the subtle programming differences between the two software packages.
To verify whether RadZIceVar is sensitive to the unchanged variables during the minimization (see Section 2b), four tests were conducted.The same noDA 4-hr forecast was used.In these tests, q r , q s , and q g 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 RadZIceVar 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.Fig. 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.
The tangent linear operator of RadZIceVar was verified through a ratio, which is given by where H and H represent the nonlinear operator and the corresponding TL version, respectively, x is the vector of the model state variables (q r , q s , q g ), 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 The meanings of all of the symbols in this equation are consistent with those in Eq. ( 30), and H T 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 7 identical digits.

Analyses in the observation space and the model space
The Exp_ref analyses (6 outer loops, 150 inner iterations, and 0.1 weight for the retrieval part) are first examined in both the reflectivity space and model space.As expected, the analyzed reflectivity agrees much more closely with the observed reflectivity than the background reflectivity for the both analysis times (00Z and 01Z) (Fig. 6).In addition, considering that the 00Z background bias is much larger than that at 01Z, 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 q r , q s , and q g , where the nonlinearity of the radar operator and the deficiency of the TL/AD operator are more pronounced.
Figure 7 shows the horizontal distributions of the observed and analyzed reflectivities at 2 km, 4 km (melting layer), and 6 km AGL.In general, they match well in observed areas with a weaker analysis at 00Z, which is likely caused by the relatively bad background.Spurious echoes appear over unobserved areas in the analyses at both 00Z and 01Z 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 one-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 Section 5.4.mixture of rainwater and snow, while a three-phase mixture of hydrometeors is present when no preprocessing is applied.
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.
In contrast, the 01Z analysis without the preprocessing step (Fig. 12c, d) closely resembles that with the preprocessing (Fig. 7h and Fig. 8e) in both the observation and model spaces, especially over the strong convective line (black dashed line in Fig. 12d).However, removing this preprocessing step results in a spurious strong echo (marked by "A" in Fig. 12c) that is much weaker 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.

Impact of the spatial correlation scale
Figure 13 shows the 00Z analysis reflectivity at 4 km AGL 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 strong spurious echo 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, 7f).Note that the background error statistics used for these 3DVAR 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 4DVAR and hybrid-3D/4DEnVar that are available in WRFDA.

Impact on the precipitation forecast
The performances of the precipitation forecasts were quantitatively evaluated using the Stage IV dataset with the fractions skill score (FSS;Roberts and Lean, 2008).Hourly precipitation forecasts between 02Z and 05Z 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 02Z to 05Z is shown in Fig. 14 for the forecasts initialized at 00Z and 01Z and for different rainfall thresholds.For the forecast from 00Z (Fig. 15a), Exp_ref obtained greater FSSs than those of noDA for all thresholds with a radius of influence smaller than 50 km.With a larger radius, this difference remains for the light rain (≤ 5 mm) but disappears for the heavier rain.Similar behavior can be observed for the forecast at 01Z (Fig. 15b) but with a more positive impact from the radar DA for heavier rainfall and for a radius of influence greater than 50 km.Further examination (not shown) indicated that the improvement in the light rain 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.

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 RadZIceVar.
RadZIceVar is effective at analyzing rainwater, snow, and graupel by directly assimilating US WSR-88D S-band radar reflectivity with WRFDA's 3DVAR.The analysis accuracy is somewhat sensitive to the numbers of outer loops and inner iterations.The results indicated that 2-3 outer loops with 50-100 iterations in each loop are needed to obtain a sufficiently accurate analysis.Two deficiencies are observed in the 3DVAR analysis.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 shortterm (2 h -5 h) precipitation forecast is improved by the direct reflectivity DA even though the inexpensive univariate 3DVAR technique is used in this first attempt of applying RadZIceVar.
A more thorough evaluation of reflectivity DA with RadZIceVar 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, RadZIceVar will be extended to include the computation of polarimetric quantities to better determine the phases of the hydrometeors, especially in the melting layers.

Appendix
This section provides details about the reorganization of all of the terms in the brackets in Eq. ( 11) in terms of the f wx power.
For convenience, the expressions in these brackets are represented by G(f wx ).
Applying Eq. ( 16) to the brackets in Eq. ( 11 Reorganizing the third term (omitting 2C) of A.
(2) in terms of the f wx power (i+j) results in The sum functions in the square brackets on the right-hand side of A.
Using the third expression of Eq. ( 20), the third term of A. ( 2) is rewritten as follows, where P Cxk has the same meaning as in Eq. ( 20).Similarly, we can rewrite the other two terms in A.
Because these three expressions (A.(4) and A. ( 5)) contain the same sum function with respect to k from 0 to 2n, A.

Figures
Fig. 1 (a) The reflectivity as (a) a function of q r (red), q s (green), and q g (blue) for pure water and dry snow/graupel and (b) a function of q r -q s (colors) and q r -q g (contours) for wet snow/graupel.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2019-67Manuscript under review for journal Geosci.Model Dev. Discussion started: 25 April 2019 c Author(s) 2019.CC BY 4.0 License.time in terms of the observed reflectivity greater than 5 dBZ reached 16 km AGL.This convective line developed further and moved to southeast Nebraska from 00Z to 04Z.During this period, a bow echo could be observed on the radar mosaic provided by the National Centers for Environmental Information (NCEI), as shown in Fig.2.By 08Z, 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 18Z on 2 June 2018.3.2Settings of the forecast modelVersion 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 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2019-67Manuscript under review for journal Geosci.Model Dev. Discussion started: 25 April 2019 c Author(s) 2019.CC BY 4.0 License.

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.
The coefficients that are multiplied by A, B, and C are functions of f w and are written as Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2019-67Manuscript under review for journal Geosci.Model Dev. Discussion started: 25 April 2019 c Author(s) 2019.CC BY 4.0 License.
the coefficients P wx and P xk are given by 5194/gmd-2019-67 Manuscript under review for journal Geosci.Model Dev. Discussion started: 25 April 2019 c Author(s) 2019.CC BY 4.0 License.