Correcting for model changes in statistical postprocessing – an approach based on response theory
- 1Institut Royal Météorologique de Belgique, Avenue Circulaire, 3, 1180 Brussels, Belgium
- 2European Meteorological Network (EUMETNET), Avenue Circulaire, 3, 1180 Brussels, Belgium
Correspondence: Jonathan Demaeyer (firstname.lastname@example.org)
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.
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.
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:
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.1 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
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
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
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
it induces a perturbation of the observable's average, which at the first order is given by2
where fτ is the flow of the system defined by Eq. (2) mapping an initial condition y0 to the system's state at time τ as y(τ)=fτ(y0), 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
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
which is the solution of the homogeneous equation . Using the chain rule, the response defined by Eq. (6) is rewritten in terms of the perturbation alone:
where the causality of the perturbation acting on the system and perturbing the averaged observable appears (Lucarini, 2008), since . 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 model3 how response theory allows us to tackle the issue.
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:
where ξx and ξy are Gaussian white-noise processes such that
These are therefore uncorrelated Ornstein–Uhlenbeck processes with noise amplitudes Qx and Qy.
We then consider a change Ψy of model y(τ), possibly improving or degrading the forecast performances as
with and . 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 . 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 ).
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: . 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 yn and observations xn, the correction of the univariate EVMOS postprocessing of variable x from a new forecast y(τ) is provided by the linear regression
The coefficients α and β are obtained by minimising the functional
and are thus given by the following equations:
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
where we note that the model is initialised with the same initial conditions as reality:
Similarly, we get the same kind of results for model , after model change Ψy:
The ratio between the parameters β is given by
For , we note that this ratio tends to
and the difference between the biases α of the two models is approximatively given by
Let us now assume that model change Ψy can be considered as a perturbation of the initial model y. Using response theory, the averages and can be estimated using the initial model y instead of the perturbed model . In turn, these new estimated averages give us the new postprocessing scheme coefficients and . 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 , 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
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 forecasts (the left-hand side) based solely on the average over the model y forecasts. The 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 observations4 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)
where ρ0 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, ρ0 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
This maps an initial condition x0 of model y to the state fτ(x0) 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.
Rearranging these two terms, we also get the following expression for the variation of the variance given by Eq. (24):
Note that the variation given by Eq. (35) corresponds exactly to the difference between the average of the two models . On the other hand the variation given by Eq. (37) lacks the term of order κ2 involving δQ that appears in the difference between and given respectively by Eqs. (27) and (24). Instead, another term of order κ2 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 expression5 (Lucarini, 2012):
Applying this to the first moment of the y models directly yields
On the other hand, integrating the stochastic integrals present in this expression for the moment 〈y(τ)2〉 gives
which corrects the κ2δK2 term in Eq. (37) and makes the response theory up to order 2 exactly match the difference between and , 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 and , the approximate averages 〈y(τ)〉+δ〈y(τ)〉y and can be used to compute and . 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 κ2, 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.
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 , with as the aspect ratio between its meridional and zonal extents Ly and Lx. 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 ( at ). The fields are expanded in Fourier modes respecting these boundary conditions:
with eigenvalues . We have thus the following decomposition
where na 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
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 kd=0.1, the friction at the bottom surface and the aspect ratio of the domain n=1.3. The β plane lies at midlatitudes (50∘) and the Coriolis parameter f0 is set accordingly.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
In the present work, the parameter hd, the Newtonian cooling coefficient is fixed to 0.3 instead of the value found in (which is hd=0.045). Two additional fields have to be specified on the domain: , 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 and hi then allow for writing these fields as sums of weighted eigenfunctions:
In the present case, we consider that the only non-zero coefficients are and h2=0.4, meaning that the radiative equilibrium profile is given by the zonally varying function 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 is larger than the one chosen in (which is ) to increase the chaotic variability in the system. Trajectories of variables θ1 and ψ2 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.
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 ψ2 and ψ3 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.
4.1 Postprocessing experiments
The model described above with 10 modes (na=10) is used, and two different postprocessing experiments are performed, one involving the Newtonian cooling parameter hd and another involving the friction parameter kd 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 hd or kd 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.
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 and 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 ,
an observable A with average 〈A(τ)〉y at the lead time τ for the system
has a first-order response of
where fτ is the flow of the unperturbed system given by Eq. (48), ρ0 is the distribution of initial conditions and δy(τ) is the solution of the equation , which can be approximated at the first order by the following linear inhomogeneous differential equation
where y(τ) is the solution of the unperturbed Eq. (48) with the initial condition y(0)=y0, 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.
The response-theory approximations of the averages of model (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 θ1. 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 δ〈θ1〉 and 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 θ1, 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.
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 θ1 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.
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.
The quasi-geostrophic model used is called QGS and was obtained by adapting the Python code of the MAOOAM ocean–atmosphere model (De Cruz et al., 2016), following the model description in Cehelsky and Tung (1987). It was recently released on Zenodo (Demaeyer and De Cruz, 2020) and is also available at https://github.com/Climdyn/qgs (last access: 18 May 2020). The additional notebooks computing the response to model changes and generating the figures are also provided in the Supplement. They have been released on Zenodo as well (Demaeyer, 2020) and are available at https://github.com/jodemaey/Postprocessing_and_response_theory_notebooks (last access: 18 May 2020).
We consider a perturbed autonomous dynamical system
with a prescribed distribution of initial conditions ρ0. For the unperturbed system
an observable A has the average at time τ
where fτ is the flow of the unperturbed system given by Eq. (A2) and where ρτ is the distribution obtained by propagating the initial distribution ρ0 with the Liouville equation (Gaspard, 2005). In this section, the variation of this average due to the presence of the perturbation is evaluated as
In other words, we compute the average of A in the system defined by Eq. (A1)
as a perturbation of the average given by Eq. (A3) for the unperturbed system defined by Eq. (A2). Here, 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):
with the operators
and define an interaction observable as
with Π0(τ)=exp (ℒ0 τ). It is easy to show that the interaction observable satisfies the differential equation:
with the interaction operator . The solution to this equation is
which can be rewritten as
Iteratively replacing the interaction observable by Eq. (A10) finally leads to the Dyson series:
for any smooth function g, we get finally a formula for the perturbations in Eq. (A4):
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
which with the change of variable can be rewritten as
where y is solution of Eq. (A2) with initial condition y0, and we have the definition
we can write the first-order variation of the average of the observable A in terms of these solutions:
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.
The ordinary differential equations of the model are given by
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 ai, j, , and ci, j are the inner products of the Fourier modes Fi defined in Sect. 4:
where the coefficients ai are given by Eq. (41) and where J is the Jacobian present in the advection terms defined as
The supplement related to this article is available online at: https://doi.org/10.5194/npg-27-307-2020-supplement.
JD and SV developed the idea of using the response theory in the context of postprocessing, together with the overall experimental setup. JD made the analytical and numerical computations. Both authors contributed to the writing of the paper.
Stéphane Vannitsem is a member of the editorial board of the journal. Jonathan Demaeyer declares that he has no conflict of interest.
This article is part of the special issue “Advances in post-processing and blending of deterministic and ensemble forecasts”. It is not associated with a conference.
The authors warmly thank Lesley De Cruz for her suggestions throughout the paper. They also thank Michaël Zamo and the anonymous reviewer for their comments and suggested improvements.
This research has been supported by EUMETNET (Postprocessing module of the NWP Cooperation Programme).
This paper was edited by Maxime Taillardat and reviewed by Michaël Zamo and one anonymous referee.
Bódai, T., Lucarini, V., and Lunkeit, F.: Can we use linear response theory to assess geoengineering strategies?, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30, 023124, https://doi.org/10.1063/1.5122255, 2020. a
Bonavita, M., Trémolet, Y., Holm, E., Lang, S. T. K., Chrust, M., Janisková, M., Lopez, P., Laloyaux, P., De Rosnay, P., Fisher, M., Hamrud, M., and English, S.: A strategy for data assimilation, European Centre for Medium Range Weather Forecasts, available at: https://www.ecmwf.int/sites/default/files/elibrary/2017/17179-strategy-data-assimilation.pdf (last access: 20 May 2020), 2017. a
Charney, J. G. and Straus, D. M.: Form-drag instability, multiple equilibria and propagating planetary waves in baroclinic, orographically forced, planetary wave systems, J. Atmos. Sci., 37, 1157–1176, 1980. a
Demaeyer, J. and Vannitsem, S.: Comparison of stochastic parameterizations in the framework of a coupled ocean–atmosphere model, Nonlin. Processes Geophys., 25, 605–631, https://doi.org/10.5194/npg-25-605-2018, 2018. a
Eyink, G., Haine, T., and Lea, D.: Ruelle's linear response formula, ensemble adjoint schemes and Lévy flights, Nonlinearity, 17, 1867, https://doi.org/10.1088/0951-7715/17/5/016, 2004. a, b, c, d, e, f
Glahn, B., Peroutka, M., Wiedenfeld, J., Wagner, J., Zylstra, G., Schuknecht, B., and Jackson, B.: MOS uncertainty estimates in an ensemble framework, Mon. Weather Rev., 137, 246–268, 2009. a
Glahn, H. R. and Lowry, D. A.: The use of model output statistics (MOS) in objective weather forecasting, J. Appl. Meteorol., 11, 1203–1211, 1972. a
Hagedorn, R., Hamill, T. M., and Whitaker, J. S.: Probabilistic forecast calibration using ECMWF and GFS ensemble reforecasts. Part I: Two-meter temperatures, Mon. Weather Rev., 136, 2608–2619, 2008. a
Hamill, T. M.: Practical aspects of statistical postprocessing, in: Statistical Postprocessing of Ensemble Forecasts, edited by: Vannitsem, S., Wilks, D. S., and Messner, J. W., 187–217, Elsevier, Amsterdam, the Netherlands, 2018. a, b
Hamill, T. M., Hagedorn, R., and Whitaker, J. S.: Probabilistic forecast calibration using ECMWF and GFS ensemble reforecasts. Part II: Precipitation, Mon. Weather Rev., 136, 2620–2632, 2008. a
Hamill, T. M., Bates, G. T., Whitaker, J. S., Murray, D. R., Fiorino, M., Galarneau Jr, T. J., Zhu, Y., and Lapenta, W.: NOAA's second-generation global medium-range ensemble reforecast dataset, B. Am. Meteorol. Soc., 94, 1553–1565, 2013. a
Johnson, C. and Bowler, N.: On the reliability and calibration of ensemble forecasts, Mon. Weather Rev., 137, 1717–1720, 2009. a
Lang, M. N., Lerch, S., Mayr, G. J., Simon, T., Stauffer, R., and Zeileis, A.: Remember the past: a comparison of time-adaptive training schemes for non-homogeneous regression, Nonlin. Processes Geophys., 27, 23–34, https://doi.org/10.5194/npg-27-23-2020, 2020. a
Lasagna, D.: Sensitivity and stability of long periodic orbits of chaotic systems, ArXiv [Preprint], arXiv:1910.06706, 2019. a
Lasagna, D., Sharma, A., and Meyers, J.: Periodic shadowing sensitivity analysis of chaotic systems, J. Comput. Phys., 391, 119–141, 2019. a
Lembo, V., Lucarini, V., and Ragone, F.: Beyond Forcing Scenarios: Predicting Climate Change through Response Operators in a Coupled General Circulation Model, ArXiv [Preprint], arXiv:1912.03996, 2019. a
Ni, A.: Hyperbolicity, shadowing directions and sensitivity analysis of a turbulent three-dimensional flow, J. Fluid Mech., 863, 644–669, 2019. a
Ni, A. and Wang, Q.: Sensitivity analysis on chaotic dynamical systems by Non-Intrusive Least Squares Shadowing (NILSS), J. Comput. Phys., 347, 56–77, 2017. a
Nicolis, C.: Dynamics of model error: Some generic features, J. Atmos. Sci., 60, 2208–2218, 2003. a
Nicolis, C., Perdigao, R. A., and Vannitsem, S.: Dynamics of prediction errors under the combined effect of initial condition and model errors, J. Atmos. Sci., 66, 766–778, 2009. a
Scheuerer, M. and Hamill, T. M.: Statistical postprocessing of ensemble precipitation forecasts by fitting censored, shifted gamma distributions, Mon. Weather Rev., 143, 4578–4596, 2015. a
Vannitsem, S.: Predictability of large-scale atmospheric motions: Lyapunov exponents and error dynamics, Chaos: An Interdisciplinary Journal of Nonlinear Science, 27, 032101, https://doi.org/10.1063/1.4979042, 2017. a
Vannitsem, S. and Nicolis, C.: Dynamical properties of model output statistics forecasts, Mon. Weather Rev., 136, 405–419, 2008. a
Vannitsem, S., Wilks, D. S., and Messner, J.: Statistical postprocessing of ensemble forecasts, Elsevier, Amsterdam, the Netherlands, 2018. a
Vissio, G. and Lucarini, V.: A proof of concept for scale-adaptive parametrizations: the case of the Lorenz'96 model, Q. J. Roy. Meteor. Soc., 144, 63–75, 2018. a
Wilks, D. S.: Statistical methods in the atmospheric sciences, vol. 100, Academic press, Oxford, UK, 2011. a
Young, L.-S.: What are SRB measures, and which dynamical systems have them?, J. Stat. Phys., 108, 733–754, 2002. a
When taking the gradient of a function A, the notation ∇yA means taking the gradient at the point y, i.e. evaluating ∇yA(y).
Here we consider that the observation are perfectly assimilated in the models and that there is no observation errors. However in operational setups, such errors are of course to be taken into account.