## Abstract

We develop and compare variants of coupled data assimilation (DA) systems based on ensemble optimal interpolation (EnOI) and ensemble transform Kalman filter (ETKF) methods. The assimilation system is first tested on a small paradigm model of the coupled tropical–extratropical climate system, then implemented for a coupled general circulation model (GCM). Strongly coupled DA was employed specifically to assess the impact of assimilating ocean observations [sea surface temperature (SST), sea surface height (SSH), and sea surface salinity (SSS), Argo, XBT, CTD, moorings] on the atmospheric state analysis update via the cross-domain error covariances from the coupled-model background ensemble. We examine the relationship between ensemble spread, analysis increments, and forecast skill in multiyear ENSO prediction experiments with a particular focus on the atmospheric response to tropical ocean perturbations. Initial forecast perturbations generated from bred vectors (BVs) project onto disturbances at and below the thermocline with similar structures to ETKF perturbations. BV error growth leads ENSO SST phasing by 6 months whereupon the dominant mechanism communicating tropical ocean variability to the extratropical atmosphere is via tropical convection modulating the Hadley circulation. We find that bred vectors specific to tropical Pacific thermocline variability were the most effective choices for ensemble initialization and ENSO forecasting.

## 1. Introduction

Near-term, or multiyear, climate prediction refers to time scales between several and a few tens of years (Meehl et al. 2014). Despite being a relatively new field, it promises to provide key information about the risks and uncertainties arising from the effect of anthropogenic forcing on the natural internal climate modes of variability. Such information is crucial for infrastructure planning, adaptation of agricultural and fisheries practices (Tommasi et al. 2017), and mitigating potential socioeconomic problems arising from imminent climate change.

Multiyear variability in the climate system is both forced by changes in the external drivers, such as incoming solar radiation and anthropogenic emissions of greenhouse gases, and internally generated due to the inherent variability of the climate system. A multiyear-to-decadal prediction system should comprise several key elements. First, a general circulation model (GCM) capable of simulating the variations and spatial patterns of the major observed climate teleconnection patterns is required. Second, reliable observations of the current climate and appropriate scenario-based boundary conditions (i.e., external radiative forcing’s relevant to the projected climate) are necessary. Third, a data assimilation and ensemble initialization/prediction system with which to provide accurate and balanced initial conditions should be included. Finally, a sufficiently long and spatially homogeneous observational record with which to constrain the model and to assess errors is needed. Ideally, all elements combine such that the forecast model tracks the future climate with some degree of fidelity.

Data assimilation (DA) is the process whereby an imperfectly observed state of the Earth system at a particular point or period in time is combined with short-term model forecasts of the same system to produce an analyzed representation of present conditions. A sequence of analyzed states represents a best-guess approximation to the evolution of the climate. The ensemble Kalman filter (EnKF) and its variants are widely utilized across a range of scales from weather to climate data assimilation applications (e.g., Houtekamer and Zhang 2016). Where observations are perturbed, insufficient sampling can lead to underestimation of the Kalman gain in the EnKF (Whitaker and Hamill 2002). To address this limitation, various deterministic square root ensemble Kalman filters have been developed (Tippett et al. 2003), including the ensemble transform Kalman filter (ETKF; Bishop et al. 2001). More computationally efficient variants that allow for massively parallel computations and favorable scaling with increasing numbers of observations, such as the local ensemble transform Kalman filter (LETKF; Ott et al. 2004; Miyoshi and Yamane 2007; Bowler et al. 2008; Sakov 2017), have also been developed.

For systems with a quadratic nonlinearity, for example turbulent geophysical flows, calculation of the error covariances requires the solution of an infinite hierarchy of moment equations (O’Kane and Frederiksen 2010). The Kalman filter is Gaussian and consequentially all odd-ordered moments are zero and by construction there are no fourth-order terms. To track regime transitions in nonlinear systems, some information about the third- and higher-order terms is required (Miller et al. 1994, 1999). This problem shares many similarities to the turbulence closure problem in statistical mechanics (O’Kane and Frederiksen 2008a, 2010). O’Kane and Frederiksen (2008a,b) showed that these higher-order processes, crucial in tracking regime transitions in 2D turbulence (Nadiga and O’Kane 2017), were also required to track the development and mature phase of coherent anticyclonic blocking in the troposphere.

The widely different spatiotemporal scales of the atmosphere and ocean represent a significant challenge for coupled DA, where the climate model provides the background covariances for and between each of the component systems. Where the ocean and atmosphere analysis increments are calculated separately, the approach is termed weakly coupled DA (Zhang et al. 2007). Whereas, when the error covariance is explicitly constructed to include cross correlations between atmospheric and oceanic state vectors, it is termed strongly coupled DA (Lu et al. 2015; Laloyaux et al. 2016). In weakly coupled DA the ocean and atmosphere communicate only via surface fluxes, whereas the state variables in strongly coupled DA are adjusted simultaneously during the analysis update. For example, using an intermediate complexity model and assimilating only atmospheric observations, Sluka et al. (2016), comparing weakly to strongly coupled DA on time scales relevant to weather prediction, found including additional information into the ocean increment based on atmosphere–ocean cross covariances reduced the ocean analysis errors, which over time improved the analyzed atmospheric state.

A common approach in attempting to establish the theoretical limits to predictability in geophysical flows is to consider the divergence of pairs of initially close states propagated forward in time using a model. Where ensembles of deterministic forecasts whose initial states are in some sense close to an initial analyzed state are employed, the problem of forecasting geophysical flows can be regarded as a statistical problem of predicting the probability density function of given states or, equivalently, of calculating the moments of the relevant prognostic variables. The inherently chaotic nature of the climate system (i.e., flow-dependent background error covariances), combined with a sparsity of observations, compounded by errors in the observed initial conditions and impacted by model biases (i.e., the ability of the model to faithfully represent climate teleconnections), leads to the failure of deterministic forecasts after a relatively short period of time.

For forecasting periods extending from a couple of weeks to seasonal time scales, it is common to initialize the coupled climate model to independently analyzed ocean and atmospheric states. Examples of this approach are the National Centers for Environmental Prediction (NCEP; Kleist and Ide 2015) and European Centre for Medium-Range Weather Forecasts (ECMWF; Bonavita et al. 2012) systems, which employ hybrid ensemble–variational DA systems for the atmosphere and the 3D variational data assimilation (3DVAR) method for the ocean (Mogensen et al. 2012; Saha et al. 2014). Another approach is to relax the coupled system toward the observed sea surface temperature (SST) field assuming that modification of the fluxes due to the ocean–atmosphere coupling, given enough time, will be sufficient to adjust the model climate state while retaining balance (Luo et al. 2008). Alternately, relaxation to independent atmospheric (winds, potential temperature, and surface pressure) and oceanic (temperature and salinity) reanalyses (Smith et al. 2007, 2013) is also employed. For near-term (1–10 yr) and decadal (10–30 yr) climate forecasting, it is often assumed that predictability resides in the ocean and that the initial state of the atmosphere can be approximated by any randomly chosen state corresponding only to the given month or season (Meehl et al. 2016).

Arguably, the most important information required in any ensemble of initial forecast perturbations corresponds to instabilities or error modes (Molteni et al. 1996; Toth and Kalnay 1997; Bowler et al. 2008; Wei et al. 2008). It has long been recognized in weather and climate forecasting that errors, such as those associated with instrumental observation errors, or inexact numerical methods including truncation error and unresolved scales, will degrade and limit the predictability of any chaotic deterministic system (Lorenz 1963). In addition, the instabilities inherent in the atmosphere–ocean system will grow, projecting onto an attractor occupying a low-dimensional subspace hence dominating error growth in dynamically active regions. Initial forecast perturbations that capture information relevant to the growth of coupled dynamical instabilities have been examined in simple models (Peña and Kalnay 2004; Norwood et al. 2013), intermediate complexity (Frederiksen et al. 2010; Sluka et al. 2016), and coupled GCMs (Yang et al. 2008; Sandery and O’Kane 2014). The addition of information regarding growing instabilities into the initial conditions effectively incorporates the most dynamically relevant information contained in the cross covariances.

As we are generally interested in applications from seasonal to near-term climate prediction and, in this study, specifically tropical processes, we will assume for convenience that predictability resides in the ocean and that atmospheric initial conditions are subdominant and henceforth we conduct data assimilation experiments on that basis. This point is first examined in an idealized simple nine-component Lorenz attractor model. Next, we examine various classes of strongly coupled DA–ensemble prediction systems (EPSs) in a GCM. The first variant consists of a 96-member ETKF assimilating ocean observations and applying atmospheric increments on the basis of ocean–atmosphere covariances. The second variant is based on ensemble optimal interpolation (EnOI) again assimilating ocean observations and applying atmospheric increments. Forecasts are then made using an ensemble of coupled ocean–atmosphere bred vectors (BVs) generated as perturbations about initial states derived by the EnOI. The coupled BVs are rescaled to the tropical ocean variability in two distinct ways, either based on a percentage of the background standard deviation (SD) or on the variance within a given spatiotemporal scale bounded by an appropriate isosurface. Here, we are interested in the analysis and background errors, biases, and the atmospheric response to upper tropical ocean perturbations. Both ETKF and EnOI variants are able to provide reanalyzed climate model initial conditions about which BV perturbations may be generated. We hypothesize that the ETKF ensemble mean should provide a less-biased analysis than EnOI but that using an ensemble of BVs specifically tuned to the dynamics relevant to the spatiotemporal scales of interest, in this case tropical coupled instabilities relevant to ENSO, would be effective as forecast initial perturbations, and particularly so where ensemble sizes are limited.

## 2. Data assimilation methodology

Data assimilation for large numerical prediction problems in the ocean and atmosphere requires calculating the full covariance matrices of the background errors. In the EnKF formalism (Evensen 2003; Houtekamer and Zhang 2016) this can be achieved either in terms of static covariances estimated as anomalies or deviations from the climatological mean calculated from a long (detrended) control simulation or, alternately, through calculating the time-evolving background covariances via propagating an ensemble of model states and constructing the forecast (background) error covariance in terms of deviations from the first moment or ensemble mean at a particular time.

### a. Ensemble Kalman filter

Following O’Kane and Frederiksen (2010), let us apply the ensemble Kalman filter methodology to a system with model dynamics given by and considering an *n*-dimensional state vector **x** at time step *t*, a *q*-dimensional state space noise process vector **w** due to disturbances and model errors, and an observation vector **d** with measurement noise (here, *n* and natural numbers). Next, let us assume the observation error and the model error **w** are white such that the analysis and forecast fields are defined as

where runs over the entire ensemble; and ^ indicate the ensemble mean and perturbation, respectively; and where , , and for , where *δ* is the Kronecker delta function. The ensemble Kalman filter propagates first and second moments of **x** recursively, where and can in general be nonlinear (∆*t* is the time step). Assuming that the perturbation field ^ runs over the entire ensemble, we can write the recursive Kalman filter equations as

where the analysis fields are calculated for each realization *i* and the covariances , where

where is termed the Kalman gain, is the positive definite state covariance error matrix, is the identity matrix, and is the linearized observational operator, mapping forecast gridpoint values onto observational points.

The EnKF methodology is based on perturbing the observations, assuming that they are not correlated with the state (Burgers et al. 1998). For unperturbed or single-realization observations the covariance update [Eq. (3)] is systematically underestimated,

### b. Ensemble transform Kalman filter

The ETKF (Bishop et al. 2001; Wei et al. 2006, 2008; Bowler et al. 2008; O’Kane et al. 2008) is based on an alternate application of the Kalman filter methodology whereby *k*-ensemble forecast and analysis anomalies **z**_{i} for , defined as

and where the state vectors are and , which are *n*-dimensional in model space, and where indicate the mean of *k*-ensemble members. The ETKF acts to choose appropriate initial forecast anomalies consistent with error covariance update equations within the vector subspace of ensemble anomalies [Eq. (2b)] formed as and . To calculate the ETKF transform matrix , we are required to calculate the matrix of ensemble perturbations in normalized observation space; that is,

where maps from into observable space and does the renormalization.

We are next required to find the eigenvectors **C** and eigenvalues of Eq. (5), which is equivalent to the matrix . The transform matrix is now defined in terms of the matrix of nonzero eigenvalues such that

where the nonzero eigenvalues are and **C** is , which corresponds to the transform matrix in spherical simplex form as described by Wang et al. (2004). The ETKF generates new analysis perturbations via the transform matrix:

where the transform matrix determines how the perturbations from different members are mixed. We directly initialize the full state vector to the analysis, whereby model biases are reduced by constraining the model to be close to the observations.

### c. Ensemble optimal interpolation

A computationally inexpensive approach to data assimilation is the EnOI method of Evensen (2003). EnOI employs a time-invariant or static background covariance matrix , typically formed from a predetermined ensemble of anomalies . For the situation where is a linear operator, EnOI is equivalent to 3DVAR, which similarly uses the time-independent background covariance from a long-term statistical mean. In practical applications, the are formed as anomalies w.r.t. the climatological seasonal cycle calculated from a long control simulation of the forecast model. An alternative to using the anomaly w.r.t. climatalogy is to use the tendency of the free-running model calculated over the same interval as the assimilation window. Our tests found tendencies of the type described were more effective than anomalies and were consequently employed in this study. The advantage of EnOI is the relatively low computational cost and the fact that it allows incorporation of anisotropic multivariate covariances. EnOI is an important and preferred assimilation method for several operational oceanography centers (Cummings et al. 2009; Martin et al. 2015; Sandery 2018). In the calculations that follow, the EnOI assimilation includes atmospheric increments due to cross covariances between the ocean and atmosphere, scaled to be no larger than the tendency of the free-running atmospheric model over a time interval determined by the assimilation window (i.e., 28 days).

### d. Covariance localization

Whether using variants of the EnKF or EnOI, the large computational costs involved in running a single climate forecast model, let alone many, imposes a finite upper bound on the dimension of our ensemble (currently 96 for the ETKF). In the ensemble formulation, unrealistically small ensemble sizes, , often result in large sampling errors manifesting as spurious background error covariances and, as a consequence, overestimation of the impact of each observation. To ameliorate this sampling error, we employ horizontal localization in the form described by Evensen (2003), Ott et al. (2004), and Hunt et al. (2007), in which local ensemble transforms are calculated for each element *i* of the state vector and the analysis, calculated independently at each grid point using only the background ensemble and observations from the previously defined surrounding subdomain or alternately in a reduced-rank subspace. For each model grid point only observations falling within a given radius are used to produce a local transform matrix. This is done using local normalized ensemble observation anomalies and local normalized innovations, obtained via a Schur product with a weighting function similar to a Gaussian function, as in the approach of Gaspari and Cohn (1999). This effectively removes any spurious long-range correlations present in the ensemble estimated background error covariance, and has also been shown to improve the relationship between ensemble spread and error as a function of latitude. Disadvantages are that the localization process can remove information due to large-scale flow-dependent inhomogeneities and create potential imbalances (Zagar et al. 2010).

In applying the Kalman filter, it is common to assimilate observations collected within a fixed window of time at each given cycle. For EnOI, this requires the observation window to be centered about the assimilation time, that is, days. In such cases, we are assuming that observations are made simultaneously at the time of assimilation. This assumption results in an averaging effect, whereby variability on time scales less than the observation window is removed. For the climate model ETKF simulations that follow, we apply the asynchronous approach of Sakov et al. (2010), accounting for the timing of when the observations were taken with the resulting analysis increment corresponding to a sum over increments due to the assimilation of observations at particular times within the window.

### e. SST bias correction

The most common manifestation of model error is as a persistent error in the background state over a sequence of assimilation cycles—a good example being the warm SST bias in the Southern Ocean during the austral summer common in many climate model simulations (Wang et al. 2014). To remove SST bias from the respective ensemble forecast members in the ETKF, we represent SST as a sum of the model SST and an unknown SST bias field. To prevent collapse of the SST bias ensemble, we follow Evensen (2003) [see also section 2.12 of Sakov (2017)] by using the first-order autoregressive [AR(1)] model of the form

applied where is the cycle length and *σ* is normalized white noise, that is, , where is the error standard deviation of **q**. Here, *λ*, that is, , , determines the persistence of the bias. Equation (8) is applied during the analysis where covariances between model SST and the model noise vector are computed, which is then updated together with the state vector.

### f. Ensemble prediction

The stability of the flow and the associated growth of initially small perturbations determines predictability. Linear error growth, where initial perturbations are infinitesimal, can be described by making a tangent linear approximation to evolving system dynamics (Nadiga and O’Kane 2017). In geophysical systems linear error growth associated with the leading Lyapunov exponent is often observed over a relatively short period in the evolution of perturbations under the full nonlinear equations. However, because of the possibility that the leading Lyapunov exponent is related to small-scale instabilities in a multiscale system where predictability originates from the larger scales (Boffetta et al. 1998), the leading Lyapunov exponent is often of limited relevance to predictability. Finite-size initial perturbations (e.g., BVs) renormalized by a suitably chosen norm over a finite time shorter than the error saturation time scale have been shown to project onto spatial structures that rapidly amplify in local regions of large forecast error (Nadiga and O’Kane 2017).

Because of their ease of construction, retention of the past, or memory and dynamical balance, we have chosen to characterize flow-dependent error growth in the model using BVs (Toth and Kalnay 1997; O’Kane and Frederiksen 2008b; Evans et al. 2004; O’Kane et al. 2011; Pazó et al. 2013; Sandery and O’Kane 2014). BVs are finite perturbations generated (or bred) by evolving the perturbed system , where is the control state under the full nonlinear governing equations, and ∆ω(*t*) is the perturbation of the model state from the control (Nadiga and O’Kane 2017). The perturbations themselves are rescaled to a given size *ε* periodically at a time interval ∆*T* as follows. The difference between the control and perturbed trajectories, , is computed at times for , whereupon the perturbation is rescaled and the perturbed system is redefined as . The perturbation is now allowed to evolve freely until the next rescaling is scheduled at time . The BV corresponds to the (finite) perturbation constructed at time *t*. Within the context of our earlier notation, the breeding method produces each new analysis perturbation via a straightforward rescaling of the forecast perturbation by a uniform factor . BVs are formulated to span some local low-dimensional subspace of a given growing dynamical mode. The local growth rate of the BVs is given by

where ∆*T* is the rescaling interval and is the BV at time *t*. More generally, we may define the relative amplification factor in terms of the vector of gridpoint values of the BV of any climate variable field as initiated at time *t* and evolved to time . We take the *L*_{2} norm as the RMS of the vector by and define the amplification factor as and the local total growth rate . A schematic of an ensemble prediction system based on BVs is depicted in Fig. 1.

The ETKF has a prominent place in operational ensemble numerical weather prediction (Bowler et al. 2008; O’Kane et al. 2008; Wei et al. 2006). In applications principally concerned with ensemble prediction and not data assimilation, the mean or first moment is unchanged between the analysis and forecast. This approach has close similarities to breeding. Comparisons of the utility of bred versus ETKF initial forecast perturbations have been performed for atmospheric forecast applications (Wang and Bishop 2003; Wei et al. 2008; Deng et al. 2012) but few if any exist for climate forecast applications. In application to weather prediction Wang and Bishop (2003) used no (or very little) localization for their comparison of ETKF versus BV performance. They found that, in the absence of localization, the ETKF produced nearly orthogonal BV-type perturbations at essentially the same cost as a similar number of BV perturbations for use in an OI-type system. A major difference in the current study is the fact that we localize the ETKF perturbations but not the BVs.

## 3. Paradigm model results and discussion

Prior to comparing the relative utility of EnOI versus ETKF assimilation systems in a high-dimensional climate model, we first apply these methods to a low-dimensional nine-component Lorenz model of the extratropical atmosphere coupled with ENSO. This model was previously described by Peña and Kalnay (2004) and in modified form by Norwood et al. (2013), who used it to investigate the instability properties of BVs, as well as other related dynamical vectors. To our knowledge, this is the first time it has been applied within a data assimilation context. In this small model, three versions of the famous Lorenz (1963) model are coupled to mimic the temporal behavior of an extratropical atmosphere [Eqs. (10)] weakly coupled to a tropical atmosphere [Eqs. (11); coupling coefficient ], which in turn is strongly coupled to a slow tropical ocean [Eqs. (12); coupling coefficient ]. This small model we use as a paradigmatic representation of the coupled climate system:

The additional parameter settings are the temporal and spatial scaling factors and , respectively, for the ocean variables, such that the ocean has a time scale an order of magnitude longer than the time scale for the tropics and extratropics but the same amplitude. Both and are uncentering parameters, while , , and are standard settings from Lorenz’s original work (Lorenz 1963) on convection. The model was integrated forward over times with a time step of 0.01 using a fourth-order Runge–Kutta scheme. These choices lead to a tropical subsystem that is completely dominated by changes in the ocean subsystem, but has an extratropical atmosphere, representative of weather noise, whose behavior is similar to the original Lorenz model. From Peña and Kalnay (2004), we note that the ocean exhibits significant deviations from its normal trajectory about every 2–8 years, which one might consider to be ENSO like. Here, 3 years corresponds to approximately 10 time units.

In the EnOI and ETKF data assimilation configurations, observations of only the slow ocean variables [Eq. (12)] are assimilated. The ocean observations were constructed by adding noise with a standard deviation of 1 to the truth simulation. We assimilate observations from the truth simulation assuming the model is perfect. The tropical and extratropical variables are unobserved but are adjusted via a cross covariance with the ocean observations. The ETKF background covariance uses a 10-member ensemble whereas the EnOI static background covariances are constructed as 50 random samples from a long control simulation. The data assimilation is full rank; that is, the ensemble size is larger than the number of variables. The ETKF uses covariance inflation with a factor of 1.01. The particulars of the data assimilation parameter settings are described in Table 1.

In Fig. 2a we show the evolution of the three states started from a set of randomly chosen initial conditions (red curve designated truth). The points in blue indicate positions along the trajectory described by the evolution of the EnOI-analyzed state, where the tropical and extratropical attractors are affected by the cross covariances due to the strongly coupled assimilation of ocean observations. Cross covariances between attractor state variables are calculated, contributing to the analysis increments of each respective prognostic variable at each assimilation cycle. The EnOI-analyzed ocean state tracks the truth trajectory. However, there are obvious systematic errors, indicating that the analysis increments are large. The tropical atmosphere, even though strongly coupled to the ocean, does not track the truth and in fact occupies a different manifold and larger volume in phase space, indicative of less certainty in the position of the analyzed state. The extratropical weather attractor is weakly coupled to the tropics and, while not tracking the truth at all, is obviously occupying a similar region of phase space.

For the ETKF case (Fig. 3), with the same parameter settings, we see that the ocean analysis closely tracks the truth, as does the tropical atmosphere after an initial period (≈500 time steps) in which it becomes progressively more constrained by the ocean observations. On the other hand, the variability of the analyzed extratropical weather attractor is quite severely restricted by the cross covariance with the ocean. This is surprising given the weak coupling to the tropical atmosphere and arises due to the cross covariances between the respective state vectors of the attractors. The reduced variability, occupying a smaller domain on the extratropical weather attractor manifold relative to the truth, arises as a consequence of a damping effect due to flow-dependent background covariances tightly constraining the ocean–tropical atmosphere while reducing variance in the extratropics. In the EnOI, the ensemble spread associated with the background covariances is time invariant and, for this small model, preserves the statistical properties of the true slow–fast system but fails to track the trajectories of either the tropical or extratropical atmospheres.

From these simple assimilation experiments, we surmise the following. First, in this system, a well-observed “ocean” (slow system) is in principle all that is necessary to constrain “tropical ocean–atmosphere” (strongly coupled slow–fast) coupled modes when flow dependency of the background error covariance is taken into account. Second, in a multiscale system, there may be unintended consequences when constraining the slow part of the system and taking into account its projection onto the fast temporal scales—even when they are weakly coupled. In near-term climate forecasting, one would ideally like to constrain the slow ocean and the intermediate-time-scale coupled climate modes, while at the same time recognizing that we cannot hope to track actual weather systems, but that we need to preserve the character of the variability that arises from them. For the purposes of the discussion that follows, we leave here the small model, noting that fundamental questions about the relationship between attractor dimension and ensemble spread, including stochastic model variants, are the subject of further investigation.

## 4. CAFE system

The modeling system we employ is the Climate Analysis Forecast Ensemble (CAFE) system developed by the CSIRO Decadal Climate Forecasting Project (https://research.csiro.au/dfp). This system comprises three key components:

A current-generation climate model; here, the model is ACCESS-D, a variant of the GFDL CM2.1 ocean (MOM4p1)–atmosphere (AM2)–land (LM2)–sea ice [Sea Ice Simulator (SIS)] system (Delworth et al. 2006) described in detail in the appendix.

Coupled data assimilation; specifically, we apply the EnKF-C software of Sakov (2017), which allows for easy implementation of EnOI and ETKF variants. We then assimilate a diverse range of surface and subsurface ocean observations described in section 4a.

An ensemble prediction system capable of generating initial perturbations specific to a given dynamical process at a given lead time, as described in section 4b.

This climate model configuration is described in more detail in the appendix. Throughout we focus on the tropical ocean–atmosphere coupling and teleconnections to the extratropical atmosphere. However, the CAFE system also has the capacity to either restore to or assimilate directly from atmospheric reanalysis datasets. A schematic of the system is shown in Fig. 4, where the components relevant to the current discussion are outlined in red.

### a. Ocean observations and data assimilation configuration

Consistent with the results from the paradigm model (section 3), in the EnOI and ETKF calculations that follow, we assimilate only ocean observations. We consider the recent period January 2002–July 2017, assimilating a comprehensive range of surface and subsurface ocean observations, as detailed in Table 2. More generally, we classify the observational products accordingly as surface satellite observations: SST, sea surface salinity (SSS), sea surface height (SSH), and in situ temperature *T* and salinity *S* data. SSH observations are in terms of derived sea level anomalies from the Radar Altimeter Database System (RADS) altimetry (http://rads.tudelft.nl/rads/rads.shtml) and in situ *T* and *S* observations [Argo, XBT, and CTD as well as TAO, PIRATA and Research Moored Array for African–Asian–Australian Monsoon Analysis and Prediction (RAMA) moorings] are downloaded from the CSIRO climatological atlas (CARS) (http://www.cmar.csiro.au/cars) database and from the WMO GTS (see http://www.wmo.int/pages/prog/www/TEM/GTS/index_en.html).

In any near-operational data assimilation system, choices are required in terms of localization length scales, weighting, or impact factors for various observations. Here, the impact factor is simply a factor for inversely scaling the magnitude of the observation error variance. For both EnOI and ETKF we use an assimilation window of 28 days and localization according to the horizontal distance only. The most pertinent data parameter settings are detailed in Table 3 with the number and type of observations assimilated in two typical data assimilation cycles (March 2012 and December 2017) shown in Table 4. In the results that follow, we have taken the approach of optimizing the parameters for the EnOI and ETKF assimilation systems independently based on minimizing the forecast mean absolute deviation (MAD) and bias at 28-day lead time.

### b. Initial forecast perturbations

In addition to comparing the relative merits of the EnOI and ETKF methods for coupled data assimilation, we further investigate the efficacy of ensemble forecasting using ensembles of BVs initialized about the respective EnOI analyzed state. Two approaches to generating the initial bred perturbations are considered. Cognizant of our aim of focusing on ENSO predictability and the extratropical response to equatorial disturbances, we construct two independent sets of 10 BV ensemble perturbations plus a control in the manner described in Fig. 1. Here, we are primarily concerned with (i) confining perturbations to the tropical regions where ENSO dynamics reside and (ii) having some control of the amplitude of those initial perturbations such that the errors grow consistent with observed thermocline disturbances (Yang et al. 2008). Note that experiments where the entire ocean state was perturbed resulted in saturated SST errors at mid- to high latitudes within a month.

The first approach adds an initial random perturbation scaled to be 1% of the background SD where the background ocean temperature SD is calculated from a 500-yr climate model control simulation (year 1990 forcing). Sampling over a long control simulation allows us to include diverse regimes of tropical variability where the simulated model ENSO has been shown to display substantially differing statistics (Wittenberg 2009). As we are interested only in the coupling between the tropical ocean and atmosphere, we next mask these perturbations such that only those within the region between 20°S and 20°N are retained.

In the second approach, we again estimate the variance in terms of SD for the ocean temperature in the upper 500 m, calculated from anomalies w.r.t. climatology from 500 yr of a control simulation. We then consider the fraction of the total variance at each spatial location partitioned into various temporal bands (in-band variance) calculated using Welch’s fast Fourier transform method (Cooley et al. 1969) applied to the temperature time series at each grid point in the ocean. This method is equivalent to the approach described in the studies of spatiotemporal variance distributions in the earth system by Monselesan et al. (2015) and O’Kane et al. (2016a). This variance is thresholded at values greater than 0.5°C such that the only nonzero values lie within the isosurface shown in Fig. 5. These values are also used to determine rescaling amplitudes (here, 1% of the SD within the isosurface) that can be applied to generate the bred vector ensemble forecast perturbations. In Fig. 5, we show respective isosurfaces for in-band variance between 1 and 2, 8 and 12, 24 and 48, and 48 and 96 months in the upper 300 m and between ±40° latitude. In the 1–2-month band the variance is mostly concentrated in the tropical band and is highly indicative of movement about the equatorial thermocline. At longer time scales the variability is shown to progressively move to regions of the subtropical oceans, as discussed in the studies of Monselesan et al. (2015) and O’Kane et al. (2016a).

Henceforth, we denote our first approach 20°S–20°N and the second as the *isosurface*. We further note that, for the isosurface BVs, due to thresholding surface perturbations are put to zero such that the tropical SST is unchanged from the control whereas, in the 20°S–20°N case the respective ensemble members have retained the SST bred perturbations. The isosurface approach can be seen as generalizing the simple approach of Frederiksen et al. (2010), who apply a mask to bred ocean perturbations or the more general application of a 3D mask based on the estimated analysis error as applied by Ma et al. (2014) to 500-hPa ETKF perturbations. We further note that Wang and Bishop (2003) also applied a masked BV approach to generating atmospheric BVs. Regardless of method, the aim is to reduce error growth (perturbations) in the less spatiotemporally relevant regions. In this way, we ensure that we are sampling from the local low-dimensional subspace of the most dynamically relevant unstable modes.

## 5. CAFE model results and discussion

### a. EnOI versus ETKF GCM assimilation metrics

In the results and discussion that follow, we define the forecast mean absolute deviation as , where are the observations and is the ensemble mean of the background forecasts. Forecast innovation bias refers to the mean magnitude of the forecast innovation; that is, . The accuracy of the analyzed and forecast ENSO (here defined in terms of Niño-3.4 index) is compared to independent metrics calculated from the NCEP–DOE AMIP-II reanalysis (NCEPv2) observational dataset (Kanamitsu et al. 2002). Forecast mean absolute deviations are calculated w.r.t. unassimilated observations at lead times of 28 days. All errors shown are w.r.t. forward independent unassimilated observations.

We compare the efficacy of EnOI and ETKF data assimilation within the CAFE framework using forecast innovation, defined by observations minus forecast (background), bias, and forecast mean absolute deviation (MAD) as metrics. Details of these experiments and their corresponding figures are given in Table 5. Figure 6 contrasts bias and MAD in the region between 20°S and 20°N and for the full ocean depth (see Fig. S2 in the online supplemental material for similar results for the globe) for EnOI and ETKF. We assimilate ocean observations over the Argo period 2002–16, where sufficient subsurface ocean observations are available to constrain the tropical upper ocean. For both DA methods, there is a large reduction in MAD during the first year of the assimilation that occurs as the observations progressively constrain the ocean. Here, we show the averaged forecast ensemble mean bias and MAD for the 96-member ETKF and EnOI.

For global SST (Fig. S2), we note a large seasonally dependent negative (positive) bias (MAD) with a maximum during the austral summer (December). This bias and large forecast error (1-month lead time) arises due to a significant warm bias in the Southern Ocean, where insufficient mixing enables warm water to be trapped in the upper few layers (≈20 m) and is exacerbated by a less prominent cold bias in the high latitudes of the North Pacific and Atlantic during the boreal winter. The ETKF assimilation has the advantage of using the bias correction previously described in section 2e. The combined effect of SST bias correction and flow-dependent covariances is found to dramatically reduce the SST forecast innovation bias to <0.5°C and the forecast MAD to <1°C globally.

For the tropical oceans SST between 20°S and 20°N (Fig. 6), excluding the high latitudes, has reduced the forecast bias to about 0.5°C post-2007 for EnOI with MAD close to that of the ETKF assimilation. For the ETKF, there is a positive mean SST innovation, which reflects a more persistent cold tongue model bias in the equatorial eastern Pacific Ocean. More importantly, the SST bias is halved with ETKF. For sea level anomaly (SLA), the EnOI forecast bias varies between 2 and 5 cm, whereas the ETKF bias is typically in the 1–2-cm range. Despite this, the EnOI and ETKF forecast MAD values are both around 5–7 cm. For the period between 2002 and 2007, we see a progressive and substantial decrease in MAD for *T* and *S* coincident with a corresponding dramatic increase in the number of subsurface observations as the Argo observing array comes on line. As for SST, there is little to distinguish between the EnOI and ETKF temperature forecast MADs, despite a negative bias of about 0.2°C for the EnOI. For salinity, both EnOI and ETKF display little bias. However, the ETKF MAD shows approximately a 0.03-psu reduction over EnOI. The global metrics for the subsurface ocean in *T* and *S*, and, not surprisingly, for SLA, which is a vertically integrated quantity, show very similar results.

A comparison of oceanic and atmospheric mean increments, temporally averaged over the boreal winters (DJF) for the period 2003–16 and along the tropics at 2°S and across the Pacific at 140°W, for both the EnOI (Figs. 7a,c) and ETKF (Figs. 7b,d) systems, reveals strong cross covariances between the tropical ocean and the extratropical atmosphere (Figs. 7a,b). Along 140°W, EnOI ocean increments are largest between 10°S and 10°N, with the spatial pattern of the mean increment centered about the equatorial Pacific thermocline with little surface (SST) expression. ETKF ocean increments in this region are in general larger, extending poleward into the equatorial countercurrent and with significant SST increments around the equator. ETKF produces large increments in the tropical atmosphere through cross covariances and adjustment of convective processes in the tropical Pacific. The EnOI increments in the tropical atmosphere are generally small, apart from the region 10°–25°N below 800 hPa, consistent with small SST corrections in the equatorial region. The increment in the NH atmosphere is more broad scale with both EnOI and ETKF systems showing very similar large positive increments (greater than 0.8°C) above 200 hPa, under which negative atmospheric increments are extending.

Despite there being no comparison to or assimilation of atmospheric observations, the implication here is that the cross covariance between the ocean data and atmosphere reduces the warm bias of the background (or the model); that is, the mean analysis increment is negative (positive), warming (cooling) the boreal winter midlatitude troposphere (stratosphere). In the SH extratropical atmosphere we again see large EnOI increments indicating a systematic shift in the position of the midlatitude jet (in the absence of SST bias correction) and compounded by reduced EnOI ensemble spread at higher latitudes relative to the ETKF (not shown).

Along the equator (2°S), both systems show remarkably similar mean ocean increments, apart from some banding in the ETKF due to asynchronous assimilation of the TAO/TRITON array subsurface observations (Ando et al. 2017). This is not apparent in the EnOI DJF mean increment as it is based on averaging single-realization increments calculated from static background covariances. The ETKF atmosphere temperature increments are clearly associated with the tropical thermocline and SST and act to increase convection (i.e., positive increments). Over the Pacific, the EnOI atmospheric temperature increments are largely confined to the lower atmosphere, with some local topographic influences (e.g., near 45°W) and do not extend above 800 hPa. Over the Indian Ocean, both the EnOI and ETKF systems display very similar patterns indicative of an adjustment of the mean maximum vertical extent of the Hadley cells and of the tropopause. If we interpret the magnitude of the increments as being representative of model and forecast error (at lead times of 28 days), then we might reasonably conclude that the EnOI performs better in the tropics where variability is large scale and ensemble spread might best be captured by time-invariant background covariances. In the tropics, the ETKF spread has large interannual variations unlike the spread in the climatological EnOI background covariances. Outside the tropics, where variability is associated with local processes requiring detailed information about the time-evolving flow, the ETKF performs significantly better, relative to EnOI.

### b. Ocean error growth

Next, we examine the relative growth rates of temperature BV perturbations using the 1–2-month isosurface mask and the latitudinally restricted 20°S–20°N mask. In generating the BV perturbation, the 28-day rescaling time scale has been chosen to allow for coupled tropical ocean–atmosphere instabilities to be given sufficient time to grow (described in section 4b). This choice of rescaling interval is consistent with previous studies of ENSO-related error growth in the tropical Pacific (Yang et al. 2008; Frederiksen et al. 2010). We show the tropical BV growth rates (month^{−1}) by level over the upper 200 m for isosurface perturbation (Fig. 8a) and 20°S–20°N perturbation ensemble prediction experiments (Fig. 8b). Here, we notice that the growth rate of the isosurface BVs at depths close to the thermocline (≈50–150 m) is a factor of 3 greater than for the 20°S–20°N case. Also, unlike the 20°S–20°N case, the threshold chosen for the isosurface results in zero SST perturbations. The initial amplitude of the isosurface BVs are up to an order of magnitude smaller than is the case for the 20°S–20°N BVs. However, despite this, the evolved day 28 perturbations have very similar amplitudes (up to ≈3°C) but with different growth rates and, as we discuss later, with some significant structural differences. In the far-right panels of Fig. 8 we show the correlation of the growth rates with the multivariate ENSO index (MEI; Wolter and Timlin 2011) by depth revealing growth rate profiles. For both correlations, we notice an inflection point at about 100-m depth, with positive correlations in the upper 50 m and significant negative correlations below 100 m. The maximum correlation of between 0.2 and 0.3 is comparable to values reported by Yang et al. (2008).

Cross-correlation sequences (Fig. 9) for the 20°S–20°N case at the surface and at 50 and 100 m reveal the maximum correlation between BV growth rate and MEI to be at 0 lag with a secondary maximum at 6 months. At 150 m, MEI lags the growth of the BV disturbance by 6 months. Similarly for the isosurface case with no SST perturbations, the maximum cross correlation occurs at zero lag at 50 and 100 m with no secondary peak and with a (minus) 6-month lag at 150 m. Given that standard indices of ENSO are based on SST, the high correlations at zero lag at the surface are to be expected. That said, the large negative lag correlation at 6 months shows that disturbances at and below the thermocline lead ENSO phasing. Initial forecast perturbations that project onto disturbances with these growth characteristics will be most effective at capturing the onset of El Niño and La Niña events.

It has been established in the literature that atmospheric ETKF perturbations share many characteristics with BVs (Wang and Bishop 2003; Bowler et al. 2008). However, few studies have examined the link for coupled systems. In Fig. 10, we compare the relationship between ETKF and EnOI spread and analysis increments with BV ensemble spread for SST. In Fig. 10 we show that a small set of BVs can add important flow-dependent information typically lacking in an EnOI analysis and that the spread of those BVs largely projects onto similar structures that are naturally captured in the ETKF ensemble. For the ETKF we plot the mean analysis increment. The BV ensemble members come from the 20°S–20°N perturbation case, whereas the ETKF spread is calculated based on the 96 members from the ETKF ensemble.

As an illustration, we consider the particular period at the end of October 2015 during the lead up to the 2016 El Niño event. In Figs. 10a and 10c we see that the pattern of the BV spread is very similar to that of the ETKF spread, which is indicative of the tropical instability wave. In Fig. 10b, we see that the BV spread correlates with the EnOI forecast error, for which the analysis increment is a proxy. This is to be hoped for, as the BVs are constructed such that they project onto growing error modes, which are then added to the EnOI-analyzed state from the previous cycle. While both the ETKF and EnOI spread project onto their respective analysis increments (Figs. 10c,d), the relationship is weaker for the EnOI case. These results make two key points. First, the spread in BV and ETKF SST perturbations both capture similar flow-dependent information about the tropical instability wave, indicating the utility of a small ensemble of BVs to add flow-dependent information to an EnOI static background covariance matrix. Second, EnOI and ETKF systems show considerable differences in background errors during a period when large disturbances are driving error growth.

More generally, we would like to examine the relationship between BVs and forecast errors in the EnOI and ETKF systems beyond a particular event and at depth. Hence, we examine the seasonal ensemble-averaged BVs compared to the EnOI (Fig. 11) and ETKF (Fig. 12) seasonally averaged analysis increments for sections along 140°W and 2°S in the Pacific Ocean. In both cases, the BVs are localized about regions where the analysis increments are largest (i.e., where the thermocline variability is high). Along 140°W, EnOI increments (Fig. 11) are generally smaller in amplitude and more localized about the maximum in , with high correlation to the BVs during MAM and JJA. However, during SON and DJF, the BVs are generally less structured with the largest amplitude in the northern equatorial region at 10°N and less related to the EnOI increment. Along the equator (2°S), the BVs are highly localized about the thermocline with a significant projection onto the EnOI analysis increments for all seasons except for SON, where the seasonal BV structure is weakest. Similar to the EnOI case, the largest values of the MAM–JJA ETKF analysis increments and BVs are collocated; however, unlike the EnOI case, the meridionally oriented BVs along 140°W for SON–DJF project onto regions where the gradients of the SON and DJF ETKF analysis increments are largest. This would indicate that meridionally oriented fluctuations of the equatorial Pacific thermocline between 5° and 10°N are highly flow dependent and are not captured by the static EnOI background covariances, whose spread is dominated by the climatological vacillations of the thermocline. Along 2°S, the ETKF analysis increments are well correlated with the seasonal BVs.

### c. Extratropical atmospheric response

The primary mechanism for tropical ocean disturbances to influence the midlatitude atmosphere is via anomalous tropical SST driving convection, which in turn modulates the zonally oriented Walker (Bjerknes 1969) and meridionally oriented Hadley (Dima and Wallace 2003) circulations. The atmospheric convective response to tropical ocean BVs in the 20°S–20°N case can be seen in the seasonal cycle of the average outgoing longwave radiation (OLR) BVs (Fig. S3). This response captures the expected asymmetry due to the boreal winter–summer cycle and in regions corresponding to the South Pacific convergence zone (SPCZ) (Vincent 1994) and intertropical convergence zone (ITCZ) (Schneider et al. 2014). In the Pacific the anomalous convective response shows the ITCZ is largely absent (north of the equator) with most anomalous convection being localized over South America, South Africa, and extending over most of northern Australia and into the region of the Pacific Ocean typically associated with the SPCZ. These processes are reflected in the oceanic BVs coupling to the tropical atmosphere via coherent modulation of convection.

### d. Tropical atmospheric response

The primary mechanism for tropical convection to influence the extratropics and the troposphere is via modulation of the Hadley circulation. We diagnose the tropospheric response to tropical convection via the mean meridional mass streamfunction defined as

where *a* is the radius of Earth, *ϕ* is the latitude, is the zonally averaged meridional velocity, *p* is the pressure, is the surface pressure, and *g* is gravity. Here, we show (10^{8} kg s^{−1}) where positive (negative) values correspond to a clockwise (counterclockwise) circulation in the NH (SH) winter cell.

We see qualitatively similar structures in the seasonal Hadley cell (HC) response to tropical disturbances in terms of the BVs for 20°S–20°N (Fig. 13a) and the isosurface (Fig. 13b). In both cases, there are significant counterclockwise anomalies during the SH winter (JJA), acting to strengthen the downward branch of the HC. During the boreal winter (DJF) it appears that the BVs are more confined to the southern HC and act to weaken it (i.e., expressed as positive BV anomalies). For the 20°S–20°N case during the boreal winter (DJF), BV anomalies act to strengthen the northern HC at the equator, which is an indication of enhanced convection, but with some indication of a contraction of the downward branch of the northern HC. The opposite occurs during the boreal spring (MAM), where the northern HC is weakened (negative anomalies) with evidence of poleward expansion of the HC and a stronger Ferrel cell. Similarly, the austral spring (SON) exhibits a weaker southern HC and some evidence of poleward expansion. From these results we conclude that the primary mechanism for tropical convection to influence the extratropics and the troposphere is via modulation of the Hadley circulation.

We next consider the midtropospheric synoptic-scale response to oceanic BV disturbances. Here, the leading eight EOFs are calculated using wintertime midtropospheric 500-hPa geopotential height anomalies constructed as differences between the control and ensemble members over a 28-day period. As each member evolves from the same initial atmospheric state (but with different BV anomalies in the tropical oceans), it is reasonable to assume that the patterns represent the anomalous response of the synoptic features in the midlatitude jets to vacillations in the Hadley and Ferrel cells (Fig. 13a) induced by perturbations in tropical convection. That is, the EOF patterns represent modulations about the simulated extratropical natural variability due to oceanic perturbations.

We observed similar responses to either class of tropical BVs and we discuss only those for the isosurface. An examination of the leading eight EOFs for wintertime midtropospheric 500-hPa geopotential height ensemble-averaged BVs in the Northern (Fig. S4) and Southern (Fig. S5) Hemispheres for the boreal wintertime reveal EOFs 1 and 2 showing features of the Arctic Oscillation (AO; Barnston and Livezey 1987), EOFs 3 and 5 project onto the Pacific–North American (PNA) pattern (Barnston and Livezey 1987) in the Pacific, EOFs 3 and 4 project onto blocking over Eurasia and Scandinavia, EOFs 6 and 8 clearly show features of the circumglobal wave train pattern (CWP; Branstator 2002), and finally EOFs 2 and 7 capture features representative of the North Atlantic Oscillation (NAO; Barnston and Livezey 1987) and blocking in the Gulf of Alaska (O’Kane and Frederiksen 2008a).

Similarly, the leading EOFs of the SH ensemble mean BVs also project onto the well-known SH teleconnection patterns. EOF 1 shows a reasonable representation of the southern annular mode (SAM) (Barnston and Livezey 1987), EOF 2 projects onto the second Pacific–South American mode (PSA2; Lau et al. 1994; O’Kane et al. 2017), and EOF 3 projects onto the wintertime polar and subtropical jets in the Pacific. The leading three EOFs explain ≈56% of the variance whereas the higher-order EOFs 4–8 explain ≈35% of the variance and project onto features consistent with the SH circumglobal wave train and regions of blocking (O’Kane et al. 2016b).

### e. Multiyear ENSO forecasts and predictability

We now examine the utility of the two BV approaches to initialization in multiyear ensemble forecasts of ENSO. The anomaly Niño-3.4 and Niño-4 indices were calculated using monthly means and climatologies over the complete model/reanalysis periods that were formed from the daily inputs. For the forecast runs, in each 2- (5) yr experiment monthly year-by-year anomalies were calculated as deviations from the NCEPv2 climatology. Niño-3.4 and Niño-4 indices were calculated from SST data over the regions 5°S–5°N and 190°–240°E and 5°S–5°N and 160°–210°E, respectively, with area weighting applied. In Fig. 14a, we compare the Niño-3.4 index calculated from the EnOI (D1) assimilation and observations (NCEPv2 reanalysis) for the period January 2002–July 2016. The D1 assimilation follows the observed index reasonably well, considering that the atmosphere is unconstrained and the assimilations have a common 28-day cycle, with the bias, here defined as the difference between D1 and the observed Niño-3.4 index (Fig. 14b), less than 1°C. The forecasts (F0 and F1) are initialized about the January 2007 analyzed states. Forecast F0 corresponds to the addition of 10 bred perturbations and control in the 20°S–20°N ocean domain whereas F1 corresponds to 10 bred perturbations and control in the isosurface. In both cases, the initial atmospheric state of each ensemble member is unperturbed. Here, we show the F1 forecast ensemble spread out to 2-yr lead time, after which errors are saturated. The F0 forecast ensemble spread is shown out to 5 yr for reference.

The ensemble plumes of the raw forecasts show the common bias (Fig. 14b) observed in many seasonal forecast systems, here manifesting as an overestimate of the magnitude of the observed 2008 La Niña. This overshoot is less due to “initialization shock” but rather due to the model equatorial Pacific cold tongue bias reestablishing coincident with the onset of the spring predictability barrier. The F1 forecasts are initialized with perturbations specific to the growing disturbances local to the tropical Pacific thermocline and the major source of forecast error associated with ENSO predictability. Interestingly, the ensemble spread in the isosurface calculation is considerably reduced relative to the 20°S–20°N case.

The relative merits of the two approaches can be better assessed by comparison of receiver operating characteristic (ROC) curves calculated from a large hindcast dataset (a total of 3696 forecast years). In Fig. 15 we compare the ROC curves of the F0 and F1 hindcasts for Niño-4 calculated over a 14-yr period between January 2002 and January 2016 (2-yr lead times, 11 members each starting every month). Here, we see better skill for the F1 (isosurface) BV initiated system relative to the F0 (20°S–20°N) at lead times of several months, indicating that the reduced spread observed in Fig. 14 is in fact an indicator of increased skill. At lead times between 155 and 217 days the F1 forecasts are markedly more skillful relative to F0. The F0 forecast skill drops close to climatology between 186 and 217 days with some reemergence of skill between 248 and 403 days. The F1 forecasts maintain skill throughout the entire forecast out to 465 days. A major difference between the F1 and F0 initial conditions is that SST remains unperturbed in the former, but is perturbed in the latter. Error growth in the isosurface (F1) ENSO forecast ensemble is slower due to a lag, as initial thermocline disturbances take time to express at the surface, whereas for the 20°S–20°N (F0) case, initial SST perturbations couple nearly instantly to the tropical atmosphere, leading to larger and more rapid error growth. The pertinent point is that, despite having no initial expression at the surface, the isosurface perturbations clearly capture the information required to project onto the relevant coupled disturbance where ENSO predictability resides.

## 6. Summary and conclusions

The coupled Lorenz paradigm model simulations point to the potential for a well-constrained ocean state to also constrain the tropical atmosphere, in large part due to the strong coupling between ocean and tropical atmosphere. Here, and within the context of the paradigm model, the additional information captured in the ocean–tropical atmosphere cross covariance is sufficient to enable the tropical atmosphere to track the observed trajectory in phase space. While flow-dependent information constrains the ocean–tropical atmosphere, despite being weakly coupled, cross-covariance information between the ocean and tropical and extratropical atmosphere attractors led to a suppression of the variance in the analyzed extratropics. In contrast, where static rather than flow-dependent cross covariances are employed, the analyzed tropical atmosphere state, despite being strongly coupled to the ocean attractor, fails to track the truth, but the variance of the analyzed extratropical attractor was largely unchanged even with increments due to the cross covariance included. These simple model experiments suggest that one might best initialize ensemble forecasts by constraining the slow modes of the ocean with only a relatively weak large-scale projection of ocean observations into the fast extratropical atmospheric circulation. It is on this basis that we then examined strongly coupled DA variants applied to the GCM, where the ocean is constrained, either with static or flow-dependent cross covariances, and where the large scales of the atmosphere are modified based on suitably scaled ocean–atmosphere cross covariances.

Our focus is on seasonal and longer time scales and, in particular, ENSO. Therefore, our premise underpinning the coupled model assimilation experiments is that predictability primarily resides in the oceans and the fast atmosphere acts as a stochastic driver on the longer-time-scale ocean variability. We considered two approaches to DA based on ETKF and EnOI, assimilating a wide range of ocean observations into a GCM. Outside of the tropics, the ETKF system produced a dramatically lower forecast bias and forecast MAD relative to the EnOI system; however, these improvements were substantially reduced in the tropics. The reason for the reduced analysis error in the EnOI system in the tropics was found to be a result of seasonally dependent fixed ensemble spread at times producing larger observation impacts relative to the tropical ETKF where interannual variations in the background covariances can lead to periods of relatively reduced spread.

BVs representative of growing coupled tropical instabilities were found to modify tropical convection, particularly in the region of the Maritime Continent, which in turn generates a coherent modulation of the Hadley circulation. A direct renormalization of thermocline disturbances was found to be most effective in communicating information from the tropical ocean to the extratropical atmosphere on time scales of a couple of weeks to a month. Comparison of ensemble forecasts based on two types of bred vectors centered about the EnOI-analyzed state reveals a substantial reduction in uncertainties (forecast spread), when disturbances not directly associated with thermocline variability are eliminated. In particular, excluding SST disturbances led to a significant reduction in forecast error and spread in multiyear ENSO predictions and noticeably increased skill at lead times out to 2 years. These results affirm the utility of using BVs explicitly constructed to project onto forecast errors entirely due to tropical subsurface ocean disturbances where the appropriate variance resides in application to ENSO prediction.

The data assimilation experiments and methods discussed form a basis for coupled DA relevant to multiyear near-term climate forecasts. The masked isosurface BV approach allows for the specific targeting of regions of large-scale variability pertinent to dynamical processes that determine predictability on seasonal-to-interannual. spatiotemporal scales. Beyond a season, strongly coupled data assimilation, where the slow ocean modes are explicitly constrained including projection onto the background atmospheric states (i.e., jets, cells, etc.) while leaving the fast atmospheric dynamics (synoptic scales) free, including targeted forecast perturbations, offers a pragmatic approach to determining the mechanisms and predictability of the key climate modes.

This work further highlights the complexity of data assimilation and forecast initialization in nonlinear multiscale systems. While we have demonstrated the advantages of flow-dependent ocean data assimilation and the usefulness of ocean observations to constrain the large scales of the atmosphere, it is apparent that the assimilation of atmospheric observations is required to guarantee the correct extratropical variability. This is a focus of our ongoing work, as are methods to identify an appropriate theoretical basis necessary for identifying causal relationships between climate modes and for determining predictability on given spatiotemporal scales of interest and as a basis for developing a generalizable approach to multiscale forecast initialization.

## Acknowledgments

The authors thank the reviewers (four anonymous and Eugenia Kalnay) for their considered comments, which led to substantial improvements to the article. The authors were supported by the Australian Commonwealth Scientific and Industrial Research Organisation (CSIRO) Decadal Forecasting Project (https://research.csiro.au/dfp). TJO and DPM acknowledge the support of the Earth Systems and Climate Change Hub through Project 2.3 “Towards an ACCESS decadal prediction system” (http://nespclimate.com.au/decadal-prediction/). NCEP Reanalysis 2 data were provided by the NOAA/OAR/ESRL/PSD, Boulder, Colorado, from their website (https://www.esrl.noaa.gov/psd/).

### APPENDIX

#### Climate Model Configuration

The CAFE climate model is based on the Geophysical Fluid Dynamics Laboratory’s (GFDL) Climate Model 2.1 (CM2.1; Delworth et al. 2006). The ocean component of the CAFE model is Modular Ocean Model, version 4.1 (MOM4p1; Griffies 2008), which is configured consistent with CM2.1, but using the 1° ocean grid as described by Bi et al. (2013). The ocean model is coupled to the land, atmospheric, and sea ice components from CM2.1, namely, LM2, AM2, and SIS, respectively. Output from this parent model, CM2.1, was included in GFDL submissions to the Coupled Model Intercomparison Project (CMIP), phases 3 and 5, and has been the basis of studies of predictability in the climate system (Song et al. 2008; Wittenberg et al. 2014; Zhang et al. 2017). CM2.1 has also been used to assimilate data and produce seasonal-to-interannual forecasts (Zhang et al. 2007; Yang et al. 2013). Additional details about the CAFE model configuration are described below:

The nominal resolution of the CAFE ocean component (Modular Ocean Model, version 4.1) is 1°, with extra latitudinal resolution in the tropics, and 0.33° at the equator, with extra horizontal resolution in the Southern Ocean, corresponding to 0.25° at 75°S. There are 50 vertical levels, with a 10-m resolution in the upper ocean, increasing to ~300 m at abyssal depth. The grid is tripolar over the Arctic, north of 65°N, to avoid the North Pole singularity.

Subgrid processes for the CAFE ocean component are adopted from CM2.1, including neutral physics (Redi diffusivity and Gent–McWilliams skew diffusion) (Griffies 2004), a Brian–Lewis vertical mixing profile, a Lagrangian friction scheme, and a

*K*-profile parameterization for the mixed layer calculation.The sea ice model (SIS) is interpolated onto the ocean grid and uses five ice categories.

The atmospheric model (AM2) has a resolution of 2° in latitude and 2.5° in longitude, and 24 hybrid (sigma-pressure or terrain-following pressure) vertical levels. The land model (LM2) is on the same horizontal grid.

Concentrations of atmospheric aerosols and radiative gases, and land cover are based on 1990 conditions.

The model’s climatology fields from a control experiment (a single 500-yr integration) of the CAFE climate model are shown in Fig. S1 in the online supplemental material. The climatology fields, which are calculated from the years 0200 to 0299, all closely resemble results from the parent CM2.1 model (Delworth et al. 2006). Sea surface temperature and salinity also show good comparisons to observations. The model does however exhibit biases in common with the parent model, including a dry bias over South America, and underestimates sea ice extent in the south, and particularly during the austral summer. In the control experiment shown (Fig. S1), bias in the mode water structure is reduced by restoring the temperature and salinity of water below 2000 m to observations in the control simulation; the justification here is that this work is primarily interested in processes on seasonal-to-interannual time scales rather than slower processes in the deep ocean. The initial ETKF states (i.e., January 2002) use randomly sampled Januaries from this long (500 yr) control simulation. The EnOI seasonal (monthly) varying background covariances are constructed from the last 100 years of the control simulation with the initial BV states constructed from random perturbations added to the analyzed EnOI state.

## REFERENCES

*Fundamentals of Ocean Climate Models*. Princeton University Press, 518 pp.

*Nonlinear and Stochastic Climate Dynamics*, C. L. E. Franzke and T. J. O’Kane, Eds., Cambridge University Press, 136–158.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0189.s1.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).