Full-field and anomaly initialization using a low-order climate model: a comparison and proposals for advanced formulations initialization using a low-order model: a comparison

. Initialization techniques for seasonal-to-decadal climate predictions fall into two main categories; namely full-ﬁeld initialization (FFI) and anomaly initialization (AI). In the FFI case the initial model state is replaced by the best possible available estimate of the real state. By doing so the initial error is efﬁciently reduced but, due to the unavoidable presence of model deﬁciencies, once the model is let free to run a prediction, its trajectory drifts away from the observations no matter how small the initial error is. This problem is partly overcome with AI where the aim is to forecast future anomalies by assimilating observed anomalies on an estimate of the model climate.


Introduction
State estimation theory in geosciences is commonly referred to as data assimilation (DA), (Daley, 1991). This term encompasses the entire sequence of operations that, starting from the observations of a system, and possibly from additional statistical or dynamical information (i.e., the model), provides the best possible estimate of its state (Kalnay, 2003). This estimate, the analysis, is then used for diagnostic purposes or as initial condition for predictions with environmental numerical models. Another notable application of DA, having paramount importance in climate science, is in the production of reanalyses, multiyear global state-of-theart gridded representations of the atmosphere/ocean generated by the same model and the same DA method (Dee et al., 2011). The ultimate goal of DA is to give a dynamically consistent reconstruction of all the elements of the climate system. By improving the initial condition, through a better use of the observations and model, DA has dramatically contributed to enhance the forecast skill in weather and ocean prediction in the last decades, and is nowadays regarded with attention from the seasonal-to-decadal (s2d) community interested in improving the initialization procedures.
At the basis of this growing interest there is the hope that optimizing the simultaneous use of model and observational information at the initialization step will help to improve the prediction skill in all those circumstances, and for those time horizons, when a nonnegligible portion of the forecast uncertainty is explained by the internal climate variability. It is known that this component accounts for a significant amount of the total prediction uncertainty at global scale for up to a decade in the future (and even longer at regional scales), and tends to monotonically decay afterward, dominated by the uncertainty associated with the model and the projected scenario (Smith et al., 2007;Hawkins and Sutton, 2009). Seasonal to decadal predictions, in contrast to continental-scale projections of climate change, are initialized using observations of the current climate state. Skillful (initialized) seasonal forecasts are nowadays run in several operational climate services worldwide; a recent review on the development and current status of the seasonal forecasting practice can be found in Doblas-Reyes et al. (2013a). Using global climate models and simulated observations, several studies have identified and assessed the benefits of initializing decadal predictions (see e.g., Latif et al., 2006;Dunstone and Smith, 2010, and references therein) and set the basis for the design of decadal prediction practice using real observations (see e.g., Smith et al., 2007;Keenlyside et al., 2008;van Oldenborgh et al., 2012;Hazeleger et al., 2013).
Current initialization techniques for seasonal-to-decadal climate predictions fall into two main categories; namely full-field initialization (FFI) and anomaly initialization (AI).
In FFI the initial model state is replaced by the best possible available estimate of the real state. By doing so the initial error is efficiently reduced but, due to the unavoidable pres-ence of model deficiencies, once the model is let free to run a prediction, its trajectory drifts away from the observations no matter how small the initial error is (e.g., Stockdale, 1997). This problem is partly overcome with the AI where the aim is to forecast future anomalies by assimilating observed climate anomalies on an estimate of the model mean climate. In this way, the initial model state is kept on (or closer to) its own attractor (e.g., Smith et al., 2007). The large variety of experimental setups, models and observational networks adopted in the studies to date makes it difficult to draw firm conclusions on the respective advantages and drawbacks of FFI and AI, let alone identifying distinctive lines for improvement. The lack of a unified mathematical framework adds an additional difficulty toward the design of adequate initialization strategies that fit the desired forecast horizon, observational network and model at hand.
Comprehensive comparisons between FFI and AI using state-of-the-art coupled climate models for seasonalto-multiyear time horizons, have recently appeared (Magnusson et al., 2012;Smith et al., 2013;Hazeleger et al., 2013). These studies represent a first attempt to analyze their impact using exactly the same observational and model setup and are therefore of central importance in guiding future development of initialized climate prediction systems. Overall, these results have indicated that initialization systematically improves over a climatology at seasonal timescale, with a slight, but clear, better skill for FFI (Magnusson et al., 2012;Smith et al., 2013). The picture is far less clear on multiyear horizons. Some advantage of AI was reported in Smith et al. (2013) at global scales using the MetOffice decadal prediction system, while using the EC-Earth v2.3 atmosphere-land-ocean-sea ice model, Hazeleger et al. (2013) have found a somewhat better skill for FFI in areas such as the North Atlantic. This behavior is thought to be related to the difficulties of generating physically consistent initial sea-ice and ocean conditions. The present study aims to contribute to the current debate on which initialization algorithm, and under which conditions, is the most suitable for seasonal-to-decadal prediction. We propose using the notation and concepts of data assimilation theory to outline a unified formalism from which FFI and AI can be derived. This furthermore allows for identifying specific features such as the initialization method sensitivity to model and observational accuracy. FFI and AI are studied in a range of different observational and model error scenarios, helping to clarify under which conditions one approach outperforms the other. Two advanced formulations are then introduced and discussed. The first, named leastsquare initialization (LSI), is aimed to improve the fit to the observations allowing for their informational content, usually restricted to the observational space in standard FFI or AI, to be propagated to the entire model domain. LSI uses a leastsquare initialization update to merge observation with the model, and the required (unknown) model error covariance is estimated using the covariance of the model anomalies. The Nonlin. Processes Geophys., 21, 521-537, 2014 www.nonlin-processes-geophys.net/21/521/2014/ estimation and reduction of the drift associated with parametric model error is the objective of the second novel formulation proposed in this study: the exploring parameter uncertainty (EPU) method provides an online correction of the drift on the basis of a linear and short-time approximation of its evolution. EPU works during the forecast run and as such can be used in combination with any initialization scheme, either FFI or AI. The numerical analysis in this work is carried out using a prototype of climate dynamics simple enough to reduce computational cost and allow for robust statistical inferences. These conditions are not easily met when dealing with realistic models and observational scenarios. In this study we aim at gaining additional insight on FFI and AI, to advise and support the big ongoing research effort in the development of initialized decadal prediction with state-of-the-art climate models. Moreover, the capabilities of the proposed LSI and EPU methods are better understood using a more controllable experimental setup before thinking of a possible implementation in a more realistic context. The paper is organized as follows. In Sect. 2 the general problem is posed and the basic assumptions are described. Section 3 describes FFI and AI using the data assimilation formalism, while LSI and EPU are presented in Sect. 4. The low-order climate model is described in Sect. 5 and the results are given in Sect. 6. Final Conclusions are drawn in Sect. 7.

Posing the problem
Let us write the prognostic climate model under the form of an autonomous dynamical system: with x and λ being the state and parameter vectors of dimensions I and P , respectively. The parameters can also have a time dependence, such as is the case of the impact of solar activity or the anthropogenic greenhouse gases, which has been dropped here to simplify the notation. The model (1) is used to describe the unknown Earth's climate system evolution; we will call our target, the real world state, "nature". We suppose that nature can be formally expressed under the form of a system similar to model (1), such as with x nat and λ nat being the unknown nature state and parameters. In this formulation, the model and the nature are assumed to span the same phase space, of dimension I . As a consequence, model error can only originate from parametric uncertainty, δλ = λ − λ nat , and from the presence of the extra term, G, accounting for all processes not properly described by the model, with | G | | F | in some proper norm (Nicolis, 2003). Note furthermore that the model is assumed to describe all scales of motion: model error coming from unresolved scales, a common drawback in geophysical modeling, is not considered in this formulation. Observations, y o i = y o (t i ), are assumed to be available at equally spaced times t i = iτ , i = 0, 1, . . .; in real applications, the vector of observations has a much smaller dimension than the model state vector, so that O I with O being the dimension of y o i , i.e., the number of observations available at each observation time. Moreover, unless specified otherwise, these observations are affected by some error represented as a white Gaussian noise, o , with zero mean and standard deviation σ o , so that o ∈ N (0, σ o ), with N standing for the Gaussian distribution and H being the observation operator mapping from nature to the observation phase space. Note finally that in this formulation the observational error, o , accounts for both the instrumental and the representativity error connected with the specification of the operator H (Cohn, 1997;Janjić and Cohn, 2006).

Initialization methods
The initialization procedures start typically from a model state obtained after a long spin-up run. By using terms and concepts borrowed from a data assimilation context (see e.g., Kalnay, 2003), the model state at the end of the transient spin-up (the control) can be interpreted as a background field, embedding all information about the real system (nature) prior to the assimilation of the observation; it will be hereafter indicated as x b . Similarly, the initial condition obtained at the end of the initialization procedure is interpreted as the analysis field, and indicated as x a in the following. Here, the additional assumption of a linear observation operator is made; this implies that the observation operator is given under the form of the O × I matrix, H.

Full-field initialization (FFI)
This approach consists of substituting the model state with the available observations at the initialization times. Using the data assimilation notation FFI can be formally written as with H T being the transpose of the linearized observation operator; note that in this case the entries of the matrix H are equal to one in correspondence with the observed variables and grid points, and zero elsewhere. By left multiplying both sides of Eq. (4) by H, we retrieve the usual FFI formula in which the control is replaced by the observations at observational locations, and left unchanged elsewhere. Alternatively to the direct replacement of the observed values, FFI can also be implemented using a nudging approach.
Nudging is an empirical data assimilation technique, consisting of adding a term to the prognostic equations that, acting like an extra-coupling term, nudges the model trajectory toward a consistent representation of the unknown system's evolution intended to be estimated (Hoke and Anthes, 1976). The coupling strength is expressed as a relaxation timescale and is usually chosen on the basis of the properties of the variable to be nudged: it has to be small enough to avoid dynamical shocks but large enough for the correction to counteract actively the error growth. Advanced formulations of the nudging approach have been proposed recently (Auroux and Blum, 2008), and the use of the nudging procedure using unstable modes of the dynamics has been also explored using simple chaotic systems (Yang et al., 2006).
The FFN can be written as follows: with N being the diagonal O×O nudging matrix, with the dimension equal to that of the observation vector and measured in units of time −1 . The matrix N contains, along its diagonal, the relaxation time rate of each nudged/observed variable. Any specific setups of a nudging scheme are equivalent to specifying the entries of the matrices H and N: the former indicates the observed variables, the latter the corresponding relaxation times. Note that the nudging observation vector y o nudg , for which the time dependence has been omitted for clarity in Eq. (5), coincides with the observation y o at observation times, t i , and with their time interpolation between two successive observation times. The initial condition, x a , is then obtained by integrating Eq. (5) up to the initialization time (i.e., the start date), t init . FFN has recently been used in Magnusson et al. (2012) to nudge SST (sea surface temperature) with the ECMWF (European Centre for Medium-Range Weather Forecasts) climate model.

Anomaly initialization (AI)
In AI (Smith et al., 2007) the model state at initial time is replaced by the observed anomaly (i.e., the difference between the current observation and its long-term average 1 ) plus the model average (i.e., the simulated climate). To write the AI equation using the data assimilation formalism, we introduce the pseudo-observation vector: with the overbar indicating the long-term average. The pseudo-observations are equal to the observations minus the difference between the observed and modeled climate. The AI equation can then be written as The observation operator H is defined as in Eq. (4) and, as above, by left multiplying both sides of Eq. (7) by H, we see that at observational locations, the background is replaced by the pseudo-observations and left unchanged elsewhere.
To get the expression for anomaly nudging, AN, we can proceed as for Eq. (5) and after introducing the pseudoobservations we have Similar to FFN, the initial condition, x a , is obtained by integrating Eq. (8) up to t init . AN has been used in Smith et al. (2007) and Smith et al. (2013) to nudge oceanic variables of the MetOffice global climate model.

FFI and AI: some properties
FFI and AI are easy to implement and do not require any hints on the relative accuracy of the estimators, the background and observations, entering the initialization procedure. These features make them very attractive in the context of seasonal-to-decadal prediction with global numerical models whose typical size and complexity limits the use of more advanced initialization strategies. FFI and AI produce different initial states, and one approach can outperform the other depending on a number of competing factors such as the accuracy and type of the observational network, the amplitude of the model biases and the desired forecast horizon. It is interesting to consider the error-scaling properties of FFI and AI, with respect to the level of accuracy of the information sources, model and observations. To this end, let us derive an expression for the initial condition (analysis) error, a = x a −x nat , by formally subtracting the unknown nature's state, x nat , from Eqs. (4) and (7). Assuming for simplicity that the full system is observed (i.e., H = I, with I being the I × I identity matrix) and after rearranging we have x nat being the background error. The expressions (Eq. 9) highlight a main feature of these methods; namely, their dependence on the accuracy of the observations and, for AI, the role of the model bias, b , as the error level in the idealized limit o → 0. By taking the average of Eq. (9), and using the property of unbiased observations, we get the expressions for the mean initial error: Nonlin. Processes Geophys., 21, 521-537, 2014 www.nonlin-processes-geophys.net/21/521/2014/ whereas for the root-mean-square error at analysis time (RMSE a ) we have: RMSE a FFI scales linearly with the observational error, while the behavior of RMSE a AI depends on the ratio b 2 σ o 2 that measures the relative accuracy of model-observations. In realistic circumstances this ratio is (much) larger than 1 ( b 2 σ o 2 , implying the observational error being much lower than model error) and RMSE a AI behaves almost independently on the accuracy of the observations. As a consequence and in contrast to RMSE a FFI , which can in principle be arbitrarily reduced by intervening on the observational accuracy, a significant reduction of RMSE a AI can only be achieved by improving the model. Moreover, when b 2 σ o 2 ≈ 1, RMSE a AI behaves almost quadratically with σ o so that for the same reduction of σ o , RMSE a AI will decrease at a slower rate than RMSE a FFI . It is worth mentioning that we have intentionally focused our discussion on the actual (unknown) error with the purpose of highlighting the response of FFI and AI to observation and model error. In real applications one usually computes errors using observations, possibly taken from a data set independent from the one used in the initialization. If we had proceeded similarly and had used pseudo-observations for the AI analysis, the two schemes would have shown similar scaling properties with respect to the observational error. Nevertheless this approach would have hidden the presence of the model bias, whose role we wanted to stress here.
Arguing on how the accuracy of the initial conditions will impact the prediction skill at seasonal-to-decadal timescales is far more difficult. In contrast to weather forecast practice (with horizons of 2 weeks) where model error is often neglected and much attention is placed toward an efficient control of chaotic error growth (see e.g., Palatella et al., 2013), in seasonal-to-decadal prediction the growth of initial error is relatively less important than the bias caused by model deficiencies. The unavoidable presence of model error causes the so called "model drift", in which a forecast initialized close to the observed state will eventually drift toward the model climate, i.e., the model attractor, following a nonlinear statedependent evolution that can appear as quasi-erratic before stabilizing. With AI, although at the price of larger errors, the initial state is kept close to the model attractor and the drift is mitigated. This makes the mean forecast error less time dependent and, as argued by Magnusson et al. (2012), the use of standard a posteriori bias correction techniques is more robust. Nevertheless, although the bias correction postprocessing for forecasts initialized with AI has the advantage of being independent from the forecast times, it still requires running long and computationally expensive transient simulations before reaching model equilibrium.
FFI is more prone to suffer initial dynamical shocks caused by the displacement of the model state to the observed values that can be out of the model attractor. These initial shocks are smaller, but still present, in the AI, since the model is forced to deviate from its mean state by only the amount of the observed anomalies. However shocks can also be generated by inconsistencies or geographical mismatches between the model climate and observed anomalies (Magnusson et al., 2012). By smoothly moving the model trajectory toward the observations, the nudging approach has the potential to reduce this problem (see e.g., Smith et al., 2013).

Advanced initialization procedures
In this section we introduce and discuss two advanced initialization methods based on data assimilation theory. The first approach, referred to as LSI, is thought to improve the fit to the observations. The second, named EPU, is designed to estimate and reduce the model drift. LSI and EPU are described here in the context of FFI, however they can equivalently be applied in the framework of AI with only minor modifications.

Least-square initialization -LSI
In FFI the observations are fit as if they were perfect and the model state is replaced by the observations in all variables and geographical locations where measurements are available, no matter how accurate these observations are with respect to the model solution. The information is not transferred from observed-to-unobserved areas and the observational impact is confined to the measurement locations, despite the possible presence of physically relevant spatial correlations. Note however that, in the practice of climate prediction with realistic models, this issue is partly mitigated by using either a full and homogeneous reanalysis field of all model variables over its entire domain, or by smoothly nudging toward the observations data set when this represents only a subset of the model variables/domain (Meehl et al., 2013;Doblas-Reyes et al., 2013b).
In data assimilation, statistical knowledge about the accuracy of each piece of information entering the analysis update is used to determine the relative weights of their contribution according to some criteria of optimality (Daley, 1991). Using a least-square framework, the background and the observations are linearly combined in order to minimize the expected analysis error variance (Jazwinski, 1970). All errors are assumed to be Gaussian and represented by the first two statistical moments only: the mean and the covariance. The background and observation error covariances are referred hereafter to as B and R, according to a standard notation in data assimilation literature (Ide el al., 1997).
The least-square analysis update, common to most data assimilation procedures, reads (Jazwinski, 1970) x When B and R are the correct error statistics, x a would coincide with the best linear unbiased estimate and the analysis solution minimizes the associated posterior error covariance. The matrix B plays the crucial role of spreading out the observational information content throughout the entire model domain according to the assumed forecast error structure; in practice the analysis correction is confined within the subspace spanned by the range of B. Classical data assimilation methods, such as optimal interpolation or 3DVar (Daley, 1991), rely on a statistical estimation of the covariance matrices, B and R. Nevertheless, while a relatively robust estimate of the latter can be easily obtained especially when it is related to the observational accuracy alone, more problems arise to estimate B. The approach that has been almost universally adopted since the 1990s is known as the NMC (National Meteorological Center)-method (Parrish and Derber, 1992). With the NMC B is estimated using the difference between two forecasts verified at the same time, let us say the 24 and 12 h forecasts; this matrix is then used as a proxy for the true (and unknown) forecast error matrix. And it is not based on observations, that are usually sparsely distributed over the globe, the NMC method provides a proper representation of the global error structure (Kalnay, 2003), and it is also consistent with the hypothesis of uncorrelated forecast and observation error used in Eq. (12). Finally, it is worth mentioning the research efforts carried out in the last decades that led to the development of the so-called ensemble-based data assimilation algorithms in which the forecast error covariance matrix, estimated on the basis of an ensemble of model trajectories, is made time dependent, a very desirable feature when dealing with chaotic systems (see e.g., Evensen, 2009). In the initialization of climate predictions, the NMC method cannot be directly applied because the background field here is a climatological (control) model solution output after a sufficiently long spin-up, (see Eqs. 4, 7). In the LSI method described here the actual B is approximated using the model climate covariance, which is then incorporated in the analysis update (12). The model covariance can be estimated rather accurately, and the robustness of this estimate is limited only by the availability of enough computational resources to make long runs feasible. Under this assumption, the proxy of the background error covariance matrix reads where x is a solution of the model (1), α a scalar tuning parameter and, as above, the overbar refers to a long-term average. The full-field least-square initialization (FF-LSI) equa-tion then becomes Equation (13) embeds the idea of interpreting the model anomalies as forecast errors (the actual B reads in fact B = b b T ) and we expect B m to reproduce at least the spatial structure of the error covariance. The parameter α compensates for the deviation of the amplitude of the simulated covariance from the error covariance and helps to adapt the weighting of the background term in Eq. (14). The occurrence of nature's statistical modifications (i.e., a climate change) that have not been tracked by the model during the reference averaging period can be a further limiting factor for the accuracy of Eq. (13). However, FF-LSI is easy to implement and has the potential to be beneficial in all the situations when only a portion of the model state vector is observed.
A very similar approach was originally introduced by Smith and Murphy (2007) to generate a global ocean analysis based on sparsely distributed observations of temperature and salinity. Using a state-of-the-art coupled climate model they demonstrate that using model-based covariance allows for successfully propagating information from datarich to data-poor areas and significantly improved the model representation of the observed variability. The LSI procedure follows a similar line, although our attention here is in the propagation of information between different model compartments, such as the atmosphere and the ocean, and the experiments described later are designed for that purpose. Another difference between the two approaches stands in the definition of the background field, x b . In LSI this is a forecast-field solution of the model equation at analysis time, while Smith and Murphy (2007) uses a climatology field. Our choice has the advantage of incorporating fresh, timedependent, information about the system, potentially propagating forward the signal of the most recent observations. However it would in principle require a corresponding timedependent background error covariance that is very difficult to get unless sophisticated data assimilation, such as the ensemble Kalman filter (Evensen, 2009), is adopted. The use of a climatological background avoids this issue but at the price of a static a priori picture of the system's state. Finally note that, instead of using the tuning coefficient α, in Smith and Murphy (2007) the misrepresentation of the actual B is treated by limiting the use of Eq. (14) to the set of observations that lead to a decrease of the posterior analysis error variance; the latter is computed using observations when available and turns the global solution (14) into a series of local analyses centered around each grid point.

Exploring parameter uncertainty -EPU
The second novel formulation discussed in this study is named full-field initialization with exploring parameter uncertainty (FFI-EPU), and is designed to reduce the model drift caused by parametric error. EPU is an online bias correction method based on a linear, and short-time, approximation of the bias evolution; the procedure is said to be online because the estimation and correction of the bias is done during the prediction run. A similar approach, in the offline mode, has been introduced in the context of sequential data assimilation by Carrassi et al. (2008), on the basis of the theory of deterministic model error dynamics given in Nicolis (2003).
An equation for the evolution of the estimation error, δx, can be obtained by linearizing the model equation, (1), along one of its solutions. By expanding Eq. (1) to the first order in δx and δλ, we obtain The solution of Eq. (15), with initial condition δx 0 , reads ∂F ∂x | x dt being the tangent linear model of Eq. (1) between t 0 and t, and δµ = ∂F ∂λ | x,λ δλ. The vector δµ embeds all information about the model error through the parametric error δλ and the functional dependence, ∂F ∂λ , of the dynamics on the uncertain parameters. Assuming that the initial error is unbiased, after taking the average of (16), we get an estimate of the model bias evolution, the drift, as Expression (17) provides the time evolution of the bias under the hypothesis of linearity aforementioned. Nevertheless, it cannot be solved in the case of realistic climate models, mainly because of the huge dimension of the systems involved, and some approximations are required.
From Eq. (17) we see that the bias evolution (the drift) is fully correlated in time. Therefore, an approximation suitable for realistic applications can be obtained by expanding it in a Taylor series up to the first order in time, so that: Equation (18) gives the short-time evolution of the bias; its accuracy is related to the linear hypothesis made above and to the duration of the short-time linear regime. The extent of duration of this regime is known to be proportional to the largest (in absolute value) Lyapunov exponent of the dynamics (Nicolis, 2003). For the large class of dissipative chaotic systems to which environmental models belong, the largest Lyapunov exponent is usually the most negative.
The basic idea behind EPU is to use Eq. (18) to estimate and remove the bias from the forecast at lead time t − t 0 . Let us define the bias correction time interval, T Bias , as the time period over which the bias evolution law (18) is used and the bias correction is applied. Every T Bias the bias is estimated and removed from the forecast at time t i according to: where the suffix "un" stands for unbiased, and the compact form x(t i ) = x i is used to simplify the notation. In Eq. (19), λ represents the value of the parameter used in the model along its entire run, while δλ i the (estimate of) parametric error at time t i . The use of Eq. (19) as a bias correction approach requires the estimation of the operator C i = ∂F ∂λ | x i−1 ,λ δλ i . The first term of C, ∂F ∂λ , gives the model functional dependence on the uncertain parameters, i.e., the model sensitivity; it depends on time through the dependence on the system's state and allows us to project the error from the parameter to the state vector space. It is a rectangular matrix having as many rows and columns as the dimension of the state vector and of the parametric error respectively, and can be computed at any time along the model integration, the only practical concern being the computational constraints arising when using large numerical models.
The second term, the parametric error, δλ, is unknown. To cope with this we introduce a guess strategy. Two hypotheses are formulated: (1) model users have identified a (possibly limited) set of relevant parameters suspected to be uncertain, and, (2) a range of possible parameter values, = [λ min , λ max ]. The parametric error at time t i , δλ i , is then sampled according to with U(a, b) being the uniform distribution in the interval (a, b), λ the mean value of the range and λ the parameter used in the model. In practice λ plays the role of the most probable parameter value and is used in Eq. (20) to discriminate between over-/underestimation of the unknown λ nat . Note that, even in the favorable situation λ nat ∈ , given that λ nat = λ in general, it is possible that the guess strategy (20) erroneously selects positive/negative parametric error (δλ i > 0/δλ i < 0, i.e., over-/underestimation) when the true parametric error is actually of the opposite sign (λ − λ nat < 0/λ − λ nat > 0).
In the experiments described in Sect. 6 we will put ourselves in this situation and study the capability of EPU to reduce the model drift of predictions initialized using FFI. FFI-EPU is thus given by the FFI update (4) in conjunction with the EPU bias correction technique (19)-(20) during the forecast run.

The low-order climate model
The low-order climate model is based on the Lorenz 3-variables model (Lorenz, 1963), and has been introduced by Peña and Kalnay (2004). To mimic the simultaneous presence of different dynamical compartments having different timescales, three similar copies of the original Lorenz model are coupled to simulate the extratropical and tropical atmosphere and the ocean. The equations read The atmospheric variables are denoted with the lower-case variables, with the subscripts e/t referring to the extratropical/tropical atmosphere; the ocean variables are denoted with capital letters. The two atmospheres are coupled through the variables x and y at a strength given by the parameter c e ; the tropical atmosphere and ocean are coupled through all variables with a strength given by the parameters c, for the x and y, and c z for the z component. The parameters, σ = 10, r = 28, and b = 8/3, are the same as in the classical 3-variable Lorenz model (Lorenz, 1963)); the "uncentering" parameters k 1 = 10 and k 2 = −11, taken from Peña and Kalnay (2004), introduce a sort of phase lag between model compartments. The parameters S and τ modulate the amplitude and timescale of the ocean: they are chosen as in Peña and Kalnay (2004) to be S = 1 and τ = 0.1 respectively, implying that the ocean variables will have the same amplitude as the atmospheres' but a slower rate by one order of magnitude.
In all subsequent experiments the nature climate evolution is represented by a solution of Eq. (21) with the coupling strength parameters set to the standard values, c e = 0.08 and c = c z = 1. This configuration implies that the tropical atmosphere and the ocean are strongly coupled, while the two atmospheres are only weakly coupled. Numerical solutions of Eq. (21) are obtained using a second-order Runge-Kutta scheme with time steps of δt = 0.01. According to Peña and Kalnay (2004), the model (21) represents an ENSO (El Niño-Southern Oscillation)-like configuration with an almost slave, small amplitude, atmosphere whose regime changes are modulated by the slow ocean component. Following their convention, a simulated year is made to correspond to an ocean regime, and the system oscillates between the normal regime, lasting between 3 and 12 years, and an El Niño regime, lasting only 1 year (equivalent to 240 time steps in the present experimental setup), (Peña and Kalnay, 2004). With the standard values for the coupling parameters given above, the model (Eq. 21) is chaotic with two positive Lyapunov exponents, γ nat 1,2 = 0.9063, 0.3150, while the third corresponds to the null one.
To simulate parametric errors originating at the level of the coupling between the different model compartments, we have modified simultaneously the tropical atmosphere/ocean coupling parameters c m and c m z , with the superscript m standing for "model", in the range 0.1-1.5 with steps of 0.1. We have assumed that the coupling between the two atmospheres is known, so that c m e = c e . To place ourselves in the situation in which the model is able to reproduce the qualitative behavior of nature, we have restricted our analysis to the 109 parameter combinations for which the model stability properties, as measured by the first three Lyapunov exponents, do not differ too much from those of nature. We further assume that the amplitude of the natural variability and its spectrum is well reproduced within each model component. This is not true and it might lead to an overestimation of the performance of AI in our conceptual model. Figure 1 shows, with shadow bars, the distributions of the model bias and of the first three Lyapunov exponents relative to these 109 configurations; the RMS bias is defined by taking the square root of the mean quadratic biases on all of the 9 model variables, each one normalized with respect to its own variance: .
The distribution of RMS bias values appears relatively homogeneous with a marked peak between 0.2 and 0.25. As expected, the range of values of the first Lyapunov exponent, γ 1 , does not depart too much from γ nat 1 = 0.9063, given that the two perturbed coupling parameters, c and c z , do not directly affect the dynamics of the extratropical atmosphere to which γ 1 is associated. In contrast, the Lyapunov exponents γ 2 and γ 3 display a much larger variability; in particular, most of the coupling configurations give rise to a less/more unstable tropical atmosphere/ocean with respect to the standard "nature" configuration. This is due to the form of the coupling between these two model compartments and to the role of the coefficients c and c z : note that in a few cases γ 3 becomes positive, and the null exponent corresponds to the fourth one.
To illustrate how the statistical properties of the model (21) respond to a change in the parameters, c and c z , the dependence of the variance versus the mean for each model variable and for all of the 109 model parameter variations is shown in Fig. 2; note that the values of the mean and the variance are normalized with respect to those of nature. The figure clearly reveals the complex interplay between the first two moments of the model probability density function, typical of nonlinear dynamics, and points out to the difficulty in estimating the climate variance of a model in a regime of bias change (i.e., a climate change). This figure illustrates the dependence of the climate variance on its mean state whereas the hypothesis of independence of those is often made in climate prediction. In particular, the classical anomaly initialization method relies on such a hypothesis. According to this figure, a scaling of the variance should be applied also in AI.

Results
In this section we describe the results of the numerical analysis. We proceed by showing first the results of the comparison between FFI and AI, while LSI and EPU are discussed in the two subsequent subsections.

FFI and AI comparison
The experimental setup is as follows. We have worked under the standard observation system simulation experiment (OSSE, Bengtsson et al., 1981) configuration in which the simulated nature evolution is sampled at discrete times to generate the series of simulated observations. The setup has then been structured following the typical hindcast format of climate prediction studies. A 40-year-long hindcast period is considered and the observations/start dates are distributed homogeneously every month during the first 30 years, making a total of 360 start dates. The effective number of independent start dates ( Van Storch and Zwiers, 2001), N eq , is around 241, 154, and 127 for the extratropics, tropics and ocean respectively; N eq has been calculated using the novel approach given in Guemas et al. (2014) (their Eq. 7). Note that in the CMIP5 (Coupled Model Intercomparison Project) initiative (Taylor et al., 2012), at most we have 52 decadal predictions, with some forecasting systems having only 10. The "nature" and the "control" runs are integrated, after a long spin-up of 60 000 time steps, during the hindcast period. The observations are generated by sampling the nature trajectory and then adding a Gaussian white noise, o , with zero mean and standard deviation, σ o , so that o ∈ N (0, σ o ). Ten-year-long predictions are initialized, at each start date, by using either FFI, Eq. (4), or AI, Eq. (7), with either observations of each model compartment independently (H with dimension 3×9) or of the whole system (H equal to the identity 9 × 9 matrix). The difference between the model and observed mean states, Hx − y o , entering the AI equation (7), is estimated using the sample of 360 observations/start dates, for the case of observing the whole system (i.e., H = I).
The distribution of this "estimated bias" is shown superimposed onto the "real bias" distribution in the top-right panel in Fig. 1: their comparison helps to visualize the difference between the two and quantify the error in the bias estimate used in AI. We see that the estimated biases reproduce relatively well the actual distribution, although some discrepancies are evident in particular in relation with the position of the respective largest peaks, slightly shifted toward higher values in the estimated bias distribution.
In most of the discussion that follows, we make use of the RMS skill score (RMSSS) to measure the skill of the initialized predictions. RMSSS is defined as: where RMSE Init/UnInit refers to the RMS error for the initialized and uninitialized predictions respectively; the latter corresponds to the control run. The RMSE is computed using the anomalies with respect to the mean error over all start dates. In Fig. 3 the RMSSS, for both FFI and AI, is shown as a function of time over a 10-year-long prediction. The RMSSS is calculated by averaging the RMSE over all start dates and is displayed after a monthly averaging of the variables; the observation error standard deviation of each variable is set to σ o = 1.5 % of the corresponding system's natural variability (i.e., the square root of the system's climate variance). Note that in a realistic prediction framework, observation error is meant to incorporate also the representativity error arising from the mismatch between the model and the observation-resolved scales. These errors are not present in our experimental setup, and this legitimates the choice of a relatively small σ o . Two cases of coupling error are considered, c m = 0.8, 0.3 and c m z = 0.9, 1.2: these configurations give rise to biases falling in the peak and in the rightmost extreme of the real-bias distribution (Fig. 1), while the spectrum of the first three Lyapunov exponents is γ 1,2,3 = 0.9036, 0.1895, 0.0, and γ 1,2,3 = 0.9032, 0.2162, 0.0153, respectively. By comparing these values with those of nature given above, we see that both configurations have a smaller second exponent and the second configuration possesses three positive exponents.
The first clear message in Fig. 3 is that the largest, longlasting RMSSSs are obtained when the whole system is simultaneously initialized (black curves). Note however that FFI performs slightly better with a clear predictive skill  (RMSSS > 20 %) until around month 30, while in AI the drop of the RMSSS after the initialization is more abrupt and skillful forecasts last less than 20 months.
The second aspect concerns the role of the ocean as the first driver of the system's predictive skill. When only the ocean is initialized, significant skill is observed (RMSSS > 15 %) until the 20th month. The behavior of FFI and AI is similar, with only very marginal advantages for FFI. The largest initial error reduction obtained with FFI does not seem to be effective in this case to systematically outperform AI: a slow drift (not shown) is present that masks the benefit of this initial error reduction. Having the slowest error growth rate, the ocean behaves like the system's memory. Initializing this compartment is thus the most efficient way to improve prediction skill and move forward the predictive horizon.
The performance when only one of the "atmospheres" is initialized is considered now. When the extratropical atmosphere is initialized with FFI, the RMSSS is large (RMSSS > 10 %) until month 20 in both model configurations. The lack of a similar behavior in AI suggests that this gain in predictability is due to the larger reduction of the actual initial error obtained with FFI. The extratropical atmosphere is in fact the most unstable model compartment, the one with the faster error growth: reducing the initial error counteracts the effects of this growth, and helps to increase the time horizon of skillful predictions. This effect is less pronounced in AI, where the initial error is not much reduced and is sized as the model bias (see Eq. 11 and related comments). Finally, when only the tropics are initialized, the prediction skill for both FFI and AI maintains very low levels.
The behavior of FFI and AI as a function of the observational accuracy at different prediction ranges is analyzed in Fig. 4. In these experiments the full system is observed, and the parametric error is set as in the first configuration (c m /c m z = 0.8/0.9); results (not shown) for the case (c m /c m z = 0.3/1.2) are qualitatively equivalent. Figure 4 has to be interpreted in the light of the discussion in Sect. 3.3 on the dependence of RMSE a FFI/AI on the observation accuracy σ o . It illustrates that this dependence is propagated and is also found at longer lead times. Note for instance that, in contrast to AI, FFI prediction skill appears clearly related to the observational accuracy up to the third year. A closer inspection of the left-bottom panel reveals that such a behavior, although at a lower extent, is present even in the longer prediction horizon (years 4-5) as long as observations are sufficiently accurate (σ o < 0.03). However in the longest lead time, and for the horizon years 4-5 when σ o > 2.5 %, AI shows slightly better skills, a behavior in qualitative agreement with results of Smith et al. (2013) with the MetOffice climate model.
In a similar way, Fig. 5 analyzes the behavior of FFI and AI in relation with the model bias. As above the full system is observed and the same averaging periods are considered and the observational error standard deviation is set to σ o = 2.5 %. In contrast to Fig. 4, we see here that while FFI is almost insensitive to the model bias, AI displays a marked dependence with a performance that systematically improves when the bias is reduced, that is to say when the model is improved. Furthermore, the two approaches converge to similar RMSSS levels in the limit of long lead times. It is worth mentioning that the independence of the FFI skill on the model accuracy contradicts what is observed in realistic weather and climate prediction practices. We believe this is due to two reasons. First, limiting the source of model error to parameter uncertainty marks a key difference with respect to realistic s2d contexts where model error is often connected to multiple simultaneous causes, including unresolved scales and numerical discretization errors. Second, the model configurations considered in this study, with errors in the coupling parameters, do not give rise to strong initial shocks when initialized; results (not shown) with error in the model forcing confirm this conjecture.
Together, Figs. 4 and 5 give an illustration of the response of FFI and AI to the variation of the two ingredients on which they are based, and the behavior of the forecast error seems to be quite consistent with the conjecture done in Sect. 3.3 in relation with the analysis error. Even if obtained with a very simple model, but with a quite large ensemble of start dates, these results point clearly to key differences between FFI and AI. FFI seems to be favorable when a good observational network is at hand and its skill is expected to improve with a refinement of the observations. However, the performance of AI is related to the model accuracy and its overall skill is expected to increase in coincidence with model upgrades.

LSI -numerical results
In this section we study the performance of FF-LSI in comparison with standard FFI. The observational error standard deviation is equal to σ o = 1.5 % and the error covariance For very small α values the LSI analysis, Eq. (14), tends to over-trust the background field, the observations are practically ignored and the RMSSS converges to the small values associated with uninitialized predictions. By increasing α the RMSSS increases monotonically and eventually exceeds the FFI skill level and stabilizes on a plateau afterward.
The fact that the best performances are obtained for large enough α values and that no improvements are observed when it is further increased is connected to the use of observations that are significantly more accurate than the background. The relevant aspect in the discussion about the potentials of LSI is that, when only a portion of the full system's state is observed, its RMSSS converging level is higher than for FFI. This indicates that the benefit of LSI comes by the spatial correlations embedded in B m that make it possible to propagate the observation's informational content, otherwise restricted to the observation subspace alone, to the entire model's domain. In fact when α is large, the observations are fitted as if they were perfect, and the model state vector is replaced by the observations at their locations. Out of these points the information is propagated according to the off-diagonal elements of B m , those that account for the correlations between different variables and model compartments. This interpretation is also supported by the fact that when a diagonal B m is used the FF-LSI's RMSSS (not shown) converges to the same values as FFI in the limit of large α, meaning that as in FFI the observations are perfectly fit but no correction is applied out of observational locations.
The adequacy of the approximation (Eq. 13) in describing the subspace spanned by the range of the actual B is diagnosed in Fig. 7, which shows the percentage of the explained variance of each of the eigenvectors of B and B m (Fig. 7a), as well as the scalar product of each pair of them (Fig. 7b), v i (B) × v i (B m ), with v i (B) and v i (B m ) being the ith eigenvector of B and B m respectively, × indicates the scalar product, and i = 1, . . . , 9. The percentage of variance explained by the eigenvector v i is calculated as Exp var (v i ) = λ i ( i λ i ) ·100[%], with λ i being the ith eigenvalue. From Fig. 7 we see that the variance distribution over the eigenmodes is accurately reproduced and that, except for the third and fourth that are almost perpendicular to each other, all remaining eigenvectors are almost fully aligned.

EPU -numerical results
Similarly to the previous section we study here the performance of the drift correction method EPU. EPU is applied during the forecast run once the prediction is initialized using FFI; these experiments are hereafter referred to as FFI-EPU. We compare its performance with standard FFI in the absence of drift correction procedure.
We have first studied FFI and FFI-EPU in a set of 109 experiments using the model configurations described in Sect. 5.1; this gives the range of parameter , while the bias correction time interval is set equal to one time step, T Bias = time step. Figure 8 shows the mean and standard deviation of the RMSE relative to the 109 configurations, for FFI (black) and FFI-EPU (red), as a function of the averaging forecast period. Results reveal the clear benefit of using EPU within the forecast year. The relative improvements over FFI are equal to 17, 10, and 8 % in the first month and in the first 6 and 12 months respectively. A certain minor advantage is also found afterward, with relative improvements equal to 4, 0.6, and 0.2 % in the forecast years 2-3, 4-5, and 6-10. The distributions (not shown) of the RMSE for FFI-EPU in the six averaging periods are all shifted toward smaller values with respect to FFI. This behavior is in part observed by looking at the corresponding standard deviations (dotted lines in Fig. 8) that are, except for the period of 4-5 years, lower than for FFI. In Fig. 8 EPU appears to be more skillful within the time horizon of 1 year, suggesting its potential usage in the context of seasonal forecasting. Errors in the sampling procedure and, most importantly, the deviation of the drift dynamics from the linear assumption on which EPU is built seems to limit its benefit on longer timescales.
The impact of the uncertainty on the width of the uniform distributions in Eq. (20) and the length of the bias correction time interval, T Bias , the two factors controlling the EPU setup, is investigated in Fig. 9a  uncertainty about the actual range of possible parameter values. To simulate this circumstance, we have used the actual parameter, λ nat , and modulated the distribution width using a scalar coefficient β, such that δλ i ∈ U 0, β λ − λ nat . As expected, the best performances are with β ≈ 2, as in this situation the mean of the sampling distribution coincides with the actual parametric error. However the important remark is that FFI-EPU outperforms FFI systematically for all values of β in both model configurations.
The sensitivity of FFI-EPU to the length of bias correction interval, T Bias , is analyzed in the right panel of Fig. 9b. The same two model configurations are considered, and are displayed with the black and red lines respectively (see corresponding labels in the figure panel); similarly the RMSE of FFI is shown with dashed lines. To focus on the impact of T Bias alone, we have set the scaling parameter β = 2. As expected, the accuracy of FFI-EPU decreases by increasing T Bias , as the limit of duration of the short-time regime of the error growth is approached. It is remarkable that, with these particular model configurations, FFI-EPU outperforms FFI in a large range of T Bias . Note moreover that the improvements over FFI are kept as long as T Bias < 40 and T Bias < 30 time steps for the first and second model configurations respectively. If these time intervals are taken as measures of the duration of the short-time regime, it is worth highlighting that the numerical result T |, in agreement with the theory of deterministic model error dynamics mentioned in Sect. 4.2 (Nicolis, 2003).
The results in Fig. 9 are quite encouraging about the robustness of EPU. Although the extent of their validity with Nonlin. Processes Geophys., 21, 521-537, 2014 www.nonlin-processes-geophys.net/21/521/2014/ more realistic models cannot be assessed with the simple dynamics used here, they at least motivate studying its performance in the situations where large uncertainties on the optimal parameter range are present, or with large climate models where reducing the frequency of its application is desirable to decrease the associated computational cost.

Conclusion and Discussion
This study provides an outlook of the initialization techniques for seasonal-to-decadal predictions using concepts and notations borrowed from the data assimilation context. Our first purpose was to give a comprehensive comparison between full-field and anomaly initialization (FFI and AI), the two main classes of initialization methods used in major climate prediction institutions, with the aim of contributing to the current debate on their respective advantages and drawbacks. The second objective of this investigation was the introduction of two methods to improve the initialization procedures; namely, the least-square initialization (LSI) and the exploring parameter uncertainty (EPU). Following a standard procedure for the analysis update in most data assimilation schemes, LSI operates a least-square fit between observations and model with weights related to the respective assumed level of accuracy. By optimally combining model and observations at the initialization step, the informational content of the observations is propagated to the entire model domain and no longer confined to the locations of the measurement, as in classical FFI or AI. The required model error covariance entering the LSI initialization is estimated using the covariance of the model anomalies taken over a long uninitialized run, in a similar manner to what was done in Smith and Murphy (2007). Using LSI, the model's initial state out of the observational locations is modified according to the structure of the model statistics. This allows to reduce the initial errors in areas otherwise totally unconstrained by the observations and has the potential to mitigate the dynamical shocks resulting from discontinuities caused by pushing the model state to observation values and leaving it unchanged elsewhere.
EPU is an online drift correction method in which the drift caused by the parametric error is estimated using a short-time evolution law and is then removed during the forecast run. It requires the computation of the model equation's first derivative with respect to the parameters assumed to be uncertain (i.e., the model parameter sensitivity), and a hypothesis about the assumed range of possible parameter values. Being based on a short-time approximation, EPU deteriorates when the length of the time interval over which it is applied is longer than the duration of the short-time regime. This duration is certainly connected to the model sensitivity to the specific set of parameters under consideration and is model dependent. The ideal situation would be to implement EPU at every time step. Nevertheless, although the computational demand of EPU is very low, reducing the frequency of its application (by increasing the drift correction interval) would reduce the computational cost.
The comparison between FFI and AI as well as the performances of LSI and EPU have been studied using an idealized coupled climate model based on the classical Lorenz 3 variables system (Lorenz, 1963)). The model, introduced by Peña and Kalnay (2004), possesses three compartments characterized by different timescales and amplitudes, and are taken as proxies for extratropical/tropical atmosphere and the ocean. The low model dimension and its low degree of complexity allowed us to run a large ensemble of trajectories (360 start dates) using a standard hindcast experimental setup over 30 years and monthly start dates.
The use of the data assimilation formalism helped us to rewrite the initialization updates in a general form and highlight specific features such as their sensitivity to observation and model accuracy. The analysis of their error-scaling properties suggests the use of FFI when a good observational network is available and reveals the direct relation of its skill with the observational accuracy. The skill of AI appears, however, mostly related to the model quality, and clear increases of skill can only be expected in coincidence with model upgrades. The numerical results confirm this behavior and reveal that, in contrast to AI, FFI prediction skill appears clearly related to the observational accuracy up to the third forecast year. As opposed to this, AI displays a marked dependence on the model bias with a performance that systematically improves when the model is improved.
We have compared FFI and AI in experiments in which either the full system or each of the model compartments was initialized independently. When the full system is observed the best performances are obtained in FFI and AI, although the former shows clearly better and longer-lasting improvements, with skillful predictions (RMSSS > 20 %) until month 30, whereas in AI the RMSSS falls below this level between months 15 and 20. In the initialization of single compartments, the best performance (with no remarkable differences between FFI and AI), is obtained when the stabler component of the model (the ocean) is initialized. Results show also that with FFI it is possible to have some predictive skill when the most unstable compartment (the extratropics) is observed. This behavior is not mirrored in AI, which supports the idea that the gain in skill is obtained by the efficient initial error reduction of FFI. However, in agreement with results from a state-of-the-art climate model (Smith et al., 2013), AI slightly outperforms FFI after the fourth year.
Finally, FFI was compared with FF-LSI and FFI-EPU. Results show that using LSI improves the performance of FFI in all the situations when only a portion of the system's state is observed. This proves that the performance of LSI is based on an efficient propagation of information from datacovered to data-uncovered areas and to some extent reduces the initial error also far from the observational locations. Estimating the background error covariance matrix using the www.nonlin-processes-geophys.net/21/521/2014/ Nonlin. Processes Geophys., 21, 521-537, 2014 statistics of the anomalies has therefore proved to be a viable choice. These results, in agreement with Smith and Murphy (2007), extend them to the case of the estimation of crosscovariances between different model compartments (atmospheres and ocean). This feature is extremely relevant for the development of coupled data assimilation systems, where the aim is to simultaneously update all coupled subsystems based on observations of one compartment alone. The use of EPU has clearly improved the skill of FFI within the first forecast year; later, the limits of accuracy of the linear and short-time assumptions at the basis of EPU are approached, and only minor advantages over FFI are recorded. Results have also demonstrated the robustness of EPU with respect to the two factors determining its implementation: the length of the drift-correction interval and the accuracy of our knowledge of the actual range of possible parameters.
This study is intended as a proof of concept in which initialization algorithms are analyzed in a simplified context in which the respective advantages and drawbacks could be easily highlighted. Results of the two proposed advanced formulations, LSI and EPU, encourage their application with models, and observational setups, of increasing complexity. The extent to which the conclusions drawn in the present work can be applied to a realistic climate prediction exercise is a central topic of research the authors are currently undertaking. Two main lines of research have been undertaken as follow-up activities. In the context of a state-of-the-art climate model we are studying the effect of initializing different areas using FFI in a multiyear prediction horizon. Similarly, we are investigating the capability of EPU to estimate and deal with the drift caused by a very limited set of coupling parameters, on seasonal timescales. A crucial issue in this type of realistic applications will be the estimate of the drift caused by a parametric error having spatial dependence. We are also committed to extend the analysis presented in this study to a larger set of uncertain parameters and to study the use of LSI in coupled climate models of intermediate complexity.