# Exploring the potential and limitations of weak-constraint 4D-Var

Patrick Laloyaux, M. Bonavita, M. Chrust, Selime Gürol

► To cite this version:

Patrick Laloyaux, M. Bonavita, M. Chrust, Selime Gürol. Exploring the potential and limitations of weak-constraint 4D-Var. Quarterly Journal of the Royal Meteorological Society, 2020, 146 (733), pp.4067-4082. 10.1002/qj.3891 . hal-04737387

HAL Id: hal-04737387

https://cnrs.hal.science/hal-04737387v1

Submitted on 24 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
DOI: 10.1002/qj.3891

Quarterly Journal of the
Royal Meteorological Society

RMetS

RESEARCH ARTICLE

# Exploring the potential and limitations of weak-constraint 4D-Var

P. Laloyaux<fnref n="1" /> | M. Bonavita<fnref n="1" /> | M. Christ<fnref n="1" /> | S. Gürol<fnref n="2" />

**Correspondence**

P. Laloyaux, European Centre for
Medium-Range Weather Forecasts,
Shinfield Park, Reading, RG2 9AX, UK.

Email: patrick.laloyaux@ecmwf.int

**Abstract**

The standard formulation of 4D-Var assumes random zero-mean errors for all sources of information used in the analysis. This assumption is usually not well verified in real-world applications. The performance of a weak-constraint 4D-Var formulation ("forcing" formulation) is studied in this paper in a simplified experimental setting using additive model errors of different length-scales and observing systems of different coverage and accuracy. A set of twin experiments is carried out and results show that weak-constraint 4D-Var can accurately estimate the actual model errors and the initial state only when background and model errors have different spatial scales and when the observations are unbiased and spatially homogeneous. We also present preliminary results from a different weak-constraint 4D-Var formulation ("state" formulation) which could in principle overcome some of these limitations, but at the cost of a substantial increase of computational and memory requirements. These findings help identify the potential but also the intrinsic limitations of the weak-constraint 4D-Var approach. They also help to clarify the experimental results seen in the operational ECMWF analysis system where the analysis and first-guess temperature bias is reduced by up to 50% in the stratosphere.

**KEYWORDS**

model error, scale separation, weak-constraint 4D-Var

# 1 | INTRODUCTION

Errors arise in numerical weather prediction (NWP) from two main sources, incorrect initial conditions and deficiencies in the dynamics and the physical parametrizations of the forecast model. To correct initial errors, 4D-Var adjusts the initial state so as to find the model trajectory that best fits (in a least-square sense) the observations over the assimilation window. However, the model errors which gradually accumulate during the model integration are not taken into account. Considerable efforts are made

to remove biases from models at operational centres, yet their effect on the quality of the analysis remains significant. The standard formulation of 4D-Var is designed to correct random, zero-mean errors in the first-guess trajectory with unbiased observations. If a model has systematic errors, then the model drift causes a bias in the resulting analysis. The size of this analysis bias depends on the accuracy as well as the frequency of the observations. This has significant consequences not only for NWP applications but also for climate reanalysis which typically runs a data assimilation system over decades. The availability of a new

<footnote n="1">¹European Centre for Medium-Range Weather Forecasts, Reading, UK</footnote>

<footnote n="2">²Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique, Toulouse, France</footnote>
observing instrument could change or reduce the analysis
bias but this might be perceived as an apparent change in
the climate (Dee, 2005).

To overcome these difficulties, an assimilation
algorithm should ideally be able to estimate the model
systematic errors. It should then use this information to
adjust subsequent model trajectories by applying a correc-
tion at each model time step throughout the assimilation
interval. Therefore, the model mean errors would rapidly
reduce in magnitude and the original hypotheses of
4D-Var would not be violated any more. As model biases
are expected to vary spatially and temporally, though on
spatial and temporal time-scales significantly longer than
random errors, the model bias estimation and correction
needs to be adaptive. Correcting the model biases at every
model time step also improves the temporal homogeneity
of the solution from one assimilation window to the next.
This feature is particularly beneficial for applications such
as reanalysis as it helps to achieve a better temporal and
spatial consistency in the climate reconstruction.

Accounting for model error in NWP data assimilation
is not a new research topic. The application of variational
analysis to the estimation of the initial conditions was
first proposed in the NWP domain by Sasaki (1970) where
the model constraint to the analysis could be applied
either as an exact (strong) constraint or as an approxi-
mate (weak) constraint. Derber (1989) developed the Vari-
tational Continuous Assimilation where the model initial
state and parameters that represent systematic errors in
the model are estimated simultaneously. The word ‘con-
tinuous’ denotes the fact that the algorithm distributes the
correction throughout the window. Dee (2005) introduced
the term bias-aware assimilation to describe any assimi-
lation algorithm which explicitly takes into account sys-
tematic errors from the model. He found that the presence
of bias can be detected on the input side of the analysis
by monitoring differences between observations and their
model-predicted equivalents, and on the output side by
examining systematic features of the analysis increments.
Separation of different bias sources requires additional
information, such as independent observations, knowl-
edge of the underlying causes, or hypotheses about the
error characteristics of possible sources. Weak-constraint
4D-Var has been introduced by several authors to denote
the same family of algorithms which relax the perfect
model assumption (Wergen, 1992; Zupanski, 1993; Ben-
nett et al., 1996; Vidard et al., 2004; Trémolet, 2006). The
model is considered in the same way as the other sources
of information, taking into account that there is a degree
of uncertainty about the information it can provide on the
evolution of the atmospheric state.

The continued reduction in sources of analysis error
since the introduction of 4D-Var means that it is no longer

possible to ignore the model itself as a source of error.
However, weak-constraint 4D-Var has never fully achieved
its objective of removing model error in an operational
forecasting system up to now. Part of the reason is that
it is intrinsically difficult to reliably diagnose model error
and its statistics (Todling, 2015b; 2015a). Observations are
the only independent source of information available to
estimate the actual model error and the model error statis-
tics. However, it is impossible to separate the estimation
of model errors from that of observation and background
errors inside a data assimilation system, unless ad hoc
assumptions are made on the magnitude, length-scales
and statistical independence of the respective error covari-
ances (Talagrand, 1999; 2003; Desroziers et al., 2005;
Todling, 2015a). Another aspect of the problem is that a
computationally affordable form of the weak-constraint
4D-Var cost function requires making strong assumptions
on the model error statistics (Trémolet, 2007). The ques-
tion then becomes one about the physical plausibility and
realism of the assumptions employed and the limitations
they impose on the type of model errors that can be esti-
mated and corrected.

New diagnostics have been recently developed to study
the biases present in the IFS atmospheric model used at
European Centre for Medium-Range Weather Forecasts
(ECMWF, 2018). Short model integrations over the assim-
ilation window have been compared to radiosonde tem-
perature measurements. While the systematic difference
between the model and these observations is small in the
troposphere, it increases significantly in the stratosphere
where a cold bias of 0.6 °C is observed in the lower/mid
stratosphere and a warm bias of 0.6 °C is observed in the
upper stratosphere in the 9 km resolution forecasts (Shep-
herd et al., 2018; Laloyaux et al., 2020). The Global Posi-
tioning System Radio Occultation (GPSRO) is another type
of instrument that can provide high vertical resolution
refractivity profiles, which depend on pressure, temper-
ature and humidity (Healy, 2007). These measurements
are particularly useful because they can be assimilated
without bias correction and they have strong sensitivity
to temperature structures in the stratosphere (Bauer et al.,
2014). Similar systematic differences between the model
and the Radio Occultation (RO) temperature retrievals
have been found, confirming the presence of the cold and
warm biases in the IFS stratosphere. While radiosondes
tend to be mainly located over land in the Northern Hemi-
sphere, RO is a homogeneous observing system which
provides a global coverage with typically 1,900 profiles per
day (September 2019). This allows the study of the spatial
structure of the model biases. It has been found that the
IFS model contains mainly large-scale biases in the strato-
sphere evolving slowly over a seasonal time-scale. On
the contrary, the background errors used in the ECMWF
<figure><img src="image_2.png" /></figure>

assimilation system have typically smaller spatial scales and evolve on hourly time-scales (Laloyaux et al., 2020). We have found that this distinct property of the model and background errors is the key element to attribute errors to their proper sources.

Following these experimental findings, the operational ECMWF weak-constraint 4D-Var formulation has been revised, with the introduction of a model error covariance matrix whose spectral density is heavily skewed towards the long length-scales (small wave numbers) relative to those characteristics of the operational background-error covariance matrix (Laloyaux et al., 2020). Preliminary results show that the bias in the stratospheric analysis has been reduced by up to 50% with the scale separation technique. The purpose of this paper is to provide a deeper understanding of those results, and illustrate the potential and the limits of applicability of the weak-constraint 4D-Var formulation used at ECMWF. This is done through experimentation in a simplified, controlled testing environment, where the true system evolution and all sources of errors are known.

This article is organised as follows. In Section 2, the weak-constraint formulation in use at ECMWF is recalled with a view to highlight the interaction between model error forcing and initial state analysis. The quasi-geostrophic model used in the numerical experiments is described in Section 3. The performance of weak-constraint 4D-Var in its forcing formulation is assessed in Section 4, where different types of model and/or background errors are introduced in the assimilation system. The state formulation of weak-constraint 4D-Var is introduced in Section 5 and its behaviour is compared to and contrasted with that of the forcing formulation under the same error forcing. Conclusions and further directions of research in the estimation and correction of model error in 4D-Var are finally presented in Section 6.

## 2 | WEAK-CONSTRAINT FORMULATION

The standard 4D-Var formulation assumes that a model state $\mathbf{x}_k$ valid at time $t_k$ in the assimilation window can be computed by a model integration from the initial state $\mathbf{x}_0$ valid at the beginning of the assimilation window

$$ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) \quad \text{for } k = 1, \dots, N, \quad (1) $$

where $\mathcal{M}_{k,k-1}$ denotes the model integration from time $t_{k-1}$ to time $t_k$ and $N$ is the number of model time steps in the assimilation window. This formulation is sometimes called strong-constraint 4D-Var as no errors

are represented in the model integration. With such an assumption, the cost function can be written as a function of only the state at the beginning of the assimilation window

$$ J_{\text{SC}}(\mathbf{x}_0) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_0^b)^{\text{T}}\mathbf{B}^{-1}(\mathbf{x}_0 - \mathbf{x}_0^b) + \frac{1}{2}\sum_{k=0}^{N}\{\mathcal{H}_k\mathcal{M}_{k,0}(\mathbf{x}_0) - \mathbf{y}_k\}^{\text{T}} \mathbf{R}_k^{-1}\{\mathcal{H}_k\mathcal{M}_{k,0}(\mathbf{x}_0) - \mathbf{y}_k\}, \quad (2) $$

where $\mathbf{x}_0^b$ is the background state, a prior estimate of the initial condition, $\mathbf{y}_k$ is the observation vector at time $t_k$ and $\mathcal{H}_k$ is the observation operator mapping the model state in the observation space. The background and observation covariance matrices are represented by $\mathbf{B}$ and $\mathbf{R}_k$ respectively. For linear model and observation operators, the initial state $\mathbf{x}_0$ that minimises this functional is found by imposing the zero-gradient condition

$$ \nabla J_{\text{SC}}(\mathbf{x}_0) = \mathbf{B}^{-1}(\mathbf{x}_0 - \mathbf{x}_0^b) + \sum_{k=0}^{N} \mathbf{M}_{k,0}^{\text{T}} \mathbf{H}_k^{\text{T}} \mathbf{R}_k^{-1} \times \{\mathbf{H}_k \mathbf{M}_{k,0}(\mathbf{x}_0 - \mathbf{x}_0^b) + \mathbf{H}_k \mathbf{M}_{k,0} \mathbf{x}_0^b - \mathbf{y}_k\} = 0, \quad (3) $$

where $\mathbf{M}_{k,0}$ and $\mathbf{H}_k$ are the tangent linear models of the nonlinear model $\mathcal{M}_{k,0}$ and the observation operator $\mathcal{H}_k$ respectively.

The solution of this equation leads to the standard Kalman Filter analysis update (Humpherys et al., 2012)

$$ \mathbf{x}_0 = \mathbf{x}_0^b + \left( \mathbf{B}^{-1} + \sum_{k=0}^{N} \mathbf{M}_{k,0}^{\text{T}} \mathbf{H}_k^{\text{T}} \mathbf{R}_k^{-1} \mathbf{H}_k \mathbf{M}_{k,0} \right)^{-1} \times \sum_{k=0}^{N} \mathbf{M}_{k,0}^{\text{T}} \mathbf{H}_k^{\text{T}} \mathbf{R}_k^{-1} (\mathbf{y}_k - \mathbf{H}_k \mathbf{M}_{k,0} \mathbf{x}_0^b). \quad (4) $$

For nonlinear models and observation operators, the solution to (2) is obtained recursively by linearising the operators around progressively more accurate guess solutions following the Gauss–Newton algorithm (Courtier et al., 1994). In the forcing formulation of weak-constraint 4D-Var (Fisher et al., 2011), the perfect model assumption (1) is relaxed and an unknown forcing $\eta$ is introduced to represent the error which gradually enters into the model trajectory

$$ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) + \eta \quad \text{for } k = 1, \dots, N. \quad (5) $$

In this paper, the model error forcing is assumed to be additive and constant within the assimilation window. We also assume that the model error $\eta$ follows a
Gaussian distribution with no cross-correlation with the background error. This set of rather stringent assumptions ensures that there is no significant interaction between background and model errors in the data assimilation system. This has been verified in the stratospheric part of the IFS model where model error is mainly large-scale evolving slowly over a seasonal time-scale and background error is typically small-scale (Laloyaux et al., 2020). From these assumptions, we can write the weak-constraint 4D-Var cost function as

$$J_{\text{WC}}(\mathbf{x}_0, \eta) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_0^b)^{\text{T}} \mathbf{B}^{-1} (\mathbf{x}_0 - \mathbf{x}_0^b) \\ + \frac{1}{2} \sum_{k=0}^{N} \{ \mathcal{H}_k(\mathbf{x}_k) - \mathbf{y}_k \}^{\text{T}} \mathbf{R}_k^{-1} \{ \mathcal{H}_k(\mathbf{x}_k) - \mathbf{y}_k \} \\ + \frac{1}{2} (\eta - \eta^b)^{\text{T}} \mathbf{Q}^{-1} (\eta - \eta^b), \quad (6)$$

where $\eta^b$ is the prior estimate of the model forcing and $\mathbf{Q}$ is the model error covariance matrix. This "forcing" formulation of weak-constraint 4D-Var estimates simultaneously the initial state $\mathbf{x}_0$ and model forcing $\eta$ that best fit the observations and the background information with respect to their error covariance matrices. To carry the information forward in time to the next assimilation window, the initial condition $\mathbf{x}_0$ is integrated using the model equation (5) and the current estimate of the model bias $\eta$ is defined as the background model error $\eta^b$ in the next assimilation window.

We have seen how the strong-constraint cost function leads to the standard Kalman Filter analysis update equation (4) in the linear case. To simplify the derivation for the weak-constraint case, but without loss of generality, we consider only one time step in the model integration from $t_0$ to $t_1$, with observations available at time $t_0$ and $t_1$. Assuming that the model and observation operators are linear, we introduce the following notations following Fisher et al., (2011)

$$\mathbf{x} = \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \end{bmatrix}, \quad \mathbf{p} = \begin{bmatrix} \mathbf{x}_0 \\ \eta \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} \mathbf{y}_0 \\ \mathbf{y}_1 \end{bmatrix}, \quad (7)$$

$$\mathbf{p} = \mathbf{L}\mathbf{x} = \begin{bmatrix} I & 0 \\ -\mathbf{M}_{1,0} & I \end{bmatrix} \mathbf{x}, \quad (8)$$

$$\mathbf{x} = \mathbf{L}^{-1}\mathbf{p} = \begin{bmatrix} I & 0 \\ \mathbf{M}_{1,0} & I \end{bmatrix} \mathbf{p}, \quad (9)$$

$$\mathbf{D} = \begin{bmatrix} \mathbf{B} & 0 \\ 0 & \mathbf{Q} \end{bmatrix}, \quad \mathbf{R} = \begin{bmatrix} \mathbf{R}_0 & 0 \\ 0 & \mathbf{R}_1 \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} \mathbf{H}_0 & 0 \\ 0 & \mathbf{H}_1 \end{bmatrix}. \quad (10)$$

The control vectors $\mathbf{x}$ or $\mathbf{p}$ can be used to fully characterise any model trajectory. The change of variable $\mathbf{L}$

expresses the relationship between the four-dimensional state and the vector consisting of the initial state and the model error corrections. The linearized weak-constraint cost function (6) can be written as

$$J_{\text{WC}}(\mathbf{p}) = \frac{1}{2}(\mathbf{p} - \mathbf{p}^b)^{\text{T}} \mathbf{D}^{-1} (\mathbf{p} - \mathbf{p}^b) \\ + \frac{1}{2}(\mathbf{H}\mathbf{L}^{-1}\mathbf{p} - \mathbf{y})^{\text{T}} \mathbf{R}^{-1}(\mathbf{H}\mathbf{L}^{-1}\mathbf{p} - \mathbf{y}), \quad (11)$$

where $\mathbf{p}^b = (\mathbf{x}_0^b, \eta^b)^{\text{T}}$ is the background of the new control space. Following a similar approach to what has been done in the strong-constraint case, the minimum of this functional can be found by setting its gradient to zero:

$$\nabla J_{\text{WC}}(\mathbf{p}) = \mathbf{D}^{-1}(\mathbf{p} - \mathbf{p}^b) + \mathbf{L}^{-\text{T}}\mathbf{H}^{\text{T}}\mathbf{R}^{-1}(\mathbf{H}\mathbf{L}^{-1}\mathbf{p} - \mathbf{y}) \\ = 0, \quad (12)$$

which is equivalent to solving the following set of coupled equations in the $(\mathbf{x}_0, \eta)^{\text{T}}$ control space:

$$\mathbf{x}_0 = \mathbf{x}_0^b + \left[ \mathbf{B}^{-1} + \sum_{k=0}^{1} \mathbf{M}_{k,0}^{\text{T}} \mathbf{H}_k^{\text{T}} \mathbf{R}_k^{-1} \mathbf{H}_k \mathbf{M}_{k,0} \right]^{-1} \\ \times [\mathbf{H}_0^{\text{T}} \mathbf{R}_0^{-1}(\mathbf{y}_0 - \mathbf{H}_0 \mathbf{x}_0^b) + \mathbf{M}_{1,0}^{\text{T}} \mathbf{H}_1^{\text{T}} \mathbf{R}_1^{-1} \\ \{\mathbf{y}_1 - \mathbf{H}_1(\mathbf{M}_{1,0} \mathbf{x}_0^b + \eta)\}], \quad (13)$$

$$\eta = \eta^b + (\mathbf{Q}^{-1} + \mathbf{H}_1^{\text{T}} \mathbf{R}_1^{-1} \mathbf{H}_1)^{-1} \mathbf{H}_1^{\text{T}} \mathbf{R}_1^{-1} \{\mathbf{y}_1 - \mathbf{H}_1(\mathbf{M}_{1,0} \mathbf{x}_0 + \eta^b)\}. \quad (14)$$

In the absence of model error ($\eta = 0$), (13) reduces to the strong-constraint analysis (4). In the state update (13), the introduction of a model error term can be seen to be analogous to applying a forcing term to the model integration. The form of the model error update (14) suggests that the model error covariance matrix $\mathbf{Q}$ performs for the model error update a similar role to the one played by the background-error covariance matrix $\mathbf{B}$ in the state update. They both determine how information is spread in the vertical and horizontal directions in space and hence the importance of the prior estimates (Bannister, 2008). If the background error and the model error have similar length-scales, the form of the coupled equations (13–14) suggests that the information from observations at the end of the assimilation window will be split into a state and a model error correction whose magnitude will be inversely proportional to the magnitude of their respective error covariances. In this situation, initial state corrections can potentially be aliased into model error corrections and vice versa. To disentangle the two, apart from accurate estimates of both error covariances $\mathbf{B}^{-1}$ and $\mathbf{Q}^{-1}$, it would be necessary to extend the weak-constraint cost function (6) to explicitly take the cross-covariance between background and model error into account. Given the intrinsic entanglement of background and model error estimates
<figure><img src="image_3.png" /></figure>

(Todling, 2015a), this appears difficult to achieve. (13)–(14) have been derived by Sakov et al., (2018) for the one-lag smoother and by Fillion et al., (2020) for the L-lag smoother where they proved that these equations are valid for any nonlinear evolution model (only the observation operators need to be linear or linearised).

A way to disentangle the effects of different error sources in the state update (13) and the model error update (14) is by enforcing scale separation. This can be achieved by restricting the model error covariance matrix to represent errors on length-scales where the spectral energy density of background errors is small. In this case the assumption of uncorrelated background and model errors becomes more realistic. In practice, this results in small-scale adjustment being done in the state update (13) and large-scale modification in the model error update (14).

## 3 | QG MODEL

The experimental results presented in this paper were obtained using a simple two-layer quasi-geostrophic channel model based on the equations of Fandry and Leslie (1984). The model captures the main dynamical processes that characterise extratropical geophysical flow, and is of sufficient size and complexity to provide a convincing demonstration of the data assimilation algorithms to which it is applied (Fisher and Gürol, 2017). Mean growth of perturbation energy is roughly comparable with current NWP models (Simmons and Hollingsworth, 2002; Snyder et al., 2003) and a relative nonlinearity measure shows that the quasi-geostrophic model is slightly less nonlinear than an NWP model, having a longer duration of the linear regime before nonlinearities dominate the evolution (Gilmour et al., 2001; Fisher and Gürol, 2017).

The equations of the model are expressed in terms of non-dimensional quasi-geostrophic potential vorticity for two well-mixed layers of constant potential temperature (Pedlosky, 1979). The layers are labelled with subscripts 1 and 2, with 1 denoting the upper layer:

$$ 
\begin{aligned}
q_1 &= \nabla^2 \psi_1 - F_1(\psi_1 - \psi_2) + \beta y, \\
q_2 &= \nabla^2 \psi_2 - F_2(\psi_2 - \psi_1) + \beta y + R_s,
\end{aligned} 
\quad (15) 
$$

where $\psi$ denotes streamfunction, $\nabla^2$ is the two-dimensional Laplacian, $\beta$ is the (non-dimensionalised) northward derivative of the Coriolis parameter, $y$ is the vertical coordinate and $R_s$ represents orography or heating. The parameters $F_1$ and $F_2$ serve to couple the layers together, and are defined below. The equations of motion correspond to conservation of

potential vorticity

$$ \frac{Dq_1}{Dt} = \frac{Dq_2}{Dt} = 0. \quad (16) $$

Denoting dimensional time, spatial coordinates and velocities with tildes, the non-dimensionalisation is defined as

$$ 
\begin{aligned}
t &= \tilde{t} \frac{\overline{U}}{L}, & x &= \tilde{x} \frac{\overline{L}}{L}, & y &= \tilde{y} \frac{\overline{L}}{L}, \\
u &= \frac{\tilde{u}}{\overline{U}}, & v &= \frac{\tilde{v}}{\overline{U}}, & \beta &= \overline{\rho_0} \frac{L^2}{\overline{U}}, \\
F_1 &= \frac{f_0^2 L^2}{D_1 g \Delta\theta / \overline{\theta}}, & F_2 &= \frac{f_0^2 L^2}{D_2 g \Delta\theta / \overline{\theta}},
\end{aligned} 
\quad (17) 
$$

where $L$ and $\overline{U}$ represent typical length and velocity, $D_1$ and $D_2$ are the depths of the two layers, $g$ is the acceleration due to gravity, $f_0$ is the Coriolis parameter at the southern boundary of the domain, and $\beta_0$ is its northward derivative. The mean potential temperature is $\overline{\theta}$, and the difference in potential temperature across the layer interface is $\Delta\theta$. For all the experiments described here, we have set $L = 10^6$ m, $\overline{U} = 10$ m·s$^{-1}$, $f_0 = 10^{-4}$ s$^{-1}$, $\beta_0 = 1.5 \times 10^{-11}$ s$^{-1}$·m$^{-1}$ and $\Delta\theta/\overline{\theta} = 0.1$. The layer depths have been set to $D_1 = 6,000$ m and $D_2 = 4,000$ m. For this choice of parameters, the Rossby number is $\epsilon = \overline{U}/(f_0 L) = 0.1$.

The solution algorithm for the model was designed for speed, stability and convenience rather than for accuracy or conservation. Time stepping consists of a simple first-order, upstream, semi-Lagrangian advection of potential vorticity. Interpolation of the upstream potential vorticity is bi-cubic, with advection from outside the domain handled by assuming the potential vorticity outside the domain to be equal to its value at the nearest boundary. Values of potential vorticity along the northern and southern boundaries are supplied by the user. The advecting wind is determined by inverting the potential vorticity equations to determine the streamfunction at each grid point, followed by centred finite-difference approximation of $u = -\partial\psi/\partial y$, $v = \partial\psi/\partial x$. A standard five-point finite-difference approximation of the Laplacian operator is used.

The equations are solved on a domain that is cyclic in the zonal direction. This allows the potential vorticity equations to be decoupled into independent equations for each zonal wavenumber, which can be solved using standard direct techniques. The meridional wind component is zero on the northern and southern boundaries, and the streamfunction on these boundaries is a user-defined constant for each layer which fixes the mean zonal wind in each layer. For the experiments described in this paper, the
<figure><img src="image_4.png" /></figure>

(dimensional) mean wind was 40 m·s⁻¹ in the upper layer and 10 m·s⁻¹ in the lower layer. The domain for the experiments described here was 12,000 km by 6,300 km, and for each layer was discretised into horizontal grids of 40 × 20 points. The dimension of the state vector of the model is thus 1,600. The symmetry of the model is broken by placing a Gaussian hill centred on the point (10, 15), with a (dimensional) height of 2,000 m and an e-folding width of 1,000 km.

## 4 | EXPERIMENTS

We have shown in Section 2 that for the linear case minimising weak-constraint 4D-Var (6) is equivalent to solving the equations for the state update (13) and the model error update (14). As these equations are coupled through the initial condition $x_0$ and the model error forcing $\eta$, there is a risk that the state increment will be aliased into a model error correction and vice versa. The purpose of this paper is to carry out numerical experiments with the quasi-geostrophic model to understand the necessary conditions for weak-constraint 4D-Var to attribute model and background errors to the correct source.

The twin experiment framework is chosen for the numerical experiments to compare the ability of each method to retrieve the true model state. In this framework, synthetic observations are assimilated and the true model state is known at any time. The model is integrated from a horizontally uniform wind in each layer with a vertical shear. As the orography is non-zero and is not symmetrically placed in the domain, this is sufficient to generate a non-trivial flow after a few days of integration. All the experiments are based on the same initial true state $x_0^t$ which has been selected after 15 days of integration. They assimilate synthetic observations that have been generated from the true trajectory adding a white Gaussian noise. The control vector consisting of the values of the streamfunction $\psi$ at every grid point is partially observed every 2 hr (except at the initial time) over the 24-hr assimilation window. At each observation time, 50 locations are randomly selected among the 1,600 grid points of the control vector to replicate the homogeneous coverage provided by satellite soundings in the stratosphere. In this configuration, the total number of observations per assimilation cycle is $N = 600$. A perturbation following a Gaussian distribution with a mean of zero and a (non-dimensional) standard deviation of 0.2 is added to the streamfunction observations to introduce observation errors in the assimilation system. The observation covariance matrix $R$ is specified using the actual value of the perturbation standard deviation. In the first set of experiments presented in Section 4.1, 4.2 and 4.3, we assume that the observations

are not biased or that the biases in the observations have been already corrected by an offline technique. A bias in the observations is introduced in Section 4.4 to illustrate the behaviour of weak-constraint 4D-Var when the observation systematic error has not been accounted for.

Different options are available to introduce an error in the model integration. The simplest one is to change the model time-step and the layer depths to mimic errors due to deficiencies in physical parametrizations, limited spatial resolution, and time-stepping errors (Fisher and Gürol, 2017). In this approach, the error is introduced implicitly by changing some model parameters and the statistical properties of the model errors cannot be explicitly specified. Another method is used in this paper to introduce model error explicitly. In the quasi-geostrophic model, error covariance matrices have isotropic, homogeneous correlations in the horizontal, with Gaussian spatial structure, and with spatially constant vertical correlation. They are fully characterised by their standard deviation $\sigma$, horizontal correlation $c_h$ and vertical correlation $c_v$ lengths. The model error $\eta$ is prescribed by multiplying the square root of an error covariance matrix with a random control vector drawn from a Gaussian distribution of zero mean and unit variance. The model error is then introduced explicitly as an additive forcing following the model equations (5). This approach is more interesting for diagnosing weak-constraint 4D-Var performance as it makes it easier to compare the estimated model error with the true model error introduced in the system.

NWP assimilation system ingesting observations over data-dense areas tends to have background errors with small length-scales. This can be observed in the 20th-century reanalysis CERA-20C (Laloyaux et al., 2018) where the variational Ensemble Data Assimilation that computes the background-error covariance matrix reacts to the observation coverage (Bonavita et al., 2016). The sparse observation coverage at the beginning of the century means that the background error contains large scales, while the denser observation network at the end of the century allows small-scale adjustments. As a relatively dense, homogeneous observing system is assimilated in the quasi-geostrophic model, the parameters of the background-error covariance matrix have been chosen to contain small scales by setting $c_h = 0.6 \times 10^6$, $c_v = 0.2$ and $\sigma_B = 0.8$. This **B** matrix is fixed in all the experiments. The background standard deviation is larger than the observation standard deviation to mimic the ECMWF operational system where the stratospheric temperature background error is larger than the observation error from the RO bending-angle measurements. Two different covariance matrices of model error have been defined, one containing small-scale patterns $Q_s$ ($c_h = 0.6 \times 10^6$ and $c_v = 0.2$) similar to **B**, and one containing large-scale patterns $Q_l$
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 1** A realisation of the model error forcing obtained with (a) $Q_s$ and (b) $Q_l$ for the streamfunction of the model upper layer</figcaption></figure>

($c_h = 1.6 \times 10^6$ and $c_v = 0.8$). The model error forcing is applied at every model time step over the 24-hr assimilation window (144 time steps of 10,min). To give the same weight to the background and model errors, the standard deviation of the model error covariance matrix is given by $\sigma_Q = \sigma_B/144 = 0.005555$. The choice of parameters and the effect on the length-scale are summarised in Table 1.

An ensemble of model error forcing of mean zero and covariance $Q_l$ or $Q_s$ is constructed by multiplying the square root of these matrices by random vectors generated from a normal distribution with a mean of zero and standard deviation of 1. Figure 1 shows one realisation of the model error forcing obtained with (a) $Q_s$ and (b) $Q_l$ for the streamfunction of the model upper layer. This illustrates the effect of the horizontal correlation parameter $c_h$ that produces small- or large-scale model errors.

In this work, the performance of weak-constraint 4D-Var is assessed using twin experiments where model errors and/or background errors are introduced with the statistics defined in Table 1. Model errors and/or background errors are applied to compute the true trajectory

used to create the observations. Weak-constraint 4D-Var is assessed by measuring how much of the actual error introduced in the system is recovered by the assimilation method. To study how the method converges towards the true trajectory over time, the assimilation system is cycled over 30 days with a 24-hr window. The error covariance matrices **B** and **Q** are kept constant over time. To obtain statistically significant results, each experiment has been repeated 100 times with different realisations of the model error and/or background error.

## 4.1 | Model error

We study in this section how weak-constraint 4D-Var can estimate a model error introduced in the system when the initial background is perfect. The true trajectory of the twin experiment is generated by the model which is initialised from $x_0^t$ and integrated over 30 days where a constant model error $\eta^t$ is added at every model time step. The streamfunction observations are generated from the true trajectory adding a random noise $\epsilon_0^t$ of mean zero and standard deviation 0.2. This is summarised in the following equations.

### True Model

$$ \begin{gathered} \mathbf{x}_0 = \mathbf{x}_0^t, \\ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) + \eta^t \quad \text{for } k = 1, \dots, N, \end{gathered} \quad (18) $$

### Observations

$$ y_k = \mathcal{H}_k \mathbf{x}_k + \epsilon_0^t \quad \text{for } k = 1, \dots, N. $$

The prescribed true model error $\eta^t$ is drawn from a normal distribution of mean zero and of either covariance $Q_l$ (containing large scales) or $Q_s$ (containing small scales). In all the experiments, the model error covariance matrix specified in weak-constraint 4D-Var (6) is consistent with the one used to produce the true model error. Parameters used to define background- and observation-error covariance matrices are summarised in Table 1. The assimilation

<table><thead><tr><th>Covariance matrix</th><th>Standard deviation (σ)</th><th>Horizontal correlation (c<sub>h</sub>)</th><th>Vertical correlation (c<sub>v</sub>)</th><th>Scales</th></tr></thead><tbody><tr><td><b>B</b></td><td>0.8</td><td>0.6 × 10<sup>6</sup></td><td>0.2</td><td>Short</td></tr><tr><td><b>Q<sub>s</sub></b></td><td>0.005555</td><td>0.6 × 10<sup>6</sup></td><td>0.2</td><td>Short</td></tr><tr><td><b>Q<sub>l</sub></b></td><td>0.005555</td><td>1.6 × 10<sup>6</sup></td><td>0.8</td><td>Long</td></tr><tr><td><b>Q<sub>k<sub>i</sub></sub></b></td><td>0.8</td><td>1.6 × 10<sup>6</sup></td><td>0.8</td><td>Long</td></tr><tr><td><b>R</b></td><td>0.2</td><td>0.0</td><td>0.0</td><td>Grid point</td></tr></tbody></table>

Note: The model error standard deviation in $Q_s$ and $Q_l$ is obtained by dividing the background standard deviation by the number of time steps in the assimilation window (144).

<caption>**TABLE 1** Standard deviation, horizontal correlation and vertical correlation parameters specified in the error covariance matrices</caption>
<figure><img src="image_7.png" /><figcaption>**FIGURE 2** Model error estimated by weak-constraint 4D-Var after 30 assimilation cycles (a) without scale separation and (b) with scale separation. The circles represent the observation locations for the last assimilation window and their colour represents the analysis departure [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

experiment starts from a perfect model state background $\mathbf{x}^b = \mathbf{x}_0^b$ and from a zero background model error $\eta^b = 0$.

A first assimilation experiment is run where the true model error is the one plotted in Figure 1a and drawn from $\mathbf{Q}_s$. In this configuration, weak-constraint 4D-Var has the same scales in $\mathbf{Q} = \mathbf{Q}_s$ and $\mathbf{B}$ since both covariance matrices used the same set of parameters to specify horizontal and vertical correlations. Figure 2a shows the model error estimated by weak-constraint 4D-Var after 30 assimilation cycles. The markers correspond to the location of the observations for the last cycle and their colour represents the analysis departure (observation minus analysis). A qualitative comparison with the true model error plotted in Figure 1a shows that the estimation of the model error is poor. The very large analysis departure values over the whole domain (RMS departure of 2.4) reveals that the state estimate is also inaccurate. Weak-constraint 4D-Var struggles to attribute correctly the source of error when $\mathbf{B}$ and $\mathbf{Q}$ have similar scales, as a significant part of the model error is corrected by changing the model initial condition. A second assimilation experiment is run where the true model error is the one plotted in Figure 1b and drawn from $\mathbf{Q}_1$. The same assimilation problem has been solved but this time-scale separation between background and model error is present as $\mathbf{Q} = \mathbf{Q}_1$ contains large scales and $\mathbf{B}$ contains small scales. The model error estimated by weak-constraint 4D-Var is plotted in Figure 2b. It captures the correct structure and amplitude of the true model error shown in Figure 1b. Analysis departures are very small (RMS departure of 0.19) and this shows that only the system with scale separation fits the observations closely.

<figure><img src="image_8.png" /><figcaption>**FIGURE 3** Timeseries of the RMSE of (a) the model state analysis and (b) the model error estimate for weak-constraint 4D-Var with scale separation (blue) and without (green) when a systematic error has been introduced in the model equations. RMSE has been averaged over 100 experiments using different realisations of the model error drawn from $\mathbf{Q}_s$ or $\mathbf{Q}_1$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

To get results that are statistically significant, these two experiments are repeated 100 times using different realisations of the model error drawn from $\mathbf{Q}_s$ or $\mathbf{Q}_1$. The performance of each configuration is summarised in Figure 3 where the evolution of the Root-Mean-Square Error (RMSE) of the model state analysis (a) and of the model error estimate (b) is averaged over the 100 experiments. Weak-constraint 4D-Var with scale separation (blue lines) converges towards the true model state and model error. The norm of the initial state increment (not shown) is kept small and decreases over time which means that weak-constraint 4D-Var correctly attributes the error to the model and not to the background. The situation is different when model and background errors have the same scales (green lines); the estimation of the model error is poor and its error increases over time. Weak-constraint 4D-Var with no scale separation generates large initial state increments which demonstrates that the two types of correction are fighting each other as the attribution has not been done properly.

The fit to the assimilated observations can be used as another metric to look at the performance of each configuration. Figure 4 shows the background (dotted) and analysis (solid) RMS departures computed at every observation time (2 hr) over the 30 assimilation windows and averaged for the 100 experiments. This plot confirms the results presented in Figure 2, as only weak-constraint 4D-Var with scale separation correctly fits the observations. The first-guess trajectory used to compute the background
<figure><img src="image_9.png" /></figure>

**FIGURE 4** Timeseries of the background (dotted) and analysis (solid) RMS departures computed at every observation time (2 hr) over the 30 assimilation windows and averaged for the 100 experiments. Weak-constraint 4D-Var with scale separation is in blue and without scale separation is in green. [Colour figure can be viewed at wileyonlinelibrary.com]

departure is corrected using the model error estimate from the previous assimilation (the model error estimate is cycled). In the scale separation case, the background departures are getting smaller over the assimilation cycles, which shows that the estimation of the model error computed by weak-constraint 4D-Var is converging towards the true model error. One can see that the background (dotted) and the analysis (solid) departures are increasing over time within each assimilation window. This behaviour is expected as the model error is additive in the model equations (5).

The observations that have been assimilated in all the experiments so far are distributed randomly in space. This tends to generate a spatially homogeneous distribution which is reminiscent of the type of observation coverage produced by current satellite sounders in the stratosphere. Observations are more heterogeneous in the troposphere, being relative sparse over ocean areas and denser over land. Large inhomogeneities in observation coverage maps also appear around airports where an enormous number of measurements come from all the aircraft landing or taking off. Therefore, it is important to understand how weak-constraint 4D-Var with scale separation can perform with heterogeneous observations distributions. To answer this question, a new set of experiments is carried out where the observing system is changed to mimic the type of dense but localised coverage we would get if an airport was located in the middle of the domain. More precisely, observations are still available every 2 hr, but they are only available in a small circle for the bottom layer. The markers

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

**FIGURE 5** Correlation between the model error estimated by weak-constraint 4D-Var with scale separation and the true model error when observations are available only in a small circle for the bottom layer, for the (a) upper and (b) lower layers. Correlation is computed after 30 days and averaged over 100 experiments using different realisations of the model error drawn from $Q_1$ [Colour figure can be viewed at wileyonlinelibrary.com]

in Figure 5 represent the assimilated observations and the contours shows the correlation between the model error estimated by weak-constraint 4D-Var with scale separation and the true model error drawn from $Q_1$. High correlations are located in the vicinity of the observations and spread horizontally and vertically thanks to the correlations specified in $Q$. It shows that weak-constraint 4D-Var with scale separation can only estimate the model error where observations are available. Although this may well appear to be an obvious conclusion, it has significant consequences for implementation in an operational NWP system. Correcting the model bias only over some areas will generate an unrealistic analysis. This is another reason for restricting weak-constraint 4D-Var in the stratosphere which is more homogeneously observed than the troposphere (Trémolet and Fisher, 2010).

## 4.2 | Background error

Diagnosing the systematic errors of a model is not straightforward and it is likely that weak-constraint 4D-Var has to estimate a bias that is present only for some physical variables, for some pressure levels or over some geographic areas. It is important to verify whether weak-constraint 4D-Var does not estimate a spurious model error forcing when the model has no bias. In this context, the performance of weak-constraint 4D-Var has been assessed when
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 6** Timeseries of the RMSE of (a) the model state analysis and (b) the model error estimate for weak-constraint 4D-Var with scale separation (blue) and without (green) when a background error has been introduced. RMSE has been averaged over 100 experiments using different realisations of the background error drawn from **B** [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

no model error is introduced but a background error $\epsilon_b^t$ is added to the perfect initial state $\mathbf{x}_0^t$ valid at the initial time of the experiment.

**True Model**

$$ \mathbf{x}_0 = \mathbf{x}_0^t + \epsilon_b^t, \\ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) \quad \text{for } k = 1, \dots, N, \qquad (19) $$

**Observations**

$$ \mathbf{y}_k = \mathcal{H}_k \mathbf{x}_k + \epsilon_0^t \quad \text{for } k = 1, \dots, N. $$

The background error $\epsilon_b^t$ is drawn from a normal distribution of mean zero and of covariance **B** which means that it contains small scales. Weak-constraint 4D-Var with scale separation ($\mathbf{Q} = \mathbf{Q}_1$) and without scale separation ($\mathbf{Q} = \mathbf{Q}_s$) was run over 100 experiments using different realisations of the background error. Figure 6 shows the evolution of the RMSE of the model state analysis (a) and of the model error estimate (b), averaged over the 100 experiments. The two configurations of weak-constraint 4D-Var fit the observations by changing the initial state and produce a small model error estimate. The scale separation formulation performs slightly better, having a model error estimate closer to zero.

## 4.3 | Model and background error

In real-world applications, model and background error are present at the same time. To study this situation,

<figure><img src="image_13.png" /><figcaption>**FIGURE 7** Timeseries of the RMSE of (a) the model state analysis and (b) the model error estimate for weak-constraint 4D-Var with scale separation (blue) and without (green) when model and background errors have been introduced. RMSE has been averaged over 100 experiments using different realisations of the background error drawn from **B** and of the model error drawn from $\mathbf{Q}_s$ or $\mathbf{Q}_1$. Results from Figure 3 are plotted as dashed lines to help the comparison [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

100 experiments have been set up where the background error $\epsilon_b^t$ is drawn from a normal distribution of mean zero and covariance **B** and where the model error $\eta^t$ is drawn from a normal distribution of mean zero and either covariance $\mathbf{Q}_1$ or $\mathbf{Q}_s$ depending whether scale separation is enforced or not. This configuration is summarised as

**True Model**

$$ \mathbf{x}_0 = \mathbf{x}_0^t + \epsilon_b^t, \\ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) + \eta^t \quad \text{for } k = 1, \dots, N, \qquad (20) $$

**Observations**

$$ \mathbf{y}_k = \mathcal{H}_k \mathbf{x}_k + \epsilon_0^t \quad \text{for } k = 1, \dots, N. $$

These experiments are essentially summing up what we found when model and background errors were introduced separately. Figure 7 shows the evolution of the RMSE of the model state analysis (a) and of the model error estimate (b), averaged over the 100 experiments. Results from Figure 3 where only model error was introduced are plotted as dashed lines to highlight the effect of adding the initial background error in the system. One can see that scale separation is the only configuration that converges towards the true model state. The initial background error is slightly increasing the RMSE over the first few cycles but has no significant impact on the asymptotic effectiveness of the method.
<figure><img src="image_14.png" /></figure>

Unbiased Observations

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

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

Biased observations

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

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

**FIGURE 8** Histograms of observation, background and analysis errors over the assimilation window have been plotted for (a, c) strong-constraint 4D-Var and (b, d) weak-constraint 4D-Var with scale separation when (a, b) unbiased observations and (c, d) biased observations are asimilated. These are computed for one assimilation experiment and one assimilation cycle where model error from a normal distribution of mean zero and of covariance *Q*₁ is introduced [Colour figure can be viewed at wileyonlinelibrary.com]

**4.4 | Biased observations**

The importance for weak-constraint 4D-Var to have unbiased observations is discussed in this section. In the previous experiments, the distributions of the observation error were unbiased, following a normal law of mean zero and of standard deviation 0.2. In that situation, it has been showed that weak-constraint 4D-Var with scale separation can detect the presence of model biases, and estimate and correct for them in the assimilation cycle. This model error estimation reduces the systematic drift in the model trajectory and produces an analysis with smaller errors. This is summarised in Figure 8a, b where the histograms of observation, background and analysis errors have been plotted for strong-constraint 4D-Var (a)

and weak-constraint 4D-Var with scale separation (b). The statistics are computed over one assimilation window of one experiment and are representative of the overall behaviour.

The background error, computed as the difference between the first-guess trajectory and the true trajectory, exhibits a bias due to the error introduced in the model equations at every time step. Scale separation correctly estimates the model systematic error and produces an analysis with a much smaller bias (–0.06). Strong-constraint 4D-Var struggles to fit the observations by correcting only the initial condition, and the final analysis is biased with a mean error of 0.33. The good performance of weak-constraint 4D-Var with scale separation is partly due the unbiased observations. This important
property helps weak-constraint 4D-Var to correct for the systematic error in the model and to produce an unbiased analysis.

In this last type of experiment, a bias $\mu_0^t$ is introduced when the observations are created from the true trajectory:

$$ \begin{gathered} \underline{\text{True Model}} \\ \mathbf{x}_0 = \mathbf{x}_0^t + \epsilon_b^t, \\ \mathbf{x}_k = \mathcal{M}_{k,k-1}(\mathbf{x}_{k-1}) + \eta^t \quad \text{for } k = 1, \dots, N, \end{gathered} \quad (21) $$

$$ \underline{\text{Observations}} \qquad \mathbf{y}_k = \mathcal{H}_k \mathbf{x}_k + \epsilon_0^t + \mu_0^t \quad \text{for } k = 1, \dots, N. $$

This observation bias is constant and equal to 5 which means that the observation error this time follows a normal distribution of mean 5 and of standard deviation 0.2. Figure 8c,d shows the results of the same experiment where a bias is added to the observations. Weak-constraint 4D-Var produces an analysis with a significant bias of 4.9 and a small standard deviation value of 0.32. The systematic error added to create the observations is absorbed in the model error estimation. The analysis from the strong-constraint 4D-Var is also biased but the standard deviation is larger with a value of 0.63. This example shows that it might be safer to use strong-constraint 4D-Var if observations contain a significant systematic error that has not been accounted for by an observation bias-correction scheme (Auligné et al., 2007). If this is not done, weak-constraint 4D-Var will produce an analysis that will tend to reinforce the observation bias.

The number of degrees of freedom is larger in weak-constraint 4D-Var since both the initial conditions and the model error forcing can be corrected to fit the observations. The method fundamentally relies on the availability of unbiased observations to correctly attribute the sources of the errors. Variational Bias Correction (VarBC) has been developed at ECMWF to remove the observation biases in NWP and reanalysis (Dee and Uppala, 2009). This method requires a certain amount of anchoring information from independent observations (currently RO and radiosondes in the IFS) to remove the systematic errors of satellite radiances and other data. Bias-corrected observations can then be used by weak-constraint 4D-Var to disentangle background and model errors. However, the model error estimates can affect the variational bias corrections in a feedback loop. The importance of anchoring data – accurate observations that do not require bias-correction – is therefore evident to ensure the good behaviour of weak-constraint 4D-Var (Eyre, 2016). How much and what kind of anchoring information is needed is still an open question which will be investigated in the future.

## 5 | 4D-STATE FORMULATION

Solving the four-dimensional weak-constraint 4D-Var problem is not practically possible without introducing some approximations. In Section 2, the model error term $\eta$ introduced in the model evolution equation (5) was assumed to be additive and constant within the assimilation window. These assumptions reduce the size of the four-dimensional problem as any model trajectory can be described from an initial state and a model error forcing. Trémolet (2006) describes another weak-constraint 4D-Var formulation where the control variable is defined as the model state at the start of several intervals in the assimilation window. Within each interval, or subwindow, the forecast model is used to define the state. This is another partial reduction of the size of the weak-constraint 4D-Var control variable where the model is used as a strong constraint between the times when the components of the control variable are defined. The control variable is defined at $m$ regularly spaced times in addition to the initial time. To fix the notations, the assimilation window is split into regular intervals of $p$ time steps each, starting at steps $\{k_i = i \times p\}_{i=0, \dots, m}$. The cost function becomes

$$ \begin{aligned} J_{WC}(\mathbf{x}_0, \mathbf{x}_1, \dots, \mathbf{x}_m) = & \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_0^b)^T \mathbf{B}^{-1} (\mathbf{x}_0 - \mathbf{x}_0^b) \\ & + \frac{1}{2} \sum_{i=0}^{m} \sum_{j=0}^{p-1} \{ \mathcal{H}_{k_i+j} \mathcal{M}_{k_i+j, k_i} (\mathbf{x}_{k_i}) - \mathbf{y}_{k_i+j} \}^T \mathbf{R}_{k_i+j}^{-1} \\ & + \frac{1}{2} \sum_{i=1}^{m} \{ \mathbf{x}_{k_i} - \mathcal{M}_{k_i, k_{i-1}} (\mathbf{x}_{k_{i-1}}) \}^T \mathbf{Q}_{k_i}^{-1} \{ \mathbf{x}_{k_i} - \mathcal{M}_{k_i, k_{i-1}} (\mathbf{x}_{k_{i-1}}) \}. \end{aligned} \quad (22) $$

This expression of the cost function is slightly more complicated than with the forcing formulation of Section 2 but can still be computed without difficulties. With one subwindow ($m=0, p=N$), the control variable reduces to the initial condition and this expression is the strong-constraint 4D-Var cost function defined in (2). When the control variable is defined at every step ($m=n, p=1$), this expression is the full 4D-state weak-constraint cost function where the control vector consists of the four-dimensional state at every model time step. Note that an algebraically equivalent forcing formulation can be achieved by using a model error forcing term $\eta$ which varies with the model time step.

This 4D-state formulation is now assessed in the quasi-geostrophic model where the true state trajectory is computed following the equations (18) where a constant model error forcing is introduced. The model error forcing is drawn from a normal distribution of mean zero and of
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 9 Timeseries of the RMSE of the model state background (dotted) and analysis (solid) for a 24-hr assimilation window divided into 12 subwindows (red) and for a 4-hr assimilation window divided into 2 subwindows (black). RMSE has been averaged over 100 experiments using different realisations of the model error drawn from $Q_1$. The RMSE of the model state analysis shown in Figure 3a (blue solid line) is plotted for comparison [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

covariance matrix $Q_1$ and the experiment is repeated 100 times using different realisations of the model error to get statistically significant results. The horizontal correlation is set to $1.6 \times 10^6$ and the vertical correlation is set to 0.8 in all the model error covariance matrices $Q_{k_i}$ that have to be specified. Since the subwindow formulation computes correction at a given time and not as a forcing, the standard deviation in $Q_{k_i}$ is set to 0.8. No background error is introduced at the initial time of the system and $B$ is therefore defined the same way as $Q_{k_i}$.

A first configuration of the 4D-state formulation is considered where the 24-hr assimilation is divided into 12 subwindows of 2 hr. Figure 9 shows in red the evolution of the RMSE of the model state background (dotted) and analysis (solid), averaged over the 100 experiments. RMSE is computed at the beginning of each subwindow using the true trajectory as a reference. Focusing on the first assimilation cycle, the assimilation method has significantly reduced the RMSE over the whole window, producing an analysis trajectory that fits the observations. However, information about model error gained during the assimilation is not cycled in the 4D-state formulation and the RMSE in the next first-guess trajectory is large because the model evolution is not corrected. This leads to growing values of RMS errors which do not appear to saturate over the test period.

One strategy to prevent a large error developing during the first-guess trajectory is to reduce the length of the

assimilation window. A new set of experiments is carried out where the assimilation window is set to 4 hr with two subwindows of 2 hr. The analysis and background RMSE are represented in black in Figure 9 and show that the 4D-state formulation performs better with this configuration. This subwindows formulation can be compared to the weak-constraint formulation with scale separation presented in Section 4. To facilitate the comparison, the RMSE of the model state analysis shown in Figure 3a (blue solid line) is plotted again in Figure 9. The 4D-state formulation reveals a slightly larger RMSE which could be explained principally by two reasons. The model error introduced in the system is constant and additive following the model equations (5). This gives an advantage to weak-constraint 4D-Var with scale separation as this assumption has been made to derive the formulation, while the 4D-state formulation could estimate a model error that is time-varying over the assimilation window. Furthermore, the choice for the background-error covariance matrix $B$ and the model error covariance matrices $Q_{k_i}$ in the 4D-state formulation are probably suboptimal. The analysis quality is sensitive to their specification, but it appears cumbersome and impractical to estimate optimal correlation length-scales and standard deviations at the different subwindow times. Several sets of parameters for the standard deviation and the horizontal correlation of the covariance matrices have been tested and we have presented the results from the one giving the best results.

Emili et al., (2016) proposed a diagnostic to retrieve the accumulated model error information by looking at the differences between the analysis at time $t_k$ and the propagated analysis from time $t_0$ to $t_k$:

$$ \tilde{\eta}^k = \mathbf{x}_k^a - \mathcal{M}_{k,0} (\mathbf{x}_0^a). \quad (23) $$

Introducing $\mathbf{e}_k^a$ as the analysis error at time $t_k$, they found the approximation

$$ \tilde{\eta}^k \approx \mathbf{e}_k^a - \prod_{\ell=0}^{k-1} \mathbf{M}_{\ell+1,\ell} \mathbf{e}_0^a + \eta^k + \sum_{j=1}^{k-1} \prod_{\ell=j+1}^{k} \mathbf{M}_{\ell,\ell-1} \eta^j. \quad (24) $$

If the difference on the analysis error is small (the first two terms), the quantity $\tilde{\eta}^k$ represents the contributions of the model error on the 4D-state accumulated through the preceding subwindows. This information can then be carried to the next cycle if the model error is assumed to be constant.

Fisher et al., (2011) presented the idea of overlapping assimilation windows as another way to advance the assimilation system forward in time. For each cycle, the analysis window is advanced by a fraction of the window length, resulting in an overlap between consecutive
analysis windows. The overlap of the assimilation windows allows the initial linearisation trajectory for each analysis cycle to be constructed from the overlapping section of the preceding analysis, with only the final part of the linearisation trajectory generated from a forecast. In the overlapping portion of the analysis window, the current and preceding analyses are likely to be similar. As a consequence, analysis increments will be small and the linearisation will be accurate throughout the analysis window. Although the overlapping technique might be considerably more computationally expensive, Fisher *et al.*, (2011) showed the potential of this approach as the 24-hr 4D-Var with 12-hr overlapping outperforms the standard 12-hr 4D-Var in a reanalysis experiment assimilating only surface pressure observations. (Note however that their implementation is not statistically consistent, as is apparent from figure 21 of Fisher *et al.*, 2011.) This cycling technique could be investigated in the subwindow formulation as a way to increase the window length.

## 6 | SUMMARY AND CONCLUSIONS

Inaccuracies in the model formulation and external forcings can produce systematic errors in the model forecast that also affect the analysis cycle. Assimilation algorithms should be able to estimate (or “learn”, as it has become more fashionable to say) these model biases and use the information to adjust subsequent model trajectories. This paper studies under which conditions weak-constraint 4D-Var is able to perform this task.

We have shown that minimising the weak-constraint 4D-Var cost function is equivalent to solving a set of coupled equations in the initial model state and the model error estimate. This fact indicates that, in order to effectively distinguish between initial state and model error corrections, accurate *a priori* statistical information is required on their error statistics **B** and **Q** and their covariance. These error statistics are notoriously difficult to estimate reliably, and their covariance even more so. For this reason other assumptions are needed to disentangle initial state and model error corrections. The numerical experiments carried out with a quasi-geostrophic model showed that the forcing formulation of weak-constraint 4D-Var can perform satisfactorily when background and model errors act on different spatial scales, and this is reflected in their respective error covariances. This situation is well verified in the IFS, where model biases are prevalent in the stratosphere while initial condition errors on smaller spatial scales and shorter temporal scales are prevalent in the troposphere (Laloyaux *et al.*, 2020). The results presented in this work explain the positive impact of this error scale separation idea seen in the operational

IFS, as reported in Laloyaux *et al.*, (2020). The number of degrees of freedom is larger in weak-constraint 4D-Var since the initial condition and the model error forcing can both be corrected to obtain a better fit to the observations. The method also crucially relies on an observing system with a component, sufficient in number, that is unbiased and spatially homogeneous. The importance of anchoring data is therefore evident to ensure the good behaviour of weak-constraint 4D-Var. The interaction of weak-constraint 4D-Var and bias-correction techniques in observation space (e.g., VarBC) is also crucial and it will be explored in future work.

The assumption that the model error is additive and constant within the assimilation window implicitly limits the applicability of the proposed methodology to model errors that evolve slowly on the assimilation window's time-scales. This is consistent with the idea of correcting model errors only on scales significantly larger than background errors ("errors of the day"), and the fact that we cycle model error estimates with the identity model. However this assumption can be relaxed both in the forcing formulation and in the 4D-state weak-constraint 4D-Var formulation. In the latter, the control variable is defined as the model state at the start of several sub-intervals in the assimilation window. We found that this approach does not perform well in the case of systematic model errors, mainly because the information on model error learned during each minimisation is not propagated to the following analysis as the forcing formulation does (albeit only with an identity model). There are at least three possible ways to improve the behaviour of the 4D-state formulation in these experimental set-ups. One is to reduce the length of the assimilation window in order to prevent a large error to develop in the first-guess trajectories. However this effectively calls into question the whole idea of weak-constraint 4D-Var. Another possibility is to advance the assimilation system forward in time using sliding, overlapping windows. We have not tried this approach here but we note that (a) this requires the specification of the model error covariances at the beginning of each subwindow, and (b) this method is significantly more expensive from a computational point of view and, at least in the case of slowly evolving model errors, not efficient. Another alternative may be to diagnose the model error through the subwindows and carry this information to the next cycle. Evaluation of these issues is left for future work.

## ACKNOWLEDGEMENTS

The authors acknowledge the work done by the three anonymous reviewers and the Associate Editor who have provided a critical reading with helpful comments on an earlier version of the manuscript.
<figure><img src="image_21.png" /></figure>

ORCID

P. Laloyaux https://orcid.org/0000-0003-2808-0463

REFERENCES

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, 631–642

Bannister, R.N. (2008) A review of forecast-error covariance statistics in atmospheric variational data assimilation. I: Characteristics and measurements of forecast-error covariances. *Quarterly Journal of the Royal Meteorological Society*, 134, 1951–1970

Bauer, P., Radnóti, G., Healy, S.B. and Cardinali, C. (2014) GNSS radio occultation constellation observing system experiments. *Monthly Weather Review*, 142(2), 555–572

Bennett, F., Chua, A. and Leslie, B. (1996) Generalized inversion of a global numerical weather prediction model. *Meteorology and Atmospheric Physics*, 60, 165–178

Bonavita, M., Hólm, E.V., Isaksen, L. and Fisher, M. (2016) The evolution of the ECMWF hybrid data assimilation system. *Quarterly Journal of the Royal Meteorological Society*, 142, 287–303

Courtier, P., Thépaut, J.-N. and Hollingsworth, A. (1994) A strategy for operational implementation of 4D-Var, using an incremental approach. *Quarterly Journal of the Royal Meteorological Society*, 120, 1367–1388

Dee, D.P. (2005) Bias and data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 131, 3323–3343

Dee, D.P. and Uppala, S. (2009) Variational bias correction of satellite radiance data in the ERA-Interim reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 135, 1830–1841

Derber, J. (1989) A variational continuous assimilation technique. *Monthly Weather Review*, 117(11), 2437–2446

Desroziers, G., Berre, L., Chapnik, B. and Poli, P. (2005) Diagnosis of observation-, background- and analysis-error statistics in observation space. *Quarterly Journal of the Royal Meteorological Society*, 131, 3385–3396

ECMWF (2018). IFS documentation, cycle 45R1. ECMWF, Reading, UK.

Emili, E., Gürol, S. and Cariolle, D. (2016) Accounting for model error in air quality forecasts: an application of 4DENVar to the assimilation of atmospheric composition using QG-Chem 1.0. *Geoscientific Model Development*, 9(11), 3933–3959. https://doi.org/10.5194/gmd-9-3933-2016

Eyre, J.R. (2016) Observation bias correction schemes in data assimilation systems: a theoretical study of some of their properties. *Quarterly Journal of the Royal Meteorological Society*, 142, 2284–2291

Fandry, C.B. and Leslie, L.M. (1984) A two-layer quasi-geostrophic model of summer trough formation in the Australian subtropical easterlies. *Journal of the Atmospheric Sciences*, 41(5), 807–818

Fillion, A., Bocquet, M., Gratton, S., Gürol, S. and Sakov, P. (2020) An iterative ensemble Kalman smoother in presence of additive model error. *SIAM Journal on Uncertainty Quantification*, 8(1), 198–228. https://doi.org/10.1137/19M1244147

Fisher, M. and Gürol, S. (2017) Parallelization in the time dimension of four-dimensional variational data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 143, 1136–1147

Fisher, M., Trémolet, Y., Auvinen, H., Tan, D.G.H. and Poli, P. (2011). Weak-constraint and long window 4DVAR. Technical Memorandum 655, ECMWF, Reading, UK.

Gilmour, I., Smith, L.A. and Buizza, R. (2001) Linear regime duration: Is 24 hours a long time in synoptic weather forecasting?. *Journal of the Atmospheric Sciences*, 58(22), 3525–3539

Healy, S.B. (2007) Operational assimilation of GPS radio occultation measurements at ECMWF. *ECMWF Newsletter*, 111, 6–11

Humpherys, J., Redd, P. and West, J. (2012) A fresh look at the Kalman filter. *SIAM Review*, 54(4), 801–823

Laloyaux, P., de Boisseson, E., Balmaseda, M.A., Bidlot, J., Broennimann, S., Buizza, R., Dalghren, P., Dee, D.P., Haimberger, L., Hersbach, H., Kosaka, Y., Martin, M., Poli, P., Rayner, N., Rustemeier, E. and Schepers, D. (2018) CERA-20C: a coupled reanalysis of the twentieth century. *Journal of Advances in Modeling Earth Systems*, 10, 1172–1195. https://doi.org/10.1029/2018MS001273

Laloyaux, P., Bonavita, M., Dahoui, M., Farnan, J., Healy, S.B., Hólm, E.V. and Lang, S.T.K. (2020) Towards an unbiased stratospheric analysis. *Quarterly Journal of the Royal Meteorological Society*, 146, 2392–2409. https://doi.org/10.1002/qj.3798

Pedlosky, J. (1979) *Geophysical Fluid Dynamics*. Springer, New York, NY.

Sakov, P., Haussaire, J.-M. and Bocquet, M. (2018) An iterative ensemble Kalman filter in the presence of additive model error. *Quarterly Journal of the Royal Meteorological Society*, 144, 1297–1309

Sasaki, Y. (1970) Some basic formalisms in numerical variational analysis. *Monthly Weather Review*, 98(12), 875–883

Shepherd, T.G., Polichtchouk, I., Hogan, R.J. and Simmons, A.J. (2018). Report on Stratosphere Task Force. Technical Memorandum 824, ECMWF, Reading, UK.

Simmons, A.J. and Hollingsworth, A. (2002) Some aspects of the improvement in skill of numerical weather prediction. *Quarterly Journal of the Royal Meteorological Society*, 128, 647–677

Snyder, C., Hamill, T.M. and Trier, S.B. (2003) Linear evolution of error covariances in a quasigeostrophic model. *Monthly Weather Review*, 131(1), 189–205

Talagrand, O. (1999). *A posteriori* verification of analysis and assimilation algorithms, in Proceedings of workshop on diagnosis of data assimilation systems, November 1998. ECMWF, Reading, UK.

Talagrand, O. (2003). Objective validation and evaluation of data assimilation, in Proceedings of seminar on recent developments in data assimilation for atmosphere and ocean, September 2003. ECMWF, Reading, UK.

Todling, R. (2015a) A complementary note to "A lag-1 smoother approach to system-error estimation: sequential method". *Quarterly Journal of the Royal Meteorological Society*, 141, 2917–2922

Todling, R. (2015b) A lag-1 smoother approach to system-error estimation: sequential method. *Quarterly Journal of the Royal Meteorological Society*, 141, 1502–1513

Trémolet, Y. (2006) Accounting for an imperfect model in 4D-Var. *Quarterly Journal of the Royal Meteorological Society*, 132, 2483–2504

Trémolet, Y. (2007) Model-error estimation in 4D-Var. *Quarterly Journal of the Royal Meteorological Society*, 133, 1267–1280

Trémolet, Y. and Fisher, M. (2010) Weak Constraint 4D-Var. *ECMWF Newsletter*, 125, 12–16
<figure><img src="image_22.png" /></figure>

Vidard, P., Piacentini, A. and Dimet, F.-X.L. (2004) Variational data analysis with control of the forecast bias. *Tellus A*, 56(3), 177–188

Wergen, W. (1992) The effect of model errors in variational assimilation. *Tellus A*, 44(4), 297–313

Zupanski, M. (1993) Regional four-dimensional variational data assimilation in a quasi-operational forecasting environment. *Monthly Weather Review*, 121, 2396–2408

**How to cite this article:** Laloyaux P., Bonavita M., Crust M., Gürol S. Exploring the potential and limitations of weak-constraint 4D-Var. *Q J R Meteorol Soc.* 2020;146:4067–4082. https://doi.org/10.1002/qj.3891