# A Square-Root, Dual-Resolution 3DEnVar for the AROME Model: Formulation and Evaluation on a Summertime Convective Period

YANN MICHELª AND PIERRE BROUSSEAUª

ªCNRM, Météo-France and CNRS, Toulouse, France

(Manuscript received 17 February 2021, in final form 21 May 2021)

**ABSTRACT:** A three-dimensional ensemble-variational (3DEnVar) data assimilation algorithm has been developed for the high-resolution AROME NWP system. Building on previous work on 3DEnVar for AROME, we describe a formulation of the 3DEnVar that is based on the traditional square-root preconditioning. The localization may be performed either in gridpoint or spectral space, and allows for cross-covariances between surface pressure and the three-dimensional variables. The scheme has capacity for dual resolution, with the ensemble running at a lower 3.2-km spatial resolution than the deterministic AROME running at 1.3 km. This formulation is compatible with the variational bias correction scheme used in AROME. Hybrid covariances are implemented with climatological covariances at 1.3-km resolution being combined with ensemble perturbations that are interpolated to high resolution on the fly in the computation of the gradient. Hybrid 3DEnVar has an increased computational cost compared to 3DVar, which is partly mitigated by the use of dual resolution and the adoption of a flexible convergence criterion in the minimization. To get the full benefit from the ensemble scheme, it is recommended 1) to increase ensemble size from 25 to 50 members and 2) to decrease the localization length scale for the benefit of high-density radar observations. With those changes, the 3DEnVar outperforms the operational AROME-France 3DVar by a significant margin on the first 12 h of forecast range, as evidenced by a 3-month summer experiment. Finally, a case study reports on the improved prediction of heavy rainfall that frequently occurs in the Mediterranean region.

**KEYWORDS:** Model initialization; Data assimilation; Convective storms/systems; Numerical weather prediction/forecasting; Ensembles

## 1. Introduction

Bayes's law aims at updating a prior distribution with recent information brought by observations. If we make the Gaussian approximation, the relevant probability density functions (pdf) are fully described by their means and covariances. However, computing the mean of the posterior pdf is challenging in numerical weather prediction (NWP). This is due to the high computational cost associated with large dimensions of the state and of the observation vectors, respectively, ~10⁹ and ~10⁵ for the Application of Research to Operations at Mesoscale (AROME) limited-area NWP system. Variational data assimilation schemes estimate the mode (rather than the mean) of the posterior distribution. This is done efficiently by minimizing a cost function. In NWP, the covariance matrices of the background and observations errors are never stored explicitly (Dee 1995); rather they are implemented as matrix-free operations (Bannister 2008). The solution is then found by typically an iterative Krylov method (Saad 2003), usually the preconditioned conjugate gradient (Fisher 1998; Gürol et al. 2014). The solution also gives some information about the posterior covariance through the minimization eigenpairs (Fisher and Courtier 1995).

When the dimension of the state vector is large, it is not possible to explicitly propagate the full forecast error covariances with time. The ensemble Kalman filter (e.g., Evensen 1994, 2003) essentially represents the background error covariance matrix by the sample covariance over an ensemble of

short-range forecasts of limited size. Advancing the pdf in time is achieved by evolving each member of the ensemble through a model integration, together with the addition of model perturbations, in order to simulate both the forecast evolution of initial errors and the contribution of model errors during each forecast step. The analysis step is also performed by updating each member with the observations, with two well-known variants. The first one is a stochastic interpretation of the Kalman filter equations where every member assimilates perturbed observations upon a perturbed background state issued from the perturbed forecast integration (e.g., Houtekamer and Mitchell 2001; Houtekamer et al. 2014). Alternatively, square-root schemes avoid randomization (Whitaker and Hamill 2002; Tippett et al. 2003) but are more prone to errors when the conditions deviate from Gaussianity (Lawson and Hansen 2004).

There is a possibility to blend the matrix-free approach of variational schemes with the sample covariance from the EnKF, yielding schemes referred to as “hybrid” (Hamill and Snyder 2000). The necessity to filter sampling noise in the EnKF was pointed out by Hamill et al. (2001). One efficient approach is to perform spatial localization, which is a separation-dependent moderation of the sampled covariances. The introduction of localization together with the use of preconditioning in hybrid variational schemes was described separately by Lorenc (2003) and Buehner (2005). The equivalence between those two formulations was later proven by Wang et al. (2007). Since those developments, ensemble-variational scheme have become more and more popular in NWP. They blend the numerical efficiency of variational schemes with the possibility to evolve the background error covariances in time through an ensemble analogous to the EnKF.

Corresponding author: Yann Michel, yann.michel@meteo.fr
It is not required to use an EnKF to generate an analysis
ensemble. The variational framework can be used with a
Monte Carlo approach by adding perturbations to observa-
tions and to the model (Zagar et al. 2005; Berre et al. 2006;
Raynaud et al. 2012). This results in an ensemble of data as-
simulations (EDA) that is favored by some operational centers
due to lower maintenance costs. The EDA has the capacity to
assimilate large sets of observations, being based on a varia-
tional algorithm. The EDA has been found to provide better
covariances to the deterministic assimilation than an EnKF
variant by Bowler et al. (2017), possibly because of better
statistical consistency in the description of errors (Berre et al.
2006). Also, the EDA is sometimes described as computa-
tionally expensive (e.g., Desroziers and Berre 2012), but new
algorithms that simultaneously solve for the ensemble of
minimizations are promising (Mercier et al. 2018a,b). Some
centers use the EnKF operationally and also recenter the en-
semble of analyses on an EnVar deterministic analysis.

The 3DEnVar algorithm can also be extended to 4D as
described by Liu et al. (2008), Desroziers et al. (2014). This
important breakthrough allows to estimate a four-dimensional
state given time-distributed observations, as in the asynchro-
nous EnKF (Hunt et al. 2007). Many variants of the 4DEnVar
scheme have been described in the literature (Bannister 2017).
For cloud resolving models with sophisticated physics, de-
veloping the tangent-linear and adjoint model may be very
challenging due to limited resources, making 4DEnVar an
appealing solution.

It is believed that the EnVar-EDA framework is flexible
enough to tackle long standing issues in NWP at convective scale,
such as the introduction of hydrometeors in the control variable
(Michel et al. 2011), coupled assimilation between various Earth
components, while still allowing for the simultaneous assimilation
of dense observations from satellites and radars. Developing ad-
vanced data assimilation algorithms at the mesoscale is desirable
and has been found to bring forecast improvements (Gustafsson
et al. 2018), but operational implementation has maybe lagged a
bit behind for mesoscale models compared to global models. We
first describe the AROME NWP system and its recent progress
over the last few years before stating in more detail our contri-
bution toward an EnVar for AROME.

## 2. Overview of the recent progress with the AROME NWP system

### a. The AROME model

The AROME model (Seity et al. 2011; Brousseau et al. 2016) is the operational kilometric-scale NWP model of Météo-France. The AROME model has a nonhydrostatic dynamical core (Bubnová et al. 1995) that solves the fully elastic Euler equations in terrain-following hybrid coordinates. AROME is based on a spectral representation of variables in the horizontal with Fourier series (Haugen and Machenhauer 1993). AROME incorporates detailed moist physics borrowed from the Méso-NH research model (Lafore et al. 1998; Lac et al. 2018). In particular, the three-class ice parameterization is described by Lascaux et al. (2006) and the turbulence scheme

by Cuxart et al. (2000). The surface model SURFEX<fnref n="1" /> is also
borrowed from the Méso-NH model. AROME has been used
in operations by Météo-France since 2008 (Seity et al. 2011). It
has undergone since a major change in spatial resolution to
1.3 km in the horizontal, yielding an improved prediction of
convective activity (Brousseau et al. 2016).

### b. The AROME data assimilation scheme

As explained by Seity et al. (2011), AROME-France was
developed in particular to improve the prediction of flash flood
events that frequently occur during autumn in the south of
France near the Mediterranean sea. Ducrocq et al. (2002) in-
dicate that estimating precisely the initial state in humidity was
crucial for the forecast of those precipitating events. This
requires a good data assimilation scheme with availability of
frequent and dense observations such as surface observations,
radar reflectivities, and Meteosat data.

AROME-France data assimilation scheme is a 3DVar with
hourly cycling (Brousseau et al. 2016) assimilating a wide range
of observations such as radial winds from Doppler radars
(Montmerle and Faccani 2009), radar reflectivities through
retrievals of humidity profiles (Caumont et al. 2010; Wattrelot
et al. 2014), radiometers and infrared sounders from geostationary
or polar-orbiting satellites (Guidard et al. 2011; Montmerle et al.
2007) and atmospheric path delays from ground-based the global
navigation satellite system (GNSS) (Mahfouf et al. 2015). A more
complete picture of the observation usage in the AROME 3DVar
is shown in Tables 1–2 and a recent comparison at the interna-
tional level is proposed by Gustafsson et al. (2018). In AROME,
radar data amount to more than a half of the total number of
assimilated observations.

An overview of the algorithmic formulation of the AROME
3DVar scheme is given by Fischer et al. (2005). The static error
covariances are formulated in spectral space, with the partic-
ularity that multivariate relationships are introduced between
humidity and the other variables (Berre 2000). The AROME-
France 3DVar has undergone steady progress in performance
in particular with higher spatial resolution and more frequent
data assimilation cycling (Brousseau et al. 2016) but is lacking a
key component of modern data assimilation schemes, the flow
dependency, and time evolution of background error covari-
ances (Ménétrier et al. 2014).

### c. Toward an EnVar scheme for AROME

Recently, a prototype 3DEnVar scheme was introduced for
the AROME model (Montmerle et al. 2018, M18 hereafter).
The formulation was based on a particular choice of pre-
conditioning that does not require the square root of the
background error covariance matrix B (Gürol et al. 2014). The
algorithm was then tested against a reference 3DVar for a 5-week
wintertime period, at a lower resolution of 3.8 km with 3-h as-
similation cycle compared to the operational version of AROME
(1.3 km with 1-h assimilation cycle). The sample covariances were
taken from an experimental version of the AROME EDA, which

<footnote n="1">¹http://www.umr-cnrm.fr/surfex/</footnote>
<table><caption>TABLE 1. Glossary of observations and instruments.</caption><thead><tr><th>Term</th><th>Description</th></tr></thead><tbody><tr><td>TEMP</td><td>Format for radiosonde temperature, humidity, and wind profiles</td></tr><tr><td>PILOT</td><td>Balloon wind profile observations</td></tr><tr><td>BUOY</td><td>Reporting format for buoy measurements</td></tr><tr><td>SYNOP</td><td>Format for synoptic land surface observations</td></tr><tr><td>SHIP</td><td>Format for marine (ship) observations</td></tr><tr><td>AIREP</td><td>Aircraft reports</td></tr><tr><td>ACARS</td><td>Aircraft Communications Addressing and Reporting System</td></tr><tr><td>AMDAR</td><td>Aircraft meteorological data relay</td></tr><tr><td>MESONET</td><td>Mesoscale surface station network</td></tr><tr><td>ZTD</td><td>Zenith total delay</td></tr><tr><td>SEVIRI</td><td>Spinning Enhanced Visible and Infrared Imager (satellite instrument)</td></tr><tr><td>ATOVS</td><td>Advanced TIROS Operational Vertical Sounder (system of sounding instruments on NOAA TIROS series of satellites)</td></tr><tr><td>IASI</td><td>Infrared Atmospheric Sounding Interferometer (satellite instrument)</td></tr><tr><td>AIRS</td><td>Atmospheric Infrared Sounder (satellite instrument)</td></tr><tr><td>CRIS</td><td>Cross-track Infrared Sounder (satellite instrument)</td></tr><tr><td>ATMS</td><td>Advanced Technology Microwave Sounder (satellite instrument)</td></tr><tr><td>SSMIS</td><td>Special Sensor Microwave Imager/Sounder (satellite instrument)</td></tr><tr><td>MHS</td><td>Microwave Humidity Sounder (satellite instrument)</td></tr><tr><td>GMI</td><td>Global Microwave Imager (satellite instrument)</td></tr><tr><td>AMV</td><td>Atmospheric motion vector</td></tr></tbody></table>

since then has undergone several scientific changes and finally went into operations in July 2018. Since then, Caron et al. (2019) have highlighted the benefit to use larger ensemble sizes than the 25 members of the current AROME EDA.

The goal of this paper is to update on some results of M18 and of Caron et al. (2019) by testing 3DEnVar at higher resolution of 1.3 km, with hourly cycling, with a more recent version of the EDA, and with all observations and assimilation components active. The more recent version of the EDA has 3.2-km spatial resolution compared to 3.8 km in M18. It has a larger ensemble size, 50 members versus 25 in M18. Finally the representation of model error is also different: in this paper we use stochastic perturbations of the physics tendencies (Bouttier et al. 2012), whereas M18 use multiplicative inflation. To achieve a fair comparison between AROME 3DVar and

3DEnVar scheme, we found practical to develop a standard preconditioning based on the square root of **B** (Fisher 1998). The online variational bias correction of (Auligné et al. 2007) is also active, contrary to M18. Also, the experiments cover an extended summer period of three months that is more representative of the high precipitating events that AROME is targeted at.

The choice of a square-root formulation is possibly not final, other algorithms such as the observation-space algorithm from Gratton and Tshimanga (2009) being possibly better candidates. In the linear framework, those algorithms result in mathematically equivalent results at every iteration of the minimization (Gürol et al. 2014). The square-root implementation requires an extension of the control variable (Lorenc 2003; Buehner 2005). With a square-root implementation, it is straightforward to implement a full-**B** preconditioning with **B** = **B**¹/²**B**T/². It may also be easier to check the correctness of the code with an adjoint test. The square-root form also allows to make random draws of the background pdf. It was also favored by ECMWF for its global 4DVar, which shares many components with AROME. Altogether, the square-root implementation of EnVar is believed to be an important added feature to the AROME NWP system.

### d. Outline

This paper is structured as follows: 3 describes the newly developed square-root formulation of 3DEnVar for AROME. Several points are emphasized, first the localization is applied to a joint 2D and 3D control variable, allowing to keep the cross-covariance structure between surface pressure and the other variables. Second, both the gridpoint and spectral localization schemes described by M18 are possible. The interpolations from the coarse resolution ensemble to the high resolution gradient are performed on the fly. Section 4 describes sensitivity tests to some aspects of the scheme, in particular to the localization and to the ensemble size. Section 5 presents an extended evaluation against a version of AROME-France 3DVar that is very close to the one operational in 2019. Section 6 reports on a case study. Conclusions and perspectives are then drawn, in particular toward 4DEnVar and the treatment of nonlinearity.

## 3. A square-root formulation of 3DEnVar for AROME

### a. General formulation

In a Bayesian framework, we solve an assimilation problem to find the mode of the posterior pdf conditioned on the

<table><caption>TABLE 2. Overview of assimilated observation types in the AROME France deterministic 3DVar. See Table 1 for the description of observation types. The percentages are among a total of 3.34 × 10⁷ observations in March 2021.</caption><thead><tr><th rowspan="2">Details</th><th colspan="4">Observation types</th></tr><tr><th>Conventional</th><th>Radar</th><th>GNSS</th><th>Satellites</th></tr></thead><tbody><tr><td></td><td>TEMP, PILOT, BUOY,SYNOP, SHIP, AIREP,ACARS, AMDAR, MESONET</td><td>Doppler winds, reflectivities</td><td>Ground-based ZTD</td><td>SEVIRI, ATOVS, IASI, AIRS,CRIS, ATMS, SSM/IS, MHS, GMI,AMV, Scatterometer winds, ...</td></tr><tr><td>Avg percentage in the assimilation</td><td>30.5</td><td>54.8</td><td>2.1</td><td>12.6</td></tr></tbody></table>
observations (Lorenc 1986). The solution, defined as the analysis $\mathbf{x}_o$, is found by minimizing a cost function. In a 3D framework, e.g., when observations are given at a given time, this cost function may be written as follows:

$$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}_o - \mathbf{H}(\mathbf{x})]^T \mathbf{R}^{-1} [\mathbf{y}_o - \mathbf{H}(\mathbf{x})], \quad (1)$$

where $\mathbf{x}$ is the state we estimate, $\mathbf{x}_b$ the background, $\mathbf{y}_o$ the observation vector, $\mathbf{H}$ the nonlinear observation operator, and $\mathbf{B}$ and $\mathbf{R}$ are the error covariance matrices of the background and observation errors. 3DVar and 3DEnVar scheme share this general description, the difference being in the way $\mathbf{B}$ is represented in (1). We do not include the variational bias correction in this equation (Dee 2004) for simplicity.

The innovation vector is introduced as the departure $\mathbf{d} = \mathbf{y}_o - \mathbf{H}(\mathbf{x}_b)$. The incremental formulation (Courtier et al. 1994) solves (1) with respect to the increment $\mathbf{w} = \mathbf{x} - \mathbf{x}_b$ by a Gauss-Newton algorithm, which results in successive linearizations of $\mathbf{H}$ against more and more accurate guesses. In AROME 3DVar so far, only one linearization $\mathbf{H}$ is performed against the background. The cost function to be minimized follows:

$$J(\mathbf{w}) = \mathbf{w}^T \mathbf{B}^{-1} \mathbf{w} + \frac{1}{2} (\mathbf{d} - \mathbf{H}\mathbf{w})^T \mathbf{R}^{-1} (\mathbf{d} - \mathbf{H}\mathbf{w}). \quad (2)$$

This cost function is quadratic. Equating the gradient to zero yields the linear system to be solved:

$$(\mathbf{B}^{-1} + \mathbf{H}^T \mathbf{R}^{-1} \mathbf{H}) \mathbf{w} = \mathbf{H}^T \mathbf{R}^{-1} \mathbf{d}. \quad (3)$$

Equation (3) involves the resolution of a linear system in model space and is referred to as the primal problem. It needs to be preconditioned for efficiency, which can be achieved in split or full forms (Saad 2003). The split, or square-root, form is employed in this paper and consists in the explicit change of control variable:

$$\mathbf{w} = \mathbf{U}\mathbf{\chi}, \quad (4)$$

where $\mathbf{U}$ is a “left square root”<fnref n="2" /> of $\mathbf{B}$, e.g., a possibly rectangular matrix verifying $\mathbf{B} = \mathbf{U}\mathbf{U}^T$. The primal problem becomes then:

$$(\mathbf{I} + \mathbf{U}^T \mathbf{H}^T \mathbf{R}^{-1} \mathbf{H} \mathbf{U}) \mathbf{\chi} = \mathbf{U}^T \mathbf{H}^T \mathbf{R}^{-1} \mathbf{d}, \quad (5)$$

which can be solved with the conjugate gradient algorithm equipped with the canonical inner product (Fisher 1998). M18 employ rather a full (non-square-root) form of preconditioning (Derber and Rosati 1989; Gürol et al. 2014) that may be, in a three-dimensional context, cheaper as the scalar products are done with shorter vectors that lie in model space. Alternatively, it may be even more advantageous to solve rather for the dual problem with the restricted preconditioned conjugate gradient

as the dimension of the observation vector is much lower than the dimension of the state or control variable spaces (Gratton and Tshimanga 2009). The size of the control variable also has implication on the memory footprint, especially if full re-orthogonalization of the conjugate gradient is employed.

All formulations are compatible with the variational bias correction, but note that it would be incorrect to refer to the “full” approach as “B”-preconditioned (Gürol et al. 2014), as the preconditioning of VarBC includes a Hessian term (Dee 2004). Those formulations can also give the same iterates in exact arithmetic for nonlinear problems (Gratton et al. 2013) but this excludes so far to our knowledge other important developments such as the second-level preconditioning (Tshimanga et al. 2008; Gratton et al. 2011), the multi-incremental strategy (Veverse and Thepaut 1998) and the quasi-static approach (Pires et al. 1996; Fillion et al. 2018).

The square-root form lies in the definition of $\mathbf{U}$. The 3DEnVar developed in this paper follows a hybrid formulation (Hamill and Snyder 2000) with $\mathbf{B}$ defined as the linear combination:

$$\mathbf{B} = \beta_s^2 \mathbf{B}_s + \beta_e^2 \mathbf{B}_e, \quad (6)$$

where

* $\mathbf{B}_s$ is the static background error covariance, whose general formulation follows Berre (2000), and the application to the AROME operational 3DVar is described by Brousseau et al. (2011);

* $\mathbf{B}_e$ is the ensemble error covariance, described for instance by Hamill et al. (2001);

* $\beta_s^2$ and $\beta_e^2$ are scalar weights.

The ensemble error covariance is the localized sample covariance:

$$\mathbf{B}_e = \mathbf{C} \circ \mathbf{X} \mathbf{X}^T, \quad (7)$$

where

* $\mathbf{C}$ is the localization, a correlation matrix with positive coefficients;

* $\circ$ is the Hadamard (element-by-element) product;

* $\mathbf{X}$ is the perturbation matrix where every column contains the deviation of an ensemble member with respect to the ensemble mean, normalized by a factor of $1/\sqrt{n_e - 1}$ where $n_e$ is the ensemble size.

Localization is a standard method to reduce the sampling noise in the estimate of the covariances (Ménétrier et al. 2015a). The sample covariances are moderated, at long separation distance, by the coefficients within $\mathbf{C}$ (Hamill et al. 2001). The estimation of the localization is a long-standing issue in NWP, but some recent theoretical and practical progress is provided by Ménétrier et al. (2015b).

Then, the change of variable in Eq. (4) can be written as follows:

$$\mathbf{w} = \beta_s \mathbf{U}_s \mathbf{x}_s + \beta_e \mathbf{U}_e \mathbf{x}_e, \quad (8)$$

where $\mathbf{U}_s$ and $\mathbf{U}_e$ are the respective square roots of $\mathbf{B}_s$ and $\mathbf{U}_e$ and the control variable $\mathbf{x} = [\mathbf{\chi}_s^T, \mathbf{\chi}_e^T]^T$ has been extended.

<footnote n="2">²This abuse of notation is frequent in the data assimilation community and should not be confused with the square root of a matrix, which would be written as $\mathbf{B} = \mathbf{U}^2$.</footnote>
The expression for the square root of the static covariance is depending on the implementation, with a review given by Bannister (2008). The expression for the square root of the (localized) ensemble covariance applied to the ensemble part of a control vector $\chi_e$ is given by Buehner (2005) and Wang et al. (2007):

$$ \mathbf{U}_e \chi_e = \sum_{e=1}^{n_e} X_e \circ [\mathbf{L} \chi_e], \qquad (9) $$

where $\mathbf{L}$ is a left square root of the localization $\mathbf{C} = \mathbf{LL}^T$, the sum is over the ensemble index and $X_e$ is the $e$ column of $\mathbf{X}$, that is the normalized deviation of member $e$. Note that the element-by-element product between matrices in Eq. (7) has been replaced by $n_e$ applications of the square root of the localization in Eq. (9), making it possible in large dimensions encountered in NWP.

### b. Factorization of the localization among the variables

When the same localization is applied to the different control variables, it can be factored and applied only once, thus reducing the computational cost (Clayton et al. 2013; Desroziers et al. 2014). In AROME, the static control variables are the 3D relative vorticity $\zeta$, divergence $\eta$, specific humidity $q$, temperature $T$ and the 2D logarithm of the surface pressure $\ln P_s$. To apply a single localization, we need to group 3D and 2D variables and also we need the spatial scales of the variables to be similar. In the static covariance term, surface pressure was grouped together with temperature in a common variable ($T$, $\ln P_s$) to represent the hydrostatic balance (Berre 2000; Derber and Bouttier 1999; Parrish and Derber 1992). A similar approach is adopted here in the sense that the localization is applied to a joint control variable gathering the sets of 2D and 3D variables. However, we control the zonal and meridional components of the horizontal wind ($u, v$) rather than $(\zeta, \eta)$ to ensure that the spatial scales are more similar.

To write this formally and explicit the form of the control variable, we need to introduce two projection operators that select, respectively, the 3D and 2D parts in a control variable:

$$ \mathbf{\Pi}_3 = \begin{pmatrix} 1 & 0 \\ \vdots & \vdots \\ 0 & 1 \end{pmatrix}, \qquad (10) $$

$$ \mathbf{\Pi}_2 = \begin{pmatrix} 1 & 0 \\ \mathbf{0} & \vdots \\ 0 & 1 \end{pmatrix}, \qquad (11) $$

where $\mathbf{\Pi}_3$ and $\mathbf{\Pi}_2$ are rectangular matrices and $\mathbf{0}$ are rectangular matrices containing 0. Those matrices apply to a joint 3D + 2D control variable $\chi_e$ having $N_t + 1$ vertical levels (with $N_t$ being the number of vertical levels of any of the four 3D variables of type $u, v, T, q$). The matrix formulation of the localization may then be written as follows:

$$ \mathbf{L} = \begin{pmatrix} \mathbf{\Pi}_3 \\ \vdots \\ \mathbf{\Pi}_3 \\ \mathbf{\Pi}_2 \end{pmatrix}, \qquad (12) $$

where there is a $\mathbf{\Pi}_3$ block for every 3D variable ($u, v, T, q$).

### c. Cross localization between the variables

To reveal the cross localization between variables, we may expand the full form of $\mathbf{C} = \mathbf{LL}^T$ (not shown). It then becomes apparent that the implied localization for a 3D variable is the block:

$$ \mathbf{\Pi}_3 \mathbf{L}_s \mathbf{L}_s^T \mathbf{\Pi}_3^T, $$

that the implied localization for a 2D variable is

$$ \mathbf{\Pi}_2 \mathbf{L}_s \mathbf{L}_s^T \mathbf{\Pi}_2^T, $$

and finally that the cross localization between the 2D variable and any 3D variable is

$$ \mathbf{\Pi}_2 \mathbf{L}_s \mathbf{L}_s^T \mathbf{\Pi}_3^T. $$

Thus, the cross covariances between surface pressure and the other variables are retained up a certain amount, which depends on the considered localization approach. In the case where the horizontal localization is specified to be the same for surface pressure and for other variables (as considered here, albeit being not optimal, since surface pressure is much larger scale than e.g., low-level temperature and humidity), the localization of these cross covariances with surface pressure will depend only on the vertical localization specified in $\mathbf{L}_s$. In practice, the vertical localization is based on a separation distance using the average pressure of each vertical level (M18). The corresponding average pressure for the 2D variable is specified to be mean sea level pressure, thus the ensemble covariances are kept between surface pressure and the boundary layer.

### d. Spatial localization

In M18, the control variable is a field increment in gridpoint space. Two alternative spatial localizations were then tested, the first one being based on a spectral convolution (denoted $\mathbf{S}$) and the second one (denoted $\mathbf{G}$) being in gridpoint space that uses recursive filters (Purser et al. 2003) with a spatial deformation on the vertical (Michel 2013). In the square-root formulation, the control variable is based on a spectral representation. However, the Schur product with ensemble perturbations is performed in gridpoint space, thus requiring an inverse Fourier transform. The square root of the spatial localization can thus be formulated either in gridpoint or in spectral spaces, simply by commuting the spatial localization and this (vertical level per vertical level) inverse Fourier transform $\mathbf{F}^{-1}$. Note that the two components of the wind are localized rather than vorticity and divergence as stated above, thus this inverse Fourier transform also incorporates the corresponding change of variable. There is no increment of the mean wind over the domain, as for 3DVar (Fischer et al. 2005). The spectral formulation of spatial localization may then be written as follows:

$$ \mathbf{L}_s = \mathbf{F}^{-1} \mathbf{S}, \qquad (13) $$

whereas the gridpoint formulation of spatial localization may be written as follows:
$$L_s = \mathbf{GF}^{-1}. \qquad (14)$$

M18 indicate that the gridpoint formulation of spatial localization was performing better, a finding that was reproduced here with our set of experiments (not shown). Thus only the formulation (14) is used in the rest of this paper. A spectral formulation of the control variable is kept even when using the gridpoint formulation of spatial localization mostly to preserve the code compatibility between AROME and ARPEGE-IFS formulations, and possibly in the future ease the implementation of multi-incremental schemes (Trémolet 2007).

### e. Dual resolution

To reduce the computational cost, it is common to implement a dual resolution strategy where the ensemble is run at coarser resolution than the deterministic (hybrid) EnVar (Clayton et al. 2013; Kleist and Ide 2015). In this paper, the implementation of dual resolution for the square-root scheme is done on the fly and in the ensemble part only of Eq. (8):

$$\mathbf{w} = \beta_s \mathbf{U}_s \mathbf{\chi}_s + \beta_e \mathbf{I}_{12h} \mathbf{U}_e \mathbf{\chi}_e, \qquad (15)$$

where $\mathbf{I}_{12h}$ is the interpolation from low to high resolution, using bilinear four-point interpolation. The ensemble part of the control vector $\mathbf{\chi}_e$ is defined at lower resolution, mitigating both memory and computational costs. Intensive operators such as the localization and the Schur product are performed at low resolution. Note, however, that the static covariance term is untouched and applies to the high-resolution control vector $\mathbf{\chi}_s$, as described below.

### f. Summary of the proposed implementation

To summarize, the proposed implementation solves for Eq. (3) with the change of variable (15). This change of variable describes the square root of the hybrid background error covariance. The static part is based on the formulation of Berre (2000). The ensemble part is summarized by Eqs. (9)–(14), with a common localization applied to all variables with both spectral-based and gridpoint-based implementation possible. Cross covariances between 2D and 3D variables are retained. We now focus on the convergence of the minimization of the scheme.

### g. A flexible convergence criterion

Compared to their 4D variants, 3DVar and 3DEnVar schemes are much cheaper computationally. Still, computing time is crucial for certain applications such as short-range prediction for nowcasting (Auger et al. 2015). It is thus useful to evaluate the convergence of the minimization in order to perform only the necessary iterations and to save computational time.

The conditioning of the linear problem (5) is affected by the observation distribution and accuracy and by the length scales in the background error covariance matrix (Haben et al. 2011). For AROME, the observation density is varying a lot in time (Brousseau et al. 2014), with for instance more numerous radar observations in rainy conditions (Montmerle and Faccani 2009; Wattrelot et al. 2014) and more aircraft measurements at daytime compared to nighttime. The standard implementation of AROME-France 3DVar uses so far a fixed number of

<figure><img src="image_1.png" /><figcaption>FIG. 1. (top) Total number of assimilated observations, (middle) the conditioning of the Hessian of the cost function, and (bottom) the number of performed iterations in the minimization (limited between 25 and 50) in the AROME hybrid 3DEnVar over a 3-month period.</figcaption></figure>

iterations set to 50, because its computational cost is very affordable. Following Mercier et al. (2018b), it may be advantageous to use a convergence criterion for the minimization that allows the number of iterations to vary.

The total number of observations assimilated by an hourly cycling AROME 3DEnVar experiment is shown in Fig. 1 (top)
over a 3-month period. This number varies by a factor of 10, with a maximum of about $8 \times 10^4$. The conditioning of the Hessian is defined the ratio of the largest to the smallest eigenvalue of the matrix on the right hand side of the problem (5). It is generally a good indicator of the number of iterations needed to convergence, even if this more generally depends on the full eigenvalue spectrum and also on the initial residual. The conditioning, shown in Fig. 1 (middle), has a maximum of about $5 \times 10^3$ and exhibits significant correlation with the observation number. This favors a flexible convergence criterion to save computational time. Following Trémolet (2007) and Mercier et al. (2018b), the convergence is reached when the background term change from one iteration to the other is below 1%. The number of iterations performed by the minimization is also bounded between 25 and 50. As shown in Fig. 1 (bottom), the number of performed iterations in our implementation of 3DEnVar varies accordingly, with more iterations usually needed in the case of worse conditioning. Yet, 25 iterations are used in 65% of the situations, bringing an improvement over a fixed number of 50 iterations. Note that a similar computational saving could have been realized in the baseline 3DVar, but this has not been tried.

### h. Computational cost

For 3DEnVar, we expect the main cost to be in the successive applications of the localization. From a computational point of view, the recursive filters are used on stripes of data, which then undergo a global transposition, filtering in the other direction followed by and inverse transposition (Fujita and Purser 2001). The only difference with the scheme implemented by M18 is the presence of shared memory parallelization within the stripes achieved by OpenMP directives (Dagum and Menon 1998).

The number of iterations has been shown to vary between 25 and 50 (Fig. 1), and the cost in the worst case scenario is now examined. On 7 July 2018 the number of dense radar observations happens to be very large causing a degradation of the conditioning and of the number of iterations needed for convergence (50).

3DEnVar is more memory-consuming mostly because of the storage of the ensemble members. As shown by Table 3, it is possible to run 25- and 50-member hybrid 3DEnVar on twice the number of nodes compared to 3DVar (e.g., 36 rather than 18). The same number of MPI-tasks is used for 3DEnVar than for 3DVar, thus the number of OpenMP threads is doubled with 3DEnVar (from 5 to 10). A 25-member hybrid 3DEnVar has an increase in cost by a factor of 1.6 over 3DVar; and moving to an ensemble size of 50 members requires an additional multiplicative factor of 1.5. Thus, the cost grows almost linearly with ensemble size. This increase in cost is partly mitigated by a reduction in the number of iterations compared to operational settings as described above. Yet running hybrid 3DEnVar with 50 members doubles the wall-time in the worst case scenario when the number of iterations cannot be reduced because of bad conditioning of the Hessian.

## 4. Description and sensitivity experiments

One important aspect of this study is the use of dual resolution, with the deterministic (ensemble) running at 1.3 (3.2)-km spatial resolution. This section first describes the deterministic AROME-France 3DVar and the AROME-France EDA. The

<table><caption>TABLE 3. Wall-time of the minimization of various assimilation schemes for the AROME model run with 144 MPI tasks and 5 or 10 Open-MP threads with 50 iterations in the minimization (on the “prolix” Bullx DLC B710 platform from Météo-France equipped with Intel Xeon E5-2697v2 12C 2.7GHz processors).</caption><thead><tr><td>Scheme</td><td>Nodes</td><td>Threads</td><td>Wall-time</td></tr></thead><tbody><tr><td>3DVar</td><td>18</td><td>5</td><td>263 s</td></tr><tr><td>25-member hybrid 3DEnVar</td><td>36</td><td>10</td><td>432 s</td></tr><tr><td>50-member hybrid 3DEnVar</td><td>36</td><td>10</td><td>644 s</td></tr></tbody></table>

sensitivity of the performance of hybrid 3DEnVar to its parameters is then evaluated. Full evaluation is postponed to section 5.

### a. Description of the AROME-France EDA

The EDA is an ensemble of assimilation schemes with perturbed inputs. The EDA mimics the approach taken in the stochastic ensemble Kalman filter (Evensen 2003) by explicitly perturbing the observations, the model and in the LAM case the boundary conditions. In the case of AROME, the perturbed boundary conditions are taken from the global EDA (Raynaud et al. 2011). The AROME EDA consists in running perturbed 3DVars in parallel. It provides an ensemble of analyses and forecasts that is used to approximate the statistics of analysis and background errors (Fisher 2003; Zagar et al. 2005; Belo-Pereira and Berre 2006). There is an inconsistency between the algorithms used in the deterministic (hybrid 3DEnVar) and in the ensemble (3DVar) assimilations, and reducing this mismatch will be the subject of future research.

Over France, it is possible to run the AROME-France EDA near convection-allowing resolution. Schwartz (2016) also investigate a dual resolution analysis scheme where the deterministic background has convection-allowing resolution (4 km), but their ensemble resolution (20 km) still requires convective parameterization because of the large domain covered. For AROME, the spatial resolution of early experimental versions of the EDA was set 4 km (Raynaud and Bouttier 2015) or 3.8 km (M18). For operational implementation in July 2018, the spatial resolution was finally set to 3.2 km (yielding moderate improvements to the scores of the AROME ensemble prediction system, not shown).

The operational EDA version has an online multiplicative inflation scheme to represent the average effect of model error (Raynaud and Bouttier 2015). The inflation factor is derived from the skill over spread ratio. In our experiments, a different version of the parameterization of model errors is used where the inflation is replaced by stochastic perturbations of physics tendencies (SPPT) following Bouttier et al. (2012). The ensemble size has also been increased to 50 members, but we will report first on experiments with a subensemble of 25 members.

### b. Data assimilation workflow

The coupling between the two systems, the deterministic EnVar and the EDA, is shown in Fig. 2. The high resolution deterministic EnVar has a hourly data assimilation cycle, while the EDA is updated only every three hours. The EDA at main time HH provides 4-h-range forecasts that are used in the
<figure><img src="image_2.png" /><figcaption>FIG. 2. Data assimilation workflow of (top) the AROME deterministic hybrid 3DEnVar and (bottom) the AROME EDA showing how the ensemble information is integrated. For the deterministic 3DEnVar, $\mathbf{x}_b$ are the backgrounds, $\mathbf{y}_o$ are the observation vectors, and $\mathbf{x}_a$ are the analysis along times HH, HH+1 (h), etc. For every member (e) of the EDA, $\mathbf{x}_{b+r}^e$ are the backgrounds at forecast range r, $\mathbf{y}_o^e$ are the (perturbed) observation vectors, $\mathbf{x}_a^e$ are the analysis along times HH, HH+3, etc. The static $\mathbf{B}_s$ and ensemble $\mathbf{B}_a$ background error covariances are defined in Eqs. (6)–(7). The ensemble short-range forecast at ranges 2, 3, and 4 h are used to compute the ensemble background error covariances in the hourly cycling deterministic hybrid 3DEnVar (which is only illustrated at HH+2 in the figure for simplicity).</figcaption></figure>

ensemble background error covariance of the EnVar at times HH+2, HH+3, HH+4 (only the case HH+2 is illustrated in Fig. 2 for simplicity). This is because in real time the 1-h forecasts of the EDA would arrive too late for the corresponding AROME EnVar analysis according to operational time constraints. The EDA is not recentered at all on the deterministic EnVar, which means that the interaction between the EnVar and the EDA is only one way and done through the background error covariances (as illustrated by the dashed lines in Fig. 2).

### c. Sensitivity experiments

The performance of the AROME hybrid 3DEnVar may depend on many parameters such as the vertical and horizontal localization, the scalar weights of hybridization, the ensemble size. With an expensive state-of-the-art numerical model such as AROME at 1.3-km resolution, it is computationally

demanding to optimize those parameters over long experiments. The sensitivity of the scheme to this parameters has been therefore evaluated over a first set of experiments limited to 15 days. A preliminary experiment with hybrid 3DEnVar experiment used the settings of M18: the *e*-folding localization length scales of the horizontal and vertical localizations are set, respectively, to 170 km and 0.3 in units of log(*P*) (Buehner 2005). The ensemble size was set to 25 members and the weights are set to $\beta_e^2 = 0.8$ (vs $\beta_s^2 = 0.2$) for the ensemble (versus static) covariances. The experiments in this section compare hybrid 3DEnVars compared to this baseline version; the comparison with 3DVar is postponed to section 4.

## 1) SENSITIVITY TO LOCALIZATION

The localization is a crucial parameter in ensemble data assimilation. Objective estimation procedures generally assume
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Relative changes in RMS of background departures for (top) retrievals of relative humidity and (bottom) radial winds from Doppler radars, for a 25-member hybrid 3DEnVar with various horizontal localizations (10–120 km) compared to the reference (170 km), over a 15-day period (14–29 May 2018), stratified along elevation (in °).</figcaption></figure>

that localization is mostly due to sampling error (Anderson 2012; Ménétrier et al. 2015b). Those procedures yield a useful estimates that may need to be refined in some cases by tuning. Figure 3 presents the RMSEs of background departures for radar observations, for a cycling 25-member hybrid 3DEnVar with various horizontal localizations (10–120 km). It is clear that reducing the localization length scale has a positive impact on this measure, reducing the error by around 10% near the radar. Radar observations are assimilated at a quite high spatial density in AROME (about 8 km), which is not uncommon for hourly cycling, kilometer-scale models (Simonin et al. 2014). Those dense observations may benefit more from a reduced localization than other observation types.

Figure 4 shows the RMSEs of background departures for GNSS atmospheric path delays. The optimal setting for the

<figure><img src="image_5.png" /><figcaption>FIG. 4. RMS of background departures for GNSS delays, for a 25-member hybrid 3DEnVar, as a function of the horizontal localization and over a 15-day period (14–29 May 2018).</figcaption></figure>

horizontal localization according to this metric is 40 km. Shorter or longer localizations produce degradations by about 4%. A similar signal has been observed for the RMSEs of background departures for AIREP winds and temperature (not shown), where shorter localizations degrades the performance in the boundary layer.

The value of 40 km for the optimal localization is inferior by a factor of ~4 to the value reported by M18. Possible explanations are the facts that our experiments cover a summertime, convective period of the year rather than wintertime cases, and that our experiments are run with hourly cycling and 1.3-km spatial resolution. Another major difference lies in the fact that we use a hybrid scheme where the larger scales can be represented by the static covariances. Also finally, the square-root formulation is different from the one of M18, thus the specified length scales are not directly comparable. To be more precise, the localization of M18 is for the gridpoint formulation $GG^T$ whereas the spectral control variable requires the localization to be $GF^{-1}F^{-T}G^T$ which is different because the inverse Fourier transform $F$ incorporates an elliptic truncation. Removing the spectral transform would require a gridpoint control variable, which is envisioned in future developments. Note also that a similar sensitivity exists with respect to the vertical localization length scale (not shown).

The tuned localization has been found to depend on the observation type, and especially needs to be shorter for dense observations. Ying et al. (2018) argue that a denser network provides more independent pieces of information, thus a smaller localization radius is beneficial because it allows the same number of observations to constrain each state variable. For AROME, this result may favor the use of scale-dependent localization (Caron et al. 2019). Taking into account both the results shown in Figs. 3 and 4, we choose the horizontal localization to be 30 km. This localization is compared to the averaged correlation of temperature errors deduced from the climatological background error matrix $B_s$ in Fig. 5. The localization is a bit sharper than the average correlation, thus ensemble increments will be of smaller spatial scales. The tuned value of the vertical localization is set to about 0.15 hPa for a 25-member ensemble.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Horizontal localization (solid line) and climatological correlation (dashed line) for temperature at model level 70 (~920 hPa) as a function of separation distance (km).</figcaption></figure>

## 2) SENSITIVITY TO ENSEMBLE SIZE

The performance of ensemble data assimilation is known to be quite strongly dependent on the ensemble size, up to about 100 members (Houtekamer and Zhang 2016). The ensemble size is generally bounded by computational resources. For AROME, there is another limitation that lies in the current implementation of the LAM EDA: the boundary conditions are taken from the global EDA that currently has 50 members. Therefore the ensemble size will be taken in the range 25–50 members. Caron et al. (2019) report significant improvement of the performances of their 3DEnVar with AROME when increasing the ensemble size with a lagged approach. It is worth reinvestigating this question with a standard (nonlagged) EDA and in our context of higher spatial resolution, hourly cycling, and summertime experiment.

Figure 6 illustrates the impact of increasing ensemble size at constant horizontal localization (40 km). The relative changes in RMS of background departures reaches about 2% (4%) with 35 (50) members, compared to 25 members, for AIREP wind data. The improvement is similar for AIREP temperature data and radar retrievals of relative humidity and radial winds, but for those observations the impact is halved (about 1% and 2% with ensemble sizes of 35 and 50 members). An improvement is also visible on background departures for GNSS delays but is more limited (0.5% and 1.2% for ensemble sizes of 35 and 50 members). In summary, improving ensemble size yields an improved fit of the background to all observations, but the magnitude depends on the observation type, with more improvement for wind data compared to temperature and humidity.

## 3) SENSITIVITY TO HYBRID WEIGHTS

There is an important parameter to be tuned in the EnVar, which is the relative weight given to the ensemble covariances ($\beta_e^2$) with respect to the static ones ($\beta_s^2$). Three settings have been considered: the equal weight setting ($\beta_e^2 = \beta_s^2 = 0.5$), more weight on the ensemble ($\beta_e^2 = 0.8, \beta_s^2 = 0.2$) or more weight on the static ($\beta_e^2 = 0.2, \beta_s^2 = 0.8$). It may also be necessary, as in Clayton et al. (2013), to have ensemble and climatological

<figure><img src="image_7.png" /><figcaption>FIG. 6. Relative changes in RMS of background departures for wind observations from aircrafts, for 35- and 50-member hybrid 3DEnVar compared to the reference with 25 members, over a 15-day period (14–29 May 2018).</figcaption></figure>

weights where the squared sum does not equal to 1.0 or equivalently to use some global scaling of the hybrid covariances. Here, we chose the latter approach with $\beta_e = \beta\tilde{\beta}_e$ and $\beta_s = \beta\tilde{\beta}_s$, with optimal value being set from again optimization of background departures with $\beta = 0.6$.

In terms of innovation departures, the signal has been found to depend on the observation considered. As shown in Fig. 7, the equal weight setting gives the lower temperature RMSE compared to unequal weights. However, when considering wind observations, the best setting puts more weight on the ensemble. When considering humidity departures from GNSS delays (not shown), the equal weight setting and the setting with more weight on the ensemble give a similar RMSE, which is lower by 2% compared to the setting with more weight on the static. Occasional failure of the forecasts due to numerical instability also occurred when using a larger weight on the ensemble. Overall, the equal weight setting has been chosen but we may expect that more weight may be put on the ensemble part in the future when the ensemble will be improved in terms of spatial resolution and possibly size.

We turn now to a more complete evaluation of 3DEnVar versus 3D-Var.

# 5. Evaluation against 3DVar

## a. Settings and evaluation period

The choice of the parameters for the hybrid, square-root-based 3DEnVar version of AROME-France is based on the finding of the previous section: the static and ensemble weights are set to 0.5. The ensemble size is 50 members, and the localization is set to 0.25 in the vertical in log-pressure coordinates and 40 km in the horizontal. For the minimization, the flexible convergence criterion is used. The reference is an AROME-France 3DVar suite, which is very similar to the (2019) operational AROME-France 3DVar suite.

The period of evaluation ranges from 12 May 2018 to 14 August 2018. In May, the synoptic situation was dominated
<figure><img src="image_8.png" /><figcaption>FIG. 7. Relative changes in RMS of background departures for (top) temperature and (bottom) wind observations from aircrafts, for hybrid weights $\tilde{\beta}_e^2 = 0.2$, $\tilde{\beta}_e^2 = 0.8$ (dashed line) and $\tilde{\beta}_e^2 = 0.8$, $\tilde{\beta}_s^2 = 0.2$ (solid line) against the equal weight setting $\tilde{\beta}_e^2 = \tilde{\beta}_s^2 = 0.5$, over a 15-day period (14–29 May 2018).</figcaption></figure>

by successive blockings over Northern Europe and low pressures over Spain and the Mediterranean basin. This situation has favored the daily development of storms during several weeks. The weather turned then anticyclonic at the mid-June. In July and August, warmer conditions than the climatology have prevailed and yielded several storms in particular over mountainous areas. This period is thus a good test for AROME, and could be complemented by a period of evaluation in autumn covering Cevenol events (e.g., Ducrocq et al. 2014; Jacques 2016).

b. RMS of analysis and background departures

Figure 8 shows the relative differences between the fit of the analyses and of the background (1-h-range forecast) to various observation types measured by different instruments, between the 3DEnVar experiment and the 3DVar reference. In terms of analysis, there is a general tendency of the 3DEnVar experiment to fit most of the observations more closely in the free troposphere. On the contrary, the fit is looser near the surface

and in the stratosphere, and varies for the satellite observation that represent some average over the vertical (not shown). Background departures are generally improved by a large margin, for instance about 5% for wind observations and 3% for temperature and humidity.

c. Evaluation of longer forecasts with standard metrics

One difficulty when evaluating the performance of high resolution NWP forecast systems lies in the lack of frequent and dense observations over the limited area domain. Radiosondes that are used to evaluate global models are not very adequate for this purpose. The evaluation of forecasts is proposed here for an extended set of mesoscale observations. This evaluation is based on the relative difference (in percent) of root mean squared errors (RMSE): the final score between an experiment e and the reference r is computed as follows:

$$ \Delta \text{RMSE} = 100 \times \left( \frac{\text{RMSE}_r - \text{RMSE}_e}{\text{RMSE}_r} \right), $$

such that positive values indicate that the experiment has lower errors than the reference. The evaluation is performed separately for surface and altitude observations. For the surface, the observations considered are the 2-m temperature, 2-m relative humidity, 10-m wind, and mean sea level pressure (MSLP) recorded from ground based stations over the domain. For the altitude, the evaluation is performed against a comprehensive set of observations: aircraft temperature and wind measurements, SEVIRI radiances in water vapor channels, ground-based GNSS delays, wind profilers, Radar Doppler winds. We propose to summarize this information using scorecards as shown in Figs. 9 and 10. The scorecard presents the mean over the time period of ΔRMSE that are averaged within similar types of observations in the case of the altitude observations. Temperature scores are averaged from aircraft data; SEVIRI water vapor channels and ground based GNSS delays are used to compute scores of humidity and finally aircraft, wind profilers and Radar Doppler winds are finally used to compute scores of wind. The ΔRMSE are computed daily. Confidence intervals for the verification scores are estimated at the 95% level with a bootstrapping method. The statistical significance is represented by filling triangles in Figs. 9 and 10. Finally we compute “NWP indices” for the altitude and the surface that are the total ΔRMSE change when further averaging between variables and time ranges.

Figures 9 and 10 thus present the relative differences of RMSE between the hybrid 3DEnVar AROME experiment and the reference AROME 3DVar. The hybrid 3DEnVar performs better than 3DVar for almost all variables and forecast ranges. The improvement is larger and statistically significant especially in the first 18 h of the forecasts, where it reaches 2%–5% depending on the variable. The impact is smaller on longer-range forecasts, as expected from a change in initial conditions only (e.g., the assimilation scheme is different but the model and the lateral boundary conditions are the same). Compared to M18, the impact is similar but the improvement is not as large for MSLP. This may come from the use of dual resolution in our case, where there is likely more
<figure><img src="image_9.png" /><figcaption>FIG. 8. Relative changes in RMS of (dashed) analysis and (solid) background departures between the 3DEnVar experiment and the 3DVar reference for (top left) wind observations from radiosondes, (top center) temperature observations from radiosondes, (top right) specific humidity observations for radiosondes, (bottom left) wind observations from aircrafts, (bottom center) temperature observations from aircrafts, and (bottom right) for observations of the Special Sensor Microwave Imager/Sounder (SSM/IS), over the whole period (14 May–14 Aug 2018).</figcaption></figure>

imbalance in the ensemble perturbations due to the change of spatial resolution. It may also be due to the shorter horizontal localization used in this study that could degrade the analysis of the larger scales of the MSLP. Finally, the patterns of the scores are quite similar between the 0000 and 1200 UTC runs, but the impact is a bit inferior at 1200 UTC maybe due to the more limited predictability of convection in the afternoon.

### d. Precipitation and wind gust Brier skill scores

For mesoscale models, the RMSE-based evaluation may be completed with probabilistic scores that avoid the double penalty problem sometimes brought by high resolution. The verification of AROME with the regional Brier skill score (BSS) is described in detail by Amodei et al. (2015).

The quantitative comparison is performed against forecast rainfall and wind gusts measured from surface stations using a regional BSS with a spatial neighborhood of 50 km. Finally, a synthetic index is defined as the combination of the scores obtained for cumulative rainfall over 6 h (averaged over thresholds 0.5, 2, and 5 mm) and the maximum wind gusts over 6 h exceeding 40 km h⁻¹.

The BSS of the AROME forecasts using either hybrid 3DEnVar or 3DVar are shown in Fig. 11 for 6-h rainfall at the 5-mm threshold. At 0000 UTC, the experiment is performing better than the reference at all forecast ranges and the difference is significant at ranges from 0 to 30 h. At 1200 UTC, the experiment is also performing better than the reference but this is significant only at the 6-h range. As
ScoreCard 6C0Q vs 5AB6
20180512-20180814: HH00

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

FIG. 9. Relatives changes in RMSE (ΔRMSE) between the reference AROME 3DVar and the experimental 3DEnVar, as a function of the variable considered (ordinate) and of the forecast range (h, abscissa), over the period 12 May–14 Aug 2018. The variables considered are upper-level temperature, humidity, and wind and 2-m temperature, 2-m relative humidity, 10-m winds, and MSLP. An up-pointing (down-pointing) triangle indicates an improvement (a degradation) in RMSE. The areas are proportional to the absolute value of the difference in RMSE, and triangles are filled when the difference is statistically significant at the 95% confidence level (estimated by bootstrapping). Forecasts start at 0000 UTC.

ScoreCard 6C0Q vs 5AB6
20180512-20180814: HH12

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

FIG. 10. As in Fig. 9, but for the evaluation of forecasts starting at 1200 UTC.

stated above, this may be due to the limited predictability of convection.

The BSS for wind gusts are shown in Fig. 12 and the interpretation is very similar than for the verification of precipitation. The 3DEnVar runs at 0000 UTC has better performance than 3DVar at all forecast ranges and this is significant up to 36 h. At 1200 UTC, the performance is significantly better for 3DEnVar up to 18 h.

Because most of the part of the computational cost lies into the ensemble size, it is interesting to evaluate the performance of 3DEnVar with more limited resources beyond the simple sensitivity experiments of section 3. For this purpose, two additional hybrid 3DEnVar experiments have been evaluated on a 2-month period. They differ from the original experiment only by the ensemble size, which has been limited to, respectively, 25 and 35 members. The performance of those experiments is assessed by the AROME synthetic index (Amodei

et al. 2015) and is given in Table 4. With 25 members, there is a decrease of performance that amount to about 0.5%, compared to about 2% difference overall between 3DVar and 50-member 3DEnVar. This confirms that larger ensemble sizes and improvements in the localization should be investigated for AROME EnVar in the future to get full benefit of the scheme (Caron et al. 2019). However, this improvement can also be considered as modest, and we think that increases of the spatial resolution of the EDA should also be considered in future experimentations.

**6. Illustration with a case study**

**a. Mediterranean events**

The Southern part of France is frequently affected by high precipitation events during the autumn period (e.g., Jacques 2016). Frontal systems as well as mesoscale convective systems pump their energy from the warm Mediterranean Sea and are later blocked by the orography, leading to intense precipitation frequently exceeding 200 mm day⁻¹. Recently, an extensive field campaign has been dedicated to study those events, and more widely heavy precipitation and flash flooding occurring in the northwestern Mediterranean basin (Ducrocq et al. 2014). The prediction of those Mediterranean events is known to be sensitive to the initial conditions in humidity (Ducrocq et al. 2002). The AROME model has been deployed in particular to
<figure><img src="image_12.png" /><figcaption>RR6 Threshold 5mm</figcaption></figure>

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

FIG. 11. BSS of the AROME forecasts starting at (top) 0000 and (bottom) 1200 UTC using 3DVar (dashed lines) or hybrid 3DEnVar (solid lines) for 6-h accumulated rainfall exceeding the 5-mm threshold, averaged over the period 12 May–14 Aug 2018 and computed with a neighborhood of about 50 km, see Amodei et al. (2015). The squares indicates when the difference is statistically significant at the 95% confidence level evaluated from a bootstrap method.

improve the prediction of those events that belong to the most devastating meteorological events occurring over France (Seité et al. 2011). Thus, the evaluation of the performance of EnVar scheme to improve the prediction of Mediterranean events is crucial.

b. A high precipitating event in late summer: Montclus

The experiments presented in this study do not cover the autumn period where Mediterranean events are the most frequent. However, in 2018, Southern France was affected by a heavy rainfall event on 9 August. A midlatitude baroclinic storm has crossed the country during the day, yielding cold air in altitude. On the contrary, warm and moist air was arising from the Mediterranean sea. This combination generated quasi-stationary thunderstorms over the Gard and Ardèche regions. Cumulated rainfall over 24 h frequently exceeded 150 mm, with a maximum of 296 mm in Montclus (Gard) and 245 mm in Saint-Martin-d'Ardèche (Ardèche), as shown in Fig. 13.

According to the French National Weather service Météo-France, several rainfall records were broken this day for a

<figure><img src="image_14.png" /><figcaption>FXI6 Threshold 5.6m/s</figcaption></figure>

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

FIG. 12. As in Fig. 11, but for 6-h average wind gusts exceeding the 20 km h⁻¹ threshold (about 5.6 m s⁻¹).

month of August.<fnref n="3" /> This Mediterranean event is rare in summer but not exceptional as there are similar former situations in the same regions on 17 August 2004 or 30 July 1991. The electrical activity was also particularly intense with 42 000 lightning impacts recorded over France this day, which is the record over the year.

c. Behavior of AROME 3DVar and 3DEnVar

A complete and detailed evaluation of this case study is beyond the scope of this article. We just want to point out that hybrid 3DEnVar provides significant improvements in terms of forecast rainfall for this high precipitating event compared to 3DVar. As shown in Fig. 14, the AROME 3DVar did forecast intense precipitation with a maximum of 201 mm, lower by 50% compared to the observed value of 296 mm in Montclus or the 304 mm from the ANTILOPE (Fig. 13). The experiment with hybrid 3DEnVar offers both a sharper localization of the event and a better total amount of precipitation with a maximum of 280 mm. Some other experiments have highlighted the good behavior of the hybrid 3DEnVar scheme in forecasting

<footnote n="3">³ http://www.meteofrance.fr/actualites/65005117-retour-sur-l-episode-mediterraneen-du-9-aout</footnote>
<table><caption>TABLE 4. Synthetic indexes [average of BSS for cumulative rainfall and wind gusts; see Amodei et al. (2015)] for hybrid 3DEnVar experiments with various ensemble sizes, computed over the period from 12 May to 10 Jul.</caption><thead><tr><th>Ensemble size</th><th>0000 UTC</th><th>1200 UTC</th></tr></thead><tbody><tr><td>50</td><td>0.699</td><td>0.698</td></tr><tr><td>35</td><td>0.698</td><td>0.697</td></tr><tr><td>25</td><td>0.695</td><td>0.693</td></tr></tbody></table>

intense precipitation, including the case of the flooding in the Aude region on 15 October 2018.

## 7. Conclusions

This study presents results of a hybrid 3DEnVar scheme developed for the AROME model at kilometer scale with hourly cycling. From an algorithmic point of view and contrary to M18, a standard square-root formulation was used, but this choice is maybe not definitive. The square-root formulation is flexible enough to incorporate both the gridpoint and spectrally formulated localizations previously developed by M18. It is known to work with advanced algorithms such as the multi-incremental approach, the quasi-static implementation and standard second level preconditioners.

Our first result is that the conditioning of the assimilation problem for AROME-France is varying a lot with the meteorological situation. Thus, it is useful to adopt a flexible convergence criterion to mitigate the computational cost. We indicate that in about 65% of the cases, 25 iterations are sufficient to achieve convergence but this number may go up to 50 for some analyses. This would have a limited impact in an operational context in terms of time to solution (where a fixed and predictable length of time for the 3DEnVar to finish is required), but still could save some resources.

<figure><img src="image_16.png" /><figcaption>FIG. 13. Accumulated rainfall on 9 Aug (mm in 24 h) from the ANTILOPE product that estimates quantitative precipitations from a combination of radar data and rain gauges (Champeaux et al. 2009). The maximum is estimated to be 304 mm.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 14. Accumulated rainfall forecast on 9 Aug (mm in 24 h) by (top) the AROME 3DVar and (bottom) the AROME 3DEnVar. The maxima are 201 mm for 3DVar and 280 mm for 3DEnVar.</figcaption></figure>

The performance of hybrid 3DEnVar is shown to be sensitive to the localization length and to the ensemble size, which was expected. There is also sensitivity to the hybrid weights. The hybrid 3DEnVar presented here benefits from rather short localization lengths and equal hybrid weights. This differs a bit from the conclusions of M18, possibly because of 1) the summertime convective period, 2) the higher spatial resolution, and 3) a different version of the EDA in our study.

With a 50-member ensemble, the hybrid 3DEnVar is shown to significantly outperform a version of AROME 3DVar that is very close to the operational setting (from 2019). Most of the impact can be seen at short forecast ranges (up to 12 h), where the difference is significant and reaches about 5% in relative difference of RMSEs. For temperature and humidity, however, the impact can be noticed over the whole forecast period up to 42 h. Precipitation and wind gust probabilistic scores are also improved, by about 2% for the
AROME indicator (Amodei et al. 2015). This typically repre-
sents about 4–5 years of combined improvement for AROME-
France. Tuning the parameters of the scheme over a longer
period may bring further improvements. The scheme has the
potential to improve the prediction of severe high precipitating
events over France, as illustrated by a case study.

The next goal will be high-resolution 4DEnVar with fre-
quent and dense radar and SEVIRI observations. One issue we
are particularly interested in is the nonlinear aspects of this
scheme. Many authors have considered 4DEnVar as a scheme
with 4D covariances, e.g., making the Gaussian assumption on
the whole background trajectories. On the contrary, schemes
such as the iterative ensemble Kalman smoother (IEnKS:
Bocquet and Sakov 2014) are fully nonlinear but they are re-
stricted so far to lower-dimensional models than those encoun-
tered in NWP. We plan to develop a nonlinear, state-space
localized version of 4DEnVar inspired by the IEnKS. This is
believed to be important in the view of nonlinear effects visible
in global 4DVar (Bonavita et al. 2018).

*Acknowledgments.* We thank L. Berre and C. Fischer for their careful reviews of this manuscript. E. Arbogast and T. Montmerle contributed to the implementation by developing the interpolation routines that perform the change of spatial resolution and its adjoint. This research has been supported by the SuNDAE project funded by the ANR, Paris, France, under Grant ANR-18-CE46-0007. This article is a contribution to the HyMeX program.

REFERENCES

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

Anderson, J. L., 2012: Localization and sampling error correction in ensemble Kalman filter data assimilation. *Mon. Wea. Rev.*, **140**, 2359–2371, https://doi.org/10.1175/MWR-D-11-00013.1.

Auger, L., O. Dupont, S. Hagelin, P. Brousseau, and P. Brovelli, 2015: AROME-NWC: A new nowcasting tool based on an operational mesoscale forecasting system. *Quart. J. Roy. Meteor. Soc.*, **141**, 1603–1611, https://doi.org/10.1002/qj.2463.

Auligné, T., A. P. McNally, and D. P. Dee, 2007: Adaptive bias correction for satellite data in a numerical weather prediction system. *Quart. J. Roy. Meteor. Soc.*, **133**, 631–642, https://doi.org/10.1002/qj.56.

Bannister, R. N., 2008: A review of forecast error covariance statistics in atmospheric variational data assimilation. II: Modelling the forecast error covariance statistics. *Quart. J. Roy. Meteor. Soc.*, **134**, 1971–1996, https://doi.org/10.1002/qj.340.

—, 2017: A review of operational methods of variational and ensemble-variational data assimilation. *Quart. J. Roy. Meteor. Soc.*, **143**, 607–633, https://doi.org/10.1002/qj.2982.

Belo-Pereira, M., and L. Berre, 2006: The use of an ensemble approach to study the background error covariances in a global NWP. *Mon. Wea. Rev.*, **134**, 2466–2489, https://doi.org/10.1175/MWR3189.1.

Berre, L., 2000: Estimation of synoptic and mesoscale forecast error covariances in a limited-area model. *Mon. Wea. Rev.*, **128**, 644–667, https://doi.org/10.1175/1520-0493(2000)128<0644:EOSAMF> 2.0.CO;2.

—, S. Ecaterina-Stefanescu, and M. Belo-Pereira, 2006: The representation of the analysis effect in three error simulation techniques. *Tellus*, **58A**, 196–209, https://doi.org/10.1111/j.1600-0870.2006.00165.x.

Bocquet, M., and P. Sakov, 2014: An iterative ensemble Kalman smoother. *Quart. J. Roy. Meteor. Soc.*, **140**, 1521–1535, https://doi.org/10.1002/qj.2236.

Bonavita, M., P. Lean, and E. Holm, 2018: Nonlinear effects in 4D-Var. *Nonlinear Processes Geophys.*, **25**, 713–729, https://doi.org/10.5194/npg-25-713-2018.

Bouttier, F., B. Vié, O. Nuissier, and L. Raynaud, 2012: Impact of stochastic physics in a convection-permitting ensemble. *Mon. Wea. Rev.*, **140**, 3706–3721, https://doi.org/10.1175/MWR-D-12-00031.1.

Bowler, N. E., and Coauthors, 2017: The effect of improved en-
semble covariances on hybrid variational data assimilation.
*Quart. J. Roy. Meteor. Soc.*, **143**, 785–797, https://doi.org/
10.1002/qj.2964.

Brousseau, P., L. Berre, F. Bouttier, and G. Desroziers, 2011: Background-error covariances for a convective-scale data-assimilation system: AROME-France 3D-Var. *Quart. J. Roy. Meteor. Soc.*, **137**, 409–422, https://doi.org/10.1002/qj.750.

—, G. Desroziers, F. Bouttier, and B. Chapnik, 2014: A posteriori diagnostics of the impact of observations on the AROME-France convective-scale data assimilation system. *Quart. J. Roy. Meteor. Soc.*, **140**, 982–994, https://doi.org/10.1002/qj.2179.

—, Y. Seity, D. Ricard, and J. Léger, 2016: Improvement of the forecast of convective activity from the AROME-France system. *Quart. J. Roy. Meteor. Soc.*, **142**, 2231–2243, https://doi.org/10.1002/qj.2822.

Bubnová, R., G. Hello, P. Bénard, and J.-F. Geleyn, 1995: Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/ Aladin NWP system. *Mon. Wea. Rev.*, **123**, 515–535, https://doi.org/10.1175/1520-0493(1995)123<0515:IOTFEE>2.0.CO;2.

Buehner, M., 2005: Ensemble-derived stationary and flow-dependent background error covariances: Evaluation in a quasi-operation NWP settings. *Quart. J. Roy. Meteor. Soc.*, **131**, 1013–1043, https://doi.org/10.1256/qj.04.15.

Caron, J.-F., Y. Michel, T. Montmerle, and E. Arbogast, 2019: Improving background error covariances in a 3D ensemble-variational data assimilation system for regional NWP. *Mon. Wea. Rev.*, **147**, 135–151, https://doi.org/10.1175/MWR-D-18-0248.1.

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

Champeaux, J.-L., P. Dupuy, O. Laurantin, I. Soulan, P. Tabary, and J.-M. Soubeyroux, 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*, **95**, 28–34, https://doi.org/10.1051/lhb/2009052.

Clayton, A. M., A. Lorenc, and D. M. Barker, 2013: Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the Met Office. *Quart. J. Roy. Meteor. Soc.*, **139**, 1445–1461, https://doi.org/10.1002/qj.2054.

Courtier, P., J.-N. Thépaut, and A. Hollingsworth, 1994: A strategy for operational implementation of 4D-Var using an incremental approach. *Quart. J. Roy. Meteor. Soc.*, **120**, 1367–1387, https://doi.org/10.1002/qj.49712051912.

Cuxart, J., P. Bougeault, and J.-L. Redelsperger, 2000: A turbulence scheme allowing for mesoscale and large-eddy simulations.
Quart. J. Roy. Meteor. Soc., **126**, 1–30, https://doi.org/10.1002/qj.49712656202.

Dagum, L., and R. Menon, 1998: OpenMP: An industry standard API for shared-memory programming. *Comput. Sci. Eng. IEEE*, **5**, 46–55, https://doi.org/10.1109/99.660313.

Dee, D. P., 1995: On-line estimation of error covariance parameters for atmospheric data assimilation. *Mon. Wea. Rev.*, **123**, 1128–1145, https://doi.org/10.1175/1520-0493(1995)123<1128:OLEOEC>2.0.CO;2.

—, 2004: Variational bias correction of radiance data in the ECMWF system. *ECMWF Workshop on Assimilation of High Spectral Resolution Sounders in NWP*, ECMWF, Reading, United Kingdom, 97–112, https://www.ecmwf.int/node/8930.

Derber, J., and A. Rosati, 1989: A global oceanic data assimilation system. *J. Phys. Oceanogr.*, **19**, 1333–1347, https://doi.org/10.1175/1520-0485(1989)019<1333:AGODAS>2.0.CO;2.

—, and F. Bouttier, 1999: A reformulation of the background error covariance in the ECMWF global data assimilation system. *Tellus*, **51A**, 195–221, https://doi.org/10.3402/tellusa.v51i2.12316.

Desroziers, G., and L. Berre, 2012: Accelerating and parallelizing minimizations in ensemble and deterministic variational assimilations. *Quart. J. Roy. Meteor. Soc.*, **138**, 1599–1610, https://doi.org/10.1002/qj.1886.

—, J.-T. Camino, and L. Berre, 2014: 4DEnVar: Link with 4D state formulation of variational assimilation and different possible implementations. *Quart. J. Roy. Meteor. Soc.*, **140**, 2097–2110, https://doi.org/10.1002/qj.2325.

Ducrocq, V., D. Ricard, J.-P. Lafore, and F. Orain, 2002: Storm-scale numerical rainfall prediction for five precipitating events over France: On the importance of the initial humidity field. *Wea. Forecasting*, **17**, 1236–1256, https://doi.org/10.1175/1520-0434(2002)017<1236:SSNRPF>2.0.CO;2.

—, and Coauthors, 2014: HyMeX-SOP1: The field campaign dedicated to heavy precipitation and flash flooding in the northwestern Mediterranean. *Bull. Amer. Meteor. Soc.*, **95**, 1083–1100, https://doi.org/10.1175/BAMS-D-12-00244.1.

Evensen, G., 1994: Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics. *J. Geophys. Res.*, **99**, 10 143–10 162, https://doi.org/10.1029/94JC00572.

—, 2003: The ensemble Kalman filter: Theoretical formulation and practical implementation. *Ocean Dyn.*, **53**, 343–367, https://doi.org/10.1007/s10236-003-0036-9.

Fillion, A., M. Bocquet, and S. Gratton, 2018: Quasi-static ensemble variational data assimilation: A theoretical and numerical study with the iterative ensemble Kalman smoother. *Nonlinear Processes Geophys.*, **25**, 315–334, https://doi.org/10.5194/npg-25-315-2018.

Fischer, C., T. Montmerle, L. Berre, L. Auger, and S. Stefanescu, 2005: An overview of the variational assimilation in the ALADIN/France numerical weather-prediction system. *Quart. J. Roy. Meteor. Soc.*, **131**, 3477–3492, https://doi.org/10.1256/qj.05.115.

Fisher, M., 1998: Minimization algorithms for variational data assimilation. *Recent Developments in Numerical Methods for Atmospheric Modelling*, Reading, United Kingdom, ECMWF, 364–385.

—, 2003: Background error covariance modelling. *ECMWF Seminar on Recent Developments in Data Assimilation for Atmosphere and Ocean*, Reading, United Kingdom, ECMWF, 45–63.

—, and P. Courtier, 1995: Estimating the covariance matrices of analysis and forecast error in variational data assimilation. Tech. Memo. 220, ECMWF, 28 pp.

Fujita, T., and R. Purser, 2001: Parallel implementation of compact numerical schemes. Tech. Rep., NCEP, 34 pp.

Gratton, S., and J. Tshimanga, 2009: An observation-space formulation of variational assimilation using a restricted pre-conditioned conjugate gradient algorithm. *Quart. J. Roy. Meteor. Soc.*, **135**, 1573–1585, https://doi.org/10.1002/qj.477.

—, A. Sartenaer, and J. Tshimanga, 2011: On a class of limited memory preconditioners for large scale linear systems with multiple right-hand sides. *SIAM J. Optim.*, **21**, 912–935, https://doi.org/10.1137/08074008.

—, S. Gürol, and P. L. Toint, 2013: Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems. *Comput. Optim. Appl.*, **54**, 1–25, https://doi.org/10.1007/s10589-012-9478-7.

Guidard, V., N. Fourrié, P. Brousseau, and F. Rabier, 2011: Impact of IASI assimilation at global and convective scales and challenges for the assimilation of cloudy scenes. *Quart. J. Roy. Meteor. Soc.*, **137**, 1975–1987, https://doi.org/10.1002/qj.928.

Gürol, S., A. T. Weaver, A. M. Moore, A. Piacentini, H. G. Arango, and S. Gratton, 2014: B-preconditioned minimization algorithms for variational data assimilation with the dual formulation. *Quart. J. Roy. Meteor. Soc.*, **140**, 539–556, https://doi.org/10.1002/qj.2150.

Gustafsson, N., and Coauthors, 2018: Survey of data assimilation methods for convective-scale numerical weather prediction at operational centres. *Quart. J. Roy. Meteor. Soc.*, **144**, 1218–1256, https://doi.org/10.1002/qj.3179.

Haben, S. A., A. S. Lawless, and N. K. Nichols, 2011: Conditioning of incremental variational data assimilation, with application to the Met Office system. *Tellus*, **63A**, 782–792, https://doi.org/10.1111/j.1600-0870.2011.00527.x.

Hamill, T. M., and C. Snyder, 2000: A hybrid ensemble Kalman filter-3D variational analysis scheme. *Mon. Wea. Rev.*, **128**, 2905–2919, https://doi.org/10.1175/1520-0493(2000)128<2905: AHEKFV>2.0.CO;2.

—, J. S. Whitaker, and C. Snyder, 2001: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter. *Mon. Wea. Rev.*, **129**, 2776–2790, https://doi.org/10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2.

Haugen, J. E., and B. Machenhauer, 1993: A spectral limited-area model formulation with time-dependent boundary conditions applied to the shallow-water equations. *Mon. Wea. Rev.*, **121**, 2618–2630, https://doi.org/10.1175/1520-0493(1993)121<2618:ASLAMF>2.0.CO;2.

Houtekamer, P. L., and H. L. Mitchell, 2001: A sequential ensemble Kalman filter for atmospheric data assimilation. *Mon. Wea. Rev.*, **129**, 123–137, https://doi.org/10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2.

—, and F. Zhang, 2016: Review of the ensemble Kalman filter for atmospheric data assimilation. *Mon. Wea. Rev.*, **144**, 4489–4532, https://doi.org/10.1175/MWR-D-15-0440.1.

—, X. Deng, H. L. Mitchell, S.-J. Baek, and N. Gagnon, 2014: Higher resolution in an operational ensemble Kalman filter. *Mon. Wea. Rev.*, **142**, 1143–1162, https://doi.org/10.1175/MWR-D-13-00138.1.

Hunt, B. R., E. J. Kostelich, and I. Szunyogh, 2007: Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. *Physica D*, **230**, 112–126, https://doi.org/10.1016/j.physd.2006.11.008.

Jacques, G., 2016: Les épisodes cévenols: Un aperçu historique. *Meteorologie*, **93**, 50–57.

Kleist, D. T., and K. Ide, 2015: An OSSE-based evaluation of hybrid variational–ensemble data assimilation for the NCEP
GFS. Part I: System description and 3D-hybrid results. *Mon. Wea. Rev.*, **143**, 433–451, https://doi.org/10.1175/MWR-D-13-00351.1.

Lac, C., and Coauthors, 2018: Overview of the Meso-NH model version 5.4 and its applications. *Geosci. Model Dev.*, **11**, 1929–1969, https://doi.org/10.5194/gmd-11-1929-2018.

Lafore, J. P., and Coauthors, 1998: The Meso-NH atmospheric simulation system. Part I: Adiabatic formulation and control simulations. *Ann. Geophys.*, **16**, 90–109, https://doi.org/10.1007/s00585-997-0090-6.

Lascaux, F., E. Richard, and J.-P. Pinty, 2006: Numerical simulations of three different MAP IOPs and the associated microphysical processes. *Quart. J. Roy. Meteor. Soc.*, **132**, 1907–1926, https://doi.org/10.1256/qj.05.197.

Lawson, W. G., and J. A. Hansen, 2004: Implications of stochastic and deterministic filters as ensemble-based data assimilation methods in varying regimes of error growth. *Mon. Wea. Rev.*, **132**, 1966–1981, https://doi.org/10.1175/1520-0493(2004)132<1966:IOSADF>2.0.CO;2.

Liu, C., Q. Xiao, and B. Wang, 2008: An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technical formulation and preliminary test. *Mon. Wea. Rev.*, **136**, 3363–3373, https://doi.org/10.1175/2008MWR2312.1.

Lorenc, A. C., 1986: Analysis methods for numerical weather prediction. *Quart. J. Roy. Meteor. Soc.*, **112**, 1177–1194, https://doi.org/10.1002/qj.49711247414.

—, 2003: The potential of the ensemble Kalman filter for NWP—A comparison with 4D-Var. *Quart. J. Roy. Meteor. Soc.*, **129**, 3183–3203, https://doi.org/10.1256/qj.02.132.

Mahfouf, J.-F., F. Ahmed, P. Moll, and F. N. Teferle, 2015: Assimilation of zenith total delays in the AROME France convective scale model: A recent assessment. *Tellus*, **67A**, 26106, https://doi.org/10.3402/tellusa.v67.26106.

Ménétrier, B., T. Montmerle, L. Berre, and Y. Michel, 2014: Estimation and diagnosis of heterogeneous flow-dependent background-error covariances at the convective scale using either large or small ensembles. *Quart. J. Roy. Meteor. Soc.*, **140**, 2050–2061, https://doi.org/10.1002/qj.2267.

—, —, Y. Michel, and L. Berre, 2015a: Linear filtering of sample covariances for ensemble-based data assimilation. Part I: Optimality criteria and application to variance filtering and covariance localization. *Mon. Wea. Rev.*, **143**, 1622–1643, https://doi.org/10.1175/MWR-D-14-00157.1.

—, —, —, and —, 2015b: Linear filtering of sample covariances for ensemble-based data assimilation. Part II: Application to a convective-scale NWP model. *Mon. Wea. Rev.*, **143**, 1644–1664, https://doi.org/10.1175/MWR-D-14-00156.1.

Mercier, F., S. Gürol, P. Jolivet, Y. Michel, and T. Montmerle, 2018a: Block Krylov methods for accelerating ensembles of variational data assimilations. *Quart. J. Roy. Meteor. Soc.*, **144**, 2463–2480, https://doi.org/10.1002/qj.3329.

—, Y. Michel, T. Montmerle, P. Jolivet, and S. Gürol, 2018b: Speeding up the ensemble data assimilation system of the limited area model of Météo-France using a block Krylov algorithm. *Quart. J. Roy. Meteor. Soc.*, **145**, 910–929, https://doi.org/10.1002/qj.3428.

Michel, Y., 2013: Estimating deformations of random processes for correlation modelling: Methodology and the one-dimensional case. *Quart. J. Roy. Meteor. Soc.*, **139**, 771–783, https://doi.org/10.1002/qj.2007.

—, T. Auligné, and T. Montmerle, 2011: Heterogeneous convective-scale background error covariances with the inclusion of

hydrometeor variables. *Mon. Wea. Rev.*, **139**, 2994–3015, https://doi.org/10.1175/2011MWR3632.1.

Montmerle, T., and C. Faccani, 2009: Mesoscale assimilation of radial velocities from Doppler radars in a preoperational framework. *Mon. Wea. Rev.*, **137**, 1939–1953, https://doi.org/10.1175/2008MWR2725.1.

—, F. Rabier, and C. Fischer, 2007: Relative impact of polar-orbiting and geostationary satellite radiances in the Aladin/ France numerical weather prediction system. *Quart. J. Roy. Meteor. Soc.*, **133**, 655–671, https://doi.org/10.1002/qj.34.

—, Y. Michel, E. Arbogast, B. Ménétrier, and P. Brousseau, 2018: A 3D ensemble variational data assimilation scheme for the limited-area AROME model: Formulation and preliminary results. *Quart. J. Roy. Meteor. Soc.*, **144**, 2196–2215, https://doi.org/10.1002/qj.3334.

Parrish, D., and J. Derber, 1992: The National Meteorological Center's spectral statistical interpolation analysis system. *Mon. Wea. Rev.*, **120**, 1747–1763, https://doi.org/10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2.

Pires, C., R. Vautard, and O. Talagrand, 1996: On extending the limits of variational assimilation in nonlinear chaotic systems. *Tellus*, **48A**, 96–121, https://doi.org/10.3402/tellusa.v48i1.11634.

Purser, R., W. Wu, D. Parrish, and N. Roberts, 2003: Numerical aspects of the application of recursive filters to variational statistical analysis. Part I: Spatially homogeneous and isotropic Gaussian covariances. *Mon. Wea. Rev.*, **131**, 1524–1535, https://doi.org/10.1175/1520-0493(2003)131<1524:NAOTAO>2.0.CO;2.

Raynaud, L., and F. Bouttier, 2015: Comparison of initial perturbation methods for ensemble prediction at convective scale. *Quart. J. Roy. Meteor. Soc.*, **142**, 854–866, https://doi.org/10.1002/qj.2686.

—, L. Berre, and G. Desroziers, 2011: An extended specification of flow-dependent background error variances in the Météo-France global 4D-Var system. *Quart. J. Roy. Meteor. Soc.*, **137**, 607–619, https://doi.org/10.1002/qj.795.

—, —, and —, 2012: Accounting for model error in the Météo-France ensemble data assimilation system. *Quart. J. Roy. Meteor. Soc.*, **138**, 249–262, https://doi.org/10.1002/qj.906.

Saad, Y., 2003: *Iterative Methods for Sparse Linear Systems*. 2nd ed. Society for Industrial and Applied Mathematics, 547 pp.

Schwartz, C. S., 2016: Improving large-domain convection-allowing forecasts with high-resolution analyses and ensemble data assimilation. *Mon. Wea. Rev.*, **144**, 1777–1803, https://doi.org/10.1175/MWR-D-15-0286.1.

Seity, Y., P. Brousseau, S. Malardel, G. Hello, P. Benard, F. Bouttier, C. Lac, and V. Masson, 2011: The AROME-France convective-scale operational model. *Mon. Wea. Rev.*, **139**, 976–991, https://doi.org/10.1175/2010MWR3425.1.

Simonin, D., S. P. Ballard, and Z. Li, 2014: Doppler radar radial wind assimilation using an hourly cycling 3D-Var with a 1.5 km resolution version of the Met Office Unified Model for nowcasting. *Quart. J. Roy. Meteor. Soc.*, **140**, 2298–2314, https://doi.org/10.1002/qj.2298.

Tippett, M. K., J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker, 2003: Ensemble square root filters. *Mon. Wea. Rev.*, **131**, 1485–1490, https://doi.org/10.1175/1520-0493(2003)131<1485:ESRF>2.0.CO;2.

Trémolet, Y., 2007: Incremental 4D-Var convergence study. *Tellus*, **59A**, 706–718, https://doi.org/10.1111/j.1600-0870.2007.00271.x.

Tshimanga, J., S. Gratton, A. Weaver, and A. Sartenaer, 2008: Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation. *Quart. J. Roy. Meteor. Soc.*, **134**, 751–769, https://doi.org/10.1002/qj.228.
Veerse, F., and J.-N. Thepaut, 1998: Multiple-truncation incre-
mental approach for four-dimensional variational data as-
similation. *Quart. J. Roy. Meteor. Soc.*, **124**, 1889–1908, https://
doi.org/10.1002/qj.49712455006.

Wang, X., C. Snyder, and T. M. Hamill, 2007: On the theoretical
equivalence of differently proposed ensemble/3D-Var hybrid
analysis schemes. *Mon. Wea. Rev.*, **135**, 222–227, https://
doi.org/10.1175/MWR3282.1.

Wattrelot, E., O. Caumont, and J.-F. Mahfouf, 2014: Operational
implementation of the 1D+3D-Var assimilation method of
radar reflectivity data in the AROME model. *Mon. Wea. Rev.*,
*142*, 1852–1873, https://doi.org/10.1175/MWR-D-13-00230.1.

Whitaker, J. S., and T. M. Hamill, 2002: Ensemble data as-
similation without perturbed observations. *Mon. Wea.
Rev.*, **130**, 1913–1924, https://doi.org/10.1175/1520-0493(2002)
130<1913:EDAWPO>2.0.CO;2.

Ying, Y., F. Zhang, and J. L. Anderson, 2018: On the selection of
localization radius in ensemble filtering for multiscale quasi-
geostrophic dynamics. *Mon. Wea. Rev.*, **146**, 543–560, https://
doi.org/10.1175/MWR-D-17-0336.1.

Zagar, N., E. Andersson, and M. Fisher, 2005: Balanced tropical
data assimilation based on a study of equatorial waves in
ECMWF short-range forecast errors. *Quart. J. Roy. Meteor.
Soc.*, **131**, 987–1011, https://doi.org/10.1256/qj.04.54.