Controlling instabilities along a 3DVar analysis cycle by assimilating in the unstable subspace: a comparison with the EnKF

A hybrid scheme obtained by combining 3DVar with the Assimilation in the Unstable Subspace (3DVar-AUS) is tested in a QG model, under perfect model conditions, with a fixed observational network, with and without observational noise. The AUS scheme, originally formulated to assimilate adaptive observations, is used here to assimilate the fixed observations that are found in the region of local maxima of BDAS vectors (Bred vectors subject to assimilation), while the remaining observations are assimilated by 3DVar. The performance of the hybrid scheme is compared with that of 3DVar and of an EnKF. The improvement gained by 3DVar-AUS and the EnKF with respect to 3DVar alone is similar in the present model and observational configuration, while 3DVar-AUS outperforms the EnKF during the forecast stage. The 3DVar-AUS algorithm is easy to implement and the results obtained in the idealized conditions of this study encourage further investigation toward an implementation in more realistic contexts.


Introduction
Data assimilation in meteorology and oceanography is experiencing a rapid phase of development, with flourishing of theoretical studies and applications. Traditionally the regular, globally available measurements network has been employed in the process of combining data with model forecasts to obtain the initial condition for model integrations. Recently there has been also a growing interest in adaptive observational systems ( [22]) and satellite data thinning ( [26]).
While the merits of the relatively new ensemble methods are being compared with those of traditional variational methods, the potential benefit of the use of adaptive observations, added to the standard network, is being tested ( [32]).
The present is a follow-up study on the testing of an assimilation scheme developed by the authors ( [38,40,39,7]) and referred to as Assimilation in the Unstable Subspace (AUS). The basic paradigm of AUS is to track and control the analysis-forecast cycle instabilities. As a consequence, AUS has found its most natural applications in an adaptive observation context. Making use of just a few additional observations, properly located at each observing time in order to detect the analysis-forecast cycle instabilities estimated by Breeding on the Data Assimilation System (BDAS) it was possible to drastically reduce the forecast error.
Previous studies have in fact proven the AUS-BDAS system effectiveness in a hierarchy of dynamical systems characterized by different features and complexity, such as the Lorenz 40-variable model ( [23]), the atmospheric quasigeostrophic model used here ( [29]) and a primitive equation ocean model ( [5]). These results support, on the one hand, that BDAS, in view of its ability of tracking the instabilities of the assimilation system, provides a feasible targeting strategy; on the other hand that AUS is an efficient, computationally affordable, relatively easy to implement method to assimilate the adaptively located observations. Moreover, they show that, even dealing with high-dimensional system, an efficient error control and an accurate state estimate can be obtained by monitoring only a reduced number of unstable directions and by properly using a limited number of observations. By allowing ad hoc observations of the estimated instabilities, the deployment of an adaptive network naturally enhances the efficiency of AUS. In fact, targeting the spatial unstable structures improves the estimate of the amplitude of the unstable components of the error, thus enhancing the stabilizing effect. One natural question is therefore how AUS-BDAS performs with a fixed, standard, observational network. As long as the spatial and temporal observational coverage of the standard network is dense enough, the instabilities can be detected as they move throughout the domain, although a reduced efficiency may be expected with respect to the adaptive observations case. [39] investigated the use of fixed "satellite" surface elevation observations in a primitive equation oceanic model and pointed out the weakness of such an observational network.
A stable error reduction could not be achieved by using only surface observations since dominating instabilities developed in deeper model layers before they could be detected at the surface. Instead, an efficient error reduction and a drastic improvement of the assimilation performance were obtained by means of vertical profiles adaptively located by BDAS and assimilated by AUS.
We investigate the best use of the AUS-BDAS method in an atmospheric context where only a network of fixed soundings is available.
Errors, in an operational data assimilation cycle, are of various origin. Even if it is assumed that errors in the description of the system dynamics (the model error) are not present, noise is introduced at different stages of the observing and assimilation procedure.
It includes, among others, the measurement noise (coming also, for instance, by deficiencies in the retrieval procedures of indirect measurements), representativity errors and errors due to imperfections in the assimilation procedure itself. Moreover, if the instabilities cannot be controlled by adequate (ad hoc) observations, nonlinearities become important.
Thus, in the presence of large unstructured errors due to observation and system noise and to strong nonlinearities that can be expected to be important in a fixed observational setting, a failure of the unstable subspace paradigm of AUS is expected. These arguments suggest and motivate the combined use of AUS and a stationary assimilation method.
In the present study we propose an algorithm composed of a least-square based data assimilation scheme, the 3DVar, in combination with AUS, hereafter referred to as 3DVar-AUS. The added value of AUS over a 3DVar analysis cycle is investigated with a fixed observational network chosen to simulate a realistic observation coverage, and observation system simulation experiments are performed with both perfect and noisy observations. Given the properties of the available observing network, a practical relevant question is how AUS compares with other assimilation methods. For comparison, an ensemble Kalman filter (EnKF) ( [12,13]), has been implemented and optimized for the model and observational network used here. The EnKF is considered among the most promising ensemble-based sequential data assimilation schemes and therefore is an ideal candidate for the comparison with the proposed 3DVar-AUS scheme.
This article is organized as follows. In Section 2 the details of the model and of the observational network are given. In Section 3 the three data assimilation schemes, objects of the comparison, are briefly described along with the specific implementation choices adopted here; particular emphasis is given to the proposed 3DVar-AUS. Results are described in Section 4, while the final summary and discussion can be found in Section 5.

2
Model and observation network configuration All the experiments described in this study are performed by making use of an atmospheric numerical model based on the quasigeostrophic equations in a periodic channel ( [29] ([40, 7]) as well as in a number of other data assimilation studies ( [25,8,9]).
An analysis of its dynamical properties can be found in [30].
The forward model integration evolves from the analysis state at a given time t k to the forecast state at time t k+1 = t k + τ where k = 0, 1, 2..., and τ is the assimilation interval of a sequential assimilation scheme: where x f and x a indicate the forecast and analysis states respectively.
The observing network consists of a fixed set of observations simulating synoptic profiles located at model gridpoints and measuring model variables, (potential vorticity, PV, at the five inner levels and potential temperature, PT, at top and bottom), at each level.
In a previous work with the same model ( [7]), AUS was used to assimilate a single adaptive observation in a large data-void area, and was combined with a 3DVar that assimilated all the fixed observations over a data-rich area.
In the present study, we intentionally work with a fixed network of observations to investigate the potential usefulness of AUS in the absence of adaptive observations that naturally enhance its efficiency. The observational network used here ( Fig. 1) alternates data-rich/data-void areas as it is typical of a real observation network: the distribution of M = 125 soundings is obtained through a random selection procedure.
Observing system simulation experiments are performed: a given trajectory, solution of the model equations, is taken to represent the true atmospheric evolution. At each analysis time, every τ = 12 hours, observations are taken by sampling the "true trajectory" at the observation locations. Following the widely used notation ( [20]), observation values are stored as components of the observation vector y o : where ε(t k ) represents the additive observation error, assumed to be white in time, Gaussian distributed with covariance R and H is the observation operator, estimating the observed variables from the model state. Since the observations are located at model gridpoints and observe model variables, the observation operator is inherently linear and the notation H is used hereafter.
Observation errors are assumed to be correlated in the vertical direction only and uncorrelated for different soundings. The observational Root-Mean-Square (RMS) error at each level is set to 10 % of the system's natural variability while, as in Morss (1999), the covariance between different levels is expressed in terms of the respective variances and of a vertical correlation function (Eq. 3.19 in [4]), depending on the vertical distance only. With these assumptions the observation error covariance matrix R takes the simple form of a diagonal block matrix, with a 7 × 7 square matrix in each of the M diagonal blocks.
Following [24] and [17], observation errors are randomly generated, and added to the true values, by sampling a Gaussian distribution N (0, R) with zero mean and consistent with the assumed observation error statistics.

Data Assimilation Algorithms
The analysis state at an arbitrary time t k is obtained as a linear combination of a forecast state, taken as a background field, with the observations: where d = y o − Hx f is the innovation vector, while the expression of the gain matrix K that minimizes the analysis error variance is the Kalman gain ( [21]): where P f represents the forecast error covariance matrix. In Eqs. (3) and (4) all vectors and matrices refer to the same time t k , and T indicates matrix transposition.

3-Dimensional Variational Assimilation
The three dimensional variational algorithm, 3DVar, used here has been developed on the basis of that described in [24]. The forecast error covariance is estimated by the stationary matrix P f 3DV ar = B (background error covariance), so that the analysis update is (see e.g [33]): corresponding to the minimum of the 3D-Var objective function for a linear observation operator H. The matrix B, that is kept constant throughout the analysis cycle, was statistically estimated by [24] and has been multiplied here by a scalar coefficient in order to optimize it for the present observational setting. The optimal value, chosen by minimizing the analysis error with the present network of 125 noisy profiles ( Fig.1) is 0.97.
The initial conditions for all the experiments described in the text is the final state of a 1-year-long 3DVar analysis cycle initialized by randomly perturbing the truth state.

Ensemble Kalman Filter
The ensemble Kalman Filter, EnKF, used here is based on that described in Descamps and Talagrand (2007) and used for the study and intercomparison of ensemble initialization techniques. It basically follows the approach given by [13,14].
At each assimilation time t k , a set of perturbed observations is generated, and then used for the ensemble analysis update, by randomly perturbing the actual observation vector: where ε i are independent realizations of the observational noise N (0, R).
At the analysis time the ensemble of forecast states is updated: where the gain matrix K EnKF is obtained by setting P f = P f EnKF in Eq. (4). The ensemble of analysis states is then used to initialize N full nonlinear model forecasts.
In all the experiments described here, N = 30 members are used and the very first analysis ensemble is generated by adding random noise to the reference initial condition.
This noise is drawn from a Gaussian distribution with zero mean and variance equal to the average analysis error variance of the 3DVar experiment.
In order to avoid ensemble collapse, prevent filter divergence and optimize the EnKF performance, following [10], forecast error covariance inflation ( [2]) and localization ( [18]) have been adopted. The covariance inflation is obtained by multiplying the matrix P f by a scalar coefficient before using it in the analysis: where α is the (tunable) inflation factor. As in Houtekamer and Mitchell (2001), the 5th order [15] function is used to localize the P f EnKF , the main parameter being the zerocorrelation distance, d 0 .
Both α and d 0 have been optimized, over a 60 days analysis cycle, for the specific model and noisy observational network used here. The best results are obtained with d 0 = 3000 Km and a 7 % of inflation (α = 0.07). Figure 2 shows the normalized RMS analysis error as a function of the inflation factor (fixed zero correlation distance d 0 = 3000 Km), and as a function of the zero correlation distance (fixed inflation factor α = 0.07) in the inset.

Assimilation in the Unstable Subspace
The mathematical formulation of AUS is briefly reported in the following, along with the description of the setup of its application in the present model and fixed observations configuration. A full description of its mathematical and theoretical foundations can be found in [38], [39] and [7].
The AUS method is basically aimed at tracking, and controlling, the instabilities of the (observationally forced) analysis-forecast cyclic system in view of their prominent role in the error evolution and their impact on the overall performance of the data assimilation scheme. To this end, when observations are available, an estimate of the analysis-forecast cycle unstable subspace is used to confine the analysis correction within this subspace, thus maximizing the effect of the assimilation where it is expected to be more necessary.
By combining Eqs. (1) and (3), the evolution equation for the analysis state, the analysis-forecast cycle, is obtained: Apart from the presence of a linearized observation operator H, Eq. (10) is the governing equation of most sequential data assimilation schemes.
The analysis-forecast cycle perturbation equations are given by: where M, the Jacobian matrix of M, represents the tangent linear model evolution. In principle, from Eq. (11), it is possible to estimate the unstable manifold of the forecastanalysis system.
After storing the unstable vectors as the columns of a I × N matrix E, and confining the analysis increment in the unstable subspace, the AUS analysis becomes: where: or, equivalently: where, as before, all terms refer to the same time t k and where Γ is a N × N positive definite matrix, representing the forecast error covariance matrix in the unstable subspace: Suppose that a single mode e is detected by a set of observations, then the AUS analysis reads: where γ 2 , the forecast error variance in the e direction, can be estimated from innovations.
From Eq. (16) we see that the analysis increment vector has the direction of e, and the amplitude that best-fits the observations. This means that, in physical space, the difference between the analysis and the forecast state has the three-dimensional structure of the unstable perturbation.
In practice, the unstable subspace of the data-assimilation system is estimated by means of an extension of the Breeding technique ( [35,36]), known as Breeding on the Data Assimilation System (BDAS) ( [38]), which naturally incorporates the observational forcing in the perturbations dynamics.
Equation (11)  Furthermore, it is impractical to estimate the whole unstable subspace; this would require either a recursive orthonormalization to be applied to the set of BDAS vectors, according to the standard technique of [3] or the use of computationally demanding techniques such as those described [37] and by [45], which allow the estimation of the Lyapunov vectors. Instead, in this as in previous applications of AUS, a refresh procedure is introduced in the breeding cycle in order to systematically explore the unstable subspace of the system. In fact, as discussed in [38] and [39], the more accurately the forecast error projection on a particular unstable direction is estimated, the more dominant will become the error projection on the complementary subspace.
Once a breeding cycle has been established, the BDAS modes and the forecast error will approximately have a structure given by a different linear combinations of independent unstable vectors of the underlying dynamics.
Local structures appearing in the global BDAS modes will be positively or negatively correlated with similar structures in the forecast error. It is therefore necessary to localize these structures before entering them in an assimilation expression of type (16). In addition, the localized structure needs to be observed, i.e. He should be non-zero and large enough to provide a reliable estimate of the analysis increment, whose sign is correctly determined by the innovation 1 . It is important to notice that, by applying only a horizontal localizing function, the vertical structure of the perturbation is preserved in the definition of the analysis increment.
While EnKF-like methods are based on an estimate of the forecast error covariance built from an ensemble of trajectories, the identification of unstable directions of the given control solution, which provide the three-dimensional structure of the analysis increment whose sign is determined by the innovation, is the distinctive feature of AUS.
EnKFs are based on the traditional Kalman filter equations and use a Monte Carlo approach to estimate the forecast error covariance matrix entering the analysis update, therefore they have to face the consequences of sampling errors with or without perturbed observations ( [43]) and filter divergence. As a consequence, several ensemble filters formulations have been proposed to overcome these difficulties (for a review see [1], [13], [34]). Similarly, the AUS approach, that relies on the estimate of the forecast error projection on the unstable subspace, faces the practical difficulties of accurately estimating the unstable directions and the amplitude of the forecast error on each of them. While the theoretical premises of the two approaches are different but clear, when it comes to the implementation all schemes are subject to several practical choices that make quantitative comparisons rather subjective.
An adaptive observations approach was followed by [39] and by [7] in an oceanic and atmospheric model applications respectively. In fact, the AUS approach of tracking and controlling the data assimilation system instabilities is most efficient when observations are indeed available at the moment and in the region where an instability develops. In other words, if the set of unstable vectors E can effectively be detected (that is to say if the components of HE are large enough), then E can be used in Eq. (12) to update the background state. When, as in the present application, only fixed synoptic observations are available, the efficiency of the AUS scheme is expected to be determined by the spacing and frequency of the observational network in relation to the typical spatial patterns and growing times of the system's instabilities.

Implementation of 3DVar-AUS
Several approaches have been proposed to combine 3DVar with ensemble-based filters.
These hybrid Ensemble/3DVar analysis schemes introduce information on flow-dependent instabilities by properly weighting the static covariance with the ensemble estimated forecast error covariance ( [16,41]). [11] investigated the robustness of hybrid schemes to model error: the model error modifies the optimal values of the weights to be given to the ensemble based and static covariance matrices, significantly reducing the weights corresponding to the former.
In the present study, we introduce a hybrid scheme that uses the flow-dependent covariance and the 3DVar covariance separately to assimilate different observations. The AUS assimilation is applied first to reduce the forecast error component in the model's unstable subspace estimated by the BDAS modes, using those observations that are able to detect them. The residual forecast error is assumed to be unstructured and a static 3DVar covariance is used to assimilate the remaining observations. Possibly, error associated with unstable structures that were not present in the BDAS modes, or whose amplitude could not be reliably estimated with the available observations may be still present in the analysis and subsequent forecast: hopefully they will be accounted for at successive assimilation times.
The 3DVar-AUS forecast-assimilation and breeding cycle is implemented through the following recursive steps: Thus, in all the experiments described below, a single N = 1 BDAS mode is used at each assimilation time while the breeding time interval is set to 6 days. This latter value was chosen by tuning to minimize the average analysis error (see also [7]).
Expression (17) 17) becomes: Once the AUS analysis is completed, the analysis field is used as the background in the 3DVar analysis to assimilate the remaining M − N o i=1 M i observations (step 8).
As mentioned above, a refresh procedure (steps 1 and 9) is used in the implementation of BDAS ( [38]). At each assimilation time, after the analysis update, the set of BDAS modes is discarded after being used in the assimilation. A new one is introduced, which starts to undergo a new breeding procedure: after completion of its breeding period it will be used in an assimilation step. Although this procedure increases the computational cost of BDAS (with respect to evolving the same set of perturbed states), it is very beneficial to efficiently span the trajectory's unstable subspace and to provide a reliable estimate of the data assimilation system instabilities. The refresh was also used in [7] in the context of the same QG model used here, while a trade-off procedure between keeping and discarding the N (=6 in that case) BDAS modes altogether was used in [39] in the context of a primitive equation ocean model.
In the experiments described here, the characteristic distance is chosen, after optimization, to be d = 1500 Km. Also, a limit to the total number of the extracted assimilating structures is set,N = 20. Furthermore, the size of the searching box l and the coefficient

Experiments with perfect observations
We first compare the performance of the algorithms in a perfect observation setting.
Clearly real observations are not perfect. Anyhow, the idealized perfect observations setting is related to the problem of observability of the assimilation system and provides insight on the ability of various methods to track the instabilities of the system. The use of perfect observations with AUS allows for an accurate evaluation of the analysis increment and potentially reduces to zero the analysis error projection along the unstable structures used in the analysis update (18); the relation between the observability condition and AUS was discussed by [38]. Furthermore, the perturbed observation case provides an upper limit of performance of the assimilation schemes under comparison.   The local character of the atmospheric instabilities are described in [28] and have been efficiently exploited in AUS as well as in most of the ensemble prediction and data assimilation applications, as for instance the Local Ensemble Kalman Filter ( [27,31,9]) and the Local Ensemble Transform Kalman Filter ( [19]).

Experiments with noisy observations
Although very instructive and useful in providing insight into possible advantages and drawbacks of a given algorithm, the use of perfect observation is not realistic: in real applications, the assimilation of noisy observations introduces unstructured random errors.
Therefore, we now turn to noisy observations experiments; observations are generated according to the procedure described in Sect. 2. A set of observation realizations, Eq. (7), is used to update the ensemble forecast in the EnKF by means of Eq. (8), while Eq. (18) is used for the AUS assimilation in the 3DVar-AUS experiment. The error maximum located between 12 -20 longitude, 6 -16 latitude, "ignored" by the assimilation, is still present in the forecast error field.
The presence of the strong positive error maximum located, as at the previous analysis time, in the eastern side of the channel, at mid latitude, suggests that the analysis correction at day 204, 00UTC, was not sufficient to eliminate it. Now, at day 204, 12UTC, although its shape is again well captured by the current BDAS mode, the AUS analysis is able only to reduce it. Notice that, a very similar pattern (with opposite sign), manifestation of the same instability, is present in the EnKF forecast error. Similar considerations apply to the last instant shown, day 205, 00UTC.
For all three algorithms, the mid-level PV increments shown in Figs give a comparable contribution to error reduction. This gives further evidence that when instabilities are present the few observations used by AUS (approximately 1% of the total available observations) are almost as effective as the remaining ones. In any case, while the 3DVar always reduces the corresponding background error of approximately the same amount, the effect of AUS is more intermittent.
Finally, Fig. 11 shows the time and domain averaged analysis and forecast errors as a function of the forecast lead time. These are all deterministic forecasts; in particular the EnKF forecast is initialized from the EnKF analysis, mean among analysis ensemble members. It is interesting to note that although the EnKF and 3DVar-AUS average analysis error is comparable, the 12 hours forecast error growth is much more rapid in the former case, while it has approximately the same rate afterward.
With the specific setup choices adopted here, the computing time required to complete the 2-year 3DVar-AUS assimilation experiment described in the text, was about 60 % of the time required by the EnKF. It is remarkable that, despite the fact that 3DVar-AUS and the EnKF attain a comparable average analysis error (6.8 % and 6.6 % of the system's natural variability, respectively) the 12 hours forecast error in the EnKF is larger than in the 3DVar-AUS (9.9 % and 7.3 %, respectively) and its 24 hours forecast skill is somewhere between the 36 hours and the 48 hours forecast skill of the 3DVar-AUS..

Summary
This result confirms the efficiency of AUS in reducing the flow instabilities responsible for the error growth even in a fixed observations setting.

Conclusions
The results of the present study lead to the following interpretation: the forecast error shows very intense and localized structures that are similar in the three assimilation experiments and in particular between the 3DVar-AUS and the EnKF: accordingly, a high correlation in the time evolution of the RMS error of the two schemes is observed. As indicated by these similarities, the local maxima of the forecast error are presumably concentrated in regions where dynamical instabilities, compatible with the observational forcing, develop along the analysis-forecast cycle.
With a fixed observational network, observations are assimilated by AUS only when they are found in correspondence to the local maxima of the unstable structures. Therefore it is necessary to use 3DVar to assimilate the remaining observations. Three dimensional error structures, provide the corresponding AUS analysis correction, whose amplitude and sign are determined from the innovation: hence its ability to reduce the forecast error also at levels other than the observed ones.
In real world applications, the employment of AUS is expected to give an improvement over a 3DVar (or Optimal Interpolation) every time an unstable structure is observed by either fixed or adaptive observations. When the latter are available and there is no need to wait for the instability to travel into observed regions, the expected improvement in using AUS is enhanced.
Also the EnKF takes advantage of the observations that happen to be (or are intentionally located) in regions where instabilities develop and the present fixed observations experiments prove that the observations are exploited by the two methods with similar efficiency. The EnKF obtains the same goal of capturing and controlling the instabilities, the different approach of the two methods being exemplified by the following idealized example where a single unstable structure is present in the forecast error. In this simple case, the EnKF would identify the unstable structure as difference of the members (two members would be sufficient for this purpose) from the ensemble mean. To estimate the amplitude of the correction, however, a sufficiently large number of representative members would still be necessary to provide a reliable estimate of the associated error variance.
In contrast, AUS estimates the direction as difference of a single unstable perturbation from the control and estimates the amplitude of the correction from innovation.
The present result, where EnKF and 3DVar-AUS have a similar analysis performance, seems to indicate that both methods were able to efficiently exploit all the available observations useful to control the flow-dependent instabilities. A point in favor of the 3DVar-AUS scheme is its better performance in the forecast stage.
Apart from the differences among the different approaches proposed in previous studies (see [42] and references therein) two results deserve to be mentioned here. Usually, in these hybrid schemes, the forecast error covariance matrix is obtained by combining a static covariance with an ensemble based one, through a tunable scalar weight. Improvements over a 3DVar conveyed by the flow-dependent covariance are less important when a dense observational network is available ( [16]). Moreover the presence of model error significantly reduces the optimal value of the weight to be given to the flow-dependent covariance ( [11]).
Based on these works and on the results of the present study we draw the following conclusion.
To hybridize an ensemble based scheme or AUS with a static covariance might turn to be more convenient when instabilities loose some of their significance, either because the observational network is particularly dense ( [44]) or if the model error is such that the model does not provide an adequate representation of the true unstable modes ( [11]). In view of the latter consideration, while the hybrid approach proposed in the present study may turn out to be useful also in the presence of model error, further work is needed to address this problem.
Because the performance of different schemes depends upon many different factors, including model error and the number and distribution of observations, it remains to be seen which method is more efficient in a particular operational setting. Therefore a quantitative comparison relevant for operational purposes is out of the scope of the present study.
However, given the ubiquity of the role of instabilities in degrading the analysis accuracy, we believe that AUS may turn out to be useful in those circumstances in which the observational network, and in particular an adaptive one, allows their detection.

Acknowledgements
where , T represents the average over an appropriate time interval T and D is a scalar coefficient to be tuned. After optimization, these values are set to D = 1.35 and T = 8 days.
Based on Eq. (18), if He approaches zero the analysis correction approaches infinity.
However because only observations in proximity of maxima in the BDAS mode are chosen and in view of the selection procedure described in Sect. 3.3 (point (b)) this circumstance is never encountered.