Evaluation of two modified Kalman gain algorithms for radar data assimilation in the WRF model

This work attempts to validate two modified Kalman gain algorithms by assimilating a single radar simulation data set into the Weather Research and Forecasting model using an Ensemble Square Root Filter. Emphasis is placed on the comparison of assimilation performance between the two modified algorithms against the classical Kalman gain algorithm when the measurement operator is non-linear. Three ideal storm-scale experiments, which are configured identically except for the different Kalman gain algorithms, are designed in parallel for this purpose. The results show that the first modified algorithm can result in a better simulation of a storm, as measured by the root mean square error (RMSE). The second algorithm can also, to some extent, reduce the RMSE of the simulation of some state vectors, but with little improvement of the estimation of storm intensity. Overall, our preliminary experiments indicate that the two modified Kalman gain algorithms can benefit the assimilation of complex numerical models when the measurement operators are non-linear, confirming the earlier theoretical analysis and the results of simple models. Further work is needed to evaluate the impact of the modified Kalman gain algorithms on the assimilation performance of ensemble-based methods.


Introduction
Mesoscale convective systems, such as thunderstorms and hurricanes, which are often characterised by sudden occurrence and violent destructive force, have a profound impact on our economy and lives. It has been one of the great challenges to simulate and predict the behaviour of mesoscale convective systems by using numerical models due to their strong non-linearity and high complexity, such as their occurrence at small, local areas during very short periods of time. In recent years, with the development of a variety of improved observation techniques, data assimilation has significantly increased our capability to simulate and predict extreme weather events and climate anomalies by using numerical models to describe these mesoscale convective weather events.
Among the different data assimilation methods, the ensemble Kalman Filter (EnKF) is considered to be one of the most efficient assimilation algorithms (Evensen, 1994;Houtekamer and Mitchell, 2001). It was developed from the Extended Kalman Filter (EKF) to avoid the linearisation of non-linear models, and the prediction (background) error covariance is approximated by using an ensemble of model forecasts. The underlying conceptualisation of the EnKF is that if the dynamical model is expressed as a stochastic differential equation, the prediction error statistics, which are described by the FokkerÁPlank equation, can be estimated by using ensemble integrations (Evensen, 1994(Evensen, , 1997Ambadan and Tang, 2009). Therefore, the error covariance matrix can be calculated by integrating the ensemble of model states. The EnKF has attracted broad attention and has been widely used in atmospheric and oceanic data assimilation (Keppenne, 2000;Evensen, 2003;Lorenc, 2003;Miyoshi and Sato, 2007;Zhang and Pu, 2010;Cui et al., 2011;Deng et al., 2012). With great effort and intensive study in the last decade, the EnKF-based data assimilation technique has become mature, leading to the development of various EnKF schemes to address the issues affecting the standard EnKF, including the ensemble *Corresponding author. email: minjz@nuist.edu.cn Tellus A 2015. # 2015 C. Yang et al. This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), allowing third parties to copy and redistribute the material in any medium or format and to remix, transform, and build upon the material for any purpose, even commercially, provided the original work is properly cited and states its license. square root filter (EnSRF; Whitaker and Hamill, 2002;Schwartz and Liu, 2014), ensemble transform Kalman filter (ETKF; Bishop et al., 2001), ensemble adjustment Kalman filter (Anderson, 2001Liu et al., 2007Liu et al., , 2012; Schwartz and Liu, 2012), maximum likelihood ensemble filter (Zupanski, 2005), and localisation and inflation schemes (Ito and Xiong, 2000;Anderson, 2001;Tippett et al., 2003;Hunt et al., 2007).
In spite of these improvements of the EnKF algorithm, some concerns still exist. One is the potential problem of the Kalman gain algorithm related to the non-linear measurement operator. Ambadan and Tang (2009) and Tang and Ambandan (2009) argued that the non-linear measurement treatment in the classical EnKF contains an implicit assumption that the forecast of the measurement function is unbiased or the mean of the forecast equals the forecast of the mean. Recently, Tang et al. (2014) further analysed the current EnKF algorithm in a statistically rigorous sense and developed two modified algorithms of the Kalman gain. They argued that, in theory, the modified algorithms could lead to better results than the current one. Furthermore, they tested these algorithms by using parameter estimation experiments with the simple threecomponent Lorenz model (Lorenz, 1963). However, they did not evaluate the performance of the modified algorithms with complex weather prediction models applied to convective cases.
This work is an extension of Tang et al. (2014) that further tests the modified algorithms by using a community Weather Research and Forecasting (WRF; Skamarock et al., 2008) model and radar observations and compares the results against that of the classical algorithm. To achieve this, four idealised parallel experiments are conducted to assimilate radar observations into the WRF model by using an EnSRF data assimilation system. Section 2 describes the general Kalman (GK) gain formula and the two modified ones. Section 3 introduces the assimilation system and experimental design. Section 4 presents the results and Section 5 concludes this paper.

The general formula of Kalman gain
An ideal case in data assimilation, when the forecast and observation forward operator are both linear and the errors are Gaussian, can be described as: where t is the time index and x t denotes the true state variables at time t. h tÁ1 is the model error, o t is the observation error and y t is the observation. M and H are linear operators of the model and observation, respectively. The solution to this case among sequential methods is provided by the Kalman filter (KF, Kalman, 1960). However, if the model and observation operators are not linear, the KF cannot be directly used. Thus, the EKF is created to linearise the two operators. However, for the EKF, updating the prediction error covariance is challenging because it requires linearisation of non-linear models, which is difficult in some cases. An alternative approximation is to use an ensemble of model forecasts to estimate the prediction error covariance, as in the EnKF.
The traditional EnKF equations can be described in two parts, as follows (Hamill, 2006): Analysis part: Forecast part: x a is the analysis field, x b is the background field, y i is the perturbed observation, P b is the background error covariance, R is the observation error covariance, n is the number of ensemble members, K is the Kalman gain matrix, h is the observation operator, H is the Jacobi matrix of h, a linearised operator, and M is the non-linear operator of the prediction model.
Thus, in the traditional EnKF, the Kalman gain is defined as 2.2. The first modified algorithm As argued in Tang and Ambandan (2009) and Tang et al. (2014), eqs. (4) and (5) contain an implicit assumption that the measurement operator is linear or weakly non-linear, namely, If the measurement operator is highly non-linear, eq. (8) does not hold. Tang et al. (2014) derived the formula without the linear constraint on the observation operator in the following way: If the error is additive, the non-linear state space estimate can be described as (Tang et al., 2014) If the forecast is unbiased and the number of ensemble members is infinite, considering the random nature of x t , we can use the ensemble mean to replace the unknown true state such that we have is the mathematical expectation. g 0 t and e 0 t are both assumed to be white here. Furthermore, P x b y is the cross-covariance between state space and observation space, and for finite ensemble size, it can be written as In the same way, P yy is the sum of the forecast error covariance in the observation space alone and the observation error covariance. It has been proven that a general form of the Kalman yy (Julie et al., 1995;Simon, 2006). Thus, the modified algorithm of the Kalman gain can be written as (without the subscript t) If the observation operator is linear, eq. (14) is equivalent to eq. (7).

The second modified algorithm
In the discussion of eq. (14), we assumed that the model forecast is unbiased and used the ensemble mean (over a large ensemble) to replace the unknown true value [these assumptions are made in eq. (10)]. However, if these assumptions are not valid, they may result in systematic errors in the calculation of the Kalman gain that may be characterised as a prediction error of measurement. This error may seriously bias the assimilation analysis in some cases (e.g., Tang et al., 2014). One solution to remove the impact of the unbiased assumption on measurement prediction is to directly use perturbed actual observations to represent the measurement prediction of eq. (10), namely, Then, the Kalman gain formula is modified as Because of the use of the perturbed observation in eq. (16), the error covariance R, which is used in eq. (14), is no longer present here.
In summary, if the observation operator is non-linear, then there are three Kalman gain formulas (Table 1). The performance of these formulas will be discussed in the next two sections.

Forecast model and data assimilation system
In this work, the WRF model is used to test the Kalman gain algorithms. The WRF model is a limited-area, nonhydrostatic, terrain following, eta-coordinate mesoscale numerical weather prediction system designed to serve both atmospheric research and operational forecasting needs. The data assimilation system is based on an EnSRF method, which has been developed for WRF at the Nanjing University of Information Science and Technology, China, referred to as WRF-EnSRF (e.g., Wang et al., 2009;Min et al., 2011). The important difference between the EnSRF and the traditional EnKF is that the former does not demand the randomly perturbed observations that are required in the latter. Instead, the analysis eq. (2) is modified in the following manner: where y o is the real observation without perturbation and x a 0 i and x b 0 i are then ensemble perturbations of analysis and forecast, respectively. The over-bars represent the ensemble mean of analysis and forecast. To maintain an unchanged analysis error covariance, a parameter a is introduced (Whitaker and Hamill, 2002;Schwartz and Liu, 2014) The WRF-EnSRF assimilation system can assimilate various types of observation data, including radar, satellite and sounding. In this study, we consider the assimilation of radar data to test the impact of using different Kalman gain algorithms.

The radar data
The radar observations applied to data assimilation are generated by using two steps: (1) using the forward observation operators and true field to calculate the 'true' observations and (2) adding random errors to these unperturbed values. The simulated radar is located at x025, y 075 and all observations have 14 elevations; their interval is 5 minutes (Xu et al., 2008).
The simulated radar data include two variables: radial wind and reflectivity. The simulated radial velocity is calculated by using (Tong et al., 2005) V r ¼ u cos a sinb þ v cos a cosb þ w sin a þ random error (20) where a and b are the elevation angle and azimuth angle of the radar beams, respectively. u, v and w are the model simulated zonal, radial and vertical wind components, respectively, interpolated to the observation points.
The simulated logarithmic reflectivity factor is calculated from the following: where Z e 0Z er 'Z es 'Z eh ; the three parts of this equation are contributions from rainwater, snow and graupel   (20) and (21), we can note that the radial velocity function is linear, while the reflectivity is nonlinear. Therefore, we focus on modifying the Kalman gain formula of reflectivity to test the impact of the three schemes on the assimilation. Figure 1a and 1b illustrates the simulation of radial velocity and reflectivity, respectively, at the 90th minute of the assimilation process. It is clear that the simulated data have 14 layers and basically reflect the true characteristics of radar data.

Experiment design
In this work, we test the modified Kalman gain algorithms by using a classic super-cell storm case that occurred in Oklahoma on 20 May 1977 (Ray et al., 1981) and has been previously used as a test bed for the WRF model. In all experiments, the model is configured in a domain of 200)200)40 grids, with a resolution of 2 km in the horizontal and 0.5 km in the vertical. The storm is simulated by the WRF model initialised from modified real sounding data, as used in Xue et al. (2001). The sounding shows a strong unstable atmospheric stratification, which favours the occurrence of convective weather. To trigger the formation of the storm in the WRF model, we superimpose a 4 K ellipsoidal thermal bubble, centred at x030, y 050 and z 01.5 km with radii of 6 km in the x and y directions and 1.5 km in the z direction, onto the background field characterised by the sounding data. This new background field will be the initial conditions as the true states, while all assimilation experiments' initial conditions are perturbed true states. The model is run for a duration of 1.5 hours with a time step of 12 seconds.
For all assimilation experiments, the ensemble size is set to 40. The ensemble is generated by randomly perturbing the initial conditions on the entire model domain.
These perturbations are drawn from a Gaussian distribution with a SD of the wind of 3 m s (1 , potential temperature of 3 K, and water vapour mixing ratio of 0.0005 kg kg (1 . The relaxation inflation (Zhang et al., 2004) and Schur operator localisation (Gaspari and Cohn, 1999) schemes are also used in all assimilation experiments.
This work focuses on the investigation of how the Kalman gain impacts the assimilation when the measurement operator is non-linear by comparing four assimilation experiments with three different Kalman gain algorithms. Because only the measurement operator of the reflectivity is non-linear, the different Kalman gain algorithms are used only for the assimilation of reflectivity in these experiments. For the assimilation of radial wind, which uses a linear measurement operator, we still use the GK algorithm in all experiments. The detailed settings of the four experiments are shown in Table 2, where the duration of the control run is 90 minutes. At the 20th minute, ensemble members are created by adding perturbations. The data assimilation begins at the 40th minute and is performed every 5 minutes.

The analysis of model states
To validate the aforementioned assimilation experiments, the model control run without assimilation was conducted first. The results show that there are two circulation centres in the horizontal wind field and two updraft centres in the vertical velocity field and that both centres are close to each other. At the 60th minute (see Fig. 2, true), the two vertical velocity centres move to the northeast and southeast, respectively. The convective clouds tend to split at the rear of the updraft centre. At the 75th minute, the two updraft centres continue to move, whereas the convective cloud begins to split to form several small centres. The emergence of two small convective centres, as indicated by high values of the graupel mixing ratio at the rear of the updraft, favours the development and maintenance of all convective clouds as a whole. At the end of the experiment, i.e., the 90th minute, the storm has split and the two initial updraft centres move 20 km northeast and southeast, respectively. A new convection centre with a 10 m s (1 vertical velocity emerges at the storm's initial position. Figure 2 shows the horizontal wind, vertical velocity and graupel mixing ratio at z 04 km for the truth, Exp1, Exp2 and Exp3 at the 60th, 75th and 90th minute, respectively. A striking feature in Fig. 2 is that the three experiments' microphysical fields converge to their actual, measured values after several assimilation cycles. In the early stage (not shown), the initial perturbation of the EnKF suppresses the development of convection, resulting in the vertical velocity being smaller than that of the real storm (i.e., the simulated truth). After additional assimilation cycles, at the 60th minute in Exp1 with the GK, the north updraft area becomes irregular and the intensity of the convection cloud beyond the updraft area is stronger than in the real storm. The intensity of the south convection centre is weaker, while the updraft centre is northerly offset. In Exp2 with MK1, both the range and intensity of the north updraft are similar to the true counterparts and more regular than those in Exp1. The intensity of the

Fig. 2.
Horizontal wind (vectors, m s (1 ), vertical velocity (solid contours at intervals of 5 m s (1 , starting from 5 m s (1 ) and graupel mixing ratio (shading at intervals of 0.0005 kg kg (1 ) at z 04 km for the truth, Exp1, Exp2 and Exp3 (from top to bottom) at the 60th (a), 75th (b) and 90th (c) minute (x-axis and y-axis are model grid numbers). convection cloud of Exp2 beyond the north and south updraft and the horizontal flow are also very similar to those of the truth. However, both Exp1 and Exp2 produce a spurious convection updraft that occurs at the south centre. At the 75th minute, Exp2 is better than Exp1 compared with the truth, when the latter has a larger north centre of the graupel mixing ratio and a smaller updraft pattern. At the end of assimilation, the north centre is smaller in all assimilation experiments than in the real storms, but there is a new convection centre with a 5 m s (1 vertical velocity at the initial storm position, as shown in the truth, Exp1 and Exp2.
Exp3 uses the MK2 scheme with perturbed observations instead of the forecasted counterpart in the Kalman gain. This choice seemingly inhibits the development of convection. However, the convection range is affected, though the north updraft is more regular than in Exp1, as shown in the image of the 60th minute. At the final stage of the assimilation experiment, the simulated intensity of the north convection is closer to the truth in Exp3 than in Exp1 and Exp2, suggesting that for the north convection, Exp3 introduces more accurate information from the observation to adjust the background field. However, using perturbed observation values in the Kalman gain introduces sampling errors, resulting in relatively poor overall assimilation performance in Exp3.
To compare the difference between the first three assimilation experiments in the vertical direction, we take a vertical cross section of the microphysical field, including vertical wind, perturbation potential temperature and graupel mixing ratio, at the 65th minute (Fig. 3). For the simulated true measurements, the top height of the non-zero graupel mixing ratio is approximately 14 km (z028), and at the centre of convection, there is a clear upward movement with a faster speed, while there is a downdraft in front of the convection. The updraft and downdraft form a secondary, vertical circulation, which is good for the development and maintenance of convection. The centre of convection is 5 K warmer than the surrounding environment, with a warm area of 5 km (z010) to 11 km (z022). With this configuration of heat and dynamics, the convection can develop further. Overall, all three assimilation experiments can successfully simulate the convective central location, intensity, vertical motion structure and central heating area. The height of the warm centre and shape of convection in Exp2 is the most accurately simulated in all experiments. In Exp3, the simulated height of the warm centre is lower than in the real one, as in Exp1. The structure of the graupel mixing ratio, especially in front of the storm, is smaller in Exp3 than in Exp1 and 2, further suggesting that the introduction of perturbed observation values in the Kalman gain suppresses the development of the whole convection.

Root mean square error analysis
The root mean square error (RMSE) of the ensemble mean against the truth is used to quantify the quality of the assimilation analysis, where the RMSE is defined as where X t j is the simulated truth of the jth observation point, X i j is the value of the jth observation point of the ith ensemble member, k is the number of ensemble members and n is the number of observation points.
As shown in Fig. 4, there are some significant differences between the RMSE for different variables among the three experiments. In general, the RMSE of all variables in these three experiments decreases with the assimilation time, and the RMSE of Exp2 is smaller than in Exp1 and Exp3 for most variables. This contrast is , perturbation geopotential height (c), perturbation potential temperature (d), ice mixing ratio (e) and snow mixing ratio (f), averaged over points at which the reflectivity is greater than 10 dBZ, for Exp1 (dotted), Exp2 (solid) and Exp3 (dashed) (x-axis is assimilation step).
especially true for the water variables (qi, qs), which is related to the reflectivity operator; the results of the modified methods (Exp2 and Exp3) are all better than the results of Exp1. Figure 5 shows the averaged RMSE of these variables as a function of height over the area where reflectivity is greater than 10 dBZ at the 65th min. Generally, the RMSE of all analysis variables in the vertical is smaller in Exp2 than in Exp1 and Exp3, although, in a few cases, Exp3 is the best for some variables (e.g., qi) at some heights.
In summary, in these assimilation experiments, Exp2 has the best analysis and generates the smallest RMSE for a majority of the analysis variables at most vertical levels. For a few analysis variables, the RMSE of Exp3 is also better than that for Exp1. These results indicate the improvement of the assimilation analysis due to the modified Kalman gain schemes, in particular, the MK1 scheme, when the measurement function is non-linear.
Another useful metric to quantify the assimilation performance is the ratio of the RMSE to the ensemble spread (SPD), here, we define the ratio as RR: To some extent, the impact of noise on data assimilation cannot be displayed by the RR. It has been argued that in a good assimilation system, ideally, the RR should be equal to 1. Figure 6 shows the RMSE and RR for two chosen variables: radar velocity (RV) and radar reflectivity (RF). As shown in this figure, from the 4th assimilation step, the RR of RV is closer to 1 in Exp2 than in Exp1, and the RR of Exp3 is also better than Exp1 from the fourth step to the ninth step. The better value of RR in Exp2 than in Exp1 is, in particular, obvious for RF, which almost approaches 1, while the RR of Exp3 is worse than that of Exp1, indicating the superiority of MK1 to the original GK scheme.

The difference between h x b ð Þ and h x b ð Þ
To explore the difference between h x b ð Þ and h x b ð Þ, which characterises the difference of the Kalman gain algorithm between MK1 and GK, we calculate h x b ð Þ and h x b ð Þ from the background field of Exp2 synchronously at every assimilation time. Figure 7 shows the horizontal component of h x b ð Þ, h x b ð Þ and their difference for reflectivity at the 90th minute. As observed, at the centre of the convection, the difference between h x b ð Þ and h x b ð Þ is small,  but the latter is larger than the former at the storm fringe, suggesting that the non-linearity of the measurement function increases in this region of the storm. Figure 8 shows the averaged difference between h x b ð Þ and h x b ð Þ over the entire assimilation domain as a function of assimilation steps. As shown in this figure, h x b ð Þ is always smaller than h x b ð Þ, indicating that the unmodified algorithm overestimates the linearisation of the non-linear measurement operator. This effect is particularly obvious at the beginning of the assimilation, as shown by a rapid growth of the difference at the initial steps. After six assimilation cycles, the difference stabilises and remains between (1.5 and (1 dBZ. This result suggests that the unmodified algorithm may lead to larger errors of linearisation not only for the simple toy model, as discussed in Tang et al. (2014), but also for complex models, such as the WRF used here.

The sensitivity test for MK2
For testing the performance of the perturbed observed values used in MK2, we designed a new experiment (Exp4) with double the perturbation variance of Exp3. Figure 9 shows the RMSE difference and RR difference for RF between Exp4 and Exp3. It can be observed that after increasing the perturbation variance, the RMSE of RF becomes larger than that of Exp3 over all times. For most model valuables, the RMSE becomes larger too (Fig. 10). This means that the performance of MK2 is sensitive to the magnitude of the observation perturbations. A large perturbation can lead to poor results.

Summary and conclusions
Several recent studies argued that when the measurement function is non-linear, the current EnKF algorithm con-tains an implicit assumption, namely, that the mean of the forecast h x b ð Þ equals the forecast of the mean h x b ð Þ. Based on this argument, two modified Kalman gain algorithms were proposed, tested and validated by using the simple Lorenz model previously developed by Tang et al. (2014). In this study, we further discussed this issue. Emphasis is placed on a comprehensive comparison of the modified Kalman gain algorithms (MK1 and MK2) against the original one (GK) by comparing four parallel assimilation experiments, using a realistic weather forecast model assimilation system (WRF-EnSRF) and radar data. This study is the first to apply the modified algorithms to a high dimensional model. The result shows that using different quantities in the Kalman gain can lead to differential impacts in both horizontal distribution and in the tendency over time. During the assimilation period, all experiments analyse the evolution of the storm reasonably well. Compared to the simulated true values, MK1 generates a better simulation than the current algorithm (GK) in terms of storm intensity and evolution tendency. In general, MK1 has the smallest RMSE and RR values. MK2 also results in some improvements for some analysis variables, but these improvements are not as significant as for MK1. The reason for this result is that different perturbed observations may bring different sampling errors into the assimilation system. In addition, how to carry out the observation perturbation, including the choice of the amplitude of the perturbation, is an open question that may impact the assimilation analysis. These aspects are being studied and will be discussed in a subsequent paper.
It should be noted that this work is an ideal experiment based on the assumption that the forecast model is perfect. A further test and validation is required using a real assimilation system with observed data. Nevertheless, the experimental results presented in this paper are the first from testing the modified Kalman gain algorithms in a  The RMSE difference (solid line) and RR difference (dash line) for RF between Exp4 and Exp3 (x-axis is assimilation step). complex model, potentially providing better assimilation schemes for realistic assimilation problems when the measurement functions are non-linear.