Coupled atmosphere–ocean data assimilation experiments with a low-order climate model

Abstract

A simple idealized atmosphere–ocean climate model and an ensemble Kalman filter are used to explore different coupled ensemble data assimilation strategies. The model is a low-dimensional analogue of the North Atlantic climate system, involving interactions between large-scale atmospheric circulation and ocean states driven by the variability of the Atlantic meridional overturning circulation (MOC). Initialization of the MOC is assessed in a range of experiments, from the simplest configuration consisting of forcing the ocean with a known atmosphere to performing fully coupled ensemble data assimilation. “Daily” assimilation (that is, at the temporal frequency of the atmospheric observations) is contrasted with less frequent assimilation of time-averaged observations. Performance is also evaluated under scenarios in which ocean observations are limited to the upper ocean or are non-existent. Results show that forcing the idealized ocean model with atmospheric analyses is inefficient at recovering the slowly evolving MOC. On the other hand, daily assimilation rapidly leads to accurate MOC analyses, provided a comprehensive set of oceanic observations is available for assimilation. In the absence of sufficient observations in the ocean, the assimilation of time-averaged atmospheric observations proves to be more effective for MOC initialization, including the case where only atmospheric observations are available.

1 Introduction

Motivated by the need to properly address near-term (i.e., interannual to interdecadal) climate prediction as an initial-value problem, the development of improved initialization methods involving advanced data assimilation (DA) applied to coupled atmosphere–ocean global climate models (AOGCMs) has recently been identified as a priority (Meehl et al. 2009; Shukla et al. 2009; Mehta et al. 2011). This is supported by the findings of Sakaguchi et al. (2012), which indicate that without proper initialization, state-of-the-art AOGCMs are only skillful at predicting the externally forced climate variability at continental and larger scales over horizons of several decades and longer. Improved initial conditions (ICs) are required for sharper predictions of internal climate variability at regional scales.

Ideally, ICs should correspond to a coherent representation of atmospheric and oceanic states characterizing the current (i.e., at initial time) phase of internal variability of the climate system. In particular, some experiments have shown that a proper description of the low-frequency variability in the ocean, where the main sources of predictability are found (Collins 2002; Boer 2011), leads to more skillful predictions (Troccoli and Palmer 2007; Smith et al. 2007; Pohlmann et al. 2009; Branstator and Teng 2012; Garcia-Serrano and Doblas-Reyes 2012). Accurate initialization of the Atlantic Meridional Overturning Circulation (MOC) is of particular interest since it is believed to be an important carrier of predictability in the North Atlantic (Rahmstorf 2002; Latif and Keenlyside 2011; Boer 2011; Srokocz et al. 2012).

A number of strategies have been tested to initialize AOGCMs. The simplest technique uses global atmospheric reanalyses (e.g., Uppala et al. 2005; Dee et al. 2011; Rienecker et al. 2011; Compo et al. 2011) to spin-up an ocean model toward states consistent with the imposed atmospheric forcing (e.g., Matei et al. 2012). A slightly more comprehensive approach consists of incorporating information about the ocean using sea surface temperature (SST) analyses (e.g., Keenlyside et al. 2008). However in both instances, spin-up of the slowly evolving deep ocean and associated large-scale ocean circulations (e.g., the MOC) solely relies on the forced solutions of an ocean model, often plagued by large errors. Next in the hierarchy of initialization strategies is an approach that combines reanalyses of the atmosphere and ocean (e.g., Troccoli and Palmer 2007; Doblas-Reyes et al. 2011; Matei et al. 2012). Some coherence between atmospheric and oceanic states is maintained by using common boundary conditions. Independent SST analyses are used as a lower boundary condition by atmospheric reanalyses and in turn analyzed surface winds or atmosphere–ocean fluxes are used as upper boundary conditions by ocean reanalyses (e.g., Carton and Giese 2008; Balmaseda et al. 2008; Köhl and Stammer 2008; Masina et al. 2011). In addition, upper ocean temperatures in reanalyses are often nudged toward the SST analyses (Saha et al. 2010). Despite such care, Munoz et al. (2011) found significant differences in MOC characteristics between different ocean reanalyses, suggesting a sensitivity to the DA methodology and model error. The most comprehensive method, and the one considered here, consists of solving the fully coupled DA problem. Zhang et al. (2007) show that joint assimilation of atmospheric and oceanic observations can provide successful reconstructions of ocean heat content over large portions of the global ocean. Despite the comprehensive character of this initialization method, important challenges remain, including initializing the low frequency component of the MOC in the absence of comprehensive oceanic observations (e.g., Zhang et al. 2010).

Although prior research has tested a range of initialization strategies, several fundamental questions remain for identifying the most appropriate method used to initialize near-term climate predictions. This investigation focuses on the following questions:

• 1. Is DA needed?

• 2. If DA is needed, is the problem primarily an ocean DA problem, or is coupled DA a fundamental requirement?

• 3. What is the most efficient approach for generating ICs suitable for climate predictions?

A significant challenge in approaching these questions concerns the fact that suitable ICs require initializing components with a wide range of characteristic timescales. The fast atmospheric component is well-observed while the slow oceanic component is the source of long-range predictability but similarly comprehensive oceanic observations are unavailable, particularly prior to the deployment of Argo floats (Roemmich et al. 2009). Moreover, model error can be a significant problem. This latter issue can be dealt with by applying empirical bias corrections (Stockdale 1997) or by analyzing anomalies rather than full model states (Schneider et al. 1999; Pierce et al. 2004; Smith et al. 2007; Keenlyside et al. 2008; Pohlmann et al. 2009; Liu et al. 2012). Although we acknowledge the importance of model error, it is not addressed in the present study so that we may focus on the questions in coupled atmosphere–ocean initialization identified above.

A significant challenge to working on these problems involves choosing a tractable approach to test concepts, and evaluate different strategies. Using complex AOGCMs is problematic in this context due to the large computing requirements associated with running multiple ensemble DA experiments over long periods of time. Also, AOGCMs have not yet converged in their ability to simulate internal climate variability (Branstator et al. 2012). Therefore, simpler low-order models remain useful alternatives for investigating fundamental issues in climate variability and predictability (e.g., Shaffer and Olsen 2001; Zickfeld et al. 2004; Longworth et al. 2005; Taboada and Lorenzo 2005; Lucarini and Stone 2005; Colin de Verdière 2010; Stone and Krasovskiy 2011; Ou 2012). Here, aspects in the initialization of the low-frequency component of the MOC are explored using a low-order coupled atmosphere–ocean climate model with an ensemble Kalman filter DA approach.

The remainder of the article is organized as follows. Section 2 presents a description of the low-order model and its variability and predictability characteristics. Section 3 explores the importance of performing comprehensive data assimilation, compared to the simpler approach of forcing an ocean model with a known atmosphere (e.g., reanalyses). The impact of assimilating time-averaged observations on analysis accuracy is investigated in Sect. 4. Results from a comprehensive set of data denial experiments are presented in Sect. 5 to characterize the performance of a set of DA strategies over a wide range of climate states. Conclusions are summarized and discussed in Sect. 6.

2 Model characteristics

2.1 Model description

We use a modified version of the low-order coupled atmosphere–ocean climate model described in Roebber (1995), characterized by a fast chaotic atmosphere and oscillating or multiple equilibria solutions in the slowly overturning ocean (Roebber 1995; van Veen et al. 2001; Taboada and Lorenzo 2005). The model is fully defined in the Appendix, and here we only provide a descriptive summary. The model, illustrated schematically in Fig. 1, consists of coupling the wave–mean-flow atmospheric circulation model of Lorenz (1984, 1990) to the MOC box model of Birchfield (1989) and Birchfield et al. (1990) as an idealization of the midlatitude northern Atlantic Ocean.

Fig. 1

Schematic representation of the low-order climate model. a Model geometry and state variables and b cross section of the 3-box ocean model

The atmospheric model is composed of five variables: the meridional gradient of tropospheric temperature, or equivalently the strength of the tropospheric zonal flow, the amplitudes of cosine and sine phases of large scale transient eddies superimposed on the mean flow, and high- and low-latitude air temperatures. Internal flow interactions are represented through the amplification of eddies from the strength of the zonal flow (i.e., baroclinic eddy growth), the displacement of the eddies by the zonal flow and mechanical damping. The ocean model consists of a three-box representation of the North Atlantic basin: two representing the subtropical and subpolar upper ocean and the third the deep ocean (Fig. 1b). Upper ocean temperature is influenced by diffusion with the atmosphere, while both temperature and salinity within each box are determined through advection by the MOC and through diffusion between the upper and deep ocean. The strength of the MOC is linearly related to the meridional gradients of upper-ocean temperature and salinity (see Stommel, 1961).

Atmosphere–ocean interactions in the model consist of a linear modulation of the atmospheric zonal flow by the meridional gradient in upper ocean temperature, while atmospheric eddy amplitudes are influenced by the temperature contrast between a mid-latitude continent and the subpolar upper ocean. The upper ocean and the atmosphere exchange heat through diffusion, and upper ocean salinity is affected by a simple representation of the hydrological cycle; e.g., evaporation over the subtropical ocean, poleward water vapor transport by atmospheric eddies and associated influx of freshwater by precipitation into the subpolar box.

2.2 Variability and predictability

Some perspective on variability and initial-value predictability in the low-order model is gained from a 100-member 5,000-year ensemble simulation. The initial ensemble is constructed around a reference model state, on which small amplitude white noise is added. Initial conditions of the reference member are set such that the model solutions are in an oscillatory mode, with an initial overtuning circulation of 15 Sv (1Sv = 10⁶m³ s ⁻¹) (Fig. 2). This value is qualitatively in accord with the modern-era MOC (Munoz et al. 2011; Rayner et al. 2011). Simulated MOC variability is characterized by a superposition of large (∼10 Sv) multi-centennial oscillations and smaller amplitude (∼2–3 Sv) interdecadal fluctuations. The lower-frequency variability has an amplitude similar to fluctuations found in simulations performed with a more complex model (Schulz et al. 2007), while the strength of interdecadal variability is consistent with results from comprehensive climate models (e.g., Menary et al. 2012; Medhaug et al. 2012; Wouters et al. 2012).

Fig. 2

Ensemble simulation from perturbed initial conditions, a over the entire 5,000 years and b zoomed over the first 500 years. The deterministic reference solution (truth) is shown by the solid black line, the ensemble mean is the solid red line and ensemble spread is represented by the pink shaded area. The climatological mean and standard deviation from the deterministic simulation are shown by the circle and gray vertical line on the extreme right of the graphs

As noted by Murphy (1988), loss of predictability occurs when the ensemble mean can no longer be distinguished from the climate mean and when the ensemble variance reaches that of randomly chosen states (i.e., climatology). The ensemble mean gradually departs from the truth solution, converging toward the model’s climatological mean as the ensemble variance increases toward climatology (Fig. 2a), with an e-folding time of about 350 years. Smaller amplitude interdecadal fluctuations are more apparent when focusing on the first 500 years, as is the lack of skill of the ensemble mean at representing these variations (Fig. 2b).

3 Initialization of climate predictions: Is data assimilation necessary?

The contribution of data assimilation in generating accurate initial conditions for improved near-term climate predictions remains to be established at a fundamental level. Such a role is assessed here from the point of view of how to effectively and efficiently constrain the low-frequency evolution of the MOC in the low-order model. Two approaches are contrasted. The first approach does not involve ocean data assimilation, and reproduces the established practice of forcing an ocean model with atmospheric reanalyses to “spin-up” the ocean toward realistic states (Gulev et al. 2003; Brodeau et al. 2010). This approach has shown some success in the limited context of short-term hindcasts of interannual MOC variability (Matei et al. 2012). Despite these encouraging results, the time interval required to spin-up the ocean to states accurate enough to act as viable initial conditions for robust longer-term climate predictions remains an open question. The second approach involves comprehensive coupled atmosphere–ocean DA.

3.1 Forced–ocean experiment

Experiments are carried out reproducing the practice of forcing an ocean model with atmospheric reanalyses. Here we define a truth solution as the reference member of the fully coupled ensemble model run discussed in Sect. 2.2. An ensemble of initial ocean states is built from one hundred random draws from this “truth” simulation. Hence initial conditions in the ocean only reflect knowledge of the model climatology, without any specific information about the true initial state. Every member of the ensemble is then integrated forward with an identical sequence of time-dependent atmospheric states as boundary conditions taken from the truth solution. Forcing the ocean with the perfectly known atmosphere is performed at every time step. Constraining the evolution of ocean variables with the true atmosphere is expected to lead to the convergence of ensemble members toward the true solution, resulting in a reduction of ensemble spread over time.

Figure 3 shows the evolution of ensemble spread for the simulated MOC, averaged over five distinct 100-member ensemble runs each covering a period of 5,000 years. The five runs differ in their initial conditions used to define truth. As a reference, results from ensembles of fully interactive simulations (i.e., freely evolving coupled atmospheric and oceanic states simulated using the complete model) initialised using concurrent atmospheric and oceanic states randomly drawn from the truth simulation show an average e-folding time of about 9,000 years. This reflects the slow convergence of the solutions toward a fixed point. In comparison to this weak intrinsic damping, ocean ensembles driven by the single true atmosphere (i.e., forced mode) exhibit a faster decrease in spread, but the e-folding time is still greater than 3,000 years. These results suggest that initializing the low-frequency component of the MOC solely by forcing it with atmospheric reanalyses is impractical.

Fig. 3

Temporal evolution of MOC ensemble spread averaged over five realizations of 100-member ensembles run in free (solid blue line) and forced (solid red line) modes. Fits by exponential functions are also shown, with the corresponding e-folding times indicated in the legend

3.2 Coupled atmosphere–ocean data assimilation

We next consider atmosphere–ocean data assimilation in the idealized scenario of a well-observed ocean, with temperature and salinity observations in all three ocean boxes. Data assimilation is performed with an ensemble Kalman filter (EnKF) (Evensen 2003), based on a formulation using perturbed observations with sequential observation processing (Houtekamer and Mitchell 1998). The EnKF method has been chosen as it provides a natural framework for investigating the central issue of covariance estimation in the coupled atmosphere–ocean system and its role toward improved initialisation of probabilistic coupled climate forecasts. To perform ensemble DA experiments under realistic conditions, where a compromise has to be reached between maintaining acceptable computational cost while limiting covariance sampling errors, all experiments are performed with 100-member ensembles. Good performance of the DA system is ensured by proper ensemble calibration, using inflation applied to the background ensemble to achieve overall correspondence between background ensemble variances and mean squared forecast errors (e.g., Houtekamer et al. 2005). Covariance localization is not performed in this simplified low-dimensional model.

A perfect model framework is used in this and subsequent DA experiments. Observations are taken as instantaneous states drawn daily from the truth simulation, plus random error whose statistics are specified as follows. For atmospheric variables, error standard deviations are 10 % of the standard deviations obtained from the 5,000-year climatological simulation. For ocean variables, error standard deviations are set at 0.5 K for temperature and 0.1 psu for salinity as in Zhang et al. (2007), representing approximately 35 % of the corresponding climatological variability in respective ocean variables. These values reflect our inability to observe the ocean as well as the atmosphere. A 50-year period is chosen from the long-range reference simulation, defining truth and corresponding set of observations. DA is performed over that period, using an initial background ensemble defined using model states randomly chosen at different times over the entire 5,000 years of the truth solution. Analyses of the MOC are shown against the true solution to assess the skill in recovering this key ocean variable, which provides a concise summary of the results (note that the MOC is not an assimilated variable in any experiment).

For the first experiment, daily observations of atmospheric variables and ocean temperature and salinity in all three boxes are assimilated. The daily assimilation of this complete set of observations leads to a rapid adjustment of analyses toward the true state. Ensemble-mean analysis errors for the assimilated atmospheric and oceanic state variables are reduced to values within observation errors after the first DA cycle (not shown). Furthermore, the information from DA is propagated to the MOC, as shown by the initial error in the ensemble mean background states reduced by 70 % over the first three cycles (Fig. 4a). Good agreement of analyses is subsequently maintained over the entire DA period (Fig. 4b). Serving as a benchmark, this experiment shows, perhaps not surprisingly, that the daily assimilation of a comprehensive set of observations in the atmosphere and ocean leads to accurate analyses of the simplified coupled atmosphere–ocean system, including the variability of the overturning circulation at the annual and decadal time scales. However, a comprehensive set of ocean observations available for assimilation is the exception rather than the norm, particularly in the case of hindcasts initialized prior to the late twentieeth century. Data denial experiments, designed to assess DA performance under scenarios of limited availability of observations, are presented in the next section.

Fig. 4

Temporal evolution of ensemble-mean MOC (a) background and analyses compared against the truth (black dots) over the first month of data assimilation, and (b) over the entire DA period, from the daily assimilation of observations of all atmospheric variables, and of ocean temperature and salinity in all three boxes

3.3 Data assimilation with limited ocean observations

A series of DA experiments are conducted in which progressively fewer observations of ocean variables are assimilated. These were designed to broadly reflect the historical evolution of the ocean observing system, from the modern era back to pre-twentieth century conditions when few ocean observations were available:

• 1. Atmospheric variables and upper ocean temperature and salinity;

• 2. Atmospheric variables and upper ocean temperature;

• 3. Atmospheric variables only.

A comparison of MOC analyses resulting from these various DA configurations is shown in Fig. 5. Accurate analyses are still obtained when both upper ocean temperature and salinity are directly constrained by DA. This is not surprising as the MOC is driven by the meridional gradients of these variables. When salinity is not assimilated, errors in the analysis ensemble mean develop during the first year and persist for the remainder of the DA period. Daily DA of atmospheric observations only is not effective at constraining the ocean, leading to increasing errors in MOC analyses. The ability of a daily-cycling ensemble DA system to produce accurate MOC analyses, including the representation of the low-frequency variability of the large-scale overturning ocean circulation, depends on the availability of observations in the ocean. These findings are consistent with those obtained by Dunstone and Smith (2010) and Zhang et al. (2010) using comprehensive AOGCMs. The results also serve to reinforce the conclusions drawn from the forced–ocean experiment discussed earlier; i.e., assimilation of oceanic observations is an important component in the initialization of the low-frequency component of the MOC.

Fig. 5

Time series of analyzed MOC states resulting from the coupled daily assimilation of atmospheric observations and various sets of oceanic variables. The truth is shown by the solid black dots

4 Assimilation of time-averaged observations

A key element in effective coupled data assimilation involves maximizing the impact of any covariability between the well-observed atmosphere and the more sparsely observed ocean. Because of the difference in characteristic time scales between the two components, covariances between the noisy, fast, atmospheric and slow oceanic variables sampled at a high frequency are expected to be weak and therefore dominated by noise from other sources such as sampling error in the EnKF. Here we explore the use of averaging over the noise as an effective way of increasing the signal-to-noise ratio for DA.

4.1 Climatological correlations between the MOC and the time-averaged atmosphere

Output from the 5,000-year truth simulation is used to estimate correlations between the MOC and the other model variables. Correlations with oceanic variables (not shown) are generally important (e.g., in excess of 0.6 for upper ocean variables) and are independent of the averaging time intervals considered. In contrast, a dependence with averaging time scales is observed for atmospheric variables (Fig. 6). Correlations are weak at the daily time scale, and undergo only a marginal increase for subannual averaging intervals. The positive correlation between the MOC and the zonal wind only reaches a maximum of 0.2 when variables are averaged over a decade. This increase reflects an enhanced link between the frequency of strong atmospheric zonal flow events associated with low-frequency phases of strong overturning circulation as discussed in Roebber (1995). The strongest correlation is found for the atmospheric eddy field and, in particular, the eddy energy (sum of the squared eddy phase amplitudes). A sharp increase in the correlation magnitude occurs for averaging intervals approaching 1 year. The physical link for this relationship in the low-order model is through the eddy energy’s influence on the strength of the atmospheric branch of the hydrological cycle (see Eq. 6 in the Appendix), bringing freshwater to the subpolar ocean and weakening the MOC. Although simplified, the link between atmospheric dynamics and the MOC as represented in the low-order model is qualitatively consistent with the observed influence of atmospheric variability on deep-water formation by the local forcing of ocean-atmosphere surface heat and freshwater fluxes (Clarke and Gascard 1983; Zaucker et al. 1994; Timmermann et al. 1998; Delworth and Greatbatch 2000).

Fig. 6

Correlation between time-averaged atmospheric variables and the MOC as a function of averaging time scale

4.2 The assimilation algorithm for time-averaged observations

The stronger correlation for the longer time scales suggests that assimilation of time-averaged observations (Dirren and Hakim 2005; Huntley and Hakim 2010), may be a useful approach for coupled atmosphere–ocean DA. Using time-averaged observations also increases the time between assimilation steps, improving efficiency.

The method consists of the same steps as for canonical ensemble Kalman filters, except that time-averaged, rather than instantaneous, variables are assimilated. This assumes that the observation operator is not a function of time, and that the deviations from the time mean covary weakly with the time-averaged observations. For a given DA cycle, observations are collected over a predetermined time window and an ensemble of numerical model runs is created over the same time interval with computing model states projected in observation space. An averaging interval of 1 year is used here as it corresponds to the shortest time scale characterized by larger atmosphere–ocean correlations. Time-mean values and deviations from the time mean are computed for observations and model ensemble estimates. Background covariances and innovations (differences between observations and the corresponding model states in observation space) are then estimated on the basis of time-averaged values. The Kalman filter update equation is then applied to obtain updated time-averaged states independently from the deviations from the time mean. These deviations, assumed uncorrelated with the time-mean observations, are simply added to the updated time-averaged states to recover the full updated states. DA is otherwise performed as described in Sect. 3.2, except for analyses being generated only once a year (yearly cycling). Appropriately scaled observation error statistics are also used, e.g., standard deviations reduced by a factor of $\sqrt{N}$, where N is the number of observations used in calculating the time average. For more details, the reader is referred to appendices 1 and 2 in Huntley and Hakim (2010).

4.3 Assimilation experiments with time-averaged observations

A similar set of data denial experiments as in Sect. 3.3 is carried out to assess the impact of assimilating yearly-averaged observations of model state variables. The resulting ensemble mean analyses of time-averaged MOC (Fig. 7) are of comparable accuracy with daily DA (Fig. 5) when time-averaged ocean observations are assimilated. Also, a similar degradation in MOC analyses is observed when salinity observations are not assimilated. However, errors in time-averaged MOC analyses are smaller (by 50 %) compared to daily DA when only atmospheric observations are assimilated, due to the enhanced covariability with the time-averaged atmosphere. In spite of this reduction in errors, MOC analyses generated from atmospheric-only DA still do not accurately track with the true MOC and its interdecadal trends.

Fig. 7

Same as Fig. 5 but for experiments with the yearly assimilation of time-averaged observations

These results show that assimilating time-averaged observations in a coupled DA framework leads to MOC analyses of accuracy similar to that with daily DA in a well-observed ocean and to marginally more accurate analyses when fewer or no oceanic observations are assimilated. This is consistent with the fact that covariability between yearly averaged atmospheric state variables and the MOC are stronger than at the daily time scale but nevertheless remain small.

4.4 Assimilation of atmospheric eddy energy

The use of alternative atmospheric variables is tested by carrying out the same set of data denial DA experiments but with the assimilation of eddy energy rather than the individual eddy phase amplitudes. As previously shown, this nonlinear variable covaries more strongly with the MOC. Pseudo-observations and observation error statistics for eddy energy are determined using the same methodology as for other variables in the earlier experiments. Not surprisingly, MOC analyses generated by including the assimilation of comprehensive ocean observations are comparable to those resulting from the assimilation of individual eddy amplitudes (not shown). However, the stronger covariability between eddy energy and the MOC is evident when DA is performed only in the atmosphere (Fig. 8). Despite the absence of oceanic observations, the enhanced covariability between the MOC and atmospheric eddy energy leads to a gradual reduction in MOC analysis errors. Good agreement with the truth is obtained after about 10 years, illustrating the cumulative effect of a stronger constraint on the MOC provided by the cycled assimilation of eddy energy observations.

Fig. 8

Time series of analyzed MOC states resulting from the coupled assimilation of atmospheric observations only, with the amplitude of eddy phases as in Fig. 7, compared to the assimilation of eddy energy. The truth is shown by the solid black dots

5 Solution convergence over a large sample of climate states

The results from a single case above serve to illustrate some of the fundamental aspects in the performance of various coupled DA configurations. A more robust assessment is obtained by characterizing the evolution of the mean absolute error (MAE) in MOC ensemble-mean analyses from 100 realizations of each DA experiment described in the previous sections. Each realization has a distinct truth time sequence (e.g., different climate regime) defined by randomly choosing a different 50-year segment within the truth simulation. Initial ensembles are again composed of 100 randomly chosen model states from the truth simulation.

MAE values from MOC analyses obtained from the different DA configurations are presented in Fig. 9. Error statistics from 100 realizations of 100-member ensemble mean forecasts (no DA performed) are also shown as a baseline estimate of MOC errors in the absence DA. Mean forecast errors are shown to have a magnitude of about 2.5 Sv over the entire 50-year forecast horizon. With ocean DA (Fig. 9a, b), MAE values are reduced by about 40–50 % after the first year of assimilation (i.e., a single cycle in the case of time-averaged DA), and subsequently continue decreasing over several cycles until a minimum is reached. Not surprisingly, the most accurate MOC analyses are obtained when both ocean temperature and salinity observations are assimilated (Fig. 9a). MOC analyses are accurate to within 0.2 Sv after approximately 20–30 years of DA when instantaneous daily and yearly time-averaged observations are assimilated. Under this scenario of a well-observed upper ocean, the assimilation of atmospheric eddy energy does not provide significant advantages.

Fig. 9

Temporal evolution of the mean absolute ensemble-mean analysis error for the MOC, from 100 realizations of DA experiments using different assimilation configurations. Results obtained from model ensemble predictions (no assimilation) are also shown for reference. The shaded areas indicate 95 % confidence intervals

When upper ocean salinity is not assimilated (Fig. 9b), slightly lower levels of accuracy are obtained. In this scenario, the least accurate analyses and slowest rates of error reduction are obtained with daily DA. The assimilation of yearly averaged observations leads to faster convergence and more accurate analyses by 40 %. Moreover, the assimilation of eddy energy leads to accurate MOC analyses more rapidly. Minimum error is reached between 10 and 15 years of cycling when eddy phase amplitudes are assimilated, compared to 5–10 years when eddy energy is assimilated. The largest errors are obtained in the absence of DA in the ocean (Fig. 9c). In particular, daily assimilation of atmospheric observations fails at generating MOC analyses converging toward the truth. After an initial decrease in errors during the first year of DA, errors increase again toward an upper limit corresponding to the errors obtained without DA altogether. The assimilation of time-averaged atmospheric observations is much more effective as evidenced by the gradually decreasing MOC errors. These results also confirm the important advantage provided by the assimilation of time-averaged atmospheric eddy energy observations over any other DA configuration. MOC errors are reduced by 60 % compared to analyses obtained as a result of the assimilation of individual eddy phase amplitudes.

6 Summary and discussion

Experiments are performed with a simple low-order coupled climate model and an ensemble Kalman filter to gain insights on fundamental issues in the initialization of the coupled atmosphere–ocean system. Specific overarching questions are investigated in a simplified idealized scenario involving the meridional overturning circulation (MOC) as the key oceanic component driving the interannual to interdecadal climate variability. Various initialization strategies are considered, ranging from a simple approach requiring no specific data assimilation (DA) to fully coupled DA. DA is also considered in its traditional form (e.g., performed at the daily time scale) and a more recently developed approach involving time-averaged observations.

Question [1] asks whether data assimilation (DA) provides value for the initialization of the MOC, compared to the simpler approach of spinning–up an ocean model with known atmospheric states (e.g., reanalyses). An idealized experiment mimicking this simpler approach reveals that convergence of MOC solutions toward the true low-frequency variability occurs over very long time scales. This approach is therefore impractical for obtaining timely analyses of the coupled atmosphere–ocean system. In contrast, daily DA rapidly leads to accurate MOC estimates provided that a comprehensive set of observations is available for assimilation.

Question [2] asks if MOC initialization is primarily an ocean DA problem or whether a fully coupled approach provides value. Data denial experiments indicate that answers lie along a continuum characterized by a strong dependence on the set of oceanic observations available for assimilation. For a well-observed ocean, assimilation of oceanic observations has a dominating effect on the accuracy of the resulting MOC analyses. When the ocean is only partially observed or not observed at all, the value of coupled DA increases significantly, particularly in a configuration assimilating time-averaged observations. Fully coupled time-averaged DA is essential when only the atmosphere is observed. The positive impact of assimilating a higher-order atmospheric variable that projects more strongly on ocean states (e.g., eddy energy in this simple model) is also demonstrated. This is especially important to consider when generating hindcasts initialized at times predating the availability of comprehensive oceanic observations.

Question [3] is concerned with the identification of efficient DA approaches for generating timely and accurate analyses of the coupled atmosphere–ocean system. Simply stated, how long into the assimilation does it take before accurate MOC analyses are obtained? Daily DA in the ocean and/or atmosphere is compared against the more efficient alternative involving less frequent assimilation of time-averaged observations. Results show that time-averaged DA is more efficient, with significant gains in convergence rates in MOC analysis errors. This advantage becomes more important as less oceanic observations are assimilated, particularly when atmospheric eddy energy is assimilated.

These results indicate that averaging over high-frequency atmospheric states is an efficient way of exploiting the stronger atmosphere–ocean covariances characterizing the longer time scales. These enhanced covariabilities originate from integration of the atmospheric white noise by the ocean (Hasselmann 1976), leading to low-frequency covariability between the ocean and atmosphere (e.g., Farnetti and Vallis 2011). In ensemble DA, stronger covariances are a distinct advantage as they can be estimated with greater confidence using a reasonable ensemble size (e.g., 100 in the present experiments). On the other hand, sampling errors associated with such an ensemble size may hinder the estimation of weak atmosphere-MOC covariability characterizing short time scales. This is likely a factor in the ineffectiveness of daily atmospheric DA at producing accurate MOC analyses.

Greater insight into the role of sampling errors is gained by estimating minimum ensemble sizes required for reliable estimation of correlations between the MOC and other model variables (Fig. 10). Very large ensembles are required to estimate with high confidence the weak daily atmosphere-MOC correlations. Sampling errors become less important when yearly-averaged variables are considered, allowing the use of smaller ensembles, particularly when eddy energy is considered. On the other hand, sampling errors are not as significant for oceanic variables as stronger covariabilities characterize the link between the MOC and oceanic temperature and salinity. Hence, very large ensembles are needed if one hopes to initialize the slow component of the MOC using daily DA without comprehensive oceanic observations. This requirement is curtailed when appropriate time-averaged observations are assimilated or if sufficient oceanic observations are available for assimilation.

Fig. 10

Minimum ensemble sizes required to reliably estimate correlations at the daily and yearly time scales between the MOC and every assimilated variable. X is the zonal wind, Y and Z are the eddy phase amplitudes and E _e is eddy energy. T and S represent temperature and salinity with subscripts 1, 2, 3 indicating high-latitude upper, low-latitude upper and deep ocean boxes respectively. Values are obtained from Monte Carlo experiments composed of 10,000 realizations of N random draws of model states. For each time scale and each variable, the minimum ensemble size corresponds to the smallest N value at which the sample correlation is equal to the corresponding value shown in Fig. 6 with a corresponding p-value <0.05 among at least 95 % of the 10,000 realizations. The ensemble size used in DA experiments is shown by the dashed gray line

This investigation has yielded a clearer characterization of some fundamental aspects of coupled atmosphere–ocean DA. Better defined guidelines for the development of initialization capabilities of next-generation climate prediction systems are emerging. It is acknowledged however that these were obtained on the basis of idealized, perfect-model experiments, which typically lead to optimistic results. Further work with comprehensive coupled models is required for assessing the generality of the results. A study aimed at this question, based on data from the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Taylor et al. 2012), will be reported elsewhere.

Appendix: Model equations

The atmospheric component of the coupled low-order model is composed of the following three ordinary differential equations (ODEs) (Roebber 1995; Lorenz 1984):

$$\frac{dX}{dt} = - (Y^2 + Z^2) - aX + aF,$$

(1a)

$$\frac{dY}{dt} = XY - bXZ - Y + G,$$

(1b)

$$\frac{dZ}{dt} = XZ + bXY - Z,$$

(1c)

where X is the zonal wind and Y and Z are the amplitudes of cosine and sine phases of large scale eddies, respectively. Interactions between the mean flow and eddies are represented by the amplification of eddies at the expense of the zonal flow intensity (first terms on the right-hand side of (1)), displacement of the eddies by the zonal flow (second terms in (1b) and (1c)), as well as mechanical damping (second term in (1a) and third terms in (1b) and (1c)). F and G terms represent diabatic heating contrasts between the low- and high-latitude ocean (F term) and seasonally varying zonal heating zonal difference between the land and ocean in the mid-latitudes (G term). Formulations for F and G are detailed later.

The evolution of ocean temperature (T, in Kelvin) and salinity (S, in practical salinity unit or psu) is calculated in the three boxes using the following set of ODEs (Roebber 1995; Birchfield et al. 1990):

$$V_1\frac{dT_1}{dt} = \frac{1}{2}q(T_2 - T_3) + K_T(T_{A1} - T_1) - K_{Z}(T_1 - T_3),$$

(2a)

$$V_2\frac{dT_2}{dt} = \frac{1}{2}q(T_3 - T_1) + K_T(T_{A2} - T_2) - K_{Z}(T_2 - T_3),$$

(2b)

$$V_3\frac{dT_3}{dt} = \frac{1}{2}q(T_1 - T_2) + K_{Z}(T_1 - T_3) + K_{Z}(T_2 - T_3),$$

(2c)

$$V_1\frac{dS_1}{dt} = \frac{1}{2}q(S_2 - S_3) - K_{Z}(S_1 - S_3) - Q_s,$$

(2d)

$$V_2\frac{dS_2}{dt} = \frac{1}{2}q(S_3 - S_1) - K_{Z}(S_2 - S_3) + Q_s,$$

(2e)

$$V_3\frac{dS_3}{dt} = \frac{1}{2}q(S_1 - S_2) + K_{Z}(S_1 - S_3) + K_{Z}(S_2 - S_3).$$

(2f)

Subscripts 1, 2 and 3 in (2) respectively denote the high-, low-latitude and deep ocean boxes of respective volumes V ₁, V ₂ and V ₃. K _Z is a vertical eddy diffusion coefficient between the upper and deep ocean while K _T is the coefficient of heat exchange between the ocean and an atmosphere with an air temperature T _A . Q _s is the volume averaged equivalent salt flux, representing the net surface evaporation over the subtropical ocean, poleward transport of water vapor and influx of freshwater in the high latitude ocean, e.g., a simplified representation of the hydrological cycle.

The MOC (q) has positive values for an overturning circulation characterized by a poleward upper ocean flow, sinking at high latitudes, an equatorward return flow in the deep ocean and upwelling at low latitudes. Following Stommel (1961) and Birchfield et al. (1990), a highly parameterized thermohaline circulation is formulated as a linear relation with the density difference between the two upper-ocean boxes. By assuming a linear relation to describe ocean water density with respect to temperature and salinity as in Stommel (1961) and Birchfield et al. (1990), a simple relation for the MOC is obtained:

$$q = \mu[\alpha(T_2 - T_1) - \beta(S_2 - S_1)],$$

(3)

where α and β are coefficients for thermal and haline expansion of seawater respectively, and μ is a proportionality constant to be determined.

Coupling between the atmosphere and ocean occurs through the meridional and zonal diabatic heating terms F and G in the atmospheric model and through the equivalent salt flux (Q _s ) term influencing ocean salinity. The diabatic heating terms are expressed as:

$$F = F_0 + F_1 \cos\omega t + F_2 \frac{(T_2 - T_1)}{T_0},$$

(4a)

$$G = G_0 + G_1 \cos\omega t + G_2 \frac{T_1}{T_0},$$

(4b)

where ω is the annual frequency and T ₀, F ₀, F ₁, F ₂, G ₀, G ₁ and G ₂ are constants. Both terms are composed of background values (F ₀ and G ₀), a seasonal cycle and a dependence on ocean state. The equivalent salt flux represents the sum of the influx of freshwater into the ocean from continental runoff (Q _runoff ) and contributions from the mean ($\overline{Q_{wv}}$) and transient eddy ($Q_{wv}^\prime$) components of the atmospheric water vapor meridional transport:

$$Q_s = Q_{runoff} + \overline{Q_{wv}} + Q_{wv}^\prime ,$$

(5)

A simple parameterization for Q _s is devised by assuming that runoff and mean transport terms are constant and by postulating that the eddy water vapor transport depends linearly on the eddy energy (Y ² + Z ²) as in Stone and Yao (1990):

$$Q_s = c_1 + c_2(Y^2 + Z^2),$$

(6)

where $c_1 = Q_{runoff} + \overline{Q_{wv}}$ and c ₂ are constants to be determined. Firstly, runoff is roughly equal to two-thirds of the total flux of freshwater in the northern part of the Atlantic (Broecker et al. 1990), secondly, the total atmospheric water vapor transport at 45°N is estimated to be between 0.3 Sv (1 Sv = 10⁶m³ s⁻¹) and 0.7 Sv (Wijffels et al. 1992; Zaucker et al. 1994), with roughly half of this transport taking place through the effect of the transient eddies (Chen 1985). Assuming a total atmospheric water vapor transport corresponding to the midpoint within the range of realistic estimates, values for c ₁ and c ₂ can be readily derived.

Finally, a linear relationship is used as an analog to the thermal wind balance, linking the large-scale meridional air temperature gradient and the strength of the zonal flow:

$$T_{A1} = T_{A2} - \gamma X,$$

(7)

where T _A2 is taken as constant (=25 ^o C).

The values of model constants are provided in Table 1. The chosen K _T value corresponds to a temperature restoring time scale of approximately 25–75 years for the upper ocean boxes, values consistent with those used by Titz et al. (2002) and the upper ocean temperature restoring time scale estimated by Liu (2012). The chosen value for K _Z corresponds to a diffusivity ≈ 0.25 cm² s⁻¹, an order of magnitude consistent with available estimates (Gregg 1987; Wüest et al. 1996). The associated diffusive time scale is O(250–850 years) for the upper-ocean boxes.

Table 1 Values of constants in the low-order model

The model is integrated numerically using a second-order Runge-Kutta scheme with a time step of 3 h.