Weakly Constrained LETKF for Estimation of Hydrometeor Variables in Convective‐Scale Data Assimilation

The ensemble Kalman filter algorithm can produce negative values for non‐negative variables. To mitigate this sign problem and to simultaneously maintain the mass conservation, a new concept of combining weak constraints on mass conservation and non‐negativity has been introduced in this work, with a focus on hydrometeor variables in convective‐scale data assimilation. We modify the local ensemble transform Kalman filter with weak constraints on mass conservation for each hydrometeor variable and adopt the assimilation of clear‐air reflectivity data as a weak constraint on non‐negativity. We examine the concept by a series of sensitivity experiments using an idealized setup. Results show that both weak constraints successfully improve the mass conservation property in analyses and both reduce the biased increase in integrated mass‐flux divergence and vorticity. Furthermore, the least biased increase is obtained by combining both constraints, and the best forecasts are also achieved by the combination.


Weakly Constrained LETKF for Estimation of Hydrometeor Variables in Convective-Scale Data Assimilation
Tijana Janjić 1 and Yuefei Zeng 1 1 Meteorologisches Institut, Ludwig-Maximilians-Universität (LMU) München, Munich, Germany Key Points: • A weakly constrained LETKF for mass conservation and non-negativity is introduced and examined in convective-scale data assimilation • Combining both constraints results in the least biases in total mass of hydrometeors and in mass-flux divergence and vorticity in analyses • Best forecasts are also achieved by the combination

Supporting Information:
Supporting Information may be found in the online version of this article.
Correspondence to: Y. Zeng  that remain unresolved in case prognostic variables of the microphysics scheme are updated with EnKF . For instance, due to sparse observations, Gaussian assumption, and other approximations in these methods, the analysis ensemble members are prone to be associated with negative hydrometeor variables, which is physically impossible. In current practice, one way to deal with this issue is to execute physical consistency checks after the analysis step, which clips the negative values to zero for each ensemble member (e.g., Gustafsson et al., 2018). However, this concept inevitably causes a biased increase in the total mass of hydrometeor variables, as shown in Janjić et al. (2014) and . An alternative approach is to transform the variables by assuming that they follow pre-specified non-Gaussian distributions (Cohn, 1997;Lauvernet et al., 2009;Simon & Bertino, 2009); for example, Lien et al. (2013) applied Gaussian anamorphosis on the assimilation of precipitation with the local ensemble transform Kalman filter (LETKF, Hunt et al., 2007). These methods may avoid the sign problem but do not account for mass conservation, as shown in Janjić et al. (2014) with a simple model. Yet another alternative approach, proposed by Janjic et al. (2014), incorporates physical constraints in the data assimilation algorithm for convective scale data assimilation. To this end, a Quadratic Programming Ensemble Kalman filter (QPEns, Janjić et al., 2014), can be used since it incorporates mass conservation and non-negativity as strong constraints (i.e., constraints are perfectly satisfied) in the analysis update, so that the analysis ensemble mean and each ensemble member meet all of these requirements. The QPEns was further evaluated by Ruckstuhl and Janjić (2018) for convective-scale data assimilation, using a modified 1D shallow water model (Würsch & Craig, 2014), and it was shown that the mass conservation helps suppress the spurious convection. Furthermore, Zeng et al. (2017) extended the QPEns to include nonlinear constraints, such as conservation of total energy and enstrophy, and examined the new algorithm with a 2D shallow water model. It was concluded that the enstrophy conservation effectively reduces the unphysical energy cascade to smaller scales for global-scale data assimilation. Currently, the QPEns is still computationally expensive since the strongly constrained optimization problem that may require numerous iterations to converge has to be solved for each ensemble member. The high expense prevents its application to high-dimensional NWP models. However, recent developments (e.g., Janjić et al., 2021;Ruckstuhl et al., 2021) show that computational aspects for QPEns might be overcome in the near future.
In practice, weak constraints (i.e., constraints that are approximately satisfied) are much more established for application with high-dimensional systems. For instance, Liu and Xue (2006) imposed a weak non-negative water vapor constraint in the cost function of the three-dimensional variational data assimilation (3DVAR) algorithm for a global model. Gao et al. (1999) incorporated an anelastic mass continuity equation as a weak constraint in the storm-scale 3DVAR system of the Advanced Regional Prediction System (ARPS; Xue et al., 2000) and showed that this constraint prevents severe error accumulation in the vertical velocity. Subsequently, Ge et al. (2012) added a diagnostic pressure equation as a weak constraint to the 3DVAR of the ARPS to achieve a more dynamically consistent analysis. Li et al. (2015) implemented a dynamic weak constraint based on the steady momentum equations in the 3DVAR system of the Weather Research and Forecasting (WRF) model, and the constraint improves the accuracy of analyses in pressure, temperature, and moisture within the typhoon.
Although weak constraints have often been adopted in variational approaches, only rarely constraints have been used in ensemble methods. Yilmaz et al. (2011) implemented a weak constraint in the EnKF to mitigate the water balance residuals in the land data assimilation, and it was demonstrated that the analyses of the constrained EnKF are more accurate than those of the unconstrained filter. In Barth et al. (2016), a local assimilation scheme that satisfies conservation of mass was tested in an idealized ocean model and compared to the traditional covariance localization with an ad-hoc step enforcing conservation, and it was shown that the inclusion of the mass conservation property reduces the analysis error. In this study, we will modify the existing LETKF algorithm, which is implemented in the Kilometer-scale ENsemble Data Assimilation (KENDA, Schraff et al., 2016) system for the model of COnsortium for Small-scale Modeling (COSMO, Baldlauf et al., 2010), by imposing weak constraints on the mass conservation of hydrometeor variables. Additionally, to preserve the non-negativity of hydrometeor variables, both reflectivity and in particular clear-air reflectivity data are assimilated. It is well recognized (e.g., Aksoy et al., 2009) that the assimilation of clear-air reflectivity data can effectively suppresses spurious convection. As experimentally shown in this work, it actually reduces the number of grid points with negative values in the analysis. Janjić et al. (2014) employed a toy model to demonstrate that when the QPEns strongly preserves both mass and non-negativity, it yields better estimates than algorithms that do not preserve or preserve only mass or non-negativity. In this study, we will explore this idea with the modified LETKF, using the idealized setup for radar data assimilation .
The study is organized as follows: Section 2 presents the concept of the weakly constrained LETKF. Section 3 gives a brief introduction to the COSMO model and a detailed description of the twin experiments setup. Section 4 provides the experimental results. Section 5 gives concluding remarks. Janjić et al. (2014) showed that ensemble-type Kalman filters are mass conservative if no localization is applied. Since localization is necessary for NWP to reduce spurious long-range correlation, the LETKF is not mass conservative by its nature (e.g., Barth et al., 2016;Zeng et al., 2017;Zeng & Janjić, 2016). In this study, we modify the LETKF algorithm by adding a term for domain-wide mass conservation to the cost function and accordingly also modify the update equation. The result is an approximately mass-conserving analysis error covariance as well as analysis ensemble mean and members. Mass conservation is imposed for each type of mixing ratios (i.e., cloud water q c , cloud ice q i , rain q r , snow q s , and graupel q g ), respectively. Details of the algorithm are given in Section S1 in Supporting Information S1.

Weak Constraint on Non-Negativity
Since the hydrometeor variables are non-negative in nature, assimilation of observations that also have the non-negative property, such as radar reflectivity and cloud products, can be utilized to weakly constrain the non-negativity in the analysis step. As shown in Cohn (1997), on a simple example, if observation and background are positive at a grid point, analysis at the grid point will be positive; however, negative values can be obtained in grid points that are updated using covariances and where there are no observations.
The availability of radar reflectivity measurements indicates non-negative values of hydrometeor variables; however, due to the logarithmic unit, reflectivity values are negative for very small reflectivities. To avoid spurious convection (Aksoy et al., 2009) and overestimated analysis increments , a non-negative threshold value is set for very small reflectivities. In this work, the threshold value is 5 dBZ, which means all reflectivity values lower than 5 dBZ are set to 5 dBZ. We call those data "clear-air reflectivity data." Radar reflectivity data depend nonlinearly on hydrometeors. Further, reflectivity data (including clear-air reflectivity data) are available at radar observation locations, therefore not in every grid point of the model; thus, the example of Cohn (1997) does not apply directly. However, as we will show later, by assimilating additional clear-air reflectivity data, we are reducing negative values that could be obtained through covariance updates and asking in the approximate weak sense that non-negativity is preserved in the analysis of hydrometeors.
Note that in the LETKF, the analysis ensemble is generated around the mean that could have zero values, which will yield negative values in some ensemble members. In stochastic EnKF formulation, when assimilating clearair reflectivity for each ensemble member, the non-negativity could be approximately satisfied, in a sense discussed above, in each ensemble member. But unless we impose non-negativity strictly as in Janjić et al. (2014), we cannot completely avoid unphysical values in the analysis ensemble either.

Description of Model and Experimental Design
The COSMO model is run with a 2-km horizontal resolution and a grid size of 200 × 200 × 65, coupled with an Efficient Modular VOlume scanning RADar Operator (EMVORADO, Zeng et al., 2014. Due to the high horizontal resolution, the deep convection is explicitly resolved, whereas the shallow convection is parametrized by the Tiedtke scheme (Tiedtke, 1989). The microphysics scheme is a Lin-type one-moment bulk microphysics scheme (Lin et al., 1983;Reinhardt & Seifert, 2006). Turbulence parameterization is performed with the prognostic turbulent kinetic energy equation of Raschendorfer (2001). Periodic boundary conditions are applied. The convection is triggered by a warm bubble approach (Weisman & Klemp, 1982).
The nature run is described as follows: the supercell is triggered by warm bubble perturbations (Zeng et al., 2020) with an analytical profile (Weisman & Klemp, 1982). The nature run starts at 12:00 UTC, warm bubble perturbations are added at 12:30 UTC, and a supercell is triggered. The supercell splits into two cells at 15:00 UTC. The left-moving cell remains stable in size and shape, while the right-moving cell varies with time and further splits. An illustration of the nature run is given in Figure S1 in Supporting Supporting Information S1.
A radar network consisting of seven stations is positioned such that their scanning area covers the propagation path of supercells and the domain of interest. The synthetic radar observations are generated by adding random Gaussian noise with a standard deviation of 5.0 dBZ and 1.0 m/s to the true reflectivity and radial wind data from the nature run, respectively. The radar observations are available in 5-min intervals and are spatially superobbed to 5 km. Both radial wind and reflectivity data are assimilated. If clear-air reflectivity data are assimilated, a threshold value of 5 dBZ is set, that is, all reflectivity values smaller than 5 dBZ are set to 5 dBZ. If clear-air reflectivity data are not assimilated, all reflectivity values smaller than 5 dBZ are set to missing values. The horizontal localization radius is 8 km, and the vertical localization radius changes from 0.0075 to 0.5 (in logarithm of pressure) with the height. The observation error covariance matrix is diagonal.
The ensemble size is 80. The initial ensemble is created by the first perturbing time, location, and intensity of the warm bubble at 12:30 UTC (more details about the perturbations can be found in  and then evolving ensemble in time freely until 15:00 UTC. An illustration of the initial ensemble is given in Figure  S2 in Supporting Supporting Information S1. Data assimilation starts at 15:00 and ends at 17:00 UTC with an assimilation window of 6 min. The multiplicative inflation (Anderson & Anderson, 1999) is employed with a coefficient of 1.1. In all experiments, after the analysis update, a physical consistency check is applied to each ensemble member, which seeks negative values of hydrometeor variables and clips them to zero. We refer to the analysis after the analysis update but before the physical consistency check as "pre-analysis." The ensemble of pre-analysis could be associated with negative values for hydrometeor variables. Note that the model states at the lateral boundaries are not updated. At 16:00 UTC and 17:00 UTC, 2-h ensemble forecasts are run (forecast data are outputted every 15 min).
To explore the performance of weak constraints on mass conservation and non-negativity, four sensitivity experiments are executed: 1. E_Control: Clear-air reflectivity data are not assimilated and the mass conservation constraints are not imposed; 2. E_M ("M" indicates mass conservation constraints): Clear-air reflectivity data are not assimilated and the mass conservation constraints are imposed; 3. E_P ("P" indicates positivity constraint): Clear-air reflectivity data are assimilated and the mass conservation constraints are not imposed; 4. E_MP: Clear-air reflectivity data are assimilated and the mass conservation constraints are imposed. All experiments update model state variables including 3D wind, temperature, pressure, and hydrometeor variables. In addition, in E_M and E_MP, the ensemble means of hydrometeor variables are corrected by Equation S6 in Supporting Supporting Information S1.

Experimental Results
To evaluate the performance during data assimilation, the root mean square errors (RMSEs) of the ensemble means are calculated for the aggregate of hydrometeors (i.e., q c + q i + q r + q s + q g ) and for the meridional wind v. Furthermore, as in many previous studies in idealized setups (e.g., Tong & Xue, 2005), to obtain a more subtle and accurate measure for assessing the storm, the RMSEs are also calculated at points at which reflectivities in the nature run are greater than 5 dBZ. The relative changes (compared to true values in the nature run) in the analysis ensembles for the total mass of each hydrometeor variable and the total amount of divergence and vorticity, as well as the number of model grid points with negative values of hydrometeor variables in pre-analysis ensembles, are provided. Furthermore, the analysis ensembles are compared by grid-point-based ensemble probabilities for the reflectivity composite (i.e., for each grid point calculating the percentage of ensemble members that exceed a given threshold value of reflectivity). To evaluate the forecast skill, the fractions skill score (FSS, Robert & Lean, 2008) for the reflectivity composite is computed. Figure 1 shows the RMSEs of the ensemble means for the aggregate of hydrometeors and for the meridional wind v. From the vertical profiles in Figure 1a, E_P and E_MP produce smaller RMSEs than E_Control and E_M for the aggregate in case of considering the whole domain, however, if considering only the area of storms, the RMSEs of E_M (E_MP) are larger than E_Control (E_P). This is because the mass conservation constraints act globally, not just trying to fit radar data inside storms. If breaking down the RMSEs into cycles as in Figure 1b, no clear differences can be seen except that the RMSEs of E_Control become considerably larger at the end of cycles. For the v, it is noticed in Figure 1c that the RMSEs in E_Control and E_M are smaller than in E_P and E_MP for both the whole domain and the area of storms. This is attributed to less accurate results of E_P and E_MP in the first half of the cycles as shown in Figure 1d. Since radial wind data are directly correlated to v, in case of not assimilating clear-air reflectivity data, radial wind data have more impacts on v and more accurate v can be achieved. Otherwise, some spin-up time is required to obtain an appropriate covariance between clear-air 6 of 11 reflectivity data and v. Overall, imposing mass conservation constraints reduces the RMSEs of hydrometeors in global but not necessarily in area of storms, and in case of assimilating clear-air reflectivity data, a longer spin-up time is required for dynamical components of storms. Figure 2a shows the relative changes in the total mass for each hydrometeor variable and for the aggregate of them in the analysis ensembles of E_Control, E_M, E_P, and E_MP, compared to those of the nature run. Therefore, positive values mean positive biases, and vice versa. It is evident that all experiments produce positive biases for all hydrometeor variables but to different extents. For instance, it can be seen for q r that E_Control yields an increase in ∼115% compared to ∼70% in E_M and ∼20% in E_P. The smallest biased increase (∼15%) is achieved in E_MP. In addition, the standard deviations of the relative increases in the analysis ensembles are greatest in E_Control, followed by E_M, E_P, and E_MP. Similar patterns can be seen for the other hydrometeor variables and the aggregate. Further, relative changes in total mass of each analysis ensemble member of E_Control are more widely distributed than those of E_M, and the averages over the ensembles in E_Control are greater than that in E_M (breakdown for each experiment is shown in Figure S3 in Supporting Supporting Information S1). Comparatively, the analysis ensembles of E_P are much more narrowly distributed, and the averages are much smaller. In E_MP, the distributions are even narrower in the subsequent half of the cycles, and the averages are even smaller. All these clearly show a decrease in uncertainty in mass through time, but it is noted that this decrease does not reduce the analysis ensemble spread of hydrometeor variables (not shown). of total mass for each hydrometeor variable (i.e., q c , q i , q r , q s , and q g ) and the aggregate of them (indicated by "all") in E_Control, E_M, E_P, and E_MP compared to those of the nature run. Bars indicate relative changes averaged over all analysis ensemble members and all cycles, and error bars indicate standard deviations. Positive bar values mean greater than the nature run and vice versa. (b) Numbers of model grid points in analysis with negative q c , q i , q r , q s , and q g before the physical consistency check. Bars indicate numbers of points averaged over all analysis ensemble members and all cycles before the physical consistency check, and error bars indicate standard deviations. (c) Relative changes [%] of total amount for integrated mass-flux divergence (sdiv) and vorticity (svor) compared to those of the nature run. Figure 2b provides the number of the model grid points with negative q c , q i , q r , q s and q g in the pre-analysis ensembles. Taking q r as an example, E_Control produces approximately 90,000 grid points with negative values (∼3.5% of total grid points), and E_M produces slightly fewer grid points. E_P is associated with much fewer negative points (∼50,000) and E_MP further reduces the number of negative points to 35,000. For the other variables, it also holds that E_Control has the largest number of negative points, followed by E_M, E_P, and E_MP. Furthermore, Figure 2c illustrates the relative changes in total amount for integrated mass-flux divergence (denoted by "sdiv") and vorticity (denoted by "svor"). We regard changes in these quantities as an indicator for dynamical imbalance for the entire domain, in particular, since sdiv is proportional to the surface pressure tendency in the case of hydrostatic pressure . Positive biases can be seen for both sdiv and svor in all experiments. For both sdiv and svor, E_Control produces the greatest biased increase by average and E_MP produces the least biased increase by average. For sdiv, E_M produces slightly more than E_P, whereas the opposite results are achieved for svor. Assimilation of clear-air reflectivity as a weak constraint on non-negativity greatly contributes towards reducing the number of grid points with negative values of hydrometeor variables in pre-analyses and towards mitigating the biased increase and uncertainties in the total mass. In this regard, imposing weak constraints on mass conservation proves to be effective (especially in reducing mass bias) but not as effective as assimilating clear-air reflectivity. The best results are obtained by a combination of both constraints. Moreover, imposing these two constraints also mitigates the biased increase in divergence and vorticity, therefore, a slow error growth in forecasts is expected as shown below. Furthermore, Figure 3 compares the analysis ensembles of E_Control, E_M, E_P, and E_MP at 16:00 UTC by the reflectivity composite at elevation 0.5° (in terms of the gridpoint-based ensemble probability, i.e., the number of ensemble members exceeding the threshold value of 25 dBZ divided by the ensemble size). After 10 cycles, E_Control captures the correct cells with very high probabilities; however, it produces spurious convection in a widespread area in the rear, with fairly high probabilities at some locations. E_M also produces spurious convection, but with areal coverage much smaller than that of E_Control. The false signals are essentially removed in E_P and even less present in E_MP. This finding is consistent with Figure 2. Figure 4a provides the FSSs of reflectivity composite at 0.5° for a threshold of 25 dBZ and a box length of 16 km as a function of the forecast lead time. The FSS values of E_P and E_MP are very close at the initial time (both are greater than 0.95), but the skills of E_P decrease faster with the lead time than those of E_MP. E_M is considerably worse than E_P and E_MP, especially at the initial time (with FSS values of ∼0.72). However, E_M is still much better than E_Control throughout the forecast lead time. Figure 4b shows the FSS values of reflectivity composite as a function of the box length (from 4 to 128 km), where the values are averaged over all members and over all forecast lead times. It is clear that the FSS values of E_MP are consistently higher (approximately 0.05) than those of E_P for all box lengths. E_P outperforms E_M by approximately 0.1; the latter outperforms E_Control by approximately 0.1 for a box length of 4 km and the advantage increases slightly with increasing box lengths. The fact that experiments with weak constraints enable improved forecast skills may be, on one hand, due to less biased hydrometeor variables in the analysis, and on the other hand, due to more dynamically balanced analysis state as seen through a smaller increase in divergence and vorticity.

Conclusion and Outlook
Despite the significant benefits of EnKF for convective scale data assimilation, updates of hydrometeors with observations can produce negative values in the analysis ensemble and add mass bias to the model. To mitigate this sign problem and to simultaneously maintain mass conservation, a new concept of combining weak constraints on mass conservation and non-negativity has been introduced in this work with a focus on the hydrometeor variables. We modify the LETKF with weak constraints on the mass conservation of each hydrometeor variable and adopt the assimilation of clear-air reflectivity data as a weak constraint on non-negativity.
We examine the concept by a series of sensitivity experiments, using an idealized setup for radar data assimilation. The results show that both weak constraints may not always reduce the errors of model state estimates; however, they both successfully improve the mass conservation property of hydrometeor variables of analyses. In detail, weak constraints on both mass conservation and non-negativity considerably reduce the number of grid points with negative values of hydrometeor variables during analysis update and mitigate the biased mass increase and uncertainties in total mass, while the latter constraint is more effective. It is seen that the reduction in biased mass increase causes less spurious convection. In addition, both constraints also reduce the biased increase in integrated mass-flux divergence and vorticity. Overall, the least biased increase is obtained by combining both constraints. Weak constraints on mass conservation and non-negativity also enhance the forecast skills, which can be attributed to the less biased hydrometeor variables, integrated mass-flux divergence, and vorticity.
How the impacts of constraints vary for different assimilation windows is worth of further investigation. Based on the results in more idealized setups (Ruckstuhl & Janjić, 2018), there are two effects to consider. On one hand, if there is more time between analysis, errors grow in mass as well. On the other hand, the noise and spurious convection are more extreme for more frequent updates (Lange et al., 2017). In both cases, it seems to be beneficial to include these constraints, but in which case it is more beneficial it is not clear. Further, this method can also be combined with the adjustment filter  which works more on the dynamical balance.
As the proof of concept illustrates, method proposed here is well suited for the assimilation of radar data in convection permitting models. It is appropriate for estimation of high dimensional state and requires only minor changes to the already existing implementation of LETKF. In the next steps, we plan to apply it to the convective-scale data assimilation with real data. The total mass is unknown in real cases; therefore, it should be treated as uncertain and estimated during data assimilation (Janjić et al., 2014). As a choice, the value of total mass of each analysis ensemble member can also be restricted to that of the background ensemble member as in Ruckstuhl and Janjić (2018). Furthermore, this concept should also be of interest for the other research fields as chemistry, ecosystem, and ocean data assimilation (Bertino et al., 2003;Bocquet et al., 2010;Buehner et al., 2013;Simon & Bertino, 2012), in which mass conservation and non-negativity could also be important.