# Impact of hydrometeor control variables in a convective-scale 3DEnVar data assimilation scheme

Mayeul Destouches, Thibaut Montmerle, Yann Michel, Jean-françois Caron

► To cite this version:

Mayeul Destouches, Thibaut Montmerle, Yann Michel, Jean-françois Caron. Impact of hydrometeor control variables in a convective-scale 3DEnVar data assimilation scheme. Quarterly Journal of the Royal Meteorological Society, 2023, 149 (752), pp.757-780. 10.1002/qj.4426. hal-04695661

HAL Id: hal-04695661

https://hal.science/hal-04695661v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Impact of hydrometeor control variables in a convective-scale 3DEnVar data assimilation scheme

Mayeul Destouches<fnref n="1" /> |

Thibaut Montmerle<fnref n="1" /> |

Yann Michel<fnref n="1" /> |

Jean-François Caron<fnref n="2" />

Correspondence

M. Destouches, CNRM, Météo-France,
Université de Toulouse, Toulouse, France.
Email: mayeul destouches@umr-cnrm.fr

Funding information

Institut National des Sciences de
l'Univers, Centre National de la
Recherche Scientifique, Grant/Award
Number: LEFE-HYDRE; Ministère de
l'Enseignement supérieur, de la
Recherche et de l'Innovation,
Grant/Award Number: AMX

Abstract

Initialization of condensate water variables for convective-scale weather forecasting has been a topic of active research for the last few years. In this article, we present the addition of hydrometeor fields as control vari- ables of a three-dimensional ensemble variational scheme (3DEnVar) for the cloud-resolving model AROME-France. Even without any direct assimilation of hydrometeor observations, analysis increments of hydrometeors can be produced via covariances with observed variables in the ensemble-derived background-error covariance matrix. Cycled forecast–analysis experiments in near-operational conditions have been performed over a three-month summer period. Three configurations are compared: (i) a control experiment without hydrometeor control variables, (ii) a test experiment with hydrometeor control variables but without cycling the resulting forecasts, and (iii) an experiment with hydrometeor variables and cycled forecasts. Compared with the control experiment, both hydrometeor experiments show a positive impact up to 9 hr of forecasts in terms of cloud cover, and in the first hour for precipitation. The spin-up period, as evidenced by a reduction of precipitation forecast skill at the beginning of the forecasts, is reduced in length and intensity in the hydrometeor experiments.

KEYWORDS

convective scale, ensemble data assimilation, hydrometeor control variables

1 | INTRODUCTION

Condensate forms of atmospheric water are usually
referred to as hydrometeors. These include liquid hydrom-
eteors (rain, drizzle, cloud droplets) and solid ones (various
types of ice crystals in clouds, snow aggregates, graupel,
hail, etc.). The prediction of hydrometeor fields is a mat-
ter of importance for a variety of applications, from

the prediction of flash flood events to forecasts of solar
energy production. Accurate hydrometeor forecasts are
also needed to describe the impact of hydrometeors on
temperature and moisture properly through a variety
of physical processes (evaporation, latent heat release,
radiation, etc.).

For computational reasons, the representation of
hydrometeors in numerical weather prediction (NWP)

<footnote n="1">¹CNRM, Météo-France, Université de Toulouse, Toulouse, France</footnote>

<footnote n="2">²Meteorological Research Division, Environment and Climate Change Canada, Dorval, Quebec, Canada</footnote>
models is simplified, ranging from very crude representation in global models to more detailed representation at higher resolutions. In convective-scale applications, they are usually separated into classes according to their physical properties, for instance cloud ice, cloud water, rain, snow, and graupel (see the review by Khain et al., 2015). In each class, the information on the distribution of sizes of the hydrometeor particles is condensed into one or two statistical moments of this distribution. Hereafter, the focus is on convective-scale NWP and on single-moment bulk microphysics schemes, which describe hydrometeor variables through their mass mixing ratio.

As prognostic variables of NWP models, hydrometeors require an appropriate initialization. In practice, research on hydrometeor initialization is also especially motivated by the potential gains of assimilating observations affected by clouds and precipitation, such as observations from weather radar or space-based microwave imagers, directly. Various initialization techniques have been proposed, either before, after, or during the assimilation step. For instance, a possible way of assimilating hydrometeor-sensitive observations relies on the two-step procedure known as “1D-Bayesian + 3D/4DVar”. In this technique, hydrometeor observations are converted before the assimilation step into profiles of relative humidity through a Bayesian procedure, and the pseudo-observations obtained are then assimilated with standard 3DVar or 4DVar DA schemes (see for instance Caumont et al., 2010, for radar reflectivities and Duruisseau et al., 2019, for all-sky assimilation of microwave observations). Though these techniques are convenient, the optimal way to exploit all of the available information is theoretically provided by direct assimilation of hydrometeor observations. Having hydrometeors in the control variables of DA schemes is an obligatory milestone on the road towards this direct assimilation (Bannister et al., 2019).

One of the main difficulties on this path is the estimation and modeling of forecast-error covariances for hydrometeors (sometimes called background-error covariances in a DA context in the rest of this article). These covariances are essential ingredients of most DA schemes, and should be described carefully. When hydrometeors are accounted for, all associated equilibria should be represented in the forecast-error covariances. The high nonlinearity of the physics involving hydrometeors yields complex equilibria, which are not easily derived analytically. Besides this, hydrometeor fields are rapidly evolving in space and time, at least compared with conventional variables. The associated forecast-error covariances vary accordingly, which limits the relevance of climatological representations.

The first study including hydrometeors in the control variables of the assimilation (i.e., with hydrometeor control variables) was an attempt to assimilate Doppler radar observations in a 4DVar DA scheme (Sun and Crook, 1997; Sun and Crook, 1998). In this context, hydrometeor forecast-error covariances between timeslots, variables, and grid points were described implicitly through applications of the model tangent linear and adjoint. Nonlinearities were discarded from these by using only a moist microphysics scheme, thus excluding nonlinear processes related to solid water. The authors retrieve the structure and intensity of a convective storm successfully, but fail to estimate the hydrometeor contents properly. In particular, they mention the need for more realistic background-error covariances.

The advent of the ensemble Kalman filter (EnKF; Evensen, 1994; Burgers et al., 1998; Evensen, 2003) provided an efficient way to deal with complex forecast-error covariances. In these schemes, an ensemble of forecasts is evolved through forecast and assimilation cycles. Sample forecast-error covariances can be estimated from this ensemble, thus providing a flow-dependent description of forecast errors for all model variables, possibly including hydrometeors. Due to the limited size of the ensemble, these sample covariances are affected by sampling noise, which can be partially removed through filtering techniques (known as localization). The description of equilibria is no longer an issue with these ensemble approaches, since sampled forecast-error covariances include all the equilibria that are described implicitly by the ensemble. Snyder and Zhang (2003) assimilated synthetic radar observations in an EnKF, using liquid hydrometeor control variables only. Real-case studies followed (Dowell et al., 2004). Data assimilation with five classes of hydrometeors (including ice crystals, snow, and graupel) was then tried by Tong and Xue (2005) in twin experiments with synthetic observations. The first real-case experiments with five hydrometeor classes were performed by Aksoy et al. (2009, 2010). Promising results were obtained in these EnKF studies, thanks to the ensemble description of forecast-error covariances. Similarly to most of the studies that first updated hydrometeor variables in the analysis, they provide little to no discussion on hydrometeor initialization. Hydrometeor control variables are instead used there as a means to assimilate radar data, which indirectly convey information about the thermodynamical structure of convective systems.

Zhang et al. (2013) added five hydrometeor variables in a maximum-likelihood ensemble filter, a specific flavour of ensemble DA where the analysis update is performed through a minimization problem in ensemble space (Zupanski, 2005; Zupanski et al., 2008). In this context, they could assimilate precipitation-affected satellite
observations in the microwave domain, which resulted
in positive impacts on the prediction of a tropical storm
and a heavy rain event. The authors mention the issue of
dealing with displacement errors in this framework, and
raise the question of variable- or observation-dependent
localization when considering hydrometeors in the
analysis.

Ensemble variational methods (EnVar) combine the best of EnKF and variational methods (Hamill and Snyder, 2000; Lorenc, 2003; Buehner, 2005). Being ensemble methods, they still require some localization, which consists of filtering the sampling noise by gradually forcing sample covariances to zero as the separation distance increases. Studies including hydrometeor control variables in EnVar systems at convective scales are relatively recent. Gao et al. (2013) show the potential of 3D ensemble variational (3DEnVar) schemes for radar DA at the convective scale, using simulated data in a model-error-free context. Though their experiments make use of hydrometeor control variables, they do not discuss it further. Wang and Wang (2017) used hydrometeor variables (rain, cloud water, and graupel) in an EnVar scheme where background-error covariances were derived fully from ensembles to assimilate radar reflectivities into the Advanced Research Weather Research and Forecasting Model (WRF-ARW). The introduction of hydrometeor control variables in the EnVar scheme with hybrid background-error covariances of the Advanced Regional Prediction System has been attempted by Kong et al. (2018) with synthetic radar data, and then by Kong et al. (2021) with real observations. The 3DEnVar WRF DA system with hybrid covariances has been used further by Meng et al. (2019) to assimilate cloud-water path products, with positive impacts on the first 12 hr of forecasts of accumulated precipitation. There again, no conclusion is drawn on the importance of hydrometeor initialization independently from hydrometeor-related observations.

Here, we report on the inclusion of hydrometeor
control variables in the prototype of 3DEnVar devel-
oped recently for the convective-scale model Applica-
tion de la Recherche à l'Opérationnel à Méso-Echelle
(AROME: Seity et al., 2011; Brousseau et al., 2016; Mont-
merle et al., 2018). This is the first time hydrometeor
variables have been explicitly added in the control vari-
able of an assimilation scheme for the AROME model.
The sensitive issue of localization for the hydrometeor
ensemble covariances is discussed following Destouches
et al. (2020). In this article, a common localization is
employed for all variables. The impact of hydrome-
teor initialization is evaluated here without any direct
assimilation of hydrometeor-sensitive observations. The
observation operators used in the analysis process have

no sensitivity to hydrometeor variables. In this context,
hydrometeor increments are only produced through
forecast-error covariances between hydrometeors and
observed variables, including pseudoprofiles of relative
humidity derived from radar reflectivities. The other prog-
nostic variables are then impacted by updated hydrom-
eteor fields as the NWP model is integrated forward.
Numerical experiments over a three-month summer
period are performed to assess the impact of updating
hydrometeor fields in the analysis. Cycled and noncycled
DA frameworks are compared, to quantify the possible
differences resulting from cycling biases or imbalances
through successive DA cycles. The forecast skill is evalu-
ated for conventional variables, cloud cover, and precipita-
tion forecasts.

In Section 2, we describe the AROME 3DEnVar DA scheme used in the study. The details of hydrometeor inclusions are reported in Section 3, including considerations on covariance localization. The experimental setup for cycled forecast-assimilation experiments is introduced in Section 4, before presenting the results in a common localization setting. Section 5 concludes the article.

2 | ENSEMBLE AND
DETERMINISTIC ASSIMILATION
SCHEMES FOR AROME-FRANCE

2.1 | AROME-France

AROME-France is the nonhydrostatic limited-area model
used by Météo-France since the end of 2008 to pro-
vide operational weather forecasts at kilometer scale over
France (Seity et al., 2011). The model was originally devel-
oped to improve the prediction skill for some typical heavy
precipitation events in southeast France, resulting from
the vicinity of the Mediterranean Sea and the steep orog-
raphy of the region (Ducrocq et al., 2014). In 2015, the
horizontal resolution was refined from 2.5 km to 1.3 km,
with 90 vertical levels, leading to improvements in rain-
fall forecasts and a better representation of convective cells
(Brousseau et al., 2016).

Being a limited-area model, AROME-France requires
lateral boundary conditions, which are provided by the
global model Action de Recherche Petite Echelle Grande
Echelle (Action de Recherche Petite Echelle Grande Echelle,
Courtier et al., 1991). Surface boundary conditions are
given by the surface model SURFEX (Masson et al., 2013),
which provides vertical fluxes of momentum, heat, and
water.

The dynamical core of the AROME model is based on
a semi-implicit semi-Lagrangian spectral scheme Bénard
et al. (2010). Physical parameterizations for unresolved
processes are unidimensional (per vertical column), and include the microphysics scheme ICE3 (Lascaux et al., 2006; Lac et al., 2018). The name of this scheme originates from its three ice-phase hydrometeor variables: ice crystals $q_i$, snow $q_s$, and graupel $q_g$. The scheme also includes two liquid hydrometeors, cloud water $q_l$ and rain $q_r$. All five hydrometeor variables interact primarily with specific humidity $q_v$. They can impact temperature as well, for instance through latent heat release or through the radiation scheme. The hydrometeor variables are all advected by the wind.

The prognostic variables, that is, the variables required to describe the state of the system fully, are temperature $T$, horizontal winds $u$ and $v$ (zonal and longitudinal components), surface pressure $P_s$, specific humidity $q_v$, the five hydrometeor variables $q_i$, $q_s$, $q_g$, $q_l$ and $q_r$, the turbulent kinetic energy, the nonhydrostatic component of pressure, and the vertical divergence of the vertical wind (see Bénard et al., 2010 on the choice of the nonhydrostatic variables).

In the current operational version of AROME-France, hourly analyses are performed by the data assimilation scheme in a Rapid Update Cycling mode. The last 1-hr forecast is used as a background state for the next analysis. The updated variables (hereafter control variables) are temperature, horizontal winds, surface pressure, and specific humidity. The resulting analysis is used as an initial state for the next 1-hr forecast. Every 3 hr, the forecasts are extended up to 48 hr ahead.

Currently, the operational DA scheme for AROME-France is a 3DVar scheme, assimilating numerous observations of various types (Brousseau et al., 2016). We do not detail its specificities, since the focus of this article is on the more recent 3DEnVar scheme for AROME, presented in Section 2.3.

## 2.2 | Ensemble generation and update: AROME EDA

The ensemble forecasts used by the 3DEnVar scheme are provided by the AROME Ensemble of Data Assimilations (EDA), which we describe in this section. The goal of an EDA is to sample the probability distribution of the errors of a given assimilation and forecasting system, as it evolves through successive forecast-analysis cycles. The purpose is twofold: (i) before the assimilation step, to provide an ensemble to describe background errors to be used in deterministic ensemble-based DA schemes (Section 2.3), and (ii) after the assimilation step, to provide initial conditions to be used in an ensemble prediction system (AROME-EPS in our case: Raynaud and Bouttier, 2015). We are only interested in point (i) hereafter. Note that this

kind of scheme can also be used offline: for example, to calibrate components of climatological background-error covariances (as done by Brousseau et al., 2011 for the operational AROME 3DVar for instance).

EDAs offer some similarities with EnKF systems. Both provide an ensemble of analyses from an ensemble of backgrounds and from observations. In an EDA, in contrast to the EnKF, an analysis is performed independently for each background member. As a consequence, each background is associated with a single analysis, on a one-to-one basis. The consistency of ensemble members in an EDA and forecast and analysis errors of the associated deterministic system was first derived by Žagar et al. (2005, their appendix), in the case of weakly nonlinear systems. It is worth noting, however, that the idea is older, and dates back at least to Houtekamer et al. (1996). More details about EDAs in general can be found in Isaksen et al. (2010).

The version of AROME EDA used for this study is described in Destouches et al. (2020). It evolves 50 members through 3-hr DA cycles. The numerical cost is reduced by using a cheaper version of AROME-France: 3.2-km spatial resolution instead of 1.3 km, 100-s time steps instead of 50 s, hydrostatic version of the dynamical core. The data assimilation scheme used for the assimilation steps is AROME 3DVar.

The sources of variability that ensure the dispersion of the 50 members are the following.

* Variability in the background members, provided by the previous ensemble.

* Independent perturbations of the assimilated observations, consistent with their observation-error statistics.

* Variability in the lateral boundary conditions, provided by the ARPEGE EDA (Berre et al., 2015).

* Variability in sea-surface temperatures for surface boundary conditions.

* Variability in the time integration of the model. At each time step, the physical tendencies provided by the physical parametrizations are multiplicatively perturbed (Stochastically Perturbed Parametrization Tendencies, SPPT: Palmer et al., 2009; Bouttier et al., 2012). The random perturbations have a prescribed correlation structure in space and time.

Note that some other ensemble generation systems perform a recentering of the ensemble members on the deterministic EnVar analysis (for instance Clayton et al., 2012; Kleist and Ide, 2015). These systems are usually based on EnKFs. Here, there is no such feedback on the ensemble, so that the AROME-EDA system receives no input from the associated 3DEnVar scheme.
<figure><img src="image_3.png" /></figure>

## 2.3 | AROME 3DEnVar

### 2.3.1 | 3DEnVar schemes in general

3DEnVar schemes (Hamill and Snyder, 2000; Lorenc, 2003; Buehner, 2005) are deterministic variational DA schemes. They correct a single background to yield a single analysis. The cost function minimized by the analysis is similar to the 3DVar cost function:

$$ J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}^b)^T \mathbf{B}^{-1} (\mathbf{x} - \mathbf{x}^b) + \frac{1}{2}(\mathbf{y} - H(\mathbf{x}))^T \mathbf{R}^{-1} (\mathbf{y} - H(\mathbf{x})), \quad (1) $$

where each term has its usual signification: $\mathbf{x} \in \mathbf{R}^n$ is the control vector, $\mathbf{x}^b \in \mathbf{R}^n$ is the background state, $\mathbf{B} \in \mathbf{R}^{n \times n}$ is an estimate of the covariance matrix of the background error (or forecast error, since the background is given by the last available forecast), $\mathbf{y}$ denotes the observation vector, $\mathbf{R}$ is the covariance matrix of the observation error, and $H: \mathbf{R}^n \to \mathbf{R}^p$ is the observation operator mapping the state space to the observation space. The only difference with a 3DVar schemes lie in the description of the $\mathbf{B}$ matrix, which is derived from an ensemble valid at the moment of the analysis, typically provided by an EDA or EnKF.

From an $N$-member ensemble of forecasts $(\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(N)})$, the standard unbiased estimator of forecast-error covariances is given by

$$ \tilde{\mathbf{B}} = \frac{1}{N-1} \sum_{k=1}^{N} (\mathbf{x}^{(k)} - \bar{\mathbf{x}}) (\mathbf{x}^{(k)} - \bar{\mathbf{x}})^T, \quad (2) $$

where $\bar{\mathbf{x}}$ denotes the ensemble mean. If the normalized perturbations $\mathbf{e}^{(k)} \equiv (1/\sqrt{N-1})(\mathbf{x}^{(k)} - \bar{\mathbf{x}})$ are aggregated in a perturbation matrix $\mathbf{X} = (\mathbf{e}^{(1)} \cdots \mathbf{e}^{(N)})$, Equation (2) can be rewritten in matrix form:

$$ \tilde{\mathbf{B}} = \mathbf{X}\mathbf{X}^T. \quad (3) $$

This ensemble approach gives a flow-dependent description of forecast errors: the ensemble covariances are anisotropic, heterogeneous in space and in time. Besides this, ensemble covariances directly provide an estimation of linear balances between different variables, possibly more accurate than approximate linearizations of microphysics relationships, for instance. These properties make EnVar schemes very appealing for hydrometeor DA.

In practice, the ensemble size $N$ is small, of order 10–100, so that the estimation suffers from sampling noise, especially for near-zero covariances (in the context of AROME-France, see for instance Ménétrier et al., 2014). Localization mitigates this behavior by forcing

long-distance covariances to zero (Hamill et al., 2001):

$$ \tilde{\mathbf{B}}_e = \tilde{\mathbf{B}} \circ \mathbf{L}, \quad (4) $$

where $\mathbf{L}$ is a localization matrix acting as a distance-dependent tapering function, and $\circ$ denotes the Schur-Hadamard (elementwise) product: $\tilde{\mathbf{B}}_{e,ij} = \tilde{\mathbf{B}}_{ij} L_{ij}$. Further details on localization and how it can be applied are given in Section 2.4.

Some further filtering can be done by combining the localized ensemble matrix $\tilde{\mathbf{B}}_e$ linearly with a climatological covariance matrix $\mathbf{B}_{clim}$ (Hamill and Snyder, 2000). This approach, very common, is that of hybrid background-error covariances, in contrast to schemes that only use an ensemble term in $\mathbf{B}$.

### 2.3.2 | Implementation of AROME 3DEnVar

The prototype of the 3DEnVar scheme for AROME has been implemented by Montmerle et al. (2018, M18 hereafter). The nonlinear cost-function (Equation 1) is minimized by approximating it as a quadratic cost function:

$$ J(\delta\mathbf{x}) = \frac{1}{2}\delta\mathbf{x}^T\mathbf{B}^{-1}\delta\mathbf{x} + \frac{1}{2}(\mathbf{d} - H\delta\mathbf{x})^T\mathbf{R}^{-1}(\mathbf{d} - H\delta\mathbf{x}), \quad (5) $$

where $\delta\mathbf{x} \equiv \mathbf{x} - \mathbf{x}^b$ is the candidate analysis increment, $\mathbf{d} \equiv \mathbf{y} - H(\mathbf{x}^b)$ is the innovation vector, and $\mathbf{H}$ is the tangent linear of the observation operator at the background point $\mathbf{x}^b$. Additional linearizations could theoretically be used, and may be useful for strongly nonlinear observation operators $H$. The cost function (Equation 5) attains its minimum when its gradient is zero:

$$ (\mathbf{B}^{-1} + \mathbf{H}^T\mathbf{R}^{-1}\mathbf{H}) \delta\mathbf{x} = \mathbf{H}^T\mathbf{R}^{-1}\mathbf{d}. \quad (6) $$

Various first-level preconditioners can be used in AROME 3DEnVar. In practice, we used a $\mathbf{B}$ preconditioning. Introducing the auxiliary variable $\delta\bar{\mathbf{x}}$ defined as $\delta\bar{\mathbf{x}} = \mathbf{B}\delta\bar{\mathbf{x}}$, Equation (6) becomes

$$ (\mathbf{I} + \mathbf{H}^T\mathbf{R}^{-1}\mathbf{H}\mathbf{B}) \delta\bar{\mathbf{x}} = \mathbf{H}^T\mathbf{R}^{-1}\mathbf{d}. \quad (7) $$

An approximate solution of this linear system can be found iteratively. For instance, an adapted version of the conjugate gradient algorithm can be used, exploiting the fact that the system matrix here is auto-adjoint with respect to the $\mathbf{B}$-inner product (Derber and Rosati, 1989; Gürol et al., 2013).

In AROME 3DEnVar scheme, the ensemble covariance matrix $\tilde{\mathbf{B}}_e$ is used as a $\mathbf{B}$ matrix when computing matrix-vector products. The ensemble members
are provided by the AROME EDA system described in Section 2.2. Note that these members are generated independently of the 3DEnVar.

In this study, a few differences from the configuration of M18 are worth noting.

1. A 3.2-km horizontal resolution is used for both the EDA and the assimilation scheme (in contrast to 3.8 km in M18).

2. The EDA has 50 members instead of 25 and uses an SPPT scheme for model error instead of multiplicative online inflation as in Raynaud et al. (2011).

3. We rely only on ensemble-derived background-error covariances, while hybrid covariances were considered additionally in M18. Using only ensemble-derived covariance estimates can be considered as a first step in the inclusion of hydrometeor control variables. Adopting a hybrid approach would have required the computation of a climatological **B** matrix that includes hydrometeors. This would have been possible, following for instance Michel et al. (2011) and Descombes et al. (2015), but we have chosen to leave it for further studies.

Note that, for consistency, the AROME EDA should ideally be an ensemble of 3DEnVars, since it is used to sample the forecast errors of the AROME 3DEnVar suite. As a first level of approximation, the errors of the 3DEnVar DA suite are supposed to be similar to those of the associated 3DVar suite.

## 2.4 | Localization scheme

Covariance localization in EnVar schemes is the key component that makes the approach viable despite the small ensemble sizes that are numerically affordable. In practice, one can get rid of the elementwise product in the expression for the localized **B** matrix (Equation 4) by using the identity

$$ \tilde{\mathbf{B}}_e = \mathbf{X}\mathbf{X}^\mathrm{T} \circ \mathbf{L} \qquad (8) $$

$$ = \sum_{k=1}^{N} \mathrm{Diag}(\mathbf{e}^{(k)}) \mathbf{L} \mathrm{Diag}(\mathbf{e}^{(k)}), \qquad (9) $$

where $\mathrm{Diag}(\mathbf{e}^{(k)})$ denotes the diagonal matrix with $\mathbf{e}^{(k)}$ on the diagonal (Buehner, 2005). This gives a general formulation to apply the localized matrix $\tilde{\mathbf{B}}_e$:

$$ \tilde{\mathbf{B}}_e \mathbf{x} = \sum_{k=1}^{N} \mathbf{e}^{(k)} \circ \{ \mathbf{L} (\mathbf{e}^{(k)} \circ \mathbf{x}) \}. \qquad (10) $$

In Equation (10), **L** is used as a matrix-vector operator. In practice, the structure of **L** makes it a convolution operator, the kernel of which defines a localization function. This kernel is generally chosen as Gaussian or quasi-Gaussian. Among other things, this property allows us to model a 3D localization by applying successive 1D convolutions along the three different directions, while conserving the isotropy of the resulting localization kernel.

Common techniques to apply localization include recursive filters (Purser et al., 2003; Michel, 2013) and application in spectral space (see M18 for detail on both techniques). Following the results of M18, spatial localization with recursive filters is used in this study.

Note that the cost of applying the localization matrix to a vector can be reduced if a common localization **L** is chosen for every variable in the state vector. Indeed, the localization matrix can then be decomposed as (where $\mathbf{I}_m$ denotes the identity matrix the size of one 3D variable)

$$ \begin{bmatrix} \mathbf{L} & \cdots & \mathbf{L} \\ \vdots & \ddots & \vdots \\ \mathbf{L} & \cdots & \mathbf{L} \end{bmatrix} = \begin{bmatrix} \mathbf{I}_m \\ \vdots \\ \mathbf{I}_m \end{bmatrix} \mathbf{L} \begin{bmatrix} \mathbf{I}_m & \cdots & \mathbf{I}_m \end{bmatrix}. \qquad (11) $$

This common localization is used throughout this article.

## 3 | ADDING HYDROMETEORS TO THE CONTROL VARIABLE

In this section, we present the modifications brought to AROME 3DEnVar to add hydrometeor control variables. Control variables are defined here as those prognostic variables for which nonzero analysis increments can be obtained after the minimization step of the DA scheme.

A first issue that should be discussed beforehand is the non-Gaussian behavior of hydrometeor variables. The variational cost function (Equation 1) is designed for unbiased and Gaussian errors. However, hydrometeor forecasts exhibit strongly non-Gaussian behaviors, because of the positivity of the variables, displacement errors (Lawson and Hansen, 2005), and nonlinear processes involved (Auligné et al., 2011). About the positivity issue, a possibility would be to use a data assimilation algorithm that includes the positivity constraint explicitly. This approach has been proposed by Janjić et al. (2014) in an EnKF setting. An extension to a variational setting was tested numerically by Janjić et al. (2021) on a toy model including a rain variable, showing accurate results at promising computational cost.

As regards non-Gaussian aspects in a more general sense, various workarounds have already been proposed
(Bocquet et al., 2010, for a review). Among them, a large number are based on control variable transforms: logarithmic transforms for lognormal errors (Cohn, 1997; Fletcher, 2007), nonlinear transforms (see for instance the recent works of Chen et al., 2021 or Sun et al., 2021 with hydrometeor control variables), Gaussian anamorphosis (Amezcua and Leeuwen, 2014), transformations to a total water variable (Xiao et al., 2007), physically based transformations (Hólm and Gong, 2010), and so on. All control variable transforms are based on the variational framework. Here, as a first order of approximation, no transformation is applied to the control variables. For now, the statistics of hydrometeor errors are thus assumed to be handled satisfactorily in this Gaussian framework.

The mandatory stage to obtain hydrometeor increments consists of adding hydrometeors in the **B** matrix, so that the background term in the cost function Equation (1) is defined correctly. Since the analysis increment is restricted in practice to the column space of the **B** matrix ($δx = Bδx̄$), this step is especially important.

As explained in Section 2.3.1, the extension of **B** is straightforward in a 3DEnVar where the **B** matrix is completely derived from an ensemble of forecasts. Since hydrometeors are prognostic variables of AROME, they are present in the ensemble forecasts of the EDA (Section 2.2), which are based on the same model. As a consequence, hydrometeor forecast-error covariances can be estimated directly from there, following Equation (2). In this context, the only remaining question is the choice of localization for hydrometeor variables.

## 3.1 | Localization for hydrometeor variables

A simple choice would consist of keeping the same localization for hydrometeor and conventional variables, using a common localization scheme. However, the spatial scales of hydrometeor variables are smaller than those of conventional variables, and this difference propagates to their forecast-error covariances, as shown for instance by Michel et al. (2011). Besides this, it is known that forecast error correlations and optimal localization are closely linked (at least for Gaussian errors, as demonstrated objectively by Ménétrier et al., 2015). As a consequence, it is likely that optimal localization lengths are smaller for hydrometeor variables than for conventional variables.

This intuitive result has been confirmed by Destouches et al. (2020). Using the statistical method developed by Ménétrier et al. (2015), they diagnose optimal localization lengths on eight case studies with 50-member ensembles from the AROME EDA. They find that the resulting diagnoses are dependent on the variable considered,

weather type, and altitude. However, certain constant behaviors have been noted. For instance, for hydrometeor autocovariances, horizontal localization lengths should be shorter by a factor of 2–6 compared with conventional variables. Vertically, the optimal lengths are found to be similar among variables, except for precipitating hydrometeors, which should be localized less severely than other variables. Indeed, the latter variables are highly correlated vertically due to the sedimentation process.

Various localization schemes are proposed in section 5.3 of Destouches et al. (2020) to comply with these characteristics. The first two proposals are common localization and variable-dependent localization. The authors also mention localization schemes with parameters that would be adjusted before each assimilation cycle or flow-dependent localization scheme. Scale-dependent localization (SDL; e.g. Buehner, 2012; Buehner and Shlyaeva, 2015) is an alternative approach that leads to implicit variable-dependent localization and is believed to be self-adaptable to the different flow regimes.

In the present article, the covariances of interest are not hydrometeor auto-covariances, as in Destouches et al. (2020), but covariances between hydrometeor and conventional variables. As explained below in Section 3.3, this is linked to the fact that, in this article, no direct observation of hydrometeors is considered in the assimilation process. In this context, optimal localization lengths may possibly differ slightly from the results of Destouches et al. (2020). Since common localization is both the simplest and the cheapest localization scheme, it is retained for all the experiments presented here.

The authors have also performed the experiments using the SDL scheme of Buehner (2012) (not shown). The conclusions were unchanged, which shows that the main results presented here are not sensitive to the localization choice.

## 3.2 | Illustrations of the B matrix with hydrometeors

The purpose of this section is to illustrate the ability of AROME-EDA to provide physically consistent ensemble-based covariances that include hydrometeor variables. All illustrations were obtained from 50-member ensembles provided by the AROME EDA described in Section 2.2.

Ensemble standard deviations for ice crystal content at a particular vertical level are shown in Figure 1 for a day during autumn 2018. They represent (a portion of) the diagonal of $\tilde{B}_e$, which is unaffected by standard localization. A significant heterogeneity is displayed in the domain, with a large portion of values being null or close
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 1** Example of diagonal of the ensemble **B** matrix with hydrometeors. Ensemble standard deviations of ice crystal fields at model level 40 (approximately 570 hPa), for May 21, 2018 at 1500 UTC. Values lower than $10^{-12} \text{ kg} \cdot \text{kg}^{-1}$ are shown in white. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

to zero. No hydrometeor increments will be produced in these zones of zero variances. Together with the zero gradient problem (see Errico et al., 2007), this may become a problem if hydrometeor observations are present in these regions. Hybridizing the ensemble covariance matrix with a static covariance matrix including hydrometeors would be a way to alleviate these difficulties. Another possibility would be to force zero variances to small positive values (Bannister et al., 2019).

Note that the variances of hydrometeor background errors are not actually used in the data assimilation step if there is no direct assimilation of hydrometeor observations, as is the case in this article. In this case, only cross-covariances between hydrometeor and conventional variables are used.

Cross-covariances between hydrometeor and conventional variables are shown in Figure 2, for vertical separations only. In order to remove the sampling noise from the visualization without having to choose specific localization parameters, we show here the average value of the covariances in regions with nonzero hydrometeor content. For a given hydrometeor variable, and following the Destouches et al. (2020) approach, these regions are defined as the set of model columns for which at least one point has a nonzero hydrometeor variance in the ensemble. The hydrometeor ensemble variance is considered as zero below the threshold $10^{-24} (\text{kg} \cdot \text{kg}^{-1})^2$.

Most entries of the averaged covariance matrix are zero. These zero covariances are associated with zero

variances, for vertical levels where all ensemble members have zero hydrometeor content. Indeed, each hydrometeor variable cannot live outside a specific zone of residence, determined by the thermodynamic environment. For this particular case, the residence zone is below 700 hPa for rain, between 950 and 650 hPa for liquid cloud water, and between 750 and 200 hPa for ice crystals.

A positive covariance is observed between specific humidity around 700 hPa and rain content below (Figure 2a). This can be related to the fact that, in the ensemble, wet conditions are more often than not associated with rainy levels at lower altitudes.

A dipole structure can be observed for rain and temperature in Figure 2c. We deduce that, in the ensemble, rain is associated with temperatures lower than average below the rainy levels and warmer than average above. This could be interpreted, for instance, as signatures of latent heat exchanges, cooling the air in nonsaturated regions where the falling rain provokes some evaporation.

Other space-averaged covariances could be computed, for instance with other hydrometeor variables. These covariances are usually positive between specific humidity and hydrometeor variables, and between different hydrometeor variables (not shown). For the weather situation of Figure 2, typical values of maximum covariance between water-related variables are all of order $10^{-8}$ to $10^{-9} (\text{kg} \cdot \text{kg}^{-1})^2$, similarly to Figure 2a,b.

Structure functions that are not averaged over the domain present structures that may depend on the flow. This can be evidenced by looking at the increment that would be obtained from a single observation of specific humidity, for instance. This increment is proportional to the corresponding column of the **B** matrix, also known as the structure function. An illustration over Corsica is shown in Figure 3, after localization and normalization, so that the specific humidity “increment” is $1 \text{ g} \cdot \text{kg}^{-1}$ at the point of pseudo-observation. The position of the observation is chosen in a zone of nonzero cloud-water variance, to enable the generation of cloud-water increments. The horizontal localization length used here (measured as a Daley length, as in M18) is fixed to 150 km, and the vertical localization length is $0.4 \ln(\text{Pa})$. Since there are nonzero cross-covariances between specific humidity and cloud water, cloud-water increments are obtained.

Overall, the cloud-water increment looks realistic, does not extend above 680 hPa, and is mostly concentrated over land. The positive specific humidity “increment” is locally associated with positive increments of cloud water, consistent with Figure 2b. Some negative increments are observed further from the pseudo-observation point.

The rain increment associated with that observation is similar to the cloud-water increment, though with a twice reduced maximum (not shown).
<figure><img src="image_6.png" /><figcaption>**FIGURE 3** Localized increment of cloud water $q_l$ for a single pseudo-observation of specific humidity $q_v$ at the starred point, on May 21, 2018, 1500 UTC zoomed over Corsica. The increment is normalized to be $1 \text{ g} \cdot \text{kg}^{-1}$ for $q_v$ at the pseudo-observation point. (a) The increment at the model level of the pseudo-observation (800 hPa at the observation position). The section in figure (b) is performed along the latitude of the pseudo-observation (42.35°N) in (a). [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_7.png" /></figure>

To illustrate that hydrometeor increments can be obtained from observations of other conventional variables, we present in Figure 4 the structure function

associated with a pseudo-observation of temperature in the center of the domain. The temperature observation is located at 900 hPa, at a grid point with zero cloud-water
<figure><img src="image_8.png" /></figure>

**FIGURE 4** Vertical section of localized hydrometeor increments for a single pseudo-observation of temperature *T* at the starred point, on May 21, 2018, 1500 UTC. The increment is normalized to be 1 K for *T* at the pseudo-observation point (latitude 46.58°N). Top: cloud water *q*₁; bottom: rain *q*₂.
[Colour figure can be viewed at wileyonlinelibrary.com]

variance. The cross-covariances of temperature at this grid point with cloud water at other grid points are not zero, however. This allows for the generation of a cloud-water increment, 100–200 hPa above the pseudo-observation point. The associated rain increment is also shown in Figure 4 for illustrative purposes.

Globally, these illustrations show that hydrometeor variables have been added correctly in the ensemble B matrix. As expected, the structure functions are flow-dependent and contain some realistic information about the equilibria in the ensemble.

**3.3 | Visualization of some analysis increments**

Adding hydrometeor control variables in a 3DEnVar without adding hydrometeor-sensitive observations does not alter the conventional increments. New increments are obtained for the (unobserved) hydrometeor part of the control vector through background-error covariances, but without modifying the increment of the (observed) conventional part. Note that this is a consequence of the minimization being unconstrained. If the extension of the control variable comes with a positivity constraint for the hydrometeor analysis state, conventional increments can be altered, even though there are no actual hydrometeor

observations. Here, spurious negative hydrometeor contents produced by the analysis are set back to zero before the time integration of the forecast model.

In order to illustrate the possibility of obtaining consistent hydrometeor increments, we present the results of an analysis performed during the passing of a rain band over Southeastern France. Among the assimilated data, radar reflectivities are assimilated indirectly as pseudoprofiles of relative humidity, following the Caumont *et al.* (2010) approach. Figure 5 shows that we do obtain hydrometeor increments, for rain in that case, even without any direct assimilation of hydrometeor observations.

The innovations of conventional variables in the observation space are projected to analysis increments in model space, for both conventional variables (such as temperature, Figure 5d) and hydrometeor variables (such as rain, Figure 5c). The rain increment is added to the rain background (Figure 5a) to yield a new rain analysis (Figure 5b). Furthermore, not only are hydrometeor increments obtained, but these increments are substantial enough to affect the precipitating system strongly by extending it southward and by drying it in its northern part. As displayed in Figure 5e, such corrections are consistent with the observed radar reflectivities, which indeed display values above 40 dBZ along the main line where $q_r$ increments are highest and smaller values where $q_r$ is removed.
**FIGURE 5** Rain and temperature fields through the assimilation, altitude 917 hPa, validity date April 25, 2019, 1800 UTC. (a) Rain background, (b) rain analysis, (c) rain increment, (d) temperature increment, and (e) reflectivity composite from the French ARAMIS Radar network. Nonrealistic negative analyzed rain contents are set back to zero in panel (b). [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

# 4 | NUMERICAL EXPERIMENTS

## 4.1 | Experimental setup

### 4.1.1 | General presentation of the experiments

Adding hydrometeor control variables in the AROME 3DEnVar is expected to improve the forecast skill, especially at the beginning of the forecasts. This would be a natural consequence of improved hydrometeor initial conditions, though previous studies always showed it in combination with the impact of new observations, usually for precipitation forecasts, and in different experimental settings (MLEF for Zhang et al., 2013; EnKF for Johnson

et al., 2015; 3DEnVar with hybrid covariances for Meng et al., 2019, for instance).

This impact of hydrometeor control variables has been assessed by performing DA experiments over a three-month period in summer 2018. The period starts with the second half of May 2018, which was marked by intense convective activity over France. 12-hour forecasts were performed every 3 hr. The first forecast starts on May 13, 2018 at 0300 UTC, and the last one on August 15, 2018 at 2100 UTC. The length of the period implies that we have a reasonably large number of forecasts (94 or 95 per base hour), which increases the chance to obtain differences that are statistically significant.

The experiments use AROME-France as a NWP model, AROME 3DEnVar (Section 2.3.2) for analysis steps, and
AROME-EDA (Section 2.2) to provide the ensemble perturbations. The AROME EDA runs with a 3.2-km horizontal resolution with a hydrostatic formulation. To be consistent with this configuration, and to avoid a possible spin-up induced by interpolation to a 1.3-km grid, the AROME-France deterministic DA and forecast system are configured with same spatial resolution and same cycling frequency (3 hr).

Common localization lengths (measured as Daley lengths) are fixed to 150 km horizontally and 0.4 ln(Pa) vertically. These values are comparable with the choices of M18 and Caron et al. (2019). They have no reason to be identical, however, since the conditions of our experiments differ substantially from these previous studies: we use a 3.2-km resolution instead of 3.8 km, 50 members instead of 25, and a summer test period instead of a winter one. The conditions of the present study are actually closer to those of Destouches et al. (2020). Our choices of localization for the conventional variables are consistent with the values they diagnose over eight events of the year 2018. The same localization is applied to hydrometeor and conventional variables, as explained in Section 3.1.

The localization parameters have not been optimized manually over the test period, so that the localization scheme may not be at its optimally tuned behavior. This lack of tuning does not bring into question the validity of the conclusions drawn here, as shown by independent experiments with the scale-dependent localization of Buehner (2012) yielding the same conclusions.

### 4.1.2 | Configurations tested for hydrometeors and cycling

Three different configurations of AROME 3DEnVar are tested, differing in their use or not of hydrometeor control variables and their cycling. The first configuration is the control experiment (CTRL). It does not include hydrometeor control variables in the assimilation and runs in a continuous data assimilation cycle, each 3-hr forecast serving as a background state for the next analysis.

The two remaining configurations include hydrometeor control variables, but differ in cycling. In continuous DA cycles, the imbalances introduced by the analysis can propagate through the successive cycles if they are not eliminated by the NWP model in one cycle length. This may degrade the forecast skill and limits the minimum cycle length (see, for instance, Brousseau et al., 2016 about the DA cycle length for AROME). Similar issues might occur for continuous DA with hydrometeor control variables, for instance if cutting out negative hydrometeor contents in the analysis introduces some bias in the forecasts. To isolate this possible effect, the second configuration

(HMNCyc) includes hydrometeor control variables in a noncycled data assimilation framework. Every 3 hr, the HMNCyc analysis starts from the same background state as CTRL. The forecasts initialized from these HMNCyc analyses do not serve as backgrounds for posterior analyses. The last configuration includes hydrometeor control variables in a continuous data assimilation cycle (HMCyc). Configurations HMCyc and CTRL start from a common point at the beginning of the three-month period and are run independently afterwards.

### 4.1.3 | Observations and minimization

Our experiments assimilate most of the observations that are considered in the operational AROME-France suite. They include conventional observations, Doppler wind measurements, pseudoprofiles of relative humidity derived from radar observations (Caumont et al., 2010; Wattrelot et al., 2014), and radiances from geostationary and polar-orbiting satellites. The only observations of the operational suite that were not assimilated are satellite observations from scatterometers and from the satellite instruments Global Precipitation Measurement (GPM) Microwave Imager (GMI), Cross-track Infrared Sounder (CrIS), and Advanced Technology Microwave Sounder (ATMS).

In the version of the AROME-3DEnVar scheme used for this study, the variational Bias Correction scheme (VarBC) was not implemented for satellite measurements from the Infrared Atmospheric Sounding Interferometer (IASI) and the Spinning Enhanced Visible and InfraRed Imager (SEVIRI: Schmetz et al., 2002). When available, VarBC usually provides an adaptive correction of systematic biases in satellite observations, using some predictors derived mainly from the background state (Auligné et al., 2007). As a workaround, the regression parameters determined by VarBC during the minimizations of the operational AROME-France suite are reused here to correct the biases in IASI and SEVIRI observations.

The global cost function was solved using one iteration of the outer loop. The inner loop is based on a B-preconditioned conjugate gradient scheme with 40 iterations, whether hydrometeors are present in the control variable or not.

## 4.2 | Verification measures

### 4.2.1 | Reference observations

The impact on cloud forecasts is evaluated by comparison with brightness temperatures observed by SEVIRI,
<figure><img src="image_14.png" /></figure>

on board the Meteosat Second-Generation satellite. Comparisons for the 10.8-μm window channel are performed. At this wavelength, thick liquid and ice clouds become opaque. This property is interesting to evaluate the positioning and intensity of convective events, with lower brightness temperatures being associated with higher clouds. The frequency of the measurements is 15 min, and the spatial resolution is around 4 km over France. To translate forecast fields into observation space, the (Radiative Transfert for TIROS Operational Vertical Sounder (RTTOV) model was used (Saunders et al., 2018). The forecasts to be evaluated are extracted at the time step closest to the observation time and interpolated to the nearest observation point.

The impact on precipitation is evaluated by comparison with gridded retrievals of precipitation rates over France, at a frequency of 15 min and a spatial resolution of 1 km. These retrievals are derived from a combination of rain gauge and radar measurements (ANalyse par spaTialisation hOraire des PrEciapiations (ANTILOPE) product: Champeaux et al., 2009). The corresponding AROME-France forecasts at 3.2 km are bilinearly interpolated to the finer resolution of ANTILOPE.

### 4.2.2 | Scores

We used neighborhood-based scores to evaluate the adequacy between forecasts and observations. These scores allow for some displacement errors, thus avoiding the double penalty problem of traditional categorical scores (Ebert, 2008). The main score we use here is the Fractions Skill Score (FSS, Roberts and Lean, 2008), frequently employed for the verification of precipitation at convective scales. Recent uses of FSS in the context of hydrometeor initialization include Duda et al. (2019) and Hawkness-Smith and Simonin (2021), for instance.

Observations and forecast fields can be binarized into fields of events or nonevents, where a given threshold is exceeded or not. Let s be a given domain size, of the order of a few grid points. Let $v_i$, $i = 1 \cdots N_v$, be all the square subgrids with size $s \times s$ in the domain (with overlapping). The Fractions Skill Score compares the frequencies $f_{\text{obs}}(v_i)$ and $f_{\text{forecast}}(v_i)$ of observed and forecast events in each spatial neighborhood $v_i$:

$$FSS = \frac{FBS_{\text{ref}} - FBS}{FBS_{\text{ref}}}, \qquad (12)$$

$$FBS = \sum_{i=1}^{N_v} [f_{\text{obs}}(v_i) - f_{\text{forecast}}(v_i)]^2, \qquad (13)$$

$$FBS_{\text{ref}} = \sum_{i=1}^{N_v} [f_{\text{obs}}(v_i) + f_{\text{forecast}}(v_i)]^2, \qquad (14)$$

where FBS stands for Fractions Brier Score. The FSS varies between 0 (no skill) and 1 (perfect skill).

To complement the FSS and provide more insight on the evaluated impact, we also computed contingency tables with error compensation (Stein and Stoop, 2019). These tables include a tolerance to displacement errors by allowing pairs of false alarms and misses within some prescribed neighborhoods to be counted as pairs of hits and correct rejections. The tables we obtain have the same entries as usual, that is, counts of hits (A), false alarms (B), misses (C), and correct rejections (D). All traditional scores derived from contingency tables can be used in this context. We restrain ourselves here to the false-alarm ratio (FAR), probability of detection (POD), and frequency bias (FBI):

$$FAR = B/(A + B), \qquad (15)$$

$$POD = A/(A + C), \qquad (16)$$

$$FBI = (A + B)/(A + C). \qquad (17)$$

Interestingly, the method of Stein and Stoop (2019) ensures that the frequency bias does not depend on the chosen neighborhood size.

The retained size for spatial neighborhoods is 52 km, which is similar to values typically used for the verification of operational AROME-France forecasts (Amodei et al., 2015).

The threshold values for the SEVIRI 10.8-μm channel extend from 230-290 K in steps of 10 K. This range of values was chosen based on frequencies of occurrence, to avoid rare events the FSS is not supposed to deal with. Note that an event is said to occur for SEVIRI observations where brightness temperatures are lower than the threshold value, which corresponds to higher clouds. For accumulated precipitation over a period of 15 min, threshold values are fixed to 0.1 mm (occurrence of precipitation), 0.5 mm (medium intensity precipitation, around 45% of rainy events), and 1 mm (heavy precipitation, around 17% of rainy events).

Following Hamill (1999), the contingency tables are summed over all forecasts before computing scores such as FAR or POD. Similarly, the FBS and FBS\_ref are summed over all forecasts considered before computation of the FSS.
### 4.2.3 | Statistical significance

The statistical significance of the results is assessed from differences of paired up forecasts. For a given experiment XP, the scores are shown in terms of relative improvements compared with the control experiment CTRL:

$$ \Delta \text{FSS} = \frac{\text{FSS}(\text{XP}) - \text{FSS}(\text{CTRL})}{\text{FSS}(\text{CTRL})} \quad (18) $$

We detail here the way the ΔFSS and its confidence interval are built.

FSS and ΔFSS are computed independently for each base hour (0000 UTC, 0300 UTC, etc.). Across the evaluation period, there are 94–95 forecasts starting from each base hour, and as many values of FBS and FBSref are computed. This gives a time series of quadruplets (FBS(XP); FBSref(XP); FBS(CTRL); FBSref(CTRL)), hereafter the FBS time series. Instead of computing a time series of FSS, the FBS time series is first aggregated by summation over the evaluation period to yield an aggregated quadruplet, and from there a pair (FSS(XP), FSS(CTRL)) (using Equation 12). The relative FSS is then computed using Equation (18).

Confidence intervals for ΔFSS can be estimated using a bootstrap method (Gilleland et al., 2018). For a couple of experiments (XP, CTRL), the FBS time series is sampled with replacement $N_B = 5,000$ times. From each one of the $N_B$ resulting time series, a ΔFSS is computed. The $N_B$ values of ΔFSS form a distribution, from which 95% confidence intervals are estimated using the “basic interval” method (see Gilleland et al., 2018).

When assessing the significance of the result “ΔFSS is positive”, unilateral confidence intervals of the kind $\Delta \text{FSS} \in [\alpha, +\infty)$ are estimated. The statistical significance is assessed in this case if $\alpha$ is positive.

### 4.3 | Verification against conventional observations

Though the focus of this article is on forecasting clouds and precipitation, we first report briefly the impact of hydrometeor variables on forecasts of conventional variables. For these variables, the forecast skill is assessed by pointwise comparison with independent observations. Ground-based weather stations provide measurements for the surface variables: 2-m temperature ($T$ 2M) and relative humidity (HU 2M), 10-m wind ($U/V$ 10M), and mean sea-level pressure (MSLP). Altitude observations (ALT) consist of various measurements: aircraft wind ($U/V$) and temperature ($T$), wind profilers, radar Doppler winds, ground-based GNSS delays, and SEVIRI observations in

the water-vapour channel (HU). Relative improvements of root-mean-square error (RMSE) against control experiments for the STD localization are visualized in Figure 6 through the mean of the “score cards” introduced by M18 and used, for instance, by Caron et al. (2019) and Michel and Brousseau (2021).

The forecasts are improved globally when adding hydrometeor control variables, especially for surface variables: surface pressure, temperature, and humidity at 2 m. For surface pressure in particular, the improvements persist up to the end of the forecast range (12 hr), almost systematically (cycled and noncycled configurations, for morning or afternoon forecasts). In contrast, 10-m wind forecasts suffer from a neutral to negative impact of hydrometeor control variables. This impact is significantly negative in the first hours of forecasts and in the cycled configuration.

For altitude variables (wind, temperature, humidity), there is no clear sign of generalized improvement or deterioration. In cycled configurations, the impact tends to be more negative.

The impact on conventional observations is globally neutral or small, of the order of 1% when RMSEs are not averaged per variable and forecast range. Hereafter we instead focus on the impact of hydrometeor control variables on precipitation and cloud-cover forecasts.

### 4.4 | Impact on cloud-cover forecasts

Scores from comparison with SEVIRI observations are shown in Figure 7 for 12 hr-long forecasts. For the sake of conciseness, only forecasts starting from 1200 UTC are shown. The choice of this base hour allows us to include convective events of the afternoon in the forecast range.

For the noncycled forecasts HMNCyc, significant improvements over CTRL are observed in the first hours of the forecasts, for most thresholds of brightness temperature. The positive impact lasts longer for warmer brightness temperatures, up to 9 hr for the 290-K threshold. The scores become globally neutral afterwards. This behavior confirms the positive impact of allowing the conventional increments to propagate to hydrometeor variables through extradiagonal blocks of the B matrix.

The impact of hydrometeor variables is slightly different in the cycled configuration, where the 3-hr forecasts are used as backgrounds for the following assimilation step. The global structure of the relative FSS for HMCyc in Figure 7 is similar to the noncycled case, with a positive impact that decays with the forecast lead time. Two main differences should be highlighted, however.
<figure><img src="image_15.png" /><figcaption>FIGURE 6 Impact on RMSEs of introducing hydrometeor control variables. Top: forecasts starting at 0000 UTC; bottom: forecasts starting at 1200 UTC. Left: noncycled configuration; right: cycled configuration. Upward-pointing triangles indicate an improvement in RMSE. Triangles are filled when the impact is statistically significant (95% confidence intervals estimated with a bootstrap method). [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

* The amplitude of the positive impact is similar to that for the noncycled configuration, but with fewer statistically significant changes.

* Some negative impact is observed in the second half of the forecasts. The deterioration is displayed in a diagonal structure in the figure, suggesting a propagation

from the higher altitudes (lower brightness temperatures) to the bottom layers of the atmosphere (higher temperatures).

The first point can be investigated by looking more precisely at the confidence intervals. This is done in Figure 8,
<figure><img src="image_16.png" /><figcaption>FIGURE 7 Relative score of FSS against SEVIRI observations compared with control. Noncycled (top) and cycled (bottom) experiments. The aggregated forecasts all start at 1200 UTC. Significant improvements or deteriorations are marked with “+” or “−” signs, respectively. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

for a 250-K threshold and for centered confidence inter-
vals, which are indeed wider in the cycled configuration.
This is related to a larger variability of the relative FSS over
the evaluation period. Compared with the noncycled con-
figuration, the cycled one yields better forecasts for some
situations but worse forecasts for some others. This results
in an overall nonsignificant impact.

For the sake of completeness, we conclude this section
by detailing the scores for forecasts starting at other base
hours (in UTC hereafter). In the noncycled configuration,
most scores are positive or neutral, with positive scores
especially in the first forecast hours or for large brightness
temperature thresholds (as for base hour 12). Base hour
18 is the exception, with positive scores concentrated in
the forecast lead times 0–3 and 9–12 hr, and neutral scores
elsewhere. Negative scores are in the minority, in the sec-
ond half of forecasts only, and for forecast starting at base
hours 3, 6, 9, 12, and 15.

In the cycled configuration, the impact of hydrome-
teor control variables is less clear, with a higher number of
negative scores, a lower number of positive scores, and a
larger uncertainty in confidence intervals. Except for base
hour 12, negative scores are concentrated at low brightness
temperature thresholds (230–250 K). This implies that the
“structure” of the scores is often changed, with the separa-
tion between the first and second half of the forecasts in the
noncycled configuration shifting to a separation between
low and high thresholds of brightness temperature. The

very first hours of forecasts usually remain positively
impacted, however. The negative scores are usually in a
clear minority compared with positive scores, except for
base hours 3 (23 negative versus 39 positive) and base
hour 6 (17 negative versus 11 positive). We finally note
that cycling the hydrometeor analyses positively enhances
the forecast skill for forecasts valid between 2300 and
0100 UTC and starting at 1500, 1800, and 2100 UTC.

4.5 | Impact on precipitation forecasts

The very short term impact of hydrometeor control vari-
ables on precipitation is shown in Figure 9, for the first
15 min of the forecasts. These FAR/POD diagrams allow
us to visualize in a single figure the false-alarm ratio, prob-
ability of detection, and frequency bias (FBI). The CTRL
experiment displays a strong negative bias for all precip-
itation thresholds, with frequency biases of the order of
0.5–0.8. This deficit of precipitation at the very beginning
of the forecasts is an issue that is also found in the oper-
ational AROME-France suite (Pierre Brousseau, personal
communication). The FBI values are clearly improved
with hydrometeor control variables, with frequency biases
in the 0.8–1.1 range.

These changes in the FBI are associated with increases
of the POD, and slight increases in the FAR for low pre-
cipitation thresholds. This deterioration of the latter is
<figure><img src="image_17.png" /><figcaption>**FIGURE 8** Confidence intervals for relative score of FSS against SEVIRI observations, compared with control, for noncycled and cycled experiments. The data are the same as in Figure 7, but only threshold 250 K is shown here, and with centered confidence intervals. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 9** FAR-POD table against precipitation observations. Left: run of 0000 UTC; right: run of 1200 UTC. Contour lines of the frequency bias FBI are displayed. Precipitation thresholds are written alongside each CTRL point. The related HMNCyc and HMCyc experiment are linked by solid lines to their control. Unbiased forecasts would be on the diagonal of the figure, perfect forecasts would be in the upper-left corner. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

statistically significant for the event “occurrence of precipitation” only. Such deterioration remains small compared with the increase in POD, so that the overall impact of adding hydrometeor control variables is still significantly positive for all thresholds. Similar impacts are observed for forecasts starting at other base hours.

This positive impact is expected to decay with forecast range, similarly to the impact on clouds. The FSS of the first three hours of forecasts are shown in Figure 10 for the 0000 and 1200 UTC forecasts (comments hereafter apply to all base hours). For the control experiment, the behavior of the FSS curves is similar across all thresholds, for both base hours. A short period of spin-up is first observed in the first 45 min of forecasts. This period is evidenced by low FSS values, which increase rapidly to a maximum FSS for the

30–45 min time window in most cases. The FSS decreases almost steadily afterwards and for the rest of the 3 hr.

Adding hydrometeor control variables results in a consistently large increase of the FSS at the beginning of the forecasts, be it in cycled or noncycled configuration. The FSS improvements for the first 15 min of forecasts vary from 0.06 points (base hour 15, thresholds 0.5 and 1 mm) to 0.17 points (base hour 0, 1 mm). In relative quantities, the increases in FSS range from 10% (base hour 12, 0.1 mm; base hour 15, 0.5 mm) to 39% (base hour 0, 1 mm).

A period of spin-up can still be seen in most situations, but with less pronounced degradations of the scores. Also, the spin-up time is shorter, with FSS maxima reached for the 15–30 min window or before, that is, at least 15 min sooner than in the control experiment.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 10 Evolution of FSS for precipitation in the first three hours of forecasts. The forecasts start at 0000 UTC (left) and 1200 UTC (right). The precipitation thresholds are 0.1 mm per 15 min (top row), 0.5 mm per 15 min (middle row), and 1 mm per 15 min (bottom row). Significant impacts compared with the CTRL experiment are marked by triangles. On the x-axes, data points are placed at the end of their associated 15-min window. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Over the 3-hr forecast range, the impact is globally positive, and similar between both cycled and noncycled configurations. However, the significance can usually be assessed for longer in the noncycled configuration, with significant improvements extending up to 2 hr 15 min on

average for the occurrence of precipitation (1 hr 15 min for the largest threshold). Several cases even display improvements extending up to 2 hr 45 min (base hours 12 and 18 for occurrence of precipitation, 15 for the largest threshold). In the cycled configuration, the improvements only
extend up to 1 hr 15 min on average for the occurrence of precipitation, and 45 min for the largest threshold.

This behavior is similar to what was observed for cloud cover, namely an improvement of scores on average, but with a larger variability that prevents us from assessing the significance of the impact in cycled experiments.

5 | DISCUSSION AND CONCLUSION

Hydrometeors have been added as control variables of the AROME 3DEnVar scheme. Through their inclusion in the ensemble background-error covariance matrix, conventional increments can be propagated to hydrometeor variables, even without any direct assimilation of hydrometeor-sensitive observations.

The impact of this change has been evaluated on short-term forecasts over a three-month summer period, by comparison with real observations related to clouds and precipitation. The impact is globally positive, with positive effects on cloud cover extending up to 9 hr. The impact on precipitation is especially satisfactory, with a clear reduction of the duration and intensity of the spin-up period. The adjustment time before the FSS reaches its maximum was already short in this setting (of the order of 30–45 min). With hydrometeor control variables, it is reduced to 15–30 min. Afterwards, the positive impact of updating hydrometeor variables is visible in the first hour of forecasts, and up to the third hour in some cases of non-cycled DA. Note that the same positive impact of adding hydrometeor control variables is observed when reproducing these experiments with the scale-dependent localization scheme of Buehner (2012). This shows that our results are not sensitive to the specific choices we made for covariance localization.

These improvements can be partly explained by a bias reduction in precipitation at the beginning of the forecasts due to an increase of the hydrometeor content in the analysis. A better positioning of precipitation patterns at the beginning of the forecasts has been assessed from a better FSS, but also from higher probabilities of detection, at the expense of only small increases in the false-alarm ratios. This means that some relevant information is indeed propagated from conventional observations to hydrometeor control variables, through the cross-covariances in the **B** matrix.

The results in the cycled configuration, where 3-hr forecasts are used as background for the next analysis cycle, exhibit a larger variability than those in the noncycled configuration. Some negative scores on cloud cover are also observed in the second half of the forecast or for the highest clouds. Another possible effect of cycling could

result from the current treatment of negative hydrometeor contents in the analysis, which consists of thresholding them to zero. This implies that the increment has a nonzero mean, and possibly introduces a slight positive bias (as anticipated in more general situations by Janjić *et al.*, 2014, for instance). This bias can then build up during successive cycles. Slightly higher frequency biases are indeed observed in some cases with the cycled configuration (see for instance Figure 9), but the behavior is not generalized, so that no significant conclusion can be drawn on this. Further research is needed to investigate the cause of these behaviors.

Comparing our findings with those of previous studies of the DA community is not straightforward. To the best of our knowledge, the studies including hydrometeor control variables in convective-scale forecasting systems and with real observations all focused on the added value of new hydrometeor-sensitive observations, rather than on hydrometeor control variables *per se*. In consequence, we are not aware of any other study assessing the impact of updating hydrometeor control variables independently of new observations, or a new DA scheme. One of the most similar studies may be the two-week experiments of continuous data assimilation performed by Meng *et al.* (2019). They assess the combined impact of cloud-water path observations and hydrometeor ensemble covariances in a 3DEnVar with hybrid covariances. Among other things, they report improvements in 12-hr forecasts of temperature and specific humidity at low altitudes, consistent with our findings. The impact they observe on precipitation forecasts is equally positive and persists over the 12-hr forecast range. This longer persistence is likely to be related to the impact of the new observational data they introduce, which lets us assume that a stronger persistence of the positive impact is to be expected in our case when assimilating hydrometeor-sensitive observations.

As the first published study about the impact of hydrometeor control variables in AROME, this study has some limitations. The impact of hydrometeor control variables has been measured only for a summer period and for a horizontal resolution of 3.2 km (against 1.3 km for the operational suite). Only 50 members were used in the ensemble of data assimilations, though more members may be needed to estimate hydrometeor-related covariances correctly. Necker *et al.* (2020), for instance, recommend using 200 members for estimation of correlations involving hydrometeors. Using time-lagging as in Caron *et al.* (2019) may also be an option to increase the ensemble size, as significant improvements have been obtained in their study, especially against conventional observations. However, the possibly positive impact of time-lagging on the estimation of hydrometeor forecast-error covariances may be limited by their higher spatial and temporal
variability and lower predictability. Finally, no specific treatment has been applied to deal with non-Gaussian error statistics, no specific optimization of localization has been made to deal with hydrometeor ensemble covariances, no hybridization with a static term has been tested in the **B** matrix, and no hydrometeor observation has been assimilated directly.

These last few points are as much possibilities of improvement for later work. They make us consider the positive results of this study as very encouraging, since most of the impact of hydrometeor control variables is to be expected from the direct assimilation of hydrometeor observations, which are currently assimilated indirectly (and suboptimally) or rejected from the DA process. Ongoing studies are being conducted at Météo-France to test the impact of direct assimilation of radar reflectivities, which are mainly sensitive to liquid rain, in the particular EnVar framework described in this article. Direct assimilation of variables deduced from dual-polarimetric radar, conducted in a 1D-EnVar framework by Thomas et al. (2020), has also shown promising results. However, as 3DEnVar schemes rely only on ensemble-derived background-error covariances, the variances can be null if all ensemble members have a zero hydrometeor content at some location. This may pose a problem when observations linked to hydrometeors are available in these areas. Adding a static climatological component to the background-error covariance matrix would be a possible way to deal with this issue. Adding hydrometeors in this climatological component would be technically feasible, as discussed by Michel et al. (2011) and Descombes et al. (2015). Another more simple way would consist of imposing a minimum value of variances over the domain. Such issues will be considered in our future research on improving the analysis of hydrometeors and the assimilation of hydrometeor-sensitive observations in AROME.

Finally, note that the changes tested here on the initialization of hydrometeor variables could be extended to any prognostic variable of NWP models that is not part of the control variable of the data assimilation. For AROME, this includes the nonhydrostatic component of pressure, the vertical divergence of the vertical wind, and the turbulent kinetic energy. The inclusion of these variables may seem less critical, since they are not associated with any direct observation that could be assimilated, in contrast to hydrometeor variables. Nevertheless, the results of this study suggest that having a greater consistency between the variables used in the forecast and in the analysis could be beneficial even for unobserved variables. We can reasonably hypothesize that including all prognostic variables in the analysis would reduce adjustment phenomena such as spin-up at the start of each forecast.

## AUTHOR CONTRIBUTIONS

**Mayeul Destouches:** formal analysis; investigation; writing – original draft; writing – review and editing.
**Thibaut Montmerle:** conceptualization; investigation; methodology; software; supervision.
**Yann Michel:** conceptualization; formal analysis; investigation; methodology; software; supervision.
**Jean-François Caron:** software.

## ACKNOWLEDGEMENTS

We are grateful to Centre de Météorologie Spatiale de Lannion for providing us the SEVIRI data used in this study, and to Fabien Stoop for providing some code to compute neighborhood-based scores. We also thank an anonymous reviewer whose remarks greatly helped to improve this article. The PhD thesis of Mayeul Destouches was funded by the French Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation, through a grant for alumni from the École polytechnique. The Les Enveloppes Fluides et l'Environnement (LEFE) programme of the Institut National des Sciences de l'Univers—Centre National de la Recherche Scientifique (INSU-CNRS) funded part of this research within the framework of the Hydre project.

## CONFLICT OF INTEREST STATEMENT

The authors declare having no conflict of interest.

## ORCID

Mayeul Destouches https://orcid.org/0000-0001-6745-9009

Thibaut Montmerle https://orcid.org/0000-0003-3999-1633

## REFERENCES

Aksoy, A., Dowell, D.C. and Snyder, C. (2009) A multcase comparative assessment of the ensemble Kalman filter for assimilation of radar observations. Part I: storm-scale analyses. *Monthly Weather Review*, 137(6), 1805–1824. https://doi.org/10.1175/2008mwr2691.1.

Aksoy, A., Dowell, D.C. and Snyder, C. (2010) A multcase comparative assessment of the ensemble Kalman filter for assimilation of radar observations. Part II: short-range ensemble forecasts. *Monthly Weather Review*, 138(4), 1273–1292. https://doi.org/10.1175/2009mwr3086.1.

Amezcua, J. and Leeuwen, P.J.V. (2014) Gaussian anamorphosis in the analysis step of the EnKF: a joint state-variable/observation approach. *Tellus A: Dynamic Meteorology and Oceanography*, 66(1), 23493. https://doi.org/10.3402/tellusa.v66.23493.

Amodei, M., Sanchez, I. and Stein, J. (2015) Verification of the French operational high-resolution model AROME with the regional brier probability score. *Meteorological Applications*, 22(4), 731–745. https://doi.org/10.1002/met.1510.

Auligné, T., McNally, A.P. and Dee, D.P. (2007) Adaptive bias correction for satellite data in a numerical weather prediction system. *Quarterly Journal of the Royal Meteorological Society*, 133(624), 631–642. https://doi.org/10.1002/qj.56.
Auligné, T., Lorenc, A.C., Michel, Y., Montmerle, T., Jones, A., Hu, M. and Dudhia, J. (2011) Toward a new cloud analysis and prediction system. *Bulletin of the American Meteorological Society*, 92(2), 207–210. https://doi.org/10.1175/2010bams2978.1.

Bannister, R.N., Chipilski, H.G. and Martinez-Alvarado, O. (2019) Techniques and challenges in the assimilation of atmospheric water observations for numerical weather prediction towards convective scales. *Quarterly Journal of the Royal Meteorological Society*, 146(726), 1–48. https://doi.org/10.1002/qj.3652.

Berre, L., Varella, H. and Desroziers, G. (2015) Modelling of flow-dependent ensemble-based background-error correlations using a wavelet formulation in 4D-Var at Météo-France. *Quarterly Journal of the Royal Meteorological Society*, 141(692), 2803–2812. https://doi.org/10.1002/qj.2565.

Bocquet, M., Pires, C.A. and Wu, L. (2010) Beyond Gaussian statistical modeling in geophysical data assimilation. *Monthly Weather Review*, 138(8), 2997–3023. https://doi.org/10.1175/2010mwr3164.1.

Bouttier, F., Vié, B., Nuissier, O. and Raynaud, L. (2012) Impact of stochastic physics in a convection-permitting ensemble. *Monthly Weather Review*, 140(11), 3706–3721. https://doi.org/10.1175/mwr-d-12-00031.1.

Brousseau, P., Berre, L., Bouttier, F. and Desroziers, G. (2011) Background-error covariances for a convective-scale data-assimilation system: AROME-France 3D-Var. *Quarterly Journal of the Royal Meteorological Society*, 137(655), 409–422. https://doi.org/10.1002/qj.750.

Brousseau, P., Seity, Y., Ricard, D. and Léger, J. (2016) Improvement of the forecast of convective activity from the AROME-France system. *Quarterly Journal of the Royal Meteorological Society*, 142(699), 2231–2243. https://doi.org/10.1002/qj.2822.

Buehner, M. (2005) Ensemble-derived stationary and flow-dependent background-error covariances: evaluation in a quasi-operational NWP setting. *Quarterly Journal of the Royal Meteorological Society*, 131(607), 1013–1043. https://doi.org/10.1256/qj.04.15.

Buehner, M. (2012) Evaluation of a spatial/spectral covariance localization approach for atmospheric data assimilation. *Monthly Weather Review*, 140(2), 617–636 https://doi.org/10.1175/mwr-d-10-05052.1.

Buehner, M. and Shlyaeva, A. (2015) Scale-dependent background-error covariance localisation. *Tellus A: Dynamic Meteorology and Oceanography*, 67(1), 28027. https://doi.org/10.3402/tellusa.v67.28027.

Burgers, G., van Leeuwen, P.J. and Evensen, G. (1998) Analysis scheme in the ensemble Kalman filter. *Monthly Weather Review*, 126(6), 1719–1724. https://doi.org/10.1175/1520-0493(1998)126<1719:asitek>2.0;co.2.

Bénard, P., Vivoda, J., Mašek, J., Smolíková, P., Yessad, K., Smith, C., Brožková, R. and Geleyn, J.-F. (2010) Dynamical kernel of the Aladin-NH spectral limited-area model: revised formulation and sensitivity experiments. *Quarterly Journal of the Royal Meteorological Society*, 136(646), 155–169. https://doi.org/10.1002/qj.522.

Caron, J.-F., Michel, Y., Montmerle, T. and Arbogast, É. (2019) Improving background error Covariances in a 3D ensemble-Variational data assimilation system for regional NWP. *Monthly Weather Review*, 147(1), 135–151. https://doi.org/10.1175/mwr-d-18-0248.1.

Caumont, O., Ducrocq, V., Wattrelot, É., Jaubert, G. and Pradier-Vabre, S. (2010) 1D + 3DVar assimilation of radar reflectivity data: a proof of concept. *Tellus A: Dynamic Meteorology and Oceanography*, 62(2), 173–187. https://doi.org/10.1111/j.1600-0870.2009.00430.x.

Champeaux, J.-L., Dupuy, P., Laurantin, O., Soulan, I., Tabary, P. and Soubeyroux, J.-M. (2009) Les mesures de précipitations et l'estimation des lames d'eau à Météo-France : état de l'art et perspectives. *La Houille Blanche*, 5, 28–34. https://doi.org/10.1051/lhb/2009052.

Chen, L., Liu, C., Xue, M., Zhao, G., Kong, R. and Jung, Y. (2021) Use of power transform mixing ratios as hydrometeor control variables for direct assimilation of radar reflectivity in GSI en3dvar and tests with five convective storm cases. *Monthly Weather Review*, 149(3), 645–659. https://doi.org/10.1175/mwr-d-20-0149.1.

Clayton, A.M., Lorenc, A.C. and Barker, D.M. (2012) Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the met Office. *Quarterly Journal of the Royal Meteorological Society*, 139(675), 1445–1461. https://doi.org/10.1002/qj.2054.

Cohn, S.E. (1997) An introduction to estimation theory. *Journal of the Meteorological Society of Japan*, 75(1B), 257–288. https://doi.org/10.2151/jmsj1965.75.1b_257.

Courtier, P., Freydier, C., Geleyn, J.-F., Rabier, F. and Rochas, M. (1991) The Arpege project at Météo-France. In: *Seminar on Numerical Methods in Atmospheric Models*, Vol. II. Shinfield Park, Reading: ECMWF, pp. 193–232, https://www.ecmwf.int/node/8798.

Derber, J. and Rosati, A. (1989) A global oceanic data assimilation system. *Journal of Physical Oceanography*, 19(9), 1333–1347. https://doi.org/10.1175/1520-0485(1989)019<1333:agodas>2.0;co.2.

Descombes, G., Auligné, T., Vandenberghe, F., Barker, D.M. and Barré, J. (2015) Generalized background error covariance matrix model (GEN_BE v2.0). *Geoscientific Model Development*, 8(3), 669–696. https://doi.org/10.5194/gmd-8-669-2015.

Destouches, M., Montmerle, T., Michel, Y. and Ménétrier, B. (2020) Estimating optimal localization for sampled background-error covariances of hydrometeor variables. *Quarterly Journal of the Royal Meteorological Society*, 147(734), 74–93. https://doi.org/10.1002/qj.3906.

Dowell, D.C., Zhang, F., Wicker, L.J., Snyder, C. and Crook, N.A. (2004) Wind and temperature retrievals in the 17 May 1981 Arcadia, Oklahoma, supercell: ensemble Kalman filter experiments. *Monthly Weather Review*, 132(8), 1982–2005. https://doi.org/10.1175/1520-0493(2004)132<1982:watrit>2.0;co.2.

Ducrocq, V., Braud, I., Davolio, S., Ferretti, R., Flamant, C., Jansa, A., Kalthoff, N., Richard, E., Taupier-Letage, I., Ayral, P.-A., Belamari, S., Berne, A., Borga, M., Boudevillain, B., Bock, O., Boichard, J.-L., Bouin, M.-N., Bousquet, O., Bouvier, C., Chiggiato, J., Cimini, D., Corsmeier, U., Coppola, L., Cocquerez, P., Defer, E., Delanoë, J., Girolamo, P.D., Doerenbecher, A., Drobin-ski, P., Dufournet, Y., Fourrié, N., Gourley, J.J., Labatut, L., Lambert, D., Coz, J.L., Marzano, F.S., Molinié, G., Montani, A., Nord, G., Nuret, M., Ramage, K., Rison, W., Roussot, O., Said, F., Schwarzenboeck, A., Testor, P., Baelen, J.V., Vincen-don, B., Aran, M. and Tamayo, J. (2014) HyMeX-SOP1: the field campaign dedicated to heavy precipitation and flash flooding
in the northwestern Mediterranean. *Bulletin of the American Meteorological Society*, 95(7), 1083-1100. https://doi.org/10.1175/bams-d-12-00244.1.

Duda, J.D., Wang, X., Wang, Y. and Carley, J.R. (2019) Comparing the assimilation of radar reflectivity using the direct GSI-based ensemble-Variational (EnVar) and indirect cloud analysis methods in convection-allowing forecasts over the continental United States. *Monthly Weather Review*, 147(5), 1655-1678. https://doi.org/10.1175/mwr-d-18-0171.1.

Duruisseau, F., Chambon, P., Wattrelot, E., Barreyat, M. and Mahfouf, J.-F. (2019) Assimilating cloudy and rainy microwave observations from SAPHIR on board Megha Tropiques within the ARPEGE global model. *Quarterly Journal of the Royal Meteorological Society*, 145(719), 620-641. https://doi.org/10.1002/qj.3456.

Ebert, E.E. (2008) Fuzzy verification of high-resolution gridded forecasts: a review and proposed framework. *Meteorological Applications*, 15(1), 51-64 https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/met.25.

Errico, R.M., Bauer, P. and Mahfouf, J.-F. (2007) Issues regarding the assimilation of cloud and precipitation data. *Journal of the Atmospheric Sciences*, 64(11), 3785-3798. https://doi.org/10.1175/2006jas2044.1.

Evensen, G. (1994) Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. *Journal of Geophysical Research*, 99(C5), 10143. https://doi.org/10.1029/94jc00572.

Evensen, G. (2003) The ensemble Kalman filter: theoretical formulation and practical implementation. *Ocean Dynamics*, 53(4), 343-367. https://doi.org/10.1007/s10236-003-0036-9.

Fletcher, S.J. (2007) Implications and impacts of transforming log-normal variables into normal variables in VAR. *Meteorologische Zeitschrift*, 16(6), 755-765.

Gao, J., Xue, M. and Stensrud, D.J. (2013) The development of a hybrid EnKF-3DVAR algorithm for storm-scale data assimilation. *Advances in Meteorology*, 1-12, 2013. https://doi.org/10.1155/2013/512656.

Gilleland, E., Hering, A.S., Fowler, T.L. and Brown, B.G. (2018) Testing the tests: what are the impacts of incorrect assumptions when applying confidence intervals or hypothesis tests to compare competing forecasts? *Monthly Weather Review*, 146(6), 1685-1703. https://doi.org/10.1175/mwr-d-17-0295.1.

Gürol, S., Weaver, A.T., Moore, A.M., Piacentini, A., Arango, H.G. and Gratton, S. (2013) B-preconditioned minimization algorithms for variational data assimilation with the dual formulation. *Quarterly Journal of the Royal Meteorological Society*, 140(679), 539-556. https://doi.org/10.1002/qj.2150.

Hamill, T.M. (1999) Hypothesis tests for evaluating numerical precipitation forecasts. *Weather and Forecasting*, 14(2), 155-167. https://doi.org/10.1175/1520-0434(1999)014<0155:htfenp>2.0.co;2.

Hamill, T.M. and Snyder, C. (2000) A hybrid ensemble Kalman filter-3D variational analysis scheme. *Monthly Weather Review*, 128(8), 2905-2919. https://doi.org/10.1175/1520-0493(2000)128<2905:ahekfv>2.0.co;2.

Hamill, T.M., Whitaker, J.S. and Snyder, C. (2001) Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter. *Monthly Weather Review*, 129(11), 2776-2790. https://doi.org/10.1175/1520-0493(2001)129<2776:ddfobe>2.0.co;2.

Hawkness-Smith, L.D. and Simonin, D. (2021) Radar reflectivity assimilation using hourly cycling 4D-Var in the met Office unified model. *Quarterly Journal of the Royal Meteorological Society*, 147(736), 1516-1538. https://doi.org/10.1002/qj.3977.

Hólm, E. and Gong, J. (2010) Use of cloud condensate in the background error formulation. In: *Proc. of ECMWF/JCSDA Workshop on Assimilating Satellite Observations of Clouds and Precipitation into NWP Models*. Shinfield Park: ECMWF, pp. 15-17 https://www.ecmwf.int/node/9276.

Houtekamer, P.L., Lefaivre, L., Derome, J., Ritchie, H. and Mitchell, H.L. (1996) A system simulation approach to ensemble prediction. *Monthly Weather Review*, 124(6), 1225-1242. https://doi.org/10.1175/1520-0493(1996)124<1225:assate>2.0.co;2.

Isaksen, L., Bonavita, M., Buizza, R., Fisher, M., Haseler, J., Leutbecher, M. and Raynaud, L. (2010) Ensemble of data assimilations at ECMWF. In: *ECMWF Technical Memoranda, ECMWF, Reading*. https://doi.org/10.21957/OBKE4K60.

Janjić, T., Ruckstuhl, Y. and Toint, P.L. (2021) A data assimilation algorithm for predicting rain. *Quarterly Journal of the Royal Meteorological Society*, 147(736), 1949-1963. https://doi.org/10.1002/qj.4004.

Janjić, T., McLaughlin, D., Cohn, S.E. and Verlaan, M. (2014) Conservation of mass and preservation of positivity with ensemble-type Kalman filter algorithms. *Monthly Weather Review*, 142(2), 755-773. https://doi.org/10.1175/mwr-d-14-00056.1.

Johnson, A., Wang, X., Carley, J.R., Wicker, L.J. and Karstens, C. (2015) A comparison of multiscale GSI-based EnKF and 3DVar data assimilation using radar and conventional observations for midlatitude convective-scale precipitation forecasts. *Monthly Weather Review*, 143(8), 3087-3108. https://doi.org/10.1175/mwr-d-14-00345.1.

Khain, A.P., Beheng, K.D., Heymsfield, A., Korolev, A., Krichak, S.O., Levin, Z., Pinsky, M., Phillips, V., Prabhakaran, T., Teller, A., van den Heever, S.C. and Yano, J.-I. (2015) Representation of microphysical processes in cloud-resolving models: spectral (bin) microphysics versus bulk parameterization. *Reviews of Geophysics*, 53(2), 247-322. https://doi.org/10.1002/2014rg000468.

Kleist, D.T. and Ide, K. (2015) An OSSE-based evaluation of hybrid variational-ensemble data assimilation for the NCEP GFS. Part I: system description and 3d-hybrid results. *Monthly Weather Review*, 143(2), 433-451. https://doi.org/10.1175/mwr-d-13-00351.1.

Kong, R., Xue, M. and Liu, C. (2018) Development of a hybrid En3DVar data assimilation system and comparisons with 3DVar and EnKF for radar data assimilation with observing system simulation experiments. *Monthly Weather Review*, 146(1), 175-198. https://doi.org/10.1175/mwr-d-17-0164.1.

Kong, R., Xue, M., Liu, C. and Jung, Y. (2021) Comparisons of hybrid En3DVar with 3DVar and EnKF for radar data assimilation: tests with the 10 may 2010 Oklahoma tornado outbreak. *Monthly Weather Review*, 149(1), 21-40. https://doi.org/10.1175/mwr-d-20-0053.1.

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Auizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.-M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.-B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.-P., Brossier, C.L., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J.,
Perraud, E., Pianezze, J., Redelsperger, J.-L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V. and Wautelet, P. (2018) Overview of the Meso-NH model version 5.4 and its applications. Geoscientific Model Development, 11(5), 1929–1969. https://doi.org/10.5194/gmd-11-1929-2018.

Lascaux, F., Richard, E. and Pinty, J.-P. (2006) Numerical simulations of three different MAP IOPs and the associated microphysical processes. Quarterly Journal of the Royal Meteorological Society, 132(619), 1907–1926. https://doi.org/10.1256/qj.05.197.

Lawson, W.G. and Hansen, J.A. (2005) Alignment error models and ensemble-based data assimilation. Monthly Weather Review, 133(6), 1687–1709. https://doi.org/10.1175/mwr2945.1.

Lorenc, A.C. (2003) The potential of the ensemble Kalman filter for NWP – a comparison with 4D-Var. Quarterly Journal of the Royal Meteorological Society, 129(595), 3183–3203. https://doi.org/10.1256/qj.02.132.

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.-F., Marguin-aud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Béatrice, V., Vionnet, V. and Voldoire, A. (2013) The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes. Geoscientific Model Development, 6(4), 929–960. https://doi.org/10.5194/gmd-6-929-2013.

Ménétrier, B., Montmerle, T., Berre, L. and Michel, Y. (2014) Estimation and diagnosis of heterogeneous flow-dependent background-error covariances at the convective scale using either large or small ensembles. Quarterly Journal of the Royal Meteorological Society, 140(683), 2050–2061. https://doi.org/10.1002/qj.2267.

Ménétrier, B., Montmerle, T., Michel, Y. and Berre, L. (2015) Linear filtering of sample covariances for ensemble-based data assimilation. Part I: optimality criteria and application to variance filtering and covariance localization. Monthly Weather Review, 143(5), 1622–1643 https://doi.org/10.1175/mwr-d-14-00157.1.

Meng, D., Chen, Y., Wang, H., Gao, Y., Potthast, R. and Wang, Y. (2019) The evaluation of EnVar method including hydrometeors analysis variables for assimilating cloud liquid/ice water path on prediction of rainfall events. Atmospheric Research, 219, 1–12. https://doi.org/10.1016/j.atmosres.2018.12.017.

Michel, Y. (2013) Estimating deformations of random processes for correlation modelling: methodology and the one-dimensional case. Quarterly Journal of the Royal Meteorological Society, 139(672), 771–783. https://doi.org/10.1002/qj.2007.

Michel, Y. and Brousseau, P. (2021) A square-root, dual-resolution, 3DEnVar for the AROME model: formulation and evaluation on a summertime convective period. Monthly Weather Review, 149(9), 3135–3153.

Michel, Y., Auligné, T. and Montmerle, T. (2011) Heterogeneous convective-scale background error covariances with the inclusion of hydrometeor variables. Monthly Weather Review, 139(9), 2994–3015. https://doi.org/10.1175/2011mwr3632.1.

Montmerle, T., Michel, Y., Arbogast, E., Ménétrier, B. and Brousseau, P. (2018) A 3D ensemble variational data assimilation

scheme for the limited-area AROME model: formulation and preliminary results. Quarterly Journal of the Royal Meteorological Society, 144(716), 2196–2215. https://doi.org/10.1002/qj.3334.

Necker, T., Geiss, S., Weissmann, M., Ruiz, J., Miyoshi, T. and Lien, G. (2020) A convective-scale 1,000-member ensemble simulation and potential applications. Quarterly Journal of the Royal Meteorological Society, 146(728), 1423–1442. https://doi.org/10.1002/qj.3744.

Palmer, T.N., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G.J., Steinheimer, M. and Weisheimer, A. (2009) Stochastic parametrization and model uncertainty. In: ECMWF Technical Memoranda, ECMWF, Reading. https://doi.org/10.21957/PS8GBWBDV.

Purser, R.J., Wu, W.-S., Parrish, D.F. and Roberts, N.M. (2003) Numerical aspects of the application of recursive filters to Variational statistical analysis. Part I: spatially homogeneous and isotropic Gaussian Covariances. Monthly Weather Review, 131(8), 1524–1535. https://doi.org/10.1175/1520-0493(2003)131<1524:naotao>2.0.co;2.

Raynaud, L. and Bouttier, F. (2015) Comparison of initial perturbation methods for ensemble prediction at convective scale. Quarterly Journal of the Royal Meteorological Society, 142(695), 854–866. https://doi.org/10.1002/qj.2686.

Raynaud, L., Berre, L. and Desroziers, G. (2011) Accounting for model error in the Météo-France ensemble data assimilation system. Quarterly Journal of the Royal Meteorological Society, 138(662), 249–262. https://doi.org/10.1002/qj.906.

Roberts, N.M. and Lean, H.W. (2008) Scale-selective verification of rainfall accumulations from high-resolution forecasts of convective events. Monthly Weather Review, 136(1), 78–97. https://doi.org/10.1175/2007mwr2123.1.

Saunders, R., Hocking, J., Turner, E., Rayer, P., Rundle, D., Brunel, P., Vidot, J., Roquet, P., Matricardi, M., Geer, A., Bormann, N. and Lupu, C. (2018) An update on the RTTOV fast radiative transfer model (currently at version 12). Geoscientific Model Development, 11(7), 2717–2737. https://doi.org/10.5194/gmd-11-2717-2018.

Schmetz, J., Pili, P., Tjemkes, S., Just, D., Kerkmann, J., Rota, S. and Ratier, A. (2002) An introduction to Meteosat second generation (MSG). Bulletin of the American Meteorological Society, 83(7), 992. https://doi.org/10.1175/1520-0477(2002)083<0977:AITMSG>2.3.CO;2.

Seity, Y., Brousseau, P., Malardel, S., Hello, G., Bénard, P., Bouttier, F., Lac, C. and Masson, V. (2011) The arome-France convective-scale operational model. Monthly Weather Review, 139(3), 976–991. https://doi.org/10.1175/2010mwr3425.1.

Snyder, C. and Zhang, F. (2003) Assimilation of simulated Doppler radar observations with an ensemble Kalman filter*. Monthly Weather Review, 131(8), 1663–1677. https://doi.org/10.1175/2555.1.

Stein, J. and Stoop, F. (2019) Neighborhood-based contingency tables including errors compensation. Monthly Weather Review, 147(1), 329–344. https://doi.org/10.1175/mwr-d-17-0288.1.

Sun, J. and Crook, N.A. (1997) Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part I: model development and simulated data experiments. Journal of the Atmospheric Sciences, 54(12), 1642–1661. https://doi.org/10.1175/1520-0469(1997)054<1642:damrfd>2.0.co;2.
<figure><img src="image_21.png" /></figure>

Sun, J. and Crook, N.A. (1998) Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part II: retrieval experiments of an observed Florida convective storm. *Journal of the Atmospheric Sciences*, 55(5), 835–852. https://doi.org/10.1175/1520-0469(1998)055<0835:damrfd>2.0.co;2.

Sun, T., Chen, Y., Meng, D. and Chen, H. (2021) Background error covariance statistics of hydrometeor control variables based on Gaussian transform. *Advances in Atmospheric Sciences*, 38(5), 831–844. https://doi.org/10.1007/s00376-021-0271-3.

Thomas, G., Mahfouf, J.-F. and Montmerle, T. (2020) Toward a variational assimilation of polarimetric radar observations in a convective-scale numerical weather prediction (NWP) model. *Atmospheric Measurement Techniques*, 13(5), 2279–2298. https://doi.org/10.5194/amt-13-2279-2020.

Tong, M. and Xue, M. (2005) ISSN 1520-0493 Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model: OSS experiments. *Monthly Weather Review*, 133, 1789–1807. https://doi.org/10.1175/mwr2898.1.

Wang, Y. and Wang, X. (2017) 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. *Monthly Weather Review*, 145(4), 1447–1471. https://doi.org/10.1175/mwr-d-16-0231.1.

Wattrelot, E., Caumont, O. and Mahfouf, J.-F. (2014) Operational implementation of the 1D + 3D-Var assimilation method of radar reflectivity data in the AROME model. *Monthly Weather Review*, 142(5), 1852–1873. https://doi.org/10.1175/mwr-d-13-00230.1.

Xiao, Q., Kuo, Y.-H., Sun, J., Lee, W.-C., Barker, D.M. and Lim, E. (2007) An approach of radar reflectivity data assimilation and

its assessment with the inland QPF of typhoon Rusa (2002) at landfall. *Journal of Applied Meteorology and Climatology*, 46(1), 14–22. https://doi.org/10.1175/jam2439.1.

Žagar, N., Andersson, E. and Fisher, M. (2005) Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors. *Quarterly Journal of the Royal Meteorological Society*, 131(607), 987–1011. https://doi.org/10.1256/qj.04.54.

Zhang, S.Q., Zupanski, M., Hou, A.Y., Lin, X. and Cheung, S.H. (2013) Assimilation of precipitation-affected radiances in a cloud-resolving WRF ensemble data assimilation system. *Monthly Weather Review*, 141(2), 754–772. https://doi.org/10.1175/mwr-d-12-00055.1.

Zupanski, M. (2005) Maximum likelihood ensemble filter: theoretical aspects. *Monthly Weather Review*, 133(6), 1710–1726. https://doi.org/10.1175/mwr2946.1.

Zupanski, M., Navon, I.M. and Zupanski, D. (2008) The maximum likelihood ensemble filter as a non-differentiable minimization algorithm. *Quarterly Journal of the Royal Meteorological Society*, 134(633), 1039–1050. https://doi.org/10.1002/qj.251.

**How to cite this article:** Destouches, M., Montmerle, T., Michel, Y. & Caron, J.-F. (2023) Impact of hydrometeor control variables in a convective-scale 3DEnVar data assimilation scheme. *Quarterly Journal of the Royal Meteorological Society*, 149(752), 757–780. Available from: https://doi.org/10.1002/qj.4426