Correcting for model changes in statistical postprocessing – an approach based on response theory

For most statistical postprocessing schemes used to correct weather forecasts, changes to the forecast model induce a considerable reforecasting effort. We present a new approach based on response theory to cope with slight model changes. In this framework, the model change is seen as a perturbation of the original forecast model. The response theory allows us then to evaluate the variation induced on the parameters involved in the statistical postprocessing, provided that the magnitude of this perturbation is not too large. This approach is studied in the context of a simple Ornstein–Uhlenbeck model and then on a more realistic, yet simple, quasi-geostrophic model. The analytical results for the former case help to pose the problem, while the application to the latter provides a proof of concept and assesses the potential performance of response theory in a chaotic system. In both cases, the parameters of the statistical postprocessing used – the Error-in-Variables Model Output Statistics (EVMOS) method – are appropriately corrected when facing a model change. The potential application in an operational environment is also discussed.

1 Introduction

A generic property of the atmospheric dynamics is its sensitivity to initial conditions. This implies that probabilistic forecasts is necessary to adequately describe this behaviour (Kalnay, 2003; Wilks, 2011). Indeed, these methods represent a way to go beyond the natural predictability barrier that the chaotic atmospheric models exhibit (Vannitsem, 2017). These forecasts are at the same time subject to the impact of the presence of structural uncertainties, also known as “model errors”. Such errors degrade the forecasts as well, and their impact needs to be mitigated.

Statistical postprocessing methods are used to correct the operational predictions of the atmospheric models. An important family of statistical techniques used to postprocess the forecasts are linear-regression techniques, possibly with multiple predictors (Glahn and Lowry, 1972; Vannitsem and Nicolis, 2008), also known as model output statistics (MOS). This rather simple but very efficient technique can be adapted to ensemble forecasts (e.g. Vannitsem, 2009; Johnson and Bowler, 2009; Glahn et al., 2009; Van Schaeybroeck and Vannitsem, 2015). One of the first approaches that was proposed is called the Error-in-Variable MOS (EVMOS) method because it takes into account the presence of errors in both the observations and model observables (Vannitsem, 2009).

Despite their simplicity, most postprocessing schemes depend on the availability of a database of past forecasts, which allows one to train the regression algorithm by comparison with the observations database. Operational models are however subject to frequent evolution cycles, which are needed to improve their representation of the atmospheric processes. Therefore, there is a continuous need to recompute forecasts starting from past initial conditions with the latest model version to avoid a degradation of the postprocessing schemes due to model change. Such a recomputation of the past forecasts are called “reforecasts” and typically requires a huge data storage and management framework, as well as many computational resources (Hamill, 2018). For instance, the European Centre for Medium-range Weather Forecast (ECMWF) and the National Weather Service (NWS) both produce hundreds of reforecasts every week (Hamill et al., 2013).

Recent research has investigated non-homogeneous regression with a time-adaptative training scheme, for which a trade-off between large training data sets for stable estimates and the benefit of a shorter training period for faster adaptation to data changes is considered (Lang et al., 2020). These results can help mitigate the impact of model change on postprocessing and may call into question the need for reforecast systems. These systems do however help to better represent rare events; they increase the size of the training data sets and greatly improve sub-seasonal forecasts (Scheuerer and Hamill, 2015; Hamill, 2018), which can justify their very high cost.

The present work investigates another research direction and considers a new technique to reduce the cost of adapting a postprocessing scheme to a model change. This method relies on the response theory for dynamical systems (Ruelle, 2009) and assumes that the model change can be written as analytical perturbations of the model tendencies. In this context, parameter modifications as well as new terms in the tendencies are potential model changes.

In Sect. 2, we start by introducing the Ruelle response theory that is used to adapt past postprocessing parameters to new model versions. A didactical example of such an adaptation is considered with a simple Ornstein–Uhlenbeck model in Sect. 3. It is used to describe the methodology and the concept involved. We show that obtaining a new postprocessing scheme after a model change requires the computation of the response of the average of the involved predictors, seen as observables of the system. In the simple case considered, exact analytical results for the response can be obtained up to any order. The correction of the model observables and the postprocessing parameters due to the model change only requires the response-theory corrections up to the second order.

In Sect. 4, a more complex case is considered with a toy model of atmospheric variability in the form of a two-layer quasi-geostrophic model with an orography. We compute the linear response of the parameters of the postprocessing scheme for two model change experiments involving a modification of the friction and the horizontal temperature gradient of the model. The response-theory approach provides an efficient correction of the postprocessing scheme up to a lead time of 4 d, which matches the lead-time window where the scheme's correction is efficient.

In the last section, we discuss the implications that this new method could have on operational forecast postprocessing systems, as well as new research avenues.

2 Response theory

The systems used to produce the weather forecasts are typically non-linear dynamical systems whose time evolution is governed by multi-dimensional ordinary differential equations:

$$\begin{array}{}\text{(1)}& \dot{\mathit{y}}=\mathit{F}(t,\mathit{y}).\end{array}$$

The generic chaotic nature of these systems for some parameter values implies that they are sensitive to the initial data used to produce the forecasts. For such chaotic dynamical systems, one can assume that a well-defined time-invariant measure exists with which the averages are performed. However, the existence of such measures has been proved for systems that are uniformly hyperbolic, and they are called Sinai–Ruelle–Bowen (SRB) measures (Young, 2002), but rigorous proofs for other systems are rather difficult to obtain. A way to proceed is then to continue as if physical systems were uniformly hyperbolic. This assumption is called the Gallavotti–Cohen hypothesis (Gallavotti and Cohen, 1995a, b). With this assumption, response theory has been successfully used in various weather and climate-related problems (Demaeyer and Vannitsem, 2018; Vissio and Lucarini, 2018; Lembo et al., 2019; Bódai et al., 2020). Indeed, the systems used to produce weather forecasts are typically not uniformly hyperbolic, but thanks to the aforementioned hypothesis, one can still use what will follow and compare with the results obtained with experiments.¹ It is the rationale behind the formal presentation of the linear response theory for general systems like Eq. (1) in Ruelle (1998a). The main concepts that will be used in this article are now introduced.

2.1 Perturbations of dynamical systems

We shall assume for simplicity that the system defined by Eq. (1) is autonomous and given by

$$\begin{array}{}\text{(2)}& \dot{\mathit{y}}=\mathit{F}\left(\mathit{y}\right).\end{array}$$

In the general setting considered, let us assume that any given probability measure converges to a unique invariant measure ρ under the time evolution given by the Liouville equation of Eq. (2). This measure is used to compute the average of an arbitrary observable A (a smooth function of the state y) of the system, which is given by

$$\begin{array}{}\text{(3)}& \langle A{\rangle }_{\mathit{y}}=\int \mathit{\rho }\left(\mathrm{d}\mathit{y}\right)A\left(\mathit{y}\right),\end{array}$$

and assuming the ergodicity of the system, a time average of the observable A along a trajectory of the system on its attractor can be equivalently performed as

$$\begin{array}{}\text{(4)}& \langle A{\rangle }_{\mathit{y}}=\underset{\mathrm{T}\to \mathrm{\infty }}{lim}{\displaystyle \frac{\mathrm{1}}{\mathrm{T}}}\underset{\mathrm{0}}{\overset{\mathrm{T}}{\int }}\mathrm{d}\mathit{\tau }A\left(\mathit{y}\left(\mathit{\tau }\right)\right),\end{array}$$

where y(τ) is a solution of Eq. (2). If a perturbation Ψ of the dynamical system is introduced in the original system at time τ=0 as

$$\begin{array}{}\text{(5)}& \dot{\mathit{y}}=\mathit{F}\left(\mathit{y}\right)+\mathbf{\Psi }\left(\mathit{y}\right),\end{array}$$

it induces a perturbation of the observable's average, which at the first order is given by²

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\int \mathit{\rho }\left(\mathrm{d}{\mathit{y}}_{\mathrm{0}}\right)\mathit{\delta }A\left({\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)\\ & =\int \mathit{\rho }\left(\mathrm{d}{\mathit{y}}_{\mathrm{0}}\right)\mathit{\delta }\mathit{y}(\mathit{\tau }{)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}A,\end{array}\end{array}$$

where f^τ is the flow of the system defined by Eq. (2) mapping an initial condition y₀ to the system's state at time τ as y(τ)=f^τ(y₀), and the capital T represents the transposition. δy is the perturbation of the trajectory of the system induced by the perturbation Ψ. This formula gives the transient response to the perturbation, and the long-time average of the integrand of Eq. (6) gives the stationary response to the perturbation, i.e. the sensitivity δ〈A〉_y of the system observables to the perturbation (Eyink et al., 2004; Wang, 2013). The higher-order corrections δ^k〈A(τ)〉_y can in principle be computed as well but are quite complicated to obtain for chaotic dynamical systems; see for instance Lucarini (2009). We will show an analytically tractable case in Sect. 3.

2.2 The tangent linear model

The linear perturbation δy of the trajectories of Eq. (2) can be computed by introducing y+δy in Eq. (5) to get at the first order

$$\begin{array}{}\text{(7)}& \dot{\mathit{\delta }\mathit{y}}={\mathbf{\nabla }}_{\mathit{y}}\mathit{F}\cdot \mathit{\delta }\mathit{y}+\mathbf{\Psi }\left(\mathit{y}\right),\end{array}$$

where y is the solution of Eq. (2) and ∇_yF is the Jacobian matrix evaluated along this solution. Therefore, both Eqs. (2) and (7) have to be integrated simultaneously. In the weather forecasting context, this latter linearised equation without the perturbation term Ψ is called the tangent linear model of Eq. (2) (Kalnay, 2003). Here, Eq. (7) is initialised with δy(0)=0 and provides the linear response of the trajectory y(τ) to the perturbation Ψ. It is thus assumed that there is no interference due to initial-condition errors in the perturbation problem. Note however that the effects on the trajectories of both the initial-condition perturbation and the Ψ perturbation can be investigated through this equation by setting δy(0)≠0, although we are not aware of any study of the response to both types of perturbations together.

The tangent model provides thus the tool through which we will evaluate the impact of the model change on the average used by statistical postprocessing schemes. In other words, the tangent model will allow us to take into account the information on the model change (viewed as a perturbation of the initial model) to modify the previous postprocessing scheme and adapt it to the new model version. The solution to Eq. (7) with δy(0)=0 is by

$$\begin{array}{}\text{(8)}& \mathit{\delta }\mathit{y}\left(\mathit{\tau }\right)=\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\mathbf{M}\left(\mathit{\tau }-{\mathit{\tau }}^{\prime },{\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathbf{0}}\right)\right)\cdot \mathbf{\Psi }\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathbf{0}}\right)\right),\end{array}$$

where M is the fundamental matrix of Eq. (7) (Gaspard, 2005; Nicolis, 2016) defined as

$$\begin{array}{}\text{(9)}& \mathbf{M}(\mathit{\tau },\mathit{y})={\mathbf{\nabla }}_{\mathit{y}}{\mathit{f}}^{\mathit{\tau }},\end{array}$$

which is the solution of the homogeneous equation $\dot{\mathbf{M}}={\mathbf{\nabla }}_{\mathit{y}}\mathit{F}\cdot \mathbf{M}$. Using the chain rule, the response defined by Eq. (6) is rewritten in terms of the perturbation alone:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\int \mathit{\rho }\left(\mathrm{d}{\mathit{y}}_{\mathrm{0}}\right)\mathbf{\Psi }{\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\\ & \cdot {\mathbf{\nabla }}_{{\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)}A\left({\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right),\end{array}\end{array}$$

where the causality of the perturbation acting on the system and perturbing the averaged observable appears (Lucarini, 2008), since ${\mathit{\tau }}^{\prime }<\mathit{\tau }$. We will also use this alternative expression throughout the article. Note that when the initial perturbation δy(0) is not equal to 0, additional terms to Eqs. (8) and (10) will appear. These will not be addressed here, but some references to this can be found in Nicolis et al. (2009) and Nicolis (2016).

2.3 Non-stationary response theory

Equation (6) gives the transient, non-stationary response to the perturbation, evaluated for averages computed with the invariant measure. However, in this work, we need to evaluate the response to perturbations for averages computed with non-stationary measures evolving in time. In that sense, it is a “non-stationary response theory”, which is performed with an arbitrary initial probability density. As such, all the formulas presented are valid if the measure being used is the measure at the time when the perturbation is introduced (τ=0), as shown in Appendix A. In this case, other usual formulas obtained through substitution, for instance to obtain an adjoint representation of Eq. (10), should be used with care, since the measure is no longer invariant and an extra Jacobian term appears in the integrand.

We will thus also assume that the measures ρ_τ being used are absolutely continuous with respect to the Lebesgue measure. In this case, we can write ρ_τ(dy)=ρ_τ(y) dy. We now present the problem of model change in the framework of postprocessing and show on a simple stochastic model³ how response theory allows us to tackle the issue.

3 A simple analytical example

In order to get a first impression of the impact of a model change on a postprocessing scheme, we consider two Ornstein–Uhlenbeck processes representing reality x(τ) and a model y(τ) of reality. These processes obey the following equations:

$$\begin{array}{}\text{(11)}& {\displaystyle }\dot{x}\left(\mathit{\tau }\right)=-{\mathit{\lambda }}_{x}x\left(\mathit{\tau }\right)+K_x+Q_x{\mathit{\xi }}_x\left(\mathit{\tau }\right),\text{(12)}& {\displaystyle }\dot{y}\left(\mathit{\tau }\right)=-{\mathit{\lambda }}_{y}y\left(\mathit{\tau }\right)+K_y+Q_y{\mathit{\xi }}_y\left(\mathit{\tau }\right),\end{array}$$

where ξ_x and ξ_y are Gaussian white-noise processes such that

$$\begin{array}{rl}\langle {\mathit{\xi }}_x\left(\mathit{\tau }\right)\rangle =\langle {\mathit{\xi }}_y\left(\mathit{\tau }\right)\rangle & =\mathrm{0}\\ \langle {\mathit{\xi }}_x\left(\mathit{\tau }\right){\mathit{\xi }}_x\left({\mathit{\tau }}^{\prime }\right)\rangle & =\mathit{\delta }(\mathit{\tau }-{\mathit{\tau }}^{\prime })\\ \langle {\mathit{\xi }}_y\left(\mathit{\tau }\right){\mathit{\xi }}_y\left({\mathit{\tau }}^{\prime }\right)\rangle & =\mathit{\delta }(\mathit{\tau }-{\mathit{\tau }}^{\prime })\\ \langle {\mathit{\xi }}_x\left(\mathit{\tau }\right){\mathit{\xi }}_y\left({\mathit{\tau }}^{\prime }\right)\rangle & =\mathrm{0}.\end{array}$$

These are therefore uncorrelated Ornstein–Uhlenbeck processes with noise amplitudes Q_x and Q_y.

We then consider a change Ψ_y of model y(τ), possibly improving or degrading the forecast performances as

$$\begin{array}{}\text{(13)}& \dot{\widehat{y}}\left(\mathit{\tau }\right)=-{\mathit{\lambda }}_y\widehat{y}\left(\mathit{\tau }\right)+K_y+Q_y{\mathit{\xi }}_y\left(\mathit{\tau }\right)+{\mathrm{\Psi }}_y\left(\mathit{\tau }\right),\end{array}$$

where

$$\begin{array}{}\text{(14)}& {\mathrm{\Psi }}_y\left(\mathit{\tau }\right)=-\mathit{\kappa }(\mathit{\delta }K+\mathit{\delta }Q{\mathit{\xi }}_y(\mathit{\tau }\left)\right)\end{array}$$

with $\mathit{\delta }K=K_y-K_x$ and $\mathit{\delta }Q=Q_y-Q_x$. It can represent, for example, a better parameterisation of subgrid-scale processes or an increase of the model resolution. Note that the best correction is obtained if κ=1.

We have thus reality x(τ) and two different models of it: y(τ) and $\widehat{y}\left(\mathit{\tau }\right)$. We now want to evaluate the difference between a postprocessing scheme constructed before the model change (with the past forecasts of model y(τ)) and one constructed after (with the past forecasts of model $\widehat{y}\left(\mathit{\tau }\right)$).

3.1 The postprocessing method

We now consider a forecast situation where model y is initialised at time τ=0 with a perfect observation of reality: $y\left(\mathrm{0}\right)=x\left(\mathrm{0}\right)=x_{\mathrm{0}}$. We use the Error-in-Variables Model Output Statistics (EVMOS) postprocessing scheme (Vannitsem, 2009) to correct the forecasts of model y based on these initial conditions. In this context, given N past forecasts y_n and observations x_n, the correction of the univariate EVMOS postprocessing of variable x from a new forecast y(τ) is provided by the linear regression

$$\begin{array}{}\text{(15)}& y_C\left(\mathit{\tau }\right)=\mathit{\alpha }\left(\mathit{\tau }\right)+\mathit{\beta }\left(\mathit{\tau }\right)y\left(\mathit{\tau }\right).\end{array}$$

The coefficients α and β are obtained by minimising the functional

$$\begin{array}{}\text{(16)}& J\left(\mathit{\tau }\right)=\sum _{n=\mathrm{1}}^N{\displaystyle \frac{{\left[\mathit{\left\{}\mathit{\alpha }\right(\mathit{\tau })+\mathit{\beta }(\mathit{\tau }\left)y_n\right(\mathit{\tau }\left)\mathit{\right\}}-x_n\left(\mathit{\tau }\right)\right]}^{\mathrm{2}}}{{\mathit{\sigma }}_x^{\mathrm{2}}\left(\mathit{\tau }\right)+{\mathit{\beta }}^{\mathrm{2}}\left(\mathit{\tau }\right){\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)}}\end{array}$$

and are thus given by the following equations:

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathit{\alpha }\left(\mathit{\tau }\right)& {\displaystyle }=\langle x\left(\mathit{\tau }\right)\rangle -\mathit{\beta }\left(\mathit{\tau }\right)\langle y\left(\mathit{\tau }\right)\rangle ,\text{(18)}& {\displaystyle }\mathit{\beta }\left(\mathit{\tau }\right)& {\displaystyle }=\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_x^{\mathrm{2}}\left(\mathit{\tau }\right)}{{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)}}},\end{array}$$

where

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{\sigma }}_x^{\mathrm{2}}\left(\mathit{\tau }\right)& {\displaystyle }=〈\left(x\right(\mathit{\tau })-\langle x(\mathit{\tau })\rangle {)}^{\mathrm{2}}〉,\text{(20)}& {\displaystyle }{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)& {\displaystyle }=〈\left(y\right(\mathit{\tau })-\langle y(\mathit{\tau })\rangle {)}^{\mathrm{2}}〉,\end{array}$$

The averages 〈⋅〉 are taken over an ensemble of past forecasts and observations. This approach has been developed to obtain a correct climatological forecast calibration. It constitutes a simple setting in which the impact of model changes can be evaluated and corrected. More sophisticated approaches can be evaluated in the future (other linear MOS schemes, ensemble MOS, etc.).

Since we are dealing with simple analytical models here, we can compute the theoretical values of the coefficient α and β with an infinite ensemble of past forecasts, and the averaged quantities involved in this computation are then given by the averages of an infinite number of realisations of the Ornstein–Uhlenbeck processes, as if we had an infinite ensemble of past forecasts.

3.2 Averaging the Ornstein–Uhlenbeck processes

For reality x and model y, we directly get the averages (Gardiner, 2009)

$$\begin{array}{}\text{(21)}& {\displaystyle }\langle x\left(\mathit{\tau }\right)\rangle =\langle x_{\mathrm{0}}\rangle e^{-{\mathit{\lambda }}_x\mathit{\tau }}+{\displaystyle \frac{K_x}{{\mathit{\lambda }}_x}}\left(\mathrm{1}-e^{-{\mathit{\lambda }}_x\mathit{\tau }}\right),\text{(22)}& {\displaystyle }{\mathit{\sigma }}_x^{\mathrm{2}}\left(\mathit{\tau }\right)={\mathit{\sigma }}_{x_{\mathrm{0}}}^{\mathrm{2}}{e}^{-\mathrm{2}{\mathit{\lambda }}_x\mathit{\tau }}+{\displaystyle \frac{Q_x^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_x}}\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_x\mathit{\tau }}\right)\end{array}$$

and

$$\begin{array}{}\text{(23)}& {\displaystyle }\langle y\left(\mathit{\tau }\right)\rangle =\langle x_{\mathrm{0}}\rangle e^{-{\mathit{\lambda }}_y\mathit{\tau }}+{\displaystyle \frac{K_y}{{\mathit{\lambda }}_y}}\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right),\text{(24)}& {\displaystyle }{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)={\mathit{\sigma }}_{x_{\mathrm{0}}}^{\mathrm{2}}{e}^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}+{\displaystyle \frac{Q_y^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_y}}\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}\right),\end{array}$$

where we note that the model is initialised with the same initial conditions as reality:

$$\begin{array}{}\text{(25)}& \langle y\left(\mathrm{0}\right)\rangle =\langle x\left(\mathrm{0}\right)\rangle =\langle x_{\mathrm{0}}\rangle ,{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathrm{0}\right)={\mathit{\sigma }}_x^{\mathrm{2}}\left(\mathrm{0}\right)={\mathit{\sigma }}_{x_{\mathrm{0}}}^{\mathrm{2}}.\end{array}$$

We get the postprocessing coefficients before the model change α(τ) and β(τ) by inserting these expressions in Eqs. (17) and (18).

Similarly, we get the same kind of results for model $\widehat{y}$, after model change Ψ_y:

$$\begin{array}{}\text{(26)}& {\displaystyle }\langle \widehat{y}\left(\mathit{\tau }\right)\rangle =\langle x_{\mathrm{0}}\rangle e^{-{\mathit{\lambda }}_y\mathit{\tau }}+{\displaystyle \frac{K_y-\mathit{\kappa }\mathit{\delta }K}{{\mathit{\lambda }}_y}}\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right),\text{(27)}& {\displaystyle }{\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}\left(\mathit{\tau }\right)={\mathit{\sigma }}_{x_{\mathrm{0}}}^{\mathrm{2}}{e}^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}+{\displaystyle \frac{(Q_y-\mathit{\kappa }\mathit{\delta }Q{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_y}}\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}\right).\end{array}$$

We also obtain the postprocessing coefficients after the model change $\widehat{\mathit{\alpha }}\left(\mathit{\tau }\right)$ and $\widehat{\mathit{\beta }}\left(\mathit{\tau }\right)$ (see also the analysis in Vannitsem, 2011). We can also compute the variation of the bias α:

$$\begin{array}{}\text{(28)}& \widehat{\mathit{\alpha }}\left(\mathit{\tau }\right)-\mathit{\alpha }\left(\mathit{\tau }\right)=\mathit{\delta }\mathit{\alpha }\left(\mathit{\tau }\right)=\mathit{\beta }\left(\mathit{\tau }\right)\langle y\left(\mathit{\tau }\right)\rangle -\widehat{\mathit{\beta }}\left(\mathit{\tau }\right)\langle \widehat{y}\left(\mathit{\tau }\right)\rangle .\end{array}$$

The ratio between the parameters β is given by

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\widehat{\mathit{\beta }}\left(\mathit{\tau }\right)}{\mathit{\beta }\left(\mathit{\tau }\right)}}=\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)}{{\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}\left(\mathit{\tau }\right)}}}.\end{array}$$

For $\mathit{\tau }\gg max(\mathrm{1}/{\mathit{\lambda }}_x,\mathrm{1}/{\mathit{\lambda }}_y)$, we note that this ratio tends to

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\widehat{\mathit{\beta }}\left(\mathit{\tau }\right)}{\mathit{\beta }\left(\mathit{\tau }\right)}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\kappa }\mathit{\delta }Q/Q_y}},\end{array}$$

and the difference between the biases α of the two models is approximatively given by

$$\begin{array}{}\text{(31)}& \mathit{\delta }\mathit{\alpha }\left(\mathit{\tau }\right)\approx -\mathit{\beta }\left(\mathit{\tau }\right){\displaystyle \frac{K_y}{{\mathit{\lambda }}_y}}\left[{\displaystyle \frac{\mathrm{1}-\mathit{\kappa }\mathit{\delta }K/K_y}{\mathrm{1}-\mathit{\kappa }\mathit{\delta }Q/Q_y}}-\mathrm{1}\right].\end{array}$$

Let us now assume that model change Ψ_y can be considered as a perturbation of the initial model y. Using response theory, the averages $\langle \widehat{y}\rangle$ and ${\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}$ can be estimated using the initial model y instead of the perturbed model $\widehat{y}$. In turn, these new estimated averages give us the new postprocessing scheme coefficients $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\beta }}$. We now detail the results obtained by using this method.

3.3 Model change and response theory

After the model change, the forecasts are provided by model $\widehat{y}$, and their time evolution is given by Eq. (13). This model can be seen as a perturbation of model y by the term Ψ_y given by Eq. (14). In such a case, given an observable A, its average after the model change can then be related to its original average by

$$\begin{array}{}\text{(32)}& \langle A\left(\mathit{\tau }\right){\rangle }_{\widehat{y}}=\langle A\left(\mathit{\tau }\right){\rangle }_y+\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_y+{\mathit{\delta }}^{\mathrm{2}}\langle A\left(\mathit{\tau }\right){\rangle }_y+\mathrm{\dots },\end{array}$$

where the averages on the right-hand side are taken over the forecasts of model y. Response theory allows us to obtain the average over the model $\widehat{y}$ forecasts (the left-hand side) based solely on the average over the model y forecasts. The $\widehat{y}$ model forecasts are therefore not required to estimate the new postprocessing scheme.

The observables depend on the lead time τ of the forecast, as do the parameters α and β which determine the postprocessing correction for every lead time. This reflects the fact that the postprocessing problem is typically a non-stationary initial-value problem, since the initial conditions of the model Eqs. (12) and (13) are typically not chosen on their respective model attractor but rather as observations⁴ of reality defined by Eq. (11). As a consequence, the model averages of Eq. (32) relax toward the stationary response in the long-time limit, and the stationary response theory (Ruelle, 2009; Wang, 2013) cannot provide us their short-time relaxation behaviour. Instead, the Ruelle time-dependent response theory should be used (Ruelle, 1998a). It follows that, if the perturbation (14) is small, then the first order is given by (see Sect. 2)

$$\begin{array}{}\text{(33)}& \mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_y=\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\int \mathrm{d}{x}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left(x_{\mathrm{0}}\right)〈{\mathrm{\Psi }}_y\left({\mathit{\tau }}^{\prime }\right){\mathrm{\nabla }}_{f^{{\mathit{\tau }}^{\prime }}\left(x_{\mathrm{0}}\right)}A\left(f^{\mathit{\tau }}\left(x_{\mathrm{0}}\right)\right)〉,\end{array}$$

where ρ₀ is the distribution of the initial conditions (observations) used to initialise the models. ∇_x is the gradient evaluated at the point x, and here it is the simple derivative. As indicated by Eq. (25), in the postprocessing framework, ρ₀ is taken as the stationary distribution of reality. As shown in Appendix A, Eq. (33) can be obtained through a Kubo-type perturbative expansion (Lucarini, 2008). We remark that this example deals with stochastic models, due to which we have to perform an additional averaging over the realisations of the stochastic processes, denoted here as 〈⋅〉 (Lucarini, 2012). Finally the mapping f^τ which appears in Eq. (33) is the stochastic flow

$$\begin{array}{}\text{(34)}& f^{\mathit{\tau }}\left(x_{\mathrm{0}}\right)=x_{\mathrm{0}}{e}^{-{\mathit{\lambda }}_y\mathit{\tau }}+\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }{e}^{-{\mathit{\lambda }}_y(\mathit{\tau }-{\mathit{\tau }}^{\prime })}\left[Q_y{\mathit{\xi }}_y\right({\mathit{\tau }}^{\prime })+K_y].\end{array}$$

This maps an initial condition x₀ of model y to the state f^τ(x₀) of a realisation of this model at the later lead time τ. The principle of causality is thus implicit in Eq. (33), which estimates the impact of the perturbation Ψ_y on the subsequent perturbed-model time evolution by developing around the unperturbed-model y trajectories.

Evaluating Eq. (33) and its stochastic integrals (Gardiner, 2009) gives us the variation of the averages 〈y(τ)〉 and 〈y(τ)²〉 to the perturbation Ψ_y:

$$\begin{array}{}\text{(35)}& \mathit{\delta }\langle y\left(\mathit{\tau }\right){\rangle }_y=-\mathit{\kappa }\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\mathit{\delta }K{e}^{-{\mathit{\lambda }}_y(\mathit{\tau }-{\mathit{\tau }}^{\prime })}=-{\displaystyle \frac{\mathit{\kappa }}{{\mathit{\lambda }}_y}}\mathit{\delta }K\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right),\end{array}$$

$$\begin{array}{}\text{(36)}& \begin{array}{rl}\mathit{\delta }\langle y(\mathit{\tau }{)}^{\mathrm{2}}{\rangle }_y& =-\mathrm{2}\mathit{\kappa }\mathit{\delta }K\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\left[\langle x_{\mathrm{0}}\rangle e^{-{\mathit{\lambda }}_y(\mathrm{2}\mathit{\tau }-{\mathit{\tau }}^{\prime })}\right.\\ & \left.+{\displaystyle \frac{K_y}{{\mathit{\lambda }}_y}}{e}^{-{\mathit{\lambda }}_y(\mathit{\tau }-{\mathit{\tau }}^{\prime })}\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right)\right]\\ & -\mathrm{2}\mathit{\kappa }\mathit{\delta }Q{Q}_y\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }{e}^{-\mathrm{2}{\mathit{\lambda }}_y(\mathit{\tau }-{\mathit{\tau }}^{\prime })}\\ & =-\mathrm{2}\mathit{\kappa }{\displaystyle \frac{\mathit{\delta }K}{{\mathit{\lambda }}_y}}〈y\left(\mathit{\tau }\right)〉\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right)\\ & -{\displaystyle \frac{\mathit{\kappa }}{{\mathit{\lambda }}_y}}\mathit{\delta }Q{Q}_y\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}\right).\end{array}\end{array}$$

Rearranging these two terms, we also get the following expression for the variation of the variance given by Eq. (24):

$$\begin{array}{}\text{(37)}& \begin{array}{rl}\mathit{\delta }{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)& =-{\displaystyle \frac{\mathit{\kappa }}{{\mathit{\lambda }}_y}}\mathit{\delta }Q{Q}_y\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}\right)\\ & -{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{{\mathit{\lambda }}_y^{\mathrm{2}}}}\mathit{\delta }{K}^{\mathrm{2}}{\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right)}^{\mathrm{2}}.\end{array}\end{array}$$

Note that the variation given by Eq. (35) corresponds exactly to the difference between the average of the two models $\langle \widehat{y}\left(\mathit{\tau }\right)\rangle -\langle y\left(\mathit{\tau }\right)\rangle$. On the other hand the variation given by Eq. (37) lacks the term of order κ² involving δQ that appears in the difference between ${\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}\left(\mathit{\tau }\right)$ and ${\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)$ given respectively by Eqs. (27) and (24). Instead, another term of order κ² and involving δK is present, indicating that higher-order terms of response theory need to be considered to correct it (Ruelle, 1998b). The second-order term is given by the expression⁵ (Lucarini, 2012):

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{2}}\langle A\left(\mathit{\tau }\right){\rangle }_y& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\underset{{\mathit{\tau }}^{\prime }}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime \prime }\int \mathrm{d}y{\mathit{\rho }}_{\mathrm{0}}\left(x_{\mathrm{0}}\right)〈{\mathrm{\Psi }}_y\left({\mathit{\tau }}^{\prime }\right){\mathrm{\nabla }}_{f^{{\mathit{\tau }}^{\prime }}\left(x_{\mathrm{0}}\right)}\\ & {\mathrm{\Psi }}_y\left({\mathit{\tau }}^{\prime \prime }\right){\mathrm{\nabla }}_{f^{{\mathit{\tau }}^{\prime \prime }}\left(x_{\mathrm{0}}\right)}A\left(f^{\mathit{\tau }}\left(x_{\mathrm{0}}\right)\right)〉.\end{array}\end{array}$$

Applying this to the first moment of the y models directly yields

$$\begin{array}{}\text{(39)}& {\mathit{\delta }}^{\mathrm{2}}\langle y\left(\mathit{\tau }\right){\rangle }_y=\mathrm{0}.\end{array}$$

On the other hand, integrating the stochastic integrals present in this expression for the moment 〈y(τ)²〉 gives

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{2}}\langle y(\mathit{\tau }{)}^{\mathrm{2}}{\rangle }_y& ={\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}\mathit{\delta }{K}^{\mathrm{2}}}{{\mathit{\lambda }}_y^{\mathrm{2}}}}{\left(\mathrm{1}-e^{-{\mathit{\lambda }}_y\mathit{\tau }}\right)}^{\mathrm{2}}\\ & +{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}\mathit{\delta }{Q}^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda }}_y}}\left(\mathrm{1}-e^{-\mathrm{2}{\mathit{\lambda }}_y\mathit{\tau }}\right),\end{array}\end{array}$$

which corrects the κ²δK² term in Eq. (37) and makes the response theory up to order 2 exactly match the difference between ${\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}\left(\mathit{\tau }\right)$ and ${\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)$, for every lead time τ. In fact, the subsequent orders of the response vanish due to the linearity of the simple Ornstein–Uhlenbeck models, which enables us to truncate the response Kubo-like expansion to the second order. Finally, this shows that the (non-stationary) response theory can be used to estimate the postprocessing parameters after the model change based on the forecasts of the initial model. Indeed, instead of the averages $\langle \widehat{y}\left(\mathit{\tau }\right)\rangle$ and ${\mathit{\sigma }}_{\widehat{y}}^{\mathrm{2}}\left(\mathit{\tau }\right)$, the approximate averages 〈y(τ)〉+δ〈y(τ)〉_y and ${\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)+\mathit{\delta }{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)+{\mathit{\delta }}^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}\left(\mathit{\tau }\right)$ can be used to compute $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\beta }}$. We emphasise that the second-order contribution had to be considered in order to obtain the exact result. Nevertheless, the difference between the first- and the second-order response is of the order κ², which implies that for a small perturbation (model change), the first order will generally be a sufficiently good approximation. A more detailed derivation of the results obtained in this section can be found in the Supplement.

In order to investigate this research avenue on a case closer to those encountered in reality, we will now consider the application of postprocessing and response theory to a low-order atmospheric model displaying chaos.

4 Application to a low-order atmospheric model

A two-layer quasi-geostrophic atmospheric system on a β plane with an orography is considered (Charney and Straus, 1980; Reinhold and Pierrehumbert, 1982). This spectral model possesses well-identified large-scale flow regimes, such as “zonal” and “blocked” regimes. The horizontal nondimensionalised coordinates are denoted as x and y, with the model's domain being defined by $(\mathrm{0}\le x\le \frac{\mathrm{2}\mathit{\pi }}{n},\mathrm{0}\le y\le \mathit{\pi })$, with $n=\mathrm{2}{L}_y/L_x$ as the aspect ratio between its meridional and zonal extents L_y and L_x. The two main fields of this model are the 500 hPa pressure anomaly and temperature, which are proportional to the barotropic streamfunction ψ(x, y) and the baroclinic streamfunction θ(x, y), respectively. Both fields are defined in a zonally periodic channel with no-flux boundary conditions in the meridional direction ($\partial \cdot /\partial x\equiv \mathrm{0}$ at $y=\mathrm{0},\mathit{\pi }$). The fields are expanded in Fourier modes respecting these boundary conditions:

$$\begin{array}{rl}{F}_{\mathrm{1}}(x,y)& =\sqrt{\mathrm{2}}\cos \left(y\right),\\ F_{\mathrm{2}}(x,y)& =\mathrm{2}\cos \left(nx\right)\sin \left(y\right),\\ F_{\mathrm{3}}(x,y)& =\mathrm{2}\sin \left(nx\right)\sin \left(y\right),\\ F_{\mathrm{4}}(x,y)& =\sqrt{\mathrm{2}}\cos \left(\mathrm{2}y\right),\\ & \mathrm{⋮}\end{array}$$

such that

$$\begin{array}{}\text{(41)}& {\mathrm{\nabla }}^{\mathrm{2}}{F}_i(x,y)=-a_i^{\mathrm{2}}{F}_i(x,y)\end{array}$$

with eigenvalues $a_{\mathrm{1}}^{\mathrm{2}}=\mathrm{1},a_{\mathrm{2}}^{\mathrm{2}}=a_{\mathrm{3}}^{\mathrm{2}}=\mathrm{1}+n^{\mathrm{2}},a_{\mathrm{4}}^{\mathrm{2}}=\mathrm{4},\mathrm{\dots }$. We have thus the following decomposition

$$\begin{array}{}\text{(42)}& {\displaystyle }\mathit{\psi }(x,y)=\sum _{i=\mathrm{1}}^{n_{\mathrm{a}}}{\mathit{\psi }}_i{F}_i(x,y),\text{(43)}& {\displaystyle }\mathit{\theta }(x,y)=\sum _{i=\mathrm{1}}^{n_{\mathrm{a}}}{\mathit{\theta }}_i{F}_i(x,y),\end{array}$$

where n_a is the number of modes of the spectral expansion. The partial differential equations controlling the time evolution of the fields ψ(x, y) and θ(x, y) can then be projected on the Fourier modes to finally give a set of ordinary differential equations for the coefficients ψ_i and θ_i

$$\begin{array}{}\text{(44)}& \dot{\mathit{x}}=\mathit{F}\left(\mathit{x}\right),\mathit{x}=({\mathit{\psi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\psi }}_{n_{\mathrm{a}}},{\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_{n_{\mathrm{a}}})\end{array}$$

that can be solved with usual numerical integrators. All variables are nondimensionalised. The ordinary differential equations of the model are detailed in Appendix B.

In the version proposed by Reinhold and Pierrehumbert using the 10 first modes beyond a certain value of the zonal temperature gradient, the system displays chaos and makes transitions between the blocked and zonal flow regimes embedded in its global attractor. Here, we use their main nondimensionalised parameters values: the friction at the interface between the two layers k_d=0.1, the friction at the bottom surface $k_d^{\prime }=\mathrm{0.01}$ and the aspect ratio of the domain n=1.3. The β plane lies at midlatitudes (50^∘) and the Coriolis parameter f₀ is set accordingly.

In the present work, the parameter h_d, the Newtonian cooling coefficient is fixed to 0.3 instead of the value found in Reinhold and Pierrehumbert (which is h_d=0.045). Two additional fields have to be specified on the domain: ${\mathit{\theta }}^{*}(x,y)$, the radiative equilibrium temperature field, and h(x, y), the topographic height field. These fields can be decomposed by projecting them onto the eigenfunctions of the Laplacian as before. The corresponding coefficients ${\mathit{\theta }}_i^{*}$ and h_i then allow for writing these fields as sums of weighted eigenfunctions:

$$\begin{array}{}\text{(45)}& {\displaystyle }{\mathit{\theta }}^{*}(x,y)=\sum _{i=\mathrm{1}}^{n_{\mathrm{a}}}{\mathit{\theta }}_i^{*}{F}_i(x,y),\text{(46)}& {\displaystyle }h(x,y)=\sum _{i=\mathrm{1}}^{n_{\mathrm{a}}}{h}_i{F}_i(x,y).\end{array}$$

In the present case, we consider that the only non-zero coefficients are ${\mathit{\theta }}_{\mathrm{1}}^{*}=\mathrm{0.2}$ and h₂=0.4, meaning that the radiative equilibrium profile is given by the zonally varying function $\sqrt{\mathrm{2}}\cos \left(y\right)$ and the orography is made of a mountain and a valley shaped by the function 2 cos (nx) sin (y). Again, the value of the temperature gradient ${\mathit{\theta }}_{\mathrm{1}}^{*}$ is larger than the one chosen in Reinhold and Pierrehumbert (which is ${\mathit{\theta }}_{\mathrm{1}}^{*}=\mathrm{0.1}$) to increase the chaotic variability in the system. Trajectories of variables θ₁ and ψ₂ are depicted in Fig. 1, for the reference system (reality) and a model version (model 0) for which the friction coefficient has been slightly modified.

Figure 1 Dynamics of the reference system and model 0 of the postprocessing experiment with a modification of the friction coefficient (see Table 1) for (a) time evolution of the variable θ₁ and (b) time evolution of the variable ψ₂.

Table 1 The main parameters used and modified in the experiments. Model 0 and model 1 are respectively the forecast model of reality before and after the model change.

These parameter changes induce slight modifications of the dynamics. In particular the system possesses two distinct weather regimes, depicted in Fig. 2a and b: one characterised by a zonal circulation (see Fig. 2c) and another characterised by a blocking situation (see Fig. 2d). In the former case, the variables ψ₂ and ψ₃ characterising the strength of the meridional anomalies are small, while in the latter case they are large, indicating indeed a blocking situation. This is different from the situation considered in Reinhold and Pierrehumbert (1982), where two different blocking regimes coexist with the zonal regime.

Figure 2 Attractors for the experiment with a modification of the friction coefficient: (a) two-dimensional isodensity of the attractors estimated with a Gaussian kernel density estimator for the variables ψ₂ and ψ₃ and (b) two-dimensional scatterplot of the attractors for the variables ψ₂ and ψ₃. The attractors of reality and model 0 are qualitatively similar, with two different parts which are indicated by ellipses. The blue and red crosses correspond respectively to equilibrium points of the reference model (reality) and of model 0, respectively. The dashed ellipse corresponds on average to a zonal circulation depicted in panel (c). The dash-dotted ellipse corresponds on average to a blocking situation depicted in panel (d). In both panels (c) and (d), the underlying colour map denotes the orography in the domain, and the contours denote the geopotential height anomaly at 500 hPa.

4.1 Postprocessing experiments

The model described above with 10 modes (n_a=10) is used, and two different postprocessing experiments are performed, one involving the Newtonian cooling parameter h_d and another involving the friction parameter k_d between the two atmospheric layers. The parameter values detailed above correspond to the long-term reference (i.e. reality). A first model is defined (model 0) which is a copy of the two-layer quasi-geostrophic model defining reality, but the parameters h_d or k_d are slightly changed; i.e. the model error of the forecasting system lies in either the Newtonian cooling or the friction parameter. Then, as in Sect. 3, a model change is imposed, leading to another forecasting model (model 1) that can either improve or degrade the model error by a factor κ. The parameter variations involved in these experiments are detailed in Table 1. Without a loss of generality, we consider model changes that improve the representation of reality in the sense that the amplitude of the model error in model 1 is smaller than in model 0. The effect of the model change is depicted in Figs. 3 and 4 for the friction parameter experiment. These figures display the mean and the standard deviation of the model forecasts and observations coming from the reference forecasts, as a function of the lead time τ. We have used a set of 1 million trajectories of each system to compute these averages.

Figure 3 Behaviour of the averages as a function of the lead time τ in reality and the forecast models before (a) and after (b) the model change, for the experiment with a modification of the friction coefficient (see Table 1). The variable considered is the temperature meridional gradient θ₁. The solid lines denote the mean, while the shaded areas denote the interval of 1 standard deviation.

Figure 4 Same as Fig. 3 but for the variable ψ₃ of the streamfunction ψ.

In the framework of the EVMOS postprocessing scheme, the predictors and the predictands are the same nominal variable, and no other predictors are used. In both experiments considered, the postprocessing parameters α and β of the EVMOS for model 0, as well as $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{\beta }}$ for model 1, are computed. The main objective here is then to estimate the difference between the former and the latter using Ruelle response theory. The approach in a multivariate setting is presented below.

4.2 Model change, response theory and the tangent linear model

Let us consider again the response theory described in Sect. 3.3 but in the general multivariate deterministic case described in Sect. 2. In the postprocessing framework, models 0 and 1 evolve in time from a set of initial conditions taken outside of their respective attractors. Response formulas found in Ruelle's work have to be adapted to take this into account. One therefore has to consider the density of initial conditions as the measure. For a system with a time-independent perturbation $\mathbf{\Psi }\left(\widehat{\mathit{y}}\right)$,

$$\begin{array}{}\text{(47)}& \dot{\widehat{\mathit{y}}}=\mathit{F}\left(\widehat{\mathit{y}}\right)+\mathbf{\Psi }\left(\widehat{\mathit{y}}\right)=\widehat{\mathit{F}}\left(\widehat{\mathit{y}}\right),\end{array}$$

an observable A with average 〈A(τ)〉_y at the lead time τ for the system

$$\begin{array}{}\text{(48)}& \dot{\mathit{y}}=\mathit{F}\left(\mathit{y}\right)\end{array}$$

has a first-order response of

$$\begin{array}{}\text{(49)}& \mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}=\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathit{\delta }\mathit{y}(\mathit{\tau }{)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}A,\end{array}$$

where f^τ is the flow of the unperturbed system given by Eq. (48), ρ₀ is the distribution of initial conditions and δy(τ) is the solution of the equation $\dot{\mathit{y}}+\dot{\mathit{\delta }\mathit{y}}=\widehat{\mathit{F}}(\mathit{y}+\mathit{\delta }\mathit{y})$, which can be approximated at the first order by the following linear inhomogeneous differential equation

$$\begin{array}{}\text{(50)}& \dot{\mathit{\delta }\mathit{y}}={\mathbf{\nabla }}_{\mathit{y}}\mathit{F}\cdot \mathit{\delta }\mathit{y}+\mathbf{\Psi }\left(\mathit{y}\right),\end{array}$$

where y(τ) is the solution of the unperturbed Eq. (48) with the initial condition y(0)=y₀, and we see that the systems of Eqs. (48) and (50) have to be integrated simultaneously (Gaspard, 2005). The homogeneous part of Eq. (50) is the well-known tangent linear model of the system, and here it has to be solved with an additional boundary term which is the perturbation itself.

Equation (49) is derived in Appendix A and can be computed in the same way as the averages depicted in Figs. 3 and 4, by averaging over multiple initial conditions of the reference system. Since we initialise the unperturbed (model 0) and perturbed systems (model 1) with the same initial conditions, the initial state of the tangent model defined by Eq. (50) is δy(0)=0. Therefore we do not estimate the impact of the observation or assimilation errors but rather the direct impact of the model errors viewed as time-independent perturbations. The formulation of the problem and Eq. (50) can be adapted to take these errors into account, as described for instance by Nicolis (2016).

In what follows, we will numerically integrate Eq. (50) to evaluate the response on the average due to the perturbation induced by the model change. This will in turn, as in Sect. 3, enable us to compute the postprocessing parameters for the new model.

4.3 Main results

For each of the two experiments detailed in Table 1, we start by obtaining 1 million observations of reality that will be used to initialise the forecast models. For each observation, this is done by starting model x (the reference) with a random initial condition and running it for a very long time (100 000 nondimensionalised time units) to achieve convergence to its global attractor. Once the observations have been obtained, we run the reference model, model 0 and model 1 over 200 time units (corresponding to roughly 22 d) to obtain reality and the forecasts. The systems have been integrated using the fourth-order Runge–Kutta integration scheme with a time step of 0.1 time units corresponding to 16.15 min. The averaging over the 1 million trajectories of reality and of the forecasts at each lead time is used to compute the postprocessing coefficients α and β of the EVMOS by using Eqs. (17) and (18). For each predictand, the corresponding model variable is used as the unique predictor.

Figure 5 Corrections of the moments of θ₁ from model 0 to model 1 using the response-theory formula of Eq. (49), for the experiment with a modification of the friction coefficient.

Figure 6 Corrections of the moments of θ₁ from model 0 to model 1 using the response-theory formula of Eq. (49), for the experiment with a modification of the Newtonian cooling coefficient.

The response-theory approximations of the averages of model $\widehat{\mathit{y}}$ (model 1) averages are obtained by integrating the linearised equations of model 0 along its trajectories with the perturbation Ψ as an inhomogeneous term. This is done by integrating Eq. (50) over a lead time of 200 time units with a zero initial condition, using the same integration scheme as before. It gives us the integrand of Eq. (49) for each trajectory, and the integral is then approximated as the average of this integrand over the whole set of trajectories. The result of this integration and averaging is shown in Figs. 5 and 6 for the first and second moment of the variable θ₁. The results for other variables are available in the Supplement. The black curve shows the moments of model 0 with the addition of their linear response δ〈θ₁〉 and $\mathit{\delta }\langle {\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\rangle$ to the perturbation Ψ. This curve agrees well with the green curves of the model 1 moments up to a lead time of 4–5 d, showing the efficiency of response theory. Note that in contrast to the calculation of the averages shown in Figs. 3 and 4 and computed with 1 million trajectories, we have here considered a limited subset of 10 000 trajectories of model 0 and its tangent to compute the corrections to these averages. The correction of the moments of model 1 are accurate until 4 d for both experiments. After this critical lead time, obtaining a good accuracy requires a huge increase in the number of forecasts and tangent model integrations to perform the averaging. This problem is well-known (Nicolis, 2003; Eyink et al., 2004) and is due to the appearance of fat tails in the distribution of the perturbations δy in the integrand of Eq. (49). As it can be seen in Fig. 7 for the perturbations on θ₁, the problem worsens with the increase of the lead time: initially the distributions are near-Gaussian, and fat tails appear progressively. Therefore, the number of samples of δy needed to converge to the correct mean up to a certain precision increases exponentially as the lead time increases. This problem has consequences on the method used to perform the average. Indeed, to avoid rare and unrealistic extreme responses of the system located far in the tails of the distributions, outliers above a certain threshold (set to three nondimensional units) have been removed from the averaging.

Figure 7 Histograms of the solutions δθ₁(τ) of Eq. (50) for the perturbation δy(τ) (with θ₁ being the 11th component of y) along the trajectories of model 0, for different values of the lead time τ. The solid orange curves are fits of a Gaussian distribution function to the different histograms. The fat-tail phenomenon described in Eyink et al. (2004) is apparent and becomes more prominent as the lead time increases.

Figure 8 Coefficients α and β of the postprocessing schemes of variable θ₁ and their correction using the response theory, for the experiment with a modification of the Newtonian cooling coefficient.

Figure 9 Performance of the corrections on the variable θ₁ for the experiment with the modification of the friction coefficient. (a) Mean square error (MSE) evolution between the different forecasts and their correction and reality. (b) Mean of the different trajectories (reality, model 0 and model 1) and corrected forecasts. (c) Variance of the different trajectories and corrected forecasts.

Figure 10 Performance of the corrections on the variable θ₁ for the experiment with the modification of the Newtonian cooling coefficient. (a) Mean square error (MSE) evolution between the different forecasts and their correction and reality. (b) Mean of the different trajectories (reality, model 0 and model 1) and corrected forecasts. (c) Variance of the different trajectories and corrected forecasts.

The moments obtained by the response-theory approach are used to compute new EVMOS postprocessing α and β coefficients, thanks to Eqs. (17) and (18). These corrected coefficients for variable θ₁ are shown in Fig. 8 for the experiment with a modification of the Newtonian cooling coefficient and in panels (c) and (d) of Fig. 11 for the experiment with a modification of the friction coefficient. In Figs. 9 and 10, we compare the performances of the four postprocessing schemes hence obtained: the postprocessing of model 0 (red curves) and model 1 (green curves) obtained by averaging over their trajectories (forecasts) and the postprocessing of model 1 obtained with the past model 0 forecasts (green + crosses) and with the response-theory approach (black × crosses). In panel (a) of Figs. 9 and 10, the mean square error (MSE) between the trajectories of the models and the reference trajectories is displayed by solid curves, while the MSE between both models correction and the reference is depicted by dash-dotted curves. The EVMOS postprocessing is able to partly correct the forecasts, reducing the MSE until a lead time of the order of a few times the model's Lyapunov time (the inverse of the leading Lyapunov exponent). After that, the MSE curves of the postprocessed and uncorrected forecasts converge toward a plateau corresponding to twice the variance of the reference solution (Vannitsem, 2009). Here, the statistical postprocessing corrections are indeed efficient until lead times of 4–5 d, with a skill of the corrections decreasing with the lead time. Thus the EVMOS schemes do not become better than the original models after roughly 4 d. Note also that even if the model change is small, the postprocessing using the past forecasts of model 0 (green + crosses) completely fails to correct model 1 forecasts, highlighting the need for an adaptation of the postprocessing to the model change. In contrast, the adaptation with the response-theory method (black × crosses) produces valid corrections up to 4 d later. In panels (b) and (c) of Figs. 9 and 10, the mean and variance of the corrected forecasts are compared with those of the original models. Again, the corrections obtained with response theory are efficient until 4 d for the postprocessing schemes.

In conclusion, the correction of model 1 using the response-theory EVMOS matches almost perfectly the score of the “exact” EVMOS obtained with the forecasts of model 1 (dash-dotted green curve), up to a 4 d lead time. After that lead time, the errors due to the fat tails in the response of the first moments of the statistics induce errors in the variance needed to compute the α and β coefficients (see Eqs. 17 and 18). These coefficients therefore degrade sharply after 4 d, as shown by the solid black curve in Fig. 8 and in Fig. 11c and d. This in turn induces a degradation of the response-theory postprocessing scheme. Nevertheless, this limitation of response theory is not a concern here, since after a lead time of 4 d, the EVMOS skill improvement vanishes anyway.

Figure 11 Comparison of the efficiency of the response-theory correction for different numbers m of trajectories used to average Eq. (49), for the experiment of varying the friction coefficient: (a) mean square error with reality, (b) absolute difference between the response-theory correction and the correction based on the forecast of model 1, and (c, d) postprocessing coefficients α and β. In panels (b), (c) and (d), the higher (100 000) and the lower (20) numbers are depicted respectively by a solid black line and a dashed red line. The other cases in between are depicted by dotted lines.

5 Discussion and conclusions

Statistical postprocessing techniques used to correct numerical weather predictions (NWP) require substantial past forecast and observation databases. In the case of a model change, which frequently occurs during the normal life cycle of an operational forecast model, one has to reforecast the entire database of past forecasts (Hagedorn et al., 2008; Hamill et al., 2008) to update the postprocessing coefficients and parameters. In the present work, we proposed a new methodology based on response theory to produce these new coefficients without having to reforecast. Instead, the database of past forecasts is reused to perform integrations in the tangent space of the model. It allows us to obtain the new postprocessing coefficients as modifications of the old ones. These new coefficients were shown to be accurate enough within the lead-time range for which the postprocessing corrections improve the forecast.

Figure 11 summarises the main results of this work, with the quasi-geostrophic system described in Sect. 4, using a different number m of trajectories of model 0 and its tangent model to compute the response-theory corrections. It shows that up to a lead time of 2 d, good postprocessing scheme coefficients are obtained even with a mere 20 integrations in the tangent space.

Note however that in the context of this conceptual model, good estimates of the postprocessing coefficients α and β can be obtained by simply using a small set of reforecasts. It is indeed enough to directly integrate the updated model 1, given by the non-linear Eq. (47), with only 20 trajectories. So the response-theory approach in the present case cannot really compete with the simple reforecasting method. How this can be improved in an operational context is an important question that should be addressed in the future. For instance, we can use a simplified tangent linear model to reduce the computational burden, as is often used in data assimilation (Bonavita et al., 2017). This approach could also be implemented for short-range forecasts, say from 1 to 3 d.

The response theory is efficient because the model changes are assumed to be small in comparison with the original parameterisation of the models. The method cannot improve a postprocessing scheme, but it can efficiently adapt it to a new model version. As such, the success of this method also depends on the quality of the past postprocessing scheme. There are situations where linear response theory is known to fail, but statistical tests which allow for the identification of its breakdown have been derived in Gottwald et al. (2016) and in Wormell and Gottwald (2018). In addition, the approach presented here applies only for models for which a tangent model is available. The model change itself has to be provided as an analytic function, which can in some circumstances limit the applicability of the approach.

To test this approach, we have focused on the EVMOS statistical postprocessing method, but other methods could be considered as well. The only requirement is that the outcome of the minimisation of the cost function uses averages of the systems being considered. For instance, member-by-member methods that correct both the mean square errors and the spread of the ensemble while preserving the spatial correlation (Van Schaeybroeck and Vannitsem, 2015) could be considered. These methods generally use the covariance between the model forecasts and the observations as an important piece. Response theory can also be applied here, since this covariance can be written as an average. This will be investigated in a future work, together with the applicability of the approach to parameters of probability distributions, as is often used in meteorology (Vannitsem et al., 2018).

The impact of initial-condition errors has not been addressed here, since the purpose was to demonstrate the applicability of the approach in a perfectly controlled environment. The main limiting issue of response theory in the present context is the presence of fat tails in the distribution of the perturbations δy in the tangent model. This implies that beyond a certain lead time, typically 2–3 d for the synoptic scale, the number of trajectories of the tangent model needed for the averages to converge increases exponentially. This renders the approach impractical at lead times beyond 2–3 d. This is a well-known problem, which is typically due to the trajectories passing close to the stable manifolds structuring the dynamics of chaotic systems (Eyink et al., 2004), generating an extreme response of the system to the perturbations Ψ. This is possibly due to the exacerbated sensitivity of these manifolds to the perturbation of the system. We see two possibilities to overcome this issue in the case where a long lead-time correction is needed.

• First, as suggested by Eyink et al. (2004), the problem should be studied in other systems. It might be resolved by itself in other systems. Indeed, in very large atmospheric systems, the encounter of such manifolds might become rare. This could be related to the chaotic hypothesis (Gallavotti and Cohen, 1995a, b) which states that large systems can be considered to behave like Axiom-A hyperbolic systems for the physical quantities of interest, and thus Ruelle response theory (Ruelle, 2009) might get better as the dimensionality of a system increases. This hypothesis would be interesting to test in current state-of-the-art NWP systems.

• Secondly, another avenue would be to adapt the techniques based on the covariant Lyapunov vectors (CLVs) or on unstable periodic orbits (UPOs) to non-stationary dynamics. These techniques were recently introduced (Wang, 2013; Ni and Wang, 2017; Ni, 2019; Lasagna, 2019; Lasagna et al., 2019) to deal with stationary responses of chaotic systems, i.e. the response of a system that lies on its attractor.

The CLVs methods mentioned focus on finding an adjoint representation (Eyink et al., 2004) of the response, while in the present work the approach is based on forward integrations (direct method). The adjoint representation allows for easily changing the perturbation function Ψ for a fixed observable A, while the direct method enables the consideration of different observables while keeping the perturbation function fixed. The adjoint representation, however, requires one to integrate the tangent model backward in time. Therefore, its accuracy depends on the absolute value of the smallest Lyapunov exponent of the system, which might render its results less well than the direct forward representation.

In conclusion, the response-theory approach developed here is an effective method to deal with the problem of the impact of model change on the postprocessing scheme. Its main advantage is to be computed on the past model version and does not require reforecasts of the full model. Its operational implementation, however, is still an open question that should be addressed in the future.

Appendix A: Non-stationary response theory

We consider a perturbed autonomous dynamical system

$$\begin{array}{}\text{(A1)}& \dot{\widehat{\mathit{y}}}=\mathit{F}\left(\widehat{\mathit{y}}\right)+\mathbf{\Psi }\left(\widehat{\mathit{y}}\right)=\widehat{\mathit{F}}\left(\widehat{\mathit{y}}\right)\end{array}$$

with a prescribed distribution of initial conditions ρ₀. For the unperturbed system

$$\begin{array}{}\text{(A2)}& \dot{\mathit{y}}=\mathit{F}\left(\mathit{y}\right),\end{array}$$

an observable A has the average at time τ

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)A\left({\mathit{f}}^{\mathit{\tau }}\right({\mathit{y}}_{\mathrm{0}}\left)\right)\\ & =\int \mathrm{d}\mathit{y}{\mathit{\rho }}_{\mathit{\tau }}\left(\mathit{y}\right)A\left(\mathit{y}\right),\end{array}\end{array}$$

where f^τ is the flow of the unperturbed system given by Eq. (A2) and where ρ_τ is the distribution obtained by propagating the initial distribution ρ₀ with the Liouville equation (Gaspard, 2005). In this section, the variation of this average due to the presence of the perturbation is evaluated as

$$\begin{array}{}\text{(A4)}& \langle A\left(\mathit{\tau }\right){\rangle }_{\widehat{\mathit{y}}}=\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}+\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}+{\mathit{\delta }}^{\mathrm{2}}\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}+\mathrm{\dots }.\end{array}$$

In other words, we compute the average of A in the system defined by Eq. (A1)

$$\begin{array}{}\text{(A5)}& \langle A\left(\mathit{\tau }\right){\rangle }_{\widehat{\mathit{y}}}=\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\right({\mathit{y}}_{\mathrm{0}}\left)\right)\end{array}$$

as a perturbation of the average given by Eq. (A3) for the unperturbed system defined by Eq. (A2). Here, ${\widehat{\mathit{f}}}^{\mathit{\tau }}$ is the flow of the perturbed system defined by Eq. (A1). In the following, we will derive these corrections thanks to a Kubo-type perturbative expansion (Lucarini, 2008) that amounts to constructing a Dyson series in the interaction picture framework, where the perturbation is seen as an interaction Hamiltonian (Wouters and Lucarini, 2012). We start by considering the time evolution of the observable A in Eq. (A1):

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{\tau }}}A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)=\left({\mathcal{L}}_{\mathrm{0}}+{\mathcal{L}}_{\mathrm{1}}\right)A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)\end{array}$$

with the operators

$$\begin{array}{}\text{(A7)}& \left\{\begin{array}{lcl}{\mathcal{L}}_{\mathrm{0}}A\left(\mathit{y}\right)& =& \mathit{F}(\mathit{y}{)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{\mathit{y}}A\\ {\mathcal{L}}_{\mathrm{1}}A\left(\mathit{y}\right)& =& \mathbf{\Psi }(\mathit{y}{)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{\mathit{y}}A\end{array}\right.\end{array}$$

and define an interaction observable as

$$\begin{array}{}\text{(A8)}& A_I(\mathit{\tau },{\mathit{y}}_{\mathrm{0}})={\mathrm{\Pi }}_{\mathrm{0}}(-\mathit{\tau })A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)\end{array}$$

with Π₀(τ)=exp (ℒ₀ τ). It is easy to show that the interaction observable satisfies the differential equation:

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{\tau }}}{A}_I(\mathit{\tau },{\mathit{y}}_{\mathrm{0}})={\mathcal{L}}_I\left(\mathit{\tau }\right)A_I(\mathit{\tau },{\mathit{y}}_{\mathrm{0}})\end{array}$$

with the interaction operator ${\mathcal{L}}_I\left(\mathit{\tau }\right)={\mathrm{\Pi }}_{\mathrm{0}}(-\mathit{\tau }){\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}\left(\mathit{\tau }\right)$. The solution to this equation is

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{A}_I(\mathit{\tau },{\mathit{y}}_{\mathrm{0}})& =A_I(\mathrm{0},{\mathit{y}}_{\mathrm{0}})+\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}{\mathcal{L}}_I\left(s_{\mathrm{1}}\right)A_I(s_{\mathrm{1}},{\mathit{y}}_{\mathrm{0}})\\ & =A\left({\mathit{y}}_{\mathrm{0}}\right)+\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}{\mathcal{L}}_I\left(s_{\mathrm{1}}\right)A_I(s_{\mathrm{1}},{\mathit{y}}_{\mathrm{0}}),\end{array}\end{array}$$

which can be rewritten as

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)& ={\mathrm{\Pi }}_{\mathrm{0}}\left(\mathit{\tau }\right)A\left({\mathit{y}}_{\mathrm{0}}\right)\\ & +\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}(\mathit{\tau }-s_{\mathrm{1}}){\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}\left(s_{\mathrm{1}}\right)A_I(s_{\mathrm{1}},{\mathit{y}}_{\mathrm{0}}).\end{array}\end{array}$$

Iteratively replacing the interaction observable by Eq. (A10) finally leads to the Dyson series:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}A\left({\widehat{\mathit{f}}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)& ={\mathrm{\Pi }}_{\mathrm{0}}\left(\mathit{\tau }\right)A\left({\mathit{y}}_{\mathrm{0}}\right)\\ & +\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}(\mathit{\tau }-s_{\mathrm{1}}){\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}\left(s_{\mathrm{1}}\right)A\left({\mathit{y}}_{\mathrm{0}}\right)\\ & +\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}\underset{\mathrm{0}}{\overset{s_{\mathrm{1}}}{\int }}\mathrm{d}{s}_{\mathrm{2}}{\mathrm{\Pi }}_{\mathrm{0}}(\mathit{\tau }-s_{\mathrm{1}}){\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}(s_{\mathrm{1}}-s_{\mathrm{2}})\\ & {\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}\left(s_{\mathrm{2}}\right)A\left({\mathit{y}}_{\mathrm{0}}\right)+\mathrm{\dots }.\end{array}\end{array}$$

Using the definitions in Eqs. (A3) and (A5), as well as the fact that

$$\begin{array}{}\text{(A13)}& g\left({\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)={\mathrm{\Pi }}_{\mathrm{0}}\left(\mathit{\tau }\right)g\left({\mathit{y}}_{\mathrm{0}}\right)\end{array}$$

for any smooth function g, we get finally a formula for the perturbations in Eq. (A4):

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}\langle A\left(\mathit{\tau }\right){\rangle }_{\widehat{\mathit{y}}}& =\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}+\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right){\mathrm{\Pi }}_{\mathrm{0}}(\mathit{\tau }-s_{\mathrm{1}})\\ & {\mathcal{L}}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{0}}\left(s_{\mathrm{1}}\right)A\left({\mathit{y}}_{\mathrm{0}}\right)+\mathrm{\dots }.\end{array}\end{array}$$

We will now focus on the first term of this expansion, but the subsequent orders of the response can be treated in the same way. We thus have

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{s}_{\mathrm{1}}\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathbf{\Psi }{\left({\mathit{f}}^{\mathit{\tau }-s_{\mathrm{1}}}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\\ & \cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }-s_{\mathrm{1}}}\left({\mathit{y}}_{\mathrm{0}}\right)}A\left({\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right),\end{array}\end{array}$$

which with the change of variable $s_{\mathrm{1}}\to t-{\mathit{\tau }}^{\prime }$ can be rewritten as

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathbf{\Psi }{\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\\ & \cdot {\mathbf{\nabla }}_{{\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)}A\left({\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)\end{array}\end{array}$$

and then

$$\begin{array}{}\text{(A17)}& {\displaystyle }\begin{array}{rl}\mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathbf{\Psi }{\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\\ & \cdot {\left({\displaystyle \frac{\partial {\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}{\partial {\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)}}\right)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}A\end{array}\text{(A18)}& {\displaystyle }\begin{array}{rl}& =\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathbf{\Psi }{\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\\ & \cdot \mathbf{M}{\left(\mathit{\tau }-{\mathit{\tau }}^{\prime },{\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathrm{0}}\right)\right)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}A.\end{array}\end{array}$$

M is the fundamental matrix (Gaspard, 2005; Nicolis, 2016) of the homogeneous part of the linear differential equation

$$\begin{array}{}\text{(A19)}& \dot{\mathit{\delta }\mathit{y}}={\mathbf{\nabla }}_{\mathit{y}}\mathit{F}\cdot \mathit{\delta }\mathit{y}+\mathbf{\Psi }\left(\mathit{y}\right),\end{array}$$

where y is solution of Eq. (A2) with initial condition y₀, and we have the definition

$$\begin{array}{}\text{(A20)}& \mathbf{M}(t,\mathit{y})={\displaystyle \frac{\partial {\mathit{f}}^t\left(\mathit{y}\right)}{\partial \mathit{y}}}.\end{array}$$

Equation (A19) is the linearised approximation of Eq. (A1):

$$\begin{array}{}\text{(A21)}& \dot{\mathit{y}}+\dot{\mathit{\delta }\mathit{y}}=\mathit{F}(\mathit{y}+\mathit{\delta }\mathit{y})+\mathbf{\Psi }(\mathit{y}+\mathit{\delta }\mathit{y}),\end{array}$$

which provides a tool to estimate Eq. (A18). Indeed, since the solution of Eq. (A19) can be written as

$$\begin{array}{}\text{(A22)}& \mathit{\delta }\mathit{y}\left(\mathit{\tau }\right)=\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{d}{\mathit{\tau }}^{\prime }\mathbf{M}\left(\mathit{\tau }-{\mathit{\tau }}^{\prime },{\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathbf{0}}\right)\right)\cdot \mathbf{\Psi }\left({\mathit{f}}^{{\mathit{\tau }}^{\prime }}\left({\mathit{y}}_{\mathbf{0}}\right)\right),\end{array}$$

we can write the first-order variation of the average of the observable A in terms of these solutions:

$$\begin{array}{}\text{(A23)}& \mathit{\delta }\langle A\left(\mathit{\tau }\right){\rangle }_{\mathit{y}}=\int \mathrm{d}{\mathit{y}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}\left({\mathit{y}}_{\mathrm{0}}\right)\mathit{\delta }\mathit{y}(\mathit{\tau }{)}^{\mathrm{T}}\cdot {\mathbf{\nabla }}_{{\mathit{f}}^{\mathit{\tau }}\left({\mathit{y}}_{\mathrm{0}}\right)}A.\end{array}$$

The interpretation of this equation is that a specific averaging of an observable over the trajectories of the linear approximation given by Eq. (A19) of the perturbed Eq. (A2) provides the first-order response of the observable. It is the main result used to compute the new postprocessing scheme in the present work. It is explained in detail in Sects. 3.3 and 4.2.

Appendix B: The quasi-geostrophic model equations

The ordinary differential equations of the model are given by

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\psi }}}_i& =-a_{i,i}^{-\mathrm{1}}\sum _{j,m=\mathrm{1}}^{n_{\mathrm{a}}}{b}_{i,j,m}\left({\mathit{\psi }}_j{\mathit{\psi }}_m+{\mathit{\theta }}_j{\mathit{\theta }}_m\right)\\ & -{\displaystyle \frac{a_{i,i}^{-\mathrm{1}}}{\mathrm{2}}}\sum _{j,m=\mathrm{1}}^{n_{\mathrm{a}}}{g}_{i,j,m}{h}_m\left({\mathit{\psi }}_j-{\mathit{\theta }}_j\right)\\ & -\mathit{\beta }{a}_{i,i}^{-\mathrm{1}}\sum _{j=\mathrm{1}}^{n_{\mathrm{a}}}{c}_{i,j}{\mathit{\psi }}_j-{\displaystyle \frac{k_d}{\mathrm{2}}}\left({\mathit{\psi }}_i-{\mathit{\theta }}_i\right),\end{array}\text{(B2)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\theta }}}_i& =-a_{i,i}^{-\mathrm{1}}\sum _{j,m=\mathrm{1}}^{n_{\mathrm{a}}}{b}_{i,j,m}({\mathit{\psi }}_j{\mathit{\theta }}_m+{\mathit{\theta }}_j{\mathit{\psi }}_m)\\ & +{\displaystyle \frac{a_{i,i}^{-\mathrm{1}}}{\mathrm{2}}}\sum _{j,m=\mathrm{1}}^{n_{\mathrm{a}}}{g}_{i,j,m}{h}_m\left({\mathit{\psi }}_j-{\mathit{\theta }}_j\right)\\ & -\mathit{\beta }{a}_{i,i}^{-\mathrm{1}}\sum _{j=\mathrm{1}}^{n_{\mathrm{a}}}{c}_{i,j}{\mathit{\theta }}_j+{\displaystyle \frac{k_d}{\mathrm{2}}}\left({\mathit{\psi }}_i-{\mathit{\theta }}_i\right)\\ & -\mathrm{2}{k}_d^{\prime }{\mathit{\theta }}_i+a_{i,i}^{-\mathrm{1}}{\mathit{\omega }}_i,\end{array}\text{(B3)}& {\displaystyle }{\dot{\mathit{\theta }}}_i=-\sum _{j,m=\mathrm{1}}^{n_{\mathrm{a}}}{g}_{i,j,m}{\mathit{\psi }}_j{\mathit{\theta }}_m+{\displaystyle \frac{\mathit{\sigma }}{\mathrm{2}}}{\mathit{\omega }}_i+h_d\left({\mathit{\theta }}_i^{*}-{\mathit{\theta }}_i\right),\end{array}$$

where nondimensional parameters values and description can be found in Table 1 and Sect. 4. β is the meridional gradient of Coriolis parameter which has the nondimensional value of 0.21 at 50^∘ of latitude (Reinhold and Pierrehumbert, 1982; Cehelsky and Tung, 1987). The vertical velocity ω_i can be eliminated, leading to Eqs. (B2) and (B3) being reduced to a single equation for θ_i. The parameter σ is the nondimensional static stability of the atmosphere set typically to 0.2. The coefficients a_i, j, $g_{i,j,m}$, $b_{i,j,m}$ and c_i, j are the inner products of the Fourier modes F_i defined in Sect. 4:

$$\begin{array}{}\text{(B4)}& {\displaystyle }\begin{array}{rl}& a_{i,j}=\\ & {\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/n}{\int }}{F}_i(x,y){\mathrm{\nabla }}^{\mathrm{2}}{F}_j(x,y)\mathrm{d}x\mathrm{d}y=-{\mathit{\delta }}_{ij}{a}_i^{\mathrm{2}},\end{array}\text{(B5)}& {\displaystyle }\begin{array}{rl}& g_{i,j,m}=\\ & {\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/n}{\int }}{F}_i(x,y)J\left(F_j(x,y),F_m(x,y)\right)\mathrm{d}x\mathrm{d}y,\end{array}\text{(B6)}& {\displaystyle }\begin{array}{rl}& b_{i,j,m}=\\ & {\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/n}{\int }}{F}_i(x,y)J\left(F_j(x,y),{\mathrm{\nabla }}^{\mathrm{2}}{F}_m(x,y)\right)\mathrm{d}x\mathrm{d}y,\end{array}\end{array}$$

$$\begin{array}{}\text{(B7)}& c_{i,j}={\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/n}{\int }}{F}_i(x,y){\displaystyle \frac{\partial }{\partial x}}{F}_j(x,y)\mathrm{d}x\mathrm{d}y,\end{array}$$

where the coefficients a_i are given by Eq. (41) and where J is the Jacobian present in the advection terms defined as

$$\begin{array}{}\text{(B8)}& J(S,G)={\displaystyle \frac{\partial S}{\partial x}}{\displaystyle \frac{\partial G}{\partial y}}-{\displaystyle \frac{\partial S}{\partial y}}{\displaystyle \frac{\partial G}{\partial x}}.\end{array}$$