## Abstract

The potential use of chaos synchronization techniques in data assimilation for numerical weather prediction models is explored by coupling a Lorenz three-variable system that represents “truth” to another that represents “the model.” By adding realistic “noise” to observations of the master system, an optimal value of the coupling strength was clearly identifiable. Coupling only the *y* variable yielded the best results for a wide range of higher coupling strengths. Coupling along dynamically chosen directions identified by either singular or bred vectors could improve upon simpler chaos synchronization schemes. Generalized synchronization (with the parameter *r* of the slave system different from that of the master) could be easily achieved, as indicated by the synchronization of two identical slave systems coupled to the same master, but the slaves only provided partial information about regime changes in the master. A comparison with a standard data assimilation technique, three-dimensional variational analysis (3DVAR), demonstrated that this scheme is slightly more effective in producing an accurate analysis than the simpler synchronization scheme. Higher growth rates of bred vectors from both the master and the slave anticipated the location and size of error spikes in both 3DVAR and synchronization. With less frequent observations, synchronization using time-interpolated observational increments was competitive with 3DVAR. Adaptive synchronization, with a coupling parameter proportional to the bred vector growth rate, was successful in reducing episodes of large error growth. These results suggest that a hybrid chaos synchronization–data assimilation approach may provide an avenue to improve and extend the period for accurate weather prediction.

## 1. Introduction

Synchronization of loosely coupled oscillators is not uncommon in nature. Examples include the synchronized flashing of Asian fireflies, the emergence of scores of cicadas within a few days of each other every 17 years, and the simultaneous firing of brain neurons to create memories (Strogatz 2003). The nonperiodic solutions of two nonlinear deterministic systems with sensitive dependence on initial conditions can be made to synchronize if the coupled system has negative conditional Lyapunov exponents; that is, the average growth of small perturbations in the distance between the two coupled systems is negative (Pecora and Carroll 1990). Two systems with state vectors **u** and **v** are synchronized if the distance *d* = ‖**u**(*t*) − **v**(*t*)‖ → 0 as *t* → ∞ (Zhong et al. 2001). The synchronization of coupled chaotic subsystems has been described in electronic circuits (Pecora et al. 1997), partial differential equations (Kocarev et al. 1997), and, recently, atmospheric models (Duane and Tribbia 2001).

While these previous investigations have examined synchronization in a pair of systems, each a part of some objective reality, we propose here a different application of synchronization. One system is taken to be real, as before, but the other system is merely a model of the first. According to Brown and Rulkov (1997): “In practice one could imagine coupling the output from a physical system to a model of the system using coupling that is guaranteed to result in stable synchronization.” Their suggestion motivates the present comparison of synchronization with conventional methods of data assimilation.

Data assimilation is the process of using all available information, including short-range model forecasts and observations, to estimate the current state of a system as accurately as possible. A major application of data assimilation has been to produce analyses of the state of the atmosphere to initialize a numerical weather prediction (NWP) model (e.g., Kalnay 2003). Like chaos synchronization, data assimilation steers a prediction model toward behaving like the atmosphere by periodically forcing it with atmospheric observations. In fact, the nudging method suggested for data assimilation by Hoke and Anthes (1976), which has been used for a number of data assimilation studies (see, e.g., Stauffer and Seaman 1990; Stauffer et al. 1994; Bao and Errico 1997), is very similar to schemes used for chaos synchronization.

Indeed, the specific examples of synchronization schemes described in this paper are all forms of nudging. However, these examples illustrate a more general type of coupling of truth to model that extends nudging in a natural way. Chaos synchronization and traditional data assimilation differ in that the data assimilation objective is the minimization of the current analysis errors, which depend on observational and forecast errors, while chaos synchronization instead focuses on the minimization of the future divergence of coupled subsystems. In data assimilation the error covariances of the forecasts and of the observations of the “master” system are used to weigh more heavily the more accurate data. In contrast, typical chaos synchronization studies modify the number of variables coupled and the overall coupling scheme since certain variables are more effective in promoting convergence than others (Brown and Rulkov 1997). The data assimilation application resembles a previously suggested application of synchronization to the tracking of a target system by a model system (So et al. 1994). No direct comparisons of data assimilation and synchronization with the same system have been reported so far.

While the synchronization approach and standard data assimilation algorithms are not a priori equivalent and are related in a complex way, there is much to be learned from a simple empirical comparison. To examine whether unidirectional coupling techniques of chaos synchronization can improve upon current data assimilation schemes in producing accurate forecasts, we carried out experiments with the Lorenz (1963) three-variable model:

The Lorenz model provides a practical test case with qualitatively realistic properties. Atmospheric behavior involving barotropic and baroclinic instabilities is considered somewhat analogous to Lorenz model behavior because of the exponential instability of the model’s trajectories and its abrupt regime changes (e.g., Miller et al. 1994). The standard parameter values for the Lorenz attractor were used: *σ* = 10, *b* = 8/3, and *r* = 28 (Lorenz 1963), except in section 3 (generalized synchronization). The model was integrated, as is customary, with a fourth-order Runge–Kutta numerical scheme and a time step of 0.01. We used two sets of the Lorenz equations starting with different initial conditions: the first representing the real atmosphere (master), from which we extracted noisy “observations,” and the second representing the model (“slave”), which we tried to synchronize to the master. In section 2, we discuss our methodology and study the identical synchronization of an identical slave model, comparing different coupling techniques including one, two, or three dimensions, synchronization along bred vectors and singular vectors, and different coupling strengths. In section 3, we explore synchronization of nonidentical systems, a phenomenon known as generalized synchronization. In section 4, we describe two data assimilation methods, three-dimensional variational analysis (3DVAR) and extended Kalman filter (EKF), and compare their effectiveness with the much simpler synchronization approach and with an adaptive generalization. Section 5 is a summary and discussion.

## 2. Synchronization experiments with the Lorenz model

In these experiments the master Lorenz model followed Eq. (1), and the slave system was nudged toward values obtained from the master run:

where *k* is the coupling constant (i.e., *k*^{−1} is the time scale of the coupling), the subscript *S* represents the slave system, (*x*, *y*, *z*) are the “observations” (which may have observational errors) obtained from the master run [Eq. (1)] and is equal to one, if the observation in the *i* direction was used, and zero otherwise. The two systems were integrated in time, and an initial transient period of 500 steps was discarded in every case. In agreement with previous authors, we found synchronization between master and slave under a wide variety of conditions. In particular, synchronization was often found to occur when only one of the Δ_{i} was nonzero. This phenomenon has attracted attention because the two chaotic systems are each effectively unpredictable, yet coupling them through a single variable results in a predictable relationship, despite the sensitive dependence on the initial conditions of all three variables. Formally, synchronization occurs because the conditional Lyapunov exponents of the uncoupled slave subsystems (measuring the growth of small perturbations in the distance between the two coupled systems) are negative (Pecora and Carroll 1990; Lui et al. 1999).

### a. Simple nudging

In the first type of experiments, we used simple nudging, also known as diffusive coupling in the synchronization literature (Pecora et al. 1997). We first used Eq. (2) with different coupling constants and different combinations of master system variables assumed to be perfectly observed, as is customary in synchronization experiments. One might measure the time required for synchronization to determine the efficacy of different observed variables or coupling strengths. However, since the synchronization time depends on the initial conditions, results are strongly subject to sampling and not satisfactory for detailed comparisons. Since in a real system the exact master observations are not available, we opted instead for using observations of the master to which we added random observational errors with zero mean and standard deviation of 2. We then used simple nudging, as in Eq. (2), to nudge the slave solution toward the observations with errors, with the goal of minimizing the root-mean-square error between slave and master. Since increasing the coupling strength *k* also increases the amount of noise that is introduced, the use of noise helps to determine the optimal value of *k*. The optimal value of *k* is roughly the minimum value that would achieve synchronization in the noise-free case.

An example of identical synchronization with random observation errors is shown in Fig. 1. The results obtained are summarized in Fig. 2, showing the rms difference between the master and the slave over the last 500 steps of a 3500-step run, computed for all the integer values of *k* between 1 and 100. Much larger values eventually led to computational instability as the Courant–Friedrichs–Lewy condition was violated (*k*Δ*t* > 1). We found that, when the difference between the master and the slave is less than the observational error, the two systems synchronize perfectly when coupled without observational errors, whereas larger differences indicate that synchronization does not take place or is delayed for a very long time in the error-free case. As with all experiments described in this paper, the integration time was chosen to be long enough so that the results were statistically stationary. Experiments were repeated with a few different sets of initial conditions, but a principle of ergodicity effectively allowed exploration of such variations simply by making the run time long enough.

Figure 2 shows that the effectiveness in synchronizing with simple nudging depends on the coupling strength and the variables chosen for observation. In agreement with Pecora and Carroll (1990), we found that observing *z* did not lead to synchronization, whereas observing *x* or *y* did, but we also found that synchronizing with *y* observations is more efficient than with *x* (see Fig. 2). The relative efficiency of synchronization depends on the coupling strength. Nudging using all three observations (*x, y, z*) or two observations (*x, y*) is more effective at low values of *k* but, surprisingly, less effective than observing only *y* for larger values of *k*. For *k* = 15 or larger, the error in nudging with *y* alone is smaller than the errors using both *y* and *z*. Similarly, adding *z* observations improves synchronization with *x* and *y* for *k* < 8 but makes it worse for larger values of *k*. This implies that, when using a low coupling strength, including information from a less effective variable will help to synchronize two systems. However, it will degrade the results with a strong coupling since additional error is needlessly introduced. A heuristic explanation for the greater effectiveness of some variables may be that variables containing information about the current regime of the Lorenz system (*x* and *y*, but not *z*) are useful for synchronization.

The fact that synchronization takes place over a wide range of coupling strengths *k*, from 1 to over 100 (much larger values led to computational instability), helps explain why nudging is in practice effective even when using nudging time scales shorter than those present in the rest of the dynamics, violating a condition suggested by Hoke and Anthes (1976). Overall, the optimal coupling value was around *k =* 10 in the results presented in Fig. 2, where noisy observations of the master were used at every time step by the slave system. However, when the observations were offered less frequently (once every 8 steps or every 25 steps), the optimal value increased to about *k* = 30–100 (not shown).

Although the minimal value required for identical synchronization was *k* = 1, observing all three variables *x, y, z*, it was larger for fewer observed variables, and no identical synchronization was achieved with *k =* 0.5. (By contrast, generalized synchronization of two different systems did take place for this value, as discussed in section 3.) The manner in which synchronization breaks down for small *k* is interesting. As *k* is lowered below the synchronization threshold, the synchronized trajectories are interrupted by increasingly frequent bursts of desynchronization. This bursting behavior is an example of the more general phenomenon of on–off intermittency, known to occur when the state space of a system contains a dynamically invariant manifold (i.e., once the system is on the manifold it remains there) that in turn contains an attractor. With on–off intermittency, the system alternates chaotically, with unpredictable transition times, between two distinct phases: In the on phase the trajectories hug the invariant manifold very closely; in the off phase they burst away. In the present case the invariant synchronization manifold is the three-dimensional hyperplane in which the two systems are in the same state. This manifold is invariant because exactly synchronized subsystems remain synchronized. The attractor within the manifold is the ordinary Lorenz attractor, but defined in terms of the identical states of both subsystems. If the attractor loses stability against transverse (i.e., desynchronizing) perturbations as parameters such as coupling strength are varied, but remains an attractor within the manifold, as it does in the present case since it is still effective for each Lorenz system, then bursting ensues (Platt et al. 1993; Ott and Sommerer 1994). In a situation where the bursting depends on the presence of noise, the behavior is known as bubbling (Ashwin et al. 1994). Bubbling (see Fig. 3a) is also familiar in the data assimilation context, as will be discussed in section 4. [Intermittency in the absence of noise resembles Fig. 3a, but with rms error vanishing almost completely except during the bursting phase, Ott and Sommerer (1994).]

### b. Nudging in dynamically chosen directions

The most general coupling of slave to master that extends the nudging form in Eq. (2), and might be useful for data assimilation/synchronization, could be written as

where we have introduced vector notation, *F* describes the dynamics of a single (master or slave) system, as given for example by the right-hand side of Eq. (1), and 𝗖 is a matrix of coupling coefficients {e.g., for (2) 𝗖 = [(*kδ _{x}* 0 0), (0

*kδ*0), (0 0

_{y}*kδ*)]}, which may vary in time. Equation (3) can only be regarded as nudging if 𝗖 is a multiple of the identity or a multiple of some matrix that projects onto a subspace; that is, 𝗖 =

_{z}*c*𝗠 for a constant

*c*and a matrix 𝗠 such that 𝗠

^{2}= 𝗠. For more general matrices 𝗖 (such as those used in 3DVAR or EKF) the assimilation scheme defined by Eq. (3) is not nudging but a nontrivial extension of the nudging method.

We compared the efficacy of nudging in a single constant direction, such as *x* or *y*, with nudging in one evolving direction **n** given by either bred vectors or singular vectors, two types of dynamical vectors used as perturbations in ensemble numerical weather forecasting (Toth and Kalnay 1993; Molteni and Palmer 1993). Breeding is a nonlinear generalization of the method used to obtain Lyapunov vectors. A nonlinear system with an initial small perturbation is integrated in time *m* time steps after which the difference between the perturbed and the unperturbed solution is rescaled to its initial size, added to the unperturbed solution, and the process is repeated. The resulting differences between the perturbed and the original system are known as bred vectors (BV), which converge to the leading local Lyapunov vectors (LV) when the size of the initial perturbation is sufficiently small. In our experiments, breeding is done with the slave system including the forcing term, since the master is assumed to be unknown. The leading singular vector (SV) is the eigenvector corresponding to the largest eigenvalue of matrix 𝗠^{T}𝗠, where 𝗠 is the linear tangent model corresponding to the nonlinear model integrated for *m* time steps, and 𝗠^{T} is its transpose or adjoint model (see the appendix). As with the BVs, the tangent linear model is constructed based on the slave system, including the linear coupling term, and likewise for the adjoint model 𝗠^{T}. Both bred vectors and singular vectors are used in ensemble NWP to represent the fast-growing errors present in the analysis. The leading SVs are the fastest-growing perturbations when starting from white noise, whereas the BVs (or leading LVs) are the fastest perturbations to which all other perturbations converge. Leading SVs describe the directions in which future error will grow most rapidly, while BVs describe the most likely directions of present error, based on past growth. SVs have been used at the European Centre for Medium-Range Weather Forecasts (ECMWF) and bred vectors at the National Centers for Environmental Prediction (NCEP) to create initial perturbations for ensemble forecasting (Molteni et al. 1996; Toth and Kalnay 1997). Therefore, both of these types of perturbations would seem to be plausible candidates as effective directions for synchronization (Junge and Parlitz 2001).

Coupling in any given direction **n** is given by Eq. (3) with the matrix 𝗖 = *k***n** × **n**, where multiplication denotes the outer product (it is easily shown that 𝗠 = **n** × **n** satisfies 𝗠^{2} = 𝗠, so 𝗠 is a projection matrix). That is, nudging is performed as

where the components of **n** may vary in time. Expanding the vector **n** = (*n*_{x }*n*_{y }*n*_{z})^{T} the slave system is governed by

The experiments in this section were performed with the specific choices, **n** = BV/||BV|| or **n** = SV/||SV||, where BV and SV are a bred vector and a singular vector, respectively. We also tested the impact of defining **n** as the horizontal projection (*x, y*) of the bred or singular vector direction. In every case, **n** was renormalized to unit length.

In Fig. 4a, we compare the results of nudging in the three variables with the BV, the leading SV, and with simple nudging. For coupling strength lower that *k* = 22, simple nudging gives lower rms error than using either BV or SV. Beyond this coupling strength, the SV provides the best nudging direction to constrain the slave solution toward the master system. However, we also orthogonalized the bred vectors (Annan 2004) and defined the nudging direction as the leading eigenvector of 𝗯𝗯^{T}, where 𝗯 denotes the matrix containing a few (three) unrescaled bred vectors in each column. The orthogonalized BV were very successful in dynamically selecting the most effective nudging direction. Since including the *z* direction was found to be harmful in the BV nudging, we repeated these experiments using only the *x* and *y* components (see Fig. 4b). We found that nudging with the BV *x* and *y* components only is better than simple nudging with *x* and *y*, and for large *k* better than the *y* component alone. Nudging with SV *x* and *y* is better than with BVs for small values of *k* but worse for higher values. Table 1 summarizes the rms error obtained with the different methods and a different number of components, averaged over coupling strengths from *k* = 10 to *k* = 100. As indicated in Fig. 4, nudging with the leading vector from the three-dimensional orthogonalized BV subspace results in the lowest rms error, for arbitrarily chosen coupling strength. Our results indicate that supplying information in an unfavorable direction like the *z* component can cause desynchronization. Using dynamic vectors in the fast-growing directions indeed has an advantage over simple nudging in achieving synchronization because simple nudging does not include adaptive information on the dynamical evolution of the system. A striking result of these experiments is that coupling with a single fixed variable (*y*) is about as effective for synchronization as any of the time-varying schemes that were tested, except at very low coupling strengths. Despite the obvious idiosyncrasies of the Lorenz model, this result implies that coupling a few key variables may actually be more effective in achieving synchronization than attempting to couple all variables. This could also be applicable to data assimilation schemes.

As noted before we only used information from the slave system to construct the BV and SV, the only approach feasible in practice since the master system is never fully known. If we construct the BV and SV with the master system (not shown), both the BV and SV give lower rms errors for either two or three directions when *k* is greater than 10.

## 3. Generalized synchronization

The above experiments involve identical synchronization (IS) in which two identical nonlinear chaotic subsystems with different initial conditions will synchronize over time if they are coupled (Rulkov et al. 1995). The slave system is identical to the master aside from a coupling term and different initial conditions [Eqs. (1) and (2)].

In the case of generalized synchronization (GS), the two nonlinear chaotic systems are not identical (Boccaletti et al. 2002): GS is claimed to occur when there is a map *ϕ* from the driving (master) space D to the response (slave) space R, which takes trajectories **u**(*t*) in one space to trajectories **v**(*t*) in the other (Pyragas 1996). In terms of the geometry of the two subsystems, the evolution collapses onto the hyperplane in the full space. While GS has attracted a lot of attention since its discovery by Rulkov et al. (1995), in practice it is quite difficult to detect. Unlike IS, which gives an identical copy of the master, with GS there is often not an obvious visual relationship between the variables of the subsystems that can be gleaned from a simple plot like Fig. 1. The problem is especially grave if the correspondence *ϕ* between master and slave is given by a function that is nowhere differentiable or is multivalued (So et al. 2002), and special tests are required. Some of these tests call for intensive computational analysis, such as estimating conditional Lyapunov exponents from observed time series of both master and slave (Pyragas 1997). Another approach has been to use random sampling such as the mutual false nearest neighbors test (MFNN), based on the concept that points that are close in the phase space of the slave system should correspond to points that are close in the phase space of the master system (Rulkov et al. 1995).

A particularly simple test, the auxiliary system test, is used here to determine whether GS has been achieved (Pyragas 1996). The test does not require calculation of statistics; it only requires a second copy of the slave system and gives a clear yes/no answer. Given a master system **u** and a slave system **v**, a second slave system **v**′ is created, identical to **v** and coupled in the same manner to **u**, but started from different initial conditions. Then **u** and **v** are considered to be generally synchronized if the distance *d* = ‖**v**′ − **v**‖ → 0 as *t* → ∞. The identical synchronization of the two slave systems indicates that there is a function *ϕ* such that **u** = *ϕ*(**v**) (Pyragas 1996). (Visual inspection is commonly taken to be adequate to verify identical synchronization of the slaves since degradation occurs via on–off intermittency, giving bursts that are large and noticeable.) Unlike other tests, this test should reveal generalized synchronization even for intractable correspondence functions.

In the GS experiments, the goal is to be able to discern whether specific information about the master can be obtained from the slave. For GS, the parameters in the slave equation differ from those that determine the master’s dynamics. Specifically, the *r* parameter varied from the master *r* = 28. Parameter values were chosen arbitrarily, but so as to avoid bifurcations into different dynamical regimes. For *r* = 27, 26, and 25 in the slave Lorenz model, GS was achieved quickly, within 1000 time steps, even with a very low coupling strength *k =* 0.5 for which identical synchronization could not be reached, whereas for *r* greater or equal to 28, it took about 5000 steps to reach GS. The results were qualitatively similar for all *r* values tested. The GS experiments were carried out without observational errors.

The results of the experiment with *r* = 28 for the master and *r* = 25 for the slave are shown in Fig. 5. As indicated above, for the Lorenz model generalized synchronization is easier to attain than identical synchronization: when the master and slave systems are different, the two identical slave systems synchronize quickly even at very low coupling strengths. Figure 5 is computed with a coupling strength *k* = 0.5, smaller than the minimum coupling strength required to synchronize a slave and master with equal parameters. One possible explanation is that for the Lorenz model a synchronization manifold not given by the identity function is more stable than the manifold defined by the identity function. Figure 6a shows the behavior of the master and one of the slaves when it is uncoupled, showing similar behavior as the master but slightly slower orbits. Figure 6b is an expanded view of Fig. 5 for the last 5000 steps. When the slave is coupled (as in Fig. 6b) it exhibits smaller amplitudes, a frequency locked to that of the master, and much fewer regime changes, suggesting a much more stable attractor. It should be pointed out that the auxiliary system test indicated that GS took place even when the master was replaced by its time average or by randomly varied values, again suggesting that the GS pushes the slaves toward a more stable attractor. Even though having two identical slave systems does not allow direct communication between the slaves, at each time step both are being forced away from their normal attractor in the same direction, with a forcing proportional to the distance of each slave from the master. If the modified slave trajectory is more stable, bursting behavior such as described in section 2 is absent and the common forcing thus induces GS.

That the two Lorenz slave models synchronize so readily may cast doubt on the usefulness of the auxiliary system method and the concept of generalized synchronization. Nevertheless, the auxiliary system method does establish that GS exists in a formal sense and even a highly intractable correspondence function may be of some use. Indeed, a more careful examination of Fig. 6b shows that the slave is able to provide some information to help identify the regime in which the master is currently orbiting. The slave responds not by changing regimes itself, but by modifying the size of its orbit. For an uncoupled system, as long as the system stays on one side of the attractor, the orbits always “spin up,” getting farther away from the center of the attractor until the regime changes and the system starts at a lower orbit in the new regime (see Fig. 6a). By contrast, Fig. 6b shows that, under generalized synchronization, the slave spins up only if the master is in the same regime as the slave, but spins down or remains with a small amplitude when it is in the opposite regime. This provides information about the master state, but there are times in which the opposite is true (see arrow on Fig. 6b). As a result, even though the slaves are in perfect identical synchronization, the information they provide about the master state via the above rule is only approximately correct. This phenomenon of a weak coupling transferring only partial information has also been observed in synchronization of chaotic electrical circuits, where the phase of the oscillation may be transferred to a coupled slave circuit, but not the amplitude (phase synchronization, Boccaletti et al. (2002). A master and two slaves with identical parameters and a low coupling strength (*k* = 0.5) did not achieve identical master–slave synchronization but the two forced slaves did synchronize with each other. This occurs after a much longer integration time than for the slaves and master with different parameters (5000 compared to 2000 time steps), presumably because the two slaves were not pushed away from their unstable attractor as strongly as in the case of different parameters. As with the GS cases described above, changes in the slaves did reveal some information on regime changes in the master.

It should be noted that the failure of GS to provide a complete description of the master state may be specific to a model with multiple regimes, like the Lorenz model, where the transition is extremely unstable. For the Rössler (1976) system, given by

with a weaker form of nonlinearity that appears only in the *z* equation, and which exhibits no readily distinguishable regimes in state space, generalized synchronization behaves differently. Generalized synchronization between coupled Rössler systems, formed from Eq. (6) in the same manner as Eq. (2) is formed from Eq. (1) with different values of the parameters *a, b*, and *c* in the master and slave, is tested by checking for IS between two slaves. Where GS is verified, for typical parameter differences, the behavior of the slave system is found to provide much more complete information about the master than GS in the Lorenz system (see Fig. 6c). The improvement is most likely due to the fact the Rössler system has only one regime. Generalized synchronization was also tested between two different chaotic systems, a Lorenz master with two Rössler slaves and a Rössler master with two Lorenz slaves. Synchronization of the two slaves was achieved in each experiment using the same low *k* = 0.5 as for the Lorenz-only system experiments above. The results were similar to those of the Lorenz system: the synchronized slaves provided only very partial information about the master.

We conclude this section by pointing out that atmospheric models are imperfect representations of the real atmosphere, so the desired relationship between reality and model (good analyses and forecasts) is that of generalized synchronization. In practice, a method akin to the auxiliary system test is already used: The agreement between two operational data assimilation systems (slaves) using the same atmospheric observations from the master has been considered an indication of their agreement with the real atmosphere. The results above indicate that this conclusion could be overly optimistic: the agreement between the models indicates some definite correspondence between model states and the state of the atmosphere, but the relationship could be rather complex.

## 4. Comparison of conventional data assimilation and synchronization

Data assimilation aims to estimate the state of a system by statistically combining a short-range forecast (a.k.a. first guess or background) with observations in the most accurate way. Observations and forecasts are assigned relative weights based on their respective error covariances so that the linear combination of all available data has the smallest possible total error covariance. As with synchronization, the ultimate goal of data assimilation is to synchronize the forecast with the actual atmospheric state, using a combination of previous forecasts and current observations.

In our experiments, we used both the 3DVAR method (widely used in operational weather forecasting centers) as well as EKF, which is too expensive to implement operationally in real weather prediction systems without approximations. Both systems produce analyses (estimates of the truth) by a weighted average of a forecast and observations. 3DVAR is optimal, in the sense of giving a result that is statistically as close to truth as possible when using time-independent weights. The EKF method is optimal in the same sense, but for weights that are allowed to vary in time, and under the assumption that a linear approximation to the dynamical equations can be used to update the weights for the next assimilation cycle. The EKF method is much more computationally expensive than 3DVAR because it reestimates the forecast error covariance matrix 𝗕, used to compute the weights, each time a new observation is taken and assimilated. Although the 3DVAR method is somewhat less accurate, it is commonly used because using a single estimate of the 𝗕 matrix for all times significantly reduces computational costs (see the appendix).

Note that the various methods of data assimilation are defined in terms of a somewhat different goal than synchronization. EKF and 3DVAR are designed so as to make the best possible estimate of current truth, while synchronization is optimized so as to promote the fastest possible convergence of model to truth in the future. It has been shown that the two goals usually, but not always, yield the same result for the optimal coupling matrix 𝗖 in Eq. (3) (Duane et al. 2006, manuscript submitted to *Nonlinear Processes Geophys.*).

The extended Kalman filter method is based on assuming local linearity of the dynamics and must be modified to produce satisfactory results in the fully nonlinear case (Miller et al. 1994; Anderson 2001). We experimented with three methods to augment the background error covariance in the extended Kalman filter method, which otherwise would underestimate the forecast errors (see the appendix). In the first, small random perturbations uniformly distributed between 0 and *μ* were added to the diagonal of the analysis error covariance during the data assimilation cycle. In the second, the background covariance was multiplied by an inflation factor (1 + Δ), and in the third variation of EKF we used both random perturbations and covariance inflation, with *μ*, Δ small tunable parameters. As indicated in the appendix, although in more complex systems the background is usually kept constant during the data assimilation integration within the linear and adjoint models, we found that evolving the background state at every Runge–Kutta sub–time step was essential in order to get accurate results.

Table 2 shows the average rms error from the 500th to 5000th cycle. Each cycle consists of eight time steps and one new set of observations that, as in section 2, have randomly distributed errors with mean equal to zero and standard deviation 2. For 3DVAR, the 𝗕 matrix was first optimized to minimize the analysis error (see the appendix). Similarly, the value of *k* was optimized for synchronization for each set of observations. Finally, for the case of *x, y,* and *z* observations, we also performed adaptive synchronization where the coupling coefficient was strengthened if there was bred vector growth (see discussion of Fig. 3). So as to properly compare synchronization to the traditional data assimilation methods, the coupling term in the synchronization method [Eq. (2)] was only turned on at the same rate as analyses were performed in the traditional methods (i.e., every 8 time steps for Table 2 and 25 time steps for Table 3).

Table 2 indicates that the rms errors in the EKF using a multiplicative covariance inflation factor of 1.05 are all larger than the observation error, showing that the EKF cannot be stabilized at this (or even larger) level of inflation. By contrast, the approach to variance inflation based on adding random noise effectively stabilizes EKF, and the combination of the two approaches is, in the case of *x* observations, even more effective in avoiding EKF divergence. The best EKF always results in lower rms error compared to 3DVAR and synchronization when using the same number of observations (except for *z*, for which all systems fail). The results obtained from EKF and 3DVAR agree with those of synchronization, in the sense that rms errors from the *y* observations are the lowest when only one observation is available and are also lower than the rms error when observing both *x* and *z*. These results clearly indicate the importance of using effective observations in the optimal estimation of the true (nature) state. Indeed, the simple nudging method for synchronization is competitive with 3DVAR for the cases that include the *y* variable among the observations, a remarkable result given the simplicity of the nudging algorithm.

The leading (conditional) Lyapunov exponent of a data assimilation system provides information about the rate of convergence to the truth and can be easily estimated by using breeding (Trevisan and Uboldi 2004). For the systems presented in Table 2 we found that the more negative the exponent was, the smaller the long-term average analysis error. One problem with the 3DVAR method of data assimilation is the occasional incidence of analysis error “catastrophes,” during which rms analysis error suddenly spikes to unacceptable levels. These spikes tend to occur in conjunction with changes of regime on the Lorenz attractor. Building on the work of Evans et al. (2004), we were able to anticipate time periods more likely to contain one of these catastrophes by estimating the local (in time) Lyapunov exponents using the growth rates of the bred vectors of the analysis (not shown). Similar results were obtained with synchronization (see Fig. 3a, red curve), showing that large error bursts in synchronization with observational errors also occur when there is a strong likelihood of regime change (see Fig. 3b). We plotted stars on the evolution of the *x* variable for the master system (see Fig. 3b) with colors according to the bred vector growth rate from the slave system (conditional Lyapunov exponent), plotted on Fig. 3c. The presence of red stars, indicating an exponential growth rate greater than 0.04 (see Fig. 3b), is a good predictor of a forthcoming regime change and large synchronization errors (see Fig. 3a, red curve). Indeed, the spiking behavior is an instance of the on–off synchronization phenomenon described in section 2 as a typical mode of degradation of synchronized chaos (Baek et al. 2004). Such behavior is prone to occur near the unstable fixed point at the center of the Lorenz attractor, which is also associated with regime transitions.

These results suggest that large analysis errors may be predictable, since catastrophes tend to occur along with unstable situations, allowing us to selectively modify or use an alternate method of data assimilation at those times more likely to contain large spikes in the rms error. For example, one could increase the forecast error covariance matrix 𝗕 in 3DVAR (or the coupling *k* in synchronization) when such a situation is observed, and return it to its normal values after the large growth ends. EKF evidently does that naturally: the trace of the background error covariance matrix 𝗕, which determines the relative weight given to the forecast, also increases substantially near regime changes (not shown). We tested a simple adaptive synchronization approach: we multiplied the coupling coefficient by twice the growth of BV over eight steps, only if this growth was larger than one (see Fig. 3a, blue curve). Even without any tuning, this simple adaptive coupling was quite successful in avoiding most (but not all) the episodes of large error, not unlike the extended Kalman filter. The rms error of this adaptive synchronization approach is shown in the last column of Table 2.

In the experiments discussed so far, we used more frequent observations than Miller et al. (1994). We now compare in Table 3 the results obtained with EKF, 3DVAR, and simple nudging when observations are available every 25 time steps, as in Miller et al. (1994). Not surprisingly, results are worse than in Table 2, with observations every eight steps. The EKF diverges due to the highly nonlinear behavior of the Lorenz system, as shown by Miller et al. (1994). We found that it was possible to overcome this problem by making the random perturbations added for variance inflation adaptive, depending on the trace of the EKF forecast error covariance (*μ*′ = *μ* × trace(*P _{f}*)/3), where

*μ*is a constant. We also noticed that the analysis error “catastrophes” observed for 3DVAR are also present in EKF when the forecast error covariance 𝗕

^{b}

_{EKF}is occasionally not well conditioned. Since the times at which ill conditioning occurs (like the 3DVAR analysis error catastrophes) can be diagnosed by the BV growth rate, we also used a hybrid scheme to improve the conditioning of the forecast error covariance at these times. We estimated that, when the BV growth rate is larger then 0.06, the hybrid scheme is needed and is replaced the EKF forecast error covariance with a linear combination with the 3DVAR (constant) background error covariance:

^{b}

_{EKF}= 0.6𝗕

^{b}

_{EKF}+ 0.4𝗕

^{b}

_{3DVAR}. Combining the adaptive random perturbation and conditional hybrid scheme was the most effective way that we found to reduce the rms error. Although these results are not necessarily optimal, they point out possible adaptive strategies to increase the accuracy of methods such as the ensemble Kalman filter, which have difficulties outperforming 3DVAR with real observations (Houtekamer et al. 2005).

We also note that with infrequent observations synchronization (nudging) had much worse performance than EKF and 3DVAR (see Table 3). To address this problem we tested time interpolating the observations, with somewhat positive results (not shown). In data assimilation the analysis update is performed using “observational increments” or innovations, which have the advantage of being relatively small, so that linear operations such as interpolation can be done more accurately on the increments rather than on the full values. Based on this experience, we tested nudging interpolating, not the observations but the observational increments defined as observations minus forecast. For each 25-step assimilation window, the initial and final innovations were linearly interpolated in time. These innovations were then added back to the forecasts to create synthetic observations used for nudging every time step within the assimilation window. Table 3 shows that with this approach, synchronization becomes very competitive with the standard (optimized) 3DVAR.

## 5. Summary and discussion

Our experiments with identical synchronization of the Lorenz (1963) model with random observational errors and lower coupling strengths confirm the naïve expectation that simple nudging with all three variables leads to faster synchronization. However, the optimal choice of coupled variables depends on overall coupling strength. Nudging with *y* alone was more effective at higher coupling strengths. Including variables that are less effective for synchronization appears only to add noise when the coupling is strong. One of the potential advantages of synchronization over conventional data assimilation is that it may be effective when only a small number of variables are coupled with more frequent observations. It could provide a framework for “targeted observations” (e.g., Kalnay 2003, p. 243) in a situation in which the total number of variables or degrees of freedom that can be observed is restricted.

Choosing an optimal coupling strength for a particular combination of coupled variables was facilitated by the addition of random observational errors, which helped to identify a required level of coupling beyond which further coupling needlessly introduced more error. A wide range of coupling strengths (1 to 100) was found to be effective in achieving synchronization without compromising the model dynamics, even when the coupling time scale was smaller that the dynamical time scale of the uncoupled systems, confirming that such schemes are a practical alternative to traditional methods (Hoke and Anthes 1976; Stauffer et al. 1994).

Suggesting an extension of the nudging method as an alternative to more sophisticated data assimilation schemes may seem odd. But extended nudging [Eq. (3)] is simply the general form that includes 3DVAR and extended Kalman filtering. Since it is already known that extended Kalman filtering must be modified to give satisfactory performance in all situations (see the appendix; Miller et al. 1994), the kind of generalization proposed here is not new. Rather, the synchronization view provides a framework in which to define optimal coefficients and subspaces for extended nudging. These should be chosen so as to make the conditional Lyapunov exponents as negative as possible (e.g., Trevisan and Uboldi 2004). Good results can usually be achieved with suboptimal, but formally simpler, coupling schemes.

We tested whether synchronization along varying directions given by either singular or bred vectors obtained from the slave system was superior to simple nudging. Bred vectors did not lower the rms errors compared to simple nudging except when using both *x* and *y* observations, and singular vectors provided the best results when coupling all three variables. The best results for arbitrary coupling strength were obtained when nudging in the direction of the leading orthogonalized bred vector. Further experiments with a more realistic system would be needed to clarify whether the reduction in error could justify a similar use of singular vectors and/or bred vectors for synchronization-based data assimilation in NWP models. However, the fact that synchronization with simple nudging occurs over a large dynamic range of coupling strengths suggests that the particular choice of coupling basis may not make a large difference, so long as the basis vectors have significant projections in important directions.

We found that generalized synchronization can be attained when the slave system has parameters different from the master system even at very low coupling strengths for which identical synchronization does not occur. While some information such as orbit frequency and amplitude is lost, the slave’s constrained dynamics provides some critical information on regime changes in the master’s Lorenz attractor. For practical purposes, tests such as the auxiliary system method are overly inclusive, since the identical synchronization of the Lorenz slaves occurs even with constant or random forcing, and may reflect the fact that forcing the slaves away from their normal chaotic attractors can push them toward a more stable attractor. The task of unraveling the exact relationship between the master and its generally synchronized slave is daunting, especially when this relationship is nondifferentiable or multivalued. But, even without knowing the exact relationship, the existence of generalized synchronization might be used to provide a probabilistic description of the master state for a given state of the slave system.

Our results confirmed the effectiveness of data assimilation schemes, especially the extended Kalman filter. Although our synchronization tests did not perform quite as well as 3DVAR data assimilation, the synchronization algorithm is much simpler, so the fact that our results are comparable (especially when time interpolating observational increments) supports the use of extended nudging in conjunction with traditional data assimilation schemes. The ability to use bred vector growth rates to anticipate unstable points like regime changes also offers a way to target these situations for a short-term synchronization scheme substitution in the operational 3DVAR data assimilation system. Furthermore, there is no need to restrict nudging to a fixed direction in state space. It is possible that a judiciously chosen dynamic coupling scheme in which the direction and magnitude is allowed to vary may improve on existing data assimilation techniques. Indeed, conventional data assimilation schemes are but one form of synchronization in its most general sense.

## Acknowledgments

This work was done under NSF Grants ATM0328402 and ATM0327929. Several of the coauthors were supported as RISE scholars by the NSF-UMD RISE program. We thank three anonymous reviewers for their very constructive suggestions. Debra Baker was supported by the NSF Graduate Fellowship Program and this research was carried out in partial fulfillment of her University of Maryland doctoral requirements.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Data Assimilation Schemes

The equations used for data assimilation are (e.g., Kalnay 2003, pp. 155 and 179)

The analysis is given by

where the optimal weight matrix is given by

and the new analysis error covariance by

Here *n* denotes the analysis step (in our case the analysis is done every 8 or 25 model time steps), *M _{n}*

_{−1}is the nonlinear model that provides the forecast or background

**x**

*at step*

^{b}_{n}*n*starting from the previous analysis

**x**

^{a}

_{n−1}, 𝗠

_{n}, 𝗠

^{T}

_{n}is the corresponding linear tangent and adjoint models, 𝗥 is the observational error covariance (assumed here to be constant and equal to twice the identity matrix), 𝗕

*is the background error covariance, 𝗤*

_{n}*is the covariance of the errors introduced by the model, and 𝗛 is the observation operator that transforms the forecast into a first guess of the observations. For example, if we observe only the first two variables*

_{n}*x*and

*y*of the Lorenz model, then

In 3DVAR, the background error covariance is estimated once and for all; 𝗕 = is a time average from data assimilation cycles 500 to 5000, where *x ^{f}* is the forecast state and

*x*the true state. To have the optimal 𝗕, we iterated this process until converged. Except for the

^{t}*z*-only observation, this required less than 10 iterations.

In extended Kalman filtering the forecast error covariance is updated with the linear tangent model and its adjoint. Because of nonlinearities, even with a perfect model (𝗤* _{n}* = 0), the EKF analysis drifts away from the real solution due to an underestimation of the true forecast error covariance (e.g., Miller et al. 1994). We have successfully tested two approaches to avoid this filter divergence. The first one is multiplicative variance inflation (Anderson 2001), in which the background error covariance is multiplied by (1 + Δ). In the second method, we enhance the analysis error covariance by adding to the diagonal elements random perturbations uniformly distributed between 0 and 1, multiplied by

*μ*, before performing the time integration (Corazza et al. 2003). Both of these methods (and especially their combination) were able to avoid the extended Kalman filter divergence, and optimal results were obtained observing

*x*,

*y*, and

*z*with Δ = 0.05 and

*μ*= 0.02, values that were used unless otherwise noted. In the results of Table 2, the integration with the tangent linear model was first made keeping the background state constant for eight time steps, as frequently done for computational economy, but we found that the extended Kalman filter results can be substantially improved if instead we use the evolving background state at every Runge–Kutta time substep. When observing

*x*,

*y*, and

*z*, the analysis errors dropped from 0.49 to 0.33 (see Table 2). In the experiments where the background used in the tangent linear model integration evolved every time step, we found that the magnitude

*μ*of the random noise needed to maintain filter stability was substantially reduced, from 0.1 to 0.02.

## Footnotes

*Corresponding author address:* Gregory S. Duane, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307-3000. Email: gduane@ucar.edu

* Most of this paper is the result of a Research Internships in Science and Engineering (RISE) 10-week summer program. The RISE program encourages undergraduate women to continue their studies and pursue graduate degrees in engineering, technology, and scientific disciplines. It is funded by the National Science Foundation, the A. James Clark School of Engineering, and the University of Maryland. The program brings together women faculty (E. Kalnay), graduate mentors (S.-C. Yang and H. Li), and undergraduate fellows (D. Baker) to work with outstanding undergraduate interns (K. Cordes, M. Huff, G. Nagpal, and E. Okereke) on current research projects.