Correcting for Model Changes in Statistical Post-Processing -- An approach based on Response Theory

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

sults for the response can be obtained at any order. The correction of the model observables and the post-processing parameters due to the model change only requires the response-theory corrections up to the second order.
In section 3, a more complex case is considered with a toy model of atmospheric variability in the form of a 2-layer quasigeostrophic model with an orography. We compute the linear response of the predictors of the post-processing 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 post-processing scheme up to a lead time of 3 days, 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 post-processing systems, as well as new research avenues to be explored.
We then consider a change Ψ y of the model y(τ ), possibly improving or degrading the forecast performances: where Ψ y (τ ) = −κ δK + δQ ξ y (τ ) with δK = K y − K x and δ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 if κ = 1, the correction is perfect.
We have thus the reality x(τ ) and two different models of it: y(τ ) andŷ(τ ). We now want to evaluate the difference between a post-processing scheme constructed before the model change (with the past forecasts of the model y(τ )), and one constructed after it (with the past forecasts of modelŷ(τ )).

The post-processing method
We now consider a forecast situation where the model y is initialised at the time τ = 0 with a perfect observation of the reality: y(0) = x(0) = x 0 . We use the Error-in-Variables Model Output Statistics (EVMOS) post-processing scheme (Vannitsem, 2009) to correct the forecasts of the model y based on these initial conditions. In this context, given K past forecasts y k and observations x k , the correction of the univariate EVMOS post-processing 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 equations: β(τ ) = σ 2 x (τ ) σ 2 y (τ ) where The averages · are taken over an ensemble of past forecasts and observations. This approach has been developed to allow for 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.
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 infinite number of realisations of the Ornstein-Uhlenbeck processes, as if we had an infinite ensemble of past forecasts.

Averaging the Ornstein-Uhlenbeck processes
For the reality x and the model y, we directly get the averages (Gardiner, 2009) and where we note that the model is initialised with the same initial conditions as the reality: We get the post-processing coefficients before the model change α(τ ) and β(τ ) by inserting these expressions in the equations (7) and (8).
Similarly, we get the same kind of results for the modelŷ, after the model change Ψ y : and we also obtain the post-processing coefficients after the model changeα(τ ) andβ(τ ) (see also the analysis in Vannitsem (2011)). We can also compute the variation of the bias α: The ratio between the parameters β is given bŷ For τ max(1/λ x , 1/λ y ), we note that this ratio tends tô 20) and the difference between the biases α of the two models is approximatively given by: Let us now assume that the model change Ψ y can be considered as a perturbation of the initial model y. Using response theory, the averages ŷ and σ 2 y can be estimated using the initial model y instead of the perturbed modelŷ. In turn, the use of these new estimated averages allows for computing the new post-processing scheme coefficientsα andβ. We now detail the results obtained by using this method.

Response theory
After the model change, the forecasts are provided by the modelŷ and their time-evolution is given by Eq. (3). This model can be seen as a perturbation of the model y by the term Ψ y given by Eq. (4). In such 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 to obtain the average over the modelŷ forecasts (the left-hand side) based solely on the average over model the y forecasts. Theŷ model forecasts are therefore not required to estimate the new post-processing scheme.
The observables depend on the lead time τ of the forecast, as do the parameters α and β which determine the post-processing correction for every lead time. This reflects the fact that the post-processing problem is typically a non-stationary initial value problem, since the initial conditions of the model Eqs. (2) and (3) are typically not chosen on their respective model attractor, but rather as observations 1 of the reality (1). As a consequence, the model averages (22) 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 (4) is small, then the first order is given by (see Appendix A): 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. (15), in the post-processing framework, it is taken as the stationary/invariant distribution of the reality. It is also assumed that there is no interference due to initial condition errors in the problem. Eq. (23) 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. (23) is the stochastic flow: This maps an initial condition x 0 of the model y to the state f τ (x 0 ) of a realisation of this model at the later lead time τ .
The principle of causality is thus implicit in Eq. (23), which estimates the impact of the perturbation Ψ y on the subsequent perturbed model time-evolution by developing around the unperturbed model y trajectories. 1 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.
Evaluating Eq. (23) and its stochastic integrals gives us the variation of the averages y(τ ) and y(τ ) 2 to the perturbation Ψ y : Rearranging these two terms, we also get the following expression for the variation of the variance (14): Note that the variation (25) corresponds to the exact difference between the average of the two models ŷ(τ ) − y(τ ) . On the other hand the variation given by Eq. (27) lacks the term of order κ 2 involving δQ that appears in the exact difference σ 2 y (τ ) − σ 2 y (τ ) given by Eqs. (14) and (17). Instead, another term of order κ 2 and involving δK is present, indicating that we need to consider the higher-order term of response theory Ruelle (1998b) to correct it. The second-order term is given by the expression (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 δK 2 term in Eq. (27) and makes the response theory up to order 2 exactly match the difference σ 2 y (τ ) − σ 2 y (τ ), for every lead time τ . In fact, the subsequent orders of the response vanish due to the linearity of the simple Ornstein-Uhlenbeck models, which allows 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 post-processing parameters after the model change based on the forecasts of the initial model. Indeed, instead of the averages ŷ(τ ) and σ 2 y (τ ) , the approximate averages y(τ ) + δ y(τ ) y and σ 2 y (τ ) + δ σ 2 y (τ ) y + δ 2 σ 2 y (τ ) 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 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 supplementary material.
In order to investigate this research avenue on a case closer to those encountered in reality, we will now consider the application of post-processing and response theory to a low-order atmospheric model displaying chaos.
A 2-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 adimensionalised coordinates are denoted x and y. 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. These fields are expanded in Fourier modes: with eigenvalues a 2 1 = 1, a 2 2 = a 2 3 = 1 + n 2 , a 2 4 = 4, . . . , where n is the aspect ratio of the x-y domain. We have thus the following decomposition 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 : that can be solved with usual numerical integrators. All variables are adimensionalized.
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 adopt their notations and main adimensionalized parameters values: the friction at the interface between the two layers k = 0.05, the friction at the bottom surface k = 0.01, and the aspect ratio of the domain n = 1.3. The β-plane lies at mid-latitude (45°) and the Coriolis parameter f 0 is set accordingly.
In the present work, the parameter h, the Newtonian cooling coefficient is fixed to 0.1 instead of the value found in Reinhold and Pierrehumbert (which is h = 0.045). Two additional fields have to be specified on the domain: θ (x, y), the radiative equilibrium temperature field, and (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 θ i and i then allow for writing these fields as sums of weighted eigenfunctions: In the present case, we consider that the only non-zero coefficients are θ 1 = 0.2 and 2 = 0.2, meaning that the radiative equilibrium profile is given by the zonally varying function √ 2 cos(y) 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 θ 1 is larger than the one chosen in Reinhold and Pierrehumbert (which is θ 1 = 0.1) to increase the chaotic variability in the system. Trajectories of variables θ 1 and ψ 4 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. 2(b): one characterised by a zonal circulation (see Fig. 1(c)), and another characterised by a blocking situation (see Fig. 1(d)). This is different from the case considered in Reinhold and Pierrehumbert (1982), where two blocking regimes coexist with the zonal regime.

Post-processing experiments
The model described above with 10 modes (n a = 10) is used and two different post-processing experiments are performed, one involving the Newtonian cooling parameter h and another involving the friction parameter k between the two atmospheric layers. The parameter values detailed above correspond to the long-term reference (i.e. the reality). A first model is defined (model 0) which is a copy of the 2-layer quasi-geostrophic model defining the reality, but the parameters h or k 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 Section 2, 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 coming from the reference forecasts, as a function of the lead time τ . We have used a set of one million trajectories of each system to compute these averages.
In the framework of the EVMOS post-processing scheme, the predictors and the predictands are the same nominal variable and no other predictors are used. In both experiments considered, the post-processing 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.

Response theory and the tangent linear model
Let us consider again the response theory described in Section 2.3, but in the multivariate deterministic case. In the postprocessing framework, models 0 and 1 evolve in time from a set of initial conditions taken outside of their respective attractors.
Response formulae found in Ruelle's work have to be adapted to take this into account. One therefore has to consider the  Table 1). The variable considered is the temperature meridional gradient θ1. The solid lines denote the mean while the shaded areas denote the one standard deviation interval. 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 systeṁ has a first order response of where f τ is the flow of the unperturbed system (38), ρ 0 is the distribution of initial conditions, and δy(τ ) is the solution of the equationẏ +δy =F (y + δy) which can be approximated at first order by the following linear inhomogeneous differential equatioṅ δy = ∇ y F · δy + Ψ(y). where y(τ ) is solution of the unperturbed equation (38) with initial condition y(0) = y 0 and we see that the systems (38) and (40) have to be integrated simultaneously (Gaspard, 2005). The homogeneous part of Eq. (40) 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 (39) 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 the perturbed systems (model 1) with the same initial conditions, the initial state of the tangent model (40) 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. (40) can be adapted to take these errors into account, as described for instance by Nicolis (2016).
In what follows, we will numerically integrate Eq. (40) to evaluate the response on the average due to the perturbation induced by the model change. This will in turn, as in Section 2, allow us to compute the post-processing parameters for the new model.

Main results
For each of the two experiments detailed in Table 1, we start by obtaining one million observations of the reality that will be used to initialise the forecast models. For each observation, this is done by starting the model x (the reference) with a random initial condition and running it for a very long time (100000 adimensionalized 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 days) to obtain the 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 unit corresponding to 16.15 minutes. The averaging over the one million trajectories of the reality and of the forecasts at each lead time allow for computing the post-processing coefficients α and β of the EVMOS by using formulas (7) and (8). For each variable, the variable itself is used as the unique predictor.
The response-theory approximations of the averages of the modelŷ (model 1) averages are obtained by integrating the linearised equations of model 0 along its trajectories with the perturbation Ψ as inhomogeneous term. This is done by integrating 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. (39). As it can be seen in Fig. 9, the problem worsens with the increase of the lead time: initially the distributions are near-Gaussian and fat-tails appears 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 response of the system located far in the tails of the distributions, we have removed outliers above a certain threshold (set to 3. adimensional unit) from the averaging.
The moments obtained by the response theory approach are used to compute new EVMOS post-processing α and β coefficients, thanks to the formulas (7) and (8). In Figs. 6 and 8, we compare the three post-processing schemes hence obtained: the post-processing of model 0 (red curves) and 1 (green curves) obtained by averaging over their trajectories (forecasts), and the post-processing of model 1 obtained with the response theory approach (black crosses). In the panel (a) of these figures, 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. In general, the corrections are efficient until a lead time of 3-4 days. In fact, the skill of the corrections decreases with the lead time, and thus the EVMOS schemes become not better than the original models after 3 days. In the panels (b) and (c) of Figs. 6 and 8, the mean and variance of the corrected forecasts is compared with those of the original models. Again, these corrections are efficient until 3 days for all the post-processing scheme considered. For all the 3 panels of these figures, the correction of model 1 using the response-theory EVMOS is depicted by black crosses and it matches almost perfectly the score of the "exact" EVMOS obtained with the forecasts of model 1 (dash-dotted green curve), up to a 3 days 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. (7) and (8)). These coefficients therefore degrade sharply after 3 days, as shown by the solid black curve in Fig. 10 (c) and (d). This in turn induces a degradation of the response theory post-processing scheme, as depicted in Figs. 6 and 8.
Nevertheless, this limitation of response theory is not a concern for our present aim, since after the critical lead time of 3 days, the EVMOS skill vanishes anyway. In conclusion, obtaining accurate EVMOS coefficients is therefore critical up to 3 days, which is precisely the timespan where the response theory is the most accurate.

Discussion and Conclusions
Statistical post-processing 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 Hamill et al., 2008) to update the post-processing 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 to obtain the new post-processing coefficients as modifications of the   Note however that in the context of this conceptual quasi-geostrophic model, it turns out one can also obtain good estimates of the post-processing coefficients α and β simply by using a small set of reforecasts. For this it suffices to directly integrate the updated model 1, given by the non-linear equation (37), for only 20 trajectories. So the response theory approach in the present case cannot really compete with the simple reforcasting 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 often used in data assimilation (Bonavita et al., 2017). This approach could also be implemented for short-range forecasts, say from 1 to 3 days.
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 post-processing scheme, but it can efficiently adapt it to a new model. As such, the success of this method also depends on the quality of the past post-processing scheme. There are also situations where linear response theory is known to fail, but statistical tests which allow to identify 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 of the applicability of the approach.
To test this approach, we have focused on the EVMOS statistical post-processing method, but other methods could be considered as well. The only requirement is that the outcome of the minimisation of the cost functions uses averages of the systems being considered. In fact, we hypothesise that any statistical method using average quantities could be concerned by our approach. 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 ingredient, instead of their marginal variance for the EVMOS. However, it does not preclude the applicability of response theory since this covariance can easily be written as an ensemble average. This will be investigated in a future work, together with the applicability of the approach to parameters of probability distributions, as often used in meteorology .
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 days 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 days. 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 due to the exacerbated sensitivity of these manifolds to the perturbation of the system. Indeed, as Fig. 11 suggests, the trajectories responsible for these extreme responses concentrate near some heteroclinic connection between the two regimes mentioned in the introduction of Section 3 and depicted in Fig. 2. 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 more 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 very effective techniques based on the Covariant Lyapunov vectors (CLVs) to non-stationary dynamics. These techniques were recently introduced (Wang, 2013;Ni and Wang, 2017;Ni, 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 to change easily the perturbation function Ψ for a fixed observable A, while the direct method allows for different observables to be considered, 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 good 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 post-processing 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.
Code availability. The quasi-geostrophic model 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). The release of the full source code as well as the notebooks generating the figures is in preparation. It will be ready by the end discussion process. It will be published on GitHub, and a link will be provided through the Zenodo service.
Appendix A: Non-stationary response theory We consider a perturbed autonomous dynamical systeṁ with a prescribed distribution of initial conditions ρ 0 . For the unperturbed systeṁ an observable A has the average at time τ where f τ is the flow of the unperturbed system (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, A(τ ) ŷ = A(τ ) y + δ A(τ ) y + δ 2 A(τ ) y + . . .
In other words, we compute the average of A in system (A1) A(τ ) ŷ = dy 0 ρ 0 (y 0 ) A(f τ (y 0 )) (A5) as a perturbation of the average (A3) in the unperturbed system (A2). Here,f τ is the flow of the perturbed system (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 the system (A1): with the operators    L 0 A(y) = F (y) T · ∇ y A L 1 A(y) = Ψ(y) T · ∇ y A and define an interaction observable as with Π 0 (τ ) = exp (L 0 τ ). It is easy to show that the interaction observable satisfies the differential equation: which can be rewritten as A f τ (y 0 ) = Π 0 (τ )A(y 0 ) + τ 0 ds 1 Π 0 (τ − s 1 ) L 1 Π 0 (s 1 ) A I (s 1 , y 0 ) .
We will now focus on the first term of this expansion, but the subsequent orders of the response can be treated alike. We thus have δ A(τ ) y = τ 0 ds 1 dy 0 ρ 0 (y 0 ) Ψ f τ −s1 (y 0 ) T · ∇ f τ −s 1 (y0) A (f τ (y 0 )) (A15) which with the change of variable s 1 → t − τ can be rewritten as δ A(τ ) y = τ 0 dτ dy 0 ρ 0 (y 0 ) Ψ f τ (y 0 ) T · ∇ f τ (y0) A (f τ (y 0 )) (A16) and then δ A(τ ) y = τ 0 dτ dy 0 ρ 0 (y 0 ) Ψ f τ (y 0 ) M is the fundamental matrix (Gaspard, 2005;Nicolis, 2016) of the homogeneous part of the linear differential equatioṅ δy = ∇ y F · δy + Ψ(y) where y is solution of Eq. (A2) with initial condition y 0 , and we have the definition M(t, y) = ∂f t (y) ∂y . (A20) Equation (A19) is the linearised approximation of equation (A1): y +δy = F (y + δy) + Ψ(y + δy) that provides a tool to estimate Eq. (A18). Indeed, since the solution of Eq. (A19) can be written as we can write the first order variation of the average of the observable A in term of these solutions: δ A(τ ) y = dy 0 ρ 0 (y 0 ) δy(τ ) T · ∇ f τ (y0) A (A23) The interpretation of this equation is that the averaging of an observable over the trajectories of the linear approximation (A21) of the perturbation equation (A1) provides the first order response of the observable. It is the main ingredient used to compute the new post-processing scheme in the present work. It is explained in detail in Sections 2.3 and 3.2.