## Abstract

A central issue for understanding past climates involves the use of sparse time-integrated data to recover the physical properties of the coupled climate system. This issue is explored in a simple model of the midlatitude climate system that has attributes consistent with the observed climate. A quasigeostrophic (QG) model thermally coupled to a slab ocean is used to approximate midlatitude coupled variability, and a variant of the ensemble Kalman filter is used to assimilate time-averaged observations. The dependence of reconstruction skill on coupling and thermal inertia is explored. Results from this model are compared with those for an even simpler two-variable linear stochastic model of midlatitude air–sea interaction, for which the assimilation problem can be solved semianalytically.

Results for the QG model show that skill decreases as the length of time over which observations are averaged increases in both the atmosphere and ocean when normalized against the time-averaged climatological variance. Skill in the ocean increases with slab depth, as expected from thermal inertia arguments, but skill in the atmosphere decreases. An explanation of this counterintuitive result derives from an analytical expression for the forecast error covariance in the two-variable stochastic model, which shows that the ratio of noise to total error increases with slab ocean depth. Essentially, noise becomes trapped in the atmosphere by a thermally stiffer ocean, which dominates the decrease in initial condition error owing to improved skill in the ocean.

Increasing coupling strength in the QG model yields higher skill in the atmosphere and lower skill in the ocean, as the atmosphere accesses the longer ocean memory and the ocean accesses more atmospheric high-frequency “noise.” The two-variable stochastic model fails to capture this effect, showing decreasing skill in both the atmosphere and ocean for increased coupling strength, due to an increase in the ratio of noise to the forecast error variance. Implications for the potential for data assimilation to improve climate reconstructions are discussed.

## 1. Introduction

Determining the state of the climate prior to the twentieth century necessarily involves paleoclimate proxy data. These data may be viewed as a relatively sparse network of noisy observations that, in many cases, are indirect records of the climate mediated by biological or chemical processes. As a result, reconstructions of past climates have relied on statistical relationships between these proxies and large-scale patterns of atmospheric variability (e.g., Mann et al. 1999). One path to potentially improving upon these statistical reconstructions involves the introduction of a dynamical model to provide independent estimates of the observations. This approach, which we refer to as dynamical climate reconstruction, links the paleoclimate reconstruction problem to weather forecasting and modern reanalysis through the theory of state estimation. Practical realizations of dynamical climate reconstruction face numerous challenges, only one of which is taken up here. We address the problem of estimating the state of coupled atmosphere–ocean climate models subject to observations involving long averaging times typical of paleo-proxy data.

One key issue for dynamical climate reconstruction concerns the predictability time horizon of the model—that is, the time for which initial condition errors approach the climatological distribution. When the (model) predicted observations are indistinguishable from a random draw from the climatological distribution, the main value of using the model is lost. The predictability time scale needed for reconstruction is set by the proxy data, which span a wide range of time scales. “High frequency” paleoclimate records provide information with at best seasonal, but more typically annual, resolution (e.g., Jones et al. 2009). Since atmospheric predictability is on the order of two weeks, predictability on the time scales of these paleoclimate data requires dynamics coupled to “slower” parts of the climate system. For example, sea surface temperatures have been shown to persist on interannual time scales in observations (Frankignoul and Hasselmann 1977; Davis 1976; see Deser et al. 2003 for a modern update) and models (Saravanan et al. 2000). Interest in this persistence has led to research on variability and predictability in highly idealized models for air–sea interaction (e.g., Barsugli and Battisti 1998; Saravanan and McWilliams 1998; Scott 2003). On decadal time scales, there is evidence that the overturning circulation in the deep ocean couples to the atmosphere (e.g., Meehl et al. 2009), which appears partly responsible for predictable atmospheric signals beyond the seasonal time scale (Latif 1998; Griffies and Bryan 1997; Grötzner et al. 1999; Boer 2000; Saravanan et al. 2000; Collins 2002; Pohlmann et al. 2004; Latif et al. 2006; Boer and Lambert 2008; Koenigk and Mikolajewicz 2009).

Paleoclimate proxy data for air temperature and precipitation in midlatitudes represent climate states integrated over all of the weather during the course of a season or annual cycle. The characteristics of weather events, and the frequency with which they visit particular locations, depend on larger-scale, slower varying properties of the atmosphere and ocean. For example, extratropical cyclones organize in storm tracks, which in turn follow the location of the jet stream. Thus, the reconstruction problem is one where the details of individual high-frequency features responsible for the proxy data are unknown, but the time average of these features provides information on slower components of the system.

To address this problem, we use two idealized models of a thermally coupled midlatitude atmosphere and ocean mixed layer that represent both the relevant fast and slow dynamics. The zero-dimensional energy balance model (referred to as the BB model) was designed as an analog to the midlatitude atmosphere–ocean system by Barsugli and Battisti (1998). The quasigeostrophic (QG) model, developed by Hakim (2000), explicitly calculates the dominant atmospheric dynamics, and the slab ocean parameters were chosen to correspond to those in the BB model. Although the models are highly idealized, they capture important dynamical mechanisms for coupled climate variability relevant to paleoproxy data.

To exploit seasonal or longer predictive time scales and dynamically reconstruct past climate states, an estimation method is needed to identify accurate initial conditions. Here we will employ an ensemble Kalman filter (EnKF) technique applied to the above models of coupled atmosphere–ocean interaction to address two hypotheses: background climate reconstruction errors become smaller when (i) the thermal heat capacity of the ocean increases and (ii) the coupling between the ocean and atmosphere strengthens. The EnKF is particularly attractive for addressing these hypotheses because it relaxes the assumption of stationary statistics made in purely statistical reconstructions, such as Mann et al. (1999) and subsequent efforts along the same vein. Therefore, information about the slow component of the system (the ocean) can affect the atmosphere through both the forecast and the assimilation. Idealized models provide clean control of the processes of interest, inexpensive repeatable solutions, and results that can be tested with more complicated models. We address initial condition error and employ perfect-model experiments that, by definition, exclude model error.

The layout of the paper is as follows. Section 2 introduces the two models and the assimilation systems built around them. The QG model is described in section 2a, and the BB model is described in section 2b. Results are given in section 3, with those determined numerically for the QG model in section 3a and those obtained semianalytically for the BB model in section 3b. Portions of the error for the BB model due to the initial conditions and accumulation of noise are described in section 3c. Conclusions are provided in section 4.

## 2. Method

### a. Quasigeostrophic atmosphere model coupled to a slab ocean

The first model that we employ resolves the dynamics of midlatitude weather systems coupled to a mixed layer ocean. To good approximation, the QG equations capture the dynamics of the midlatitude storm tracks. Moreover, since potential vorticity (PV) is well mixed in the troposphere, the dynamics may be further approximated by a uniform layer of PV bounded by rigid boundaries representing the surface and tropopause, which reduces the calculation in three dimensions to a two-dimensional boundary value problem (Hakim 2000). This model, when coupled to a slab ocean model, is referred to as the QG model.

The atmosphere component of the QG model applies to an *f* plane that is thermally relaxed to a specified jet (Hoskins and West 1979) on a time scale of 8.6 days. The model is periodic in *y* and *x* (latitude and longitude), has an Ekman layer at the bottom boundary, and solves exactly for the full three-dimensional fields of temperature and wind given the surface potential temperature. The model was originally described in Hakim (2000), who provides additional details on the solution methods. An uncoupled version of the atmosphere model was used by Huntley and Hakim (2009) to explore the time-averaged data assimilation algorithm for paleoclimate reconstruction proposed by Dirren and Hakim (2005). Huntley and Hakim found a predictability horizon beyond which the model no longer provided information distinct from climatology (i.e., a fixed error covariance matrix). However, this model lacked coupling to components of the climate system that have slower time scales, which we address here by thermally coupling the QG atmosphere to a slab ocean.

The slab ocean does not move and only exchanges heat with the atmosphere. To preserve the uniform-PV assumption of the atmospheric component of the QG model, heat is exchanged uniformly at the surface and tropopause (i.e., through the entire depth of the atmosphere), which is essentially an equivalent barotropic forcing of the atmosphere. This can be thought of as each ocean grid box exchanging heat instantaneously with the entire depth of the atmosphere directly above it. The heat flux is determined by the temperature anomaly difference between the atmosphere’s surface and the slab ocean, modulated by a coupling constant. More realistic coupling coefficients would also involve the wind speed magnitude or its square, but these were sacrificed here for simplicity and to provide a link to the analytical model described below. Control values for the constants in the model are listed in Table 1.

This simple model provides a framework for addressing the two hypotheses in terms of two parameters: the coupling constant and the depth of the slab ocean. While this coupled model provides a test bed to study how the atmosphere, modulated by the slow thermal inertia of a slab ocean, may be reconstructed on long time scales, it should be clear that this model does not include nonlinear aspects of heat flux exchange (e.g., due to wind speed), ocean dynamics, changes in the mixed layer and its heat capacity, the seasonal cycle, or any precipitation–evaporation-related processes. However, it does capture the dynamics of heat exchange between a rapidly evolving atmosphere and a mixed layer with greater thermal inertia, which is fundamental in midlatitudes. We proceed now to describe the assimilation system used to reconstruct the model state from noisy observations of time averages of the atmosphere and ocean temperature.

The EnKF data assimilation system for the QG model is implemented in the Jamaica release of the Data Assimilation Research Testbed (DART) (Anderson et al. 2009), a community-oriented ensemble data assimilation system. The ensemble adjustment Kalman filter (Anderson 2001) is used without covariance inflation and, to deal with sampling error, with covariance “localization” using the Gaspari–Cohn function (Gaspari and Cohn 1999) with a half-width of half the domain. Observations are drawn from a “truth” simulation at every other grid point (for a total of 32 in *x* and 16 in *y*), each assumed to have an error of 1/10 the domain-mean climatological variance for the averaging time of the observations. The domain is 28 000 km in *x* and about 11 000 km in *y*; about four midlatitude storms are typically fit in the domain in *x*. At each observing location, temperature observations are created for the slab ocean, the surface of the atmosphere, and the tropopause. The purpose of this dense network of observations is to create a situation in which the density of observations is not a limiting factor for reconstruction skill—an issue that was thoroughly explored by Huntley and Hakim (2009). Ensembles of 48 members are initialized with random draws of instantaneous states from a 2500-day integration of the model. Each experiment consists of one ensemble integration of 100 assimilation cycles without assimilation (the control case) and another integration with assimilation. For all experiments greater than an averaging time of one day, the background ensemble spread and error of the ensemble mean are of the same magnitude; that is, the ensemble is well calibrated. For an averaging time of one day, this ensemble-mean error is too large for the atmosphere, which implies that the ensemble is not well calibrated in this case; however, our interest lies with the long, not short, time limit.

To investigate the accuracy of the reconstruction method, we examine the background ensemble error variance, which is the error variance at the end of a forecast cycle immediately before new observations are incorporated in the assimilation cycle. For a control, we use an ensemble that is integrated without assimilating observations and define “skill” to be the error of the background ensemble normalized by the error of the control ensemble. Specifically, error is defined by

the error variance is given by

and skill by

where is the time-averaged potential temperature at a point in space for one assimilation cycle, is the analogous quantity from the truth simulation, is the ensemble mean value, and *ɛ*(*x*, *t*) is the ensemble-mean error; *R* denotes the error variance calculated from each of *M* times for each grid point, and *R _{c}* is analogous quantity for the control. We will present time- and space-averaged values of

*s*, taken over all assimilation cycles and the center half of grid points in the

*y*direction for a particular experiment. If skill is greater than zero, then useful information remains in the background from the previously assimilated observations. We note that

*R*is analogous to the ensemble-mean error variance, a commonly used metric in data assimilation studies in which background forecasts apply to times much shorter than the limit of predictability. Normalization by the control

*R*is useful for both identifying when the forecast has lost skill and for accounting for the decrease in variance that occurs as the averaging time decreases. Longer averaging times leave more variance in the unobserved high frequencies, which are not included in the reconstruction. Without this normalization, the change in variance with averaging time dominates the signal.

_{c}### b. Stochastic energy balance model

A simpler zero-dimensional model of the atmosphere–ocean system provides insights into the results of the more complex system described above. The model was first introduced by Barsugli and Battisti (1998) and will be referred to as the BB model. The model utilizes the large time-scale separation in the variance of the atmosphere and ocean to parameterize atmospheric nonlinearity as a stochastic process. Stochastic climate models, such as the BB model, have a long history of providing a null hypothesis for problems in climate dynamics, dating from Hasselmann (1976) who used a stochastic climate model to explore how persistent, large-scale sea surface temperature anomalies could arise in midlatitude oceans. The stochastic noise term in the BB model corresponds to the weather variations in the QG model, which can be accurately described by random noise for sufficiently long sampling periods (Hasselmann 1976; Frankignoul and Hasselmann 1977; reviewed by Kushnir et al. 2002).

The BB model has thermal exchange between an atmosphere and ocean, damping, and random forcing applied to the atmosphere:

and

where *T _{a}* is the air temperature,

*T*is the mixed layer temperature,

_{o}*d*is the ratio of heat capacities between the atmosphere and mixed layer,

*c*is the coupling strength,

*τ*is the radiative damping time scale,

*b*is a constant controlling the deterministic portion of the atmospheric forcing, and

*N*is the white noise forcing. Equations (4) and (5) can be written compactly as

where

Here, is a vector of the temperatures and **N** represents a Gaussian white noise process with amplitude *N*. The constants *b*, *c*, *d*, *τ*, and *N* are chosen to be representative of climate-scale atmosphere–ocean interaction in the midlatitudes (see Table 1 and Barsugli and Battisti 1998).

In the absence of noise forcing (when *N* = 0), this system supports two normal-mode solutions that are damped in time, which provide a convenient basis for understanding the system dynamics with noise. The first mode damps quickly, with an *e*-folding time scale of 3 days for the parameters in Table 1, and consists of a mainly atmospheric temperature anomaly with a weaker ocean temperature anomaly of opposite sign, as noted by Schopf (1985) for a similar deterministic model. The second mode has an *e*-folding time scale of 89 days and consists of ocean and atmosphere temperature anomalies of the same sign, with larger magnitude in the ocean by roughly a factor of 2. This second mode contains the “memory” in the system, which arises because the ocean anomaly is large and of the same sign as the atmosphere so that heat exchange is minimized between the two, resulting in a slow decay.

Because the BB model contains only two variables, the data assimilation reconstruction problem can be solved semianalytically. Since we are interested in errors, we need only solve for the error covariance matrix. The expected error covariance over infinite experiments, derived in the appendix, is

where *ɛ** is the instantaneous error of an ensemble member, is the time-averaged error of an ensemble member, *L* is the averaging time, is the matrix of model parameters, and is a 2 × 2 matrix with only one nonzero element (see the appendix). The first term in (7) represents the time evolution of the initial error covariance, which for this damped system leads to a decay of error covariance. The second term involves an integral over the noise, which gives exponential growth of error covariance. We will refer to this second effect as the “accumulation” of noise. The skill metric for assimilation with the BB system is analogous to the metric used for the QG system.

We note that the semianalytic solution of (7), described in the appendix, is a solution to the Kalman filter, not the EnKF, and is therefore not subject to sampling error. Furthermore, according to Fitzgerald (1971), systems for which all system eigenvectors are observed and driven by noise, such as the BB system, have solutions that converge exponentially in time; therefore, filter divergence is impossible.

## 3. Results

Before presenting results, we review the hypotheses outlined in section 1. The first is that errors in climate reconstruction will be smaller when the thermal heat capacity of the ocean is larger. The crux of this idea is that increasing the thermal heat capacity lengthens the persistence time scale of anomalies in the ocean, as discussed in section 2b, which increases the information carried forward by the model to the time of the analysis. The second hypothesis is that stronger coupling between the atmosphere and ocean should also reduce errors in climate reconstruction, especially for the atmosphere. Again, since the ocean is the source of memory, stronger coupling should imprint that memory onto the atmosphere, enhancing predictability. For the ocean, the role of coupling is less clear. Stronger coupling imparts more high-frequency atmospheric “noise” into the ocean, reducing predictability, but some of the ocean’s persistence is imparted to the atmosphere, which may mitigate this effect.

The data assimilation algorithm used here is designed for dealing with paleoclimate proxy data, which integrates the climate signal over a period of time. Because interpreting the effect of different averaging times is not as straightforward as for instantaneous observations, we consider first the role of averaging time on the error. For instantaneous data assimilation, analysis errors should decrease as observation frequency increases since, all other things being equal, more frequent observations provide more information about the instantaneous state of the system if the model forecast is well within the predictability time horizon. For time-averaged data assimilation, the climatological variance of the states also decreases with increasing averaging time. Therefore, we scale the observation error with climatological variance for the time scale of the observations to maintain constant error relative to climatology. We expect that skill will decrease with increasing averaging time.

### a. Results for the coupled QG model

First, we consider the dependence of background skill on slab ocean depth and averaging time, holding the coupling strength fixed at *c* = 0.25 (Fig. 1). Each pixel represents the skill of the background time-averaged state for a data assimilation experiment composed of 100 assimilation cycles, area-averaged over the meridional center of the domain (the region of basic-state meridional temperature gradient). We show skill for the surface temperature in the atmosphere; skill at the tropopause is nearly identical. The slab ocean depth *d* is represented by the ratio of the slab ocean heat capacity to the atmospheric heat capacity. Results show that for both atmosphere and ocean, and for all depths of the slab ocean, skill decreases as the averaging time increases, consistent with prior expectation. There is a sharp drop of skill between 5 and 10 days, which we suspect corresponds to the model damping time scale of 8.6 days.

As slab ocean depth increases, results generally show increasing skill in the ocean, as expected. In contrast, the pattern observed in the atmosphere is generally decreasing skill with increasing slab ocean depth, despite an ocean temperature reconstruction that is more skillful. We will revisit these conflicting results after reviewing those for the BB model.

Figure 2 shows the analogous experiment for a fixed nondimensional slab ocean depth ratio of *d* = 4 as a function of the coupling parameter. For the atmosphere, skill increases slightly with increasing coupling coefficient for averaging times from 1 to 10 days. For averaging times greater than 1 day, skill in the ocean decreases with increasing coupling coefficient. These results are consistent with expectations, in the sense that the atmosphere is a source of noise for the ocean and the ocean is a source of thermal memory for the atmosphere.

The skill for 1-day averaging behaves in the opposite way to the broad pattern in the ocean both as *c* increases and as *d* increases. This period is short enough so that most ocean grid cells have encountered only one weather system, and the dense observing network allows the ensemble to resolve these individual systems well. So, the effect of the ocean’s thermal inertia is different when there has been one weather system compared to many. We include this case because it approximates the familiar weather data assimilation problem; however, the relevant case for paleoclimate is the longer averaging times over multiple weather systems. Moreover, the scenario in the QG model for 1-day averaging is not valid for the BB model because the stochastic approximation does not resolve this scenario.

In both Figs. 1 and 2 there is a striking difference between skill in the ocean and atmosphere. This is due to two factors. The first is simply the longer thermal inertia of the ocean. The process time scales in the ocean are each a factor of *d* longer than in the atmosphere, so persistence and predictability in the ocean should be longer. But even for *d* = 1, when the ocean and atmosphere have the same heat capacity, the ocean has longer predictability than the atmosphere. For *d* = 1 (first column of Fig. 1), the ocean has greater than 50% more skill than the atmosphere for averaging times of 5 and 10 days. One possible explanation is that advection occurs in the atmosphere but not in the ocean. This additional process leads to faster dynamics in the atmosphere that degrade predictability. The ocean essentially averages over this fast-moving variability, leading to higher skill at longer time averages compared to the atmosphere.

### b. Results for the BB (stochastic energy balance) model

Before reporting the results from the BB model, we will compare its equilibrium statistical properties with those of the QG model, which it is meant to approximate. Cospectral analysis of the surface air and slab ocean temperatures is shown in Fig. 3, calculated from the last 8000 days of an 8500-day model integration for the QG model and analytically for the BB model. For the QG model, the results are shown by plotting the distribution of values taken from grid points in the center half of the domain in the *y* direction. Ten-day blocks of unnormalized time series are averaged together, and then cospectral analysis is carried out with a Hanning window of length 800 days using the Matlab spectrum function. Power spectra (*P _{a}* and

*P*), cospectra (

_{o}*F*), and phase (

_{ao}*φ*) from the BB model are calculated from Eqs. (14) and (15) in Barsugli and Battisti (1998):

and

where *σ* is the frequency, the asterisk indicates the complex conjugate, and *α* = *bc*. Just one line is shown for the BB model, calculated with parameters *b* = 0.5 and *N* = 5, chosen by trial and error to correspond well to the QG model spectra; other parameters are the same as those for the QG model.

As with the BB model, the spectra of the QG ocean and atmosphere increase toward lower frequencies. The cospectral power increases as frequency decreases in both BB and QG models, and there is a low-frequency peak in the quadrature spectrum of both models. The atmosphere and ocean show high coherence at frequencies below 2π/100 days^{−1}, which then degrades owing to nonlinear effects at high frequencies. The phase relationship of the QG model follows very closely the BB phase at frequencies for which it is well defined. The close correspondence of the spectral characteristics of the two models lends confidence that the BB model provides a good approximation to the QG model for further analysis in a simpler framework.

The BB model data assimilation system discussed in section 2b provides the expected error statistics over an infinite number of experiments and ensemble members. The numerical evaluation of Eq. (7), described more completely in the appendix, begins with zero error and cycles until equilibrium is reached. Observations of both the air and slab ocean temperatures are assimilated with observational error taken to be 10% of the equilibrated error variance in the absence of assimilation. The evolution of the solution shows a rapid convergence to equilibrium (Fig. 4).

Skill as a function of averaging time and slab ocean depth for fixed coupling strength is shown in Fig. 5 (cf. with QG model results in Fig. 1). As in the QG experiments, skill decreases with increasing averaging time in both the atmosphere and ocean. Skill in the atmosphere decreases with increasing depth, and skill in the ocean increases with increasing depth, which is also consistent with results from the QG model. A notable difference between the QG and BB model results is lower skill in the atmosphere of the BB model for all slab ocean depths and averaging times.

Skill as a function of averaging time and coupling strength, for fixed slab ocean depth, is shown in Fig. 6 (cf. with QG model results in Fig. 2). Skill for air temperature decreases as coupling increases, which is opposite to the behavior found in QG experiments. Skill for ocean temperature decreases as coupling increases, which is consistent with results for the QG model. In summary, compared to the QG results, the BB system shows a similar response to varying slab ocean depth, the same response to coupling for the ocean, but an opposite response to coupling for the atmosphere.

### c. Decomposing error: Accumulation of noise and initial condition error

Further insight into the dependence of reconstruction skill on coupling strength and slab ocean depth may be derived from the error evolution equation (7) for the BB system. Recall that BB is a linear stochastic model, where nonlinearity has been parameterized by a noise process, and the system is damped in the absence of the noise. As discussed in section 2b, the error covariance evolution depends on two factors: the initial error at the start of an assimilation cycle and the accumulation of stochastic noise. Although the initial error depends on the noise error, observation error, and data assimilation system, the decomposition discussed here provides a useful perspective on the relative importance of initial error compared to accumulating noise in the forecast. Figure 7 shows the total error broken into these two terms as a function of slab ocean depth. Assimilation is performed as for the previous assimilation experiments; the converged solution for the total error and its two component parts are plotted separately for the atmosphere and ocean. For the atmosphere, the fraction of the error variance attributable to initial error decreases slightly with increasing depth, while the part proportional to noise accumulation increases. The increase due to noise accumulation dominates so that the total error increases with depth. For the ocean, the fraction of error due to noise decreases strongly with increasing slab ocean depth, whereas the fraction due to initial error peaks at small depth and decreases gradually with depth so that it dominates the total error for large depth.

Results for sensitivity to the coupling strength are shown in Fig. 8. For the atmosphere, noise error dominates over initial condition error, as is the case for the varying slab ocean depth experiments. For the ocean, initial condition error dominates, as for the depth experiments. A notable difference from the depth experiments is that the ocean noise error increases slowly with coupling, whereas it decreases sharply with increasing depth. A physical interpretation of these sensitivity results is that deeper slab oceans are less sensitive to atmospheric noise due to increased thermal inertia, whereas stronger coupling increases the amount of noise introduced from the atmosphere—but this latter effect is weak compared to initial condition error.

## 4. Conclusions

The problem of reconstructing past climates from time-averaged observations has been explored within an idealized framework for midlatitude atmosphere–ocean interaction. A QG model coupled to a slab ocean was used to explore the hypotheses that skill in reconstructions will increase when (i) the depth of the slab ocean increases and (ii) the strength of the coupling increases. Results show that for increasing slab ocean depth, ocean skill increases as expected, but skill in the atmosphere does not. Skill increases slightly in the atmosphere for stronger coupling, but it decreases in the ocean. For both atmosphere and ocean, skill decreases with increasing observation averaging time.

The average freely evolving behavior of the QG system was shown to correspond well to that of a simple two-variable linear stochastic model for midlatitude air–sea interaction (the BB system). Results for assimilation experiments for this system are similar to the QG system for sensitivity to slab ocean depth, but not for coupling. An analytical expression for the BB error covariance forecast evolution elucidates these results. The atmosphere loses state-dependent skill faster than the ocean because the stochastic component of error accumulates more quickly than in the ocean. Moreover, the fraction of the total atmospheric forecast error variance due to stochastic noise increases with both slab ocean depth and coupling strength; that is, in the BB system, noise controls the solution more than initial condition error. A physical interpretation of these results is that the ocean becomes thermally stiffer for increasing depth and coupling so that noise “accumulates” in the atmosphere. For changes in coupling, responses are not consistent between the models. Increasing the coupling coefficient yields more spectral power in slab ocean temperature at all frequencies for both QG and BB models; however, for surface air temperature, power increases in the BB model and decreases in the QG model (not shown). Therefore, one explanation for the different response to changes in coupling coefficient between the models is that the fixed-amplitude white noise approximation of the BB model does not properly represent the forcing in the QG model as the coupling coefficient changes (see Fig. 6.3 of Barsugli 1995). We speculate that, by increasing the coupling coefficient, the ocean acts to damp the growth of air temperature anomalies in the QG model, which results in increased skill in the atmosphere.

Our objective in this work is to illuminate the basic impact of assimilating proxy information for midlatitude climate reconstruction using time-averaged data. The models used here are assumed perfect with respect to truth, which generates the observations, so that our results may be considered an upper bound on reconstruction accuracy. Adding structural errors to the “true” model integration would be an ad hoc exercise in these idealized models, and we believe that the basic reconstruction limitations involving atmosphere–ocean interaction discussed here are insensitive to this issue.

What do these results imply for climate reconstructions from data assimilation? The simple mixed layer models considered here lack skill beyond a few months, which implies that to the extent that this model is applicable to the real world, dynamical reconstruction may not be more skillful in midlatitudes than statistical reconstruction on annual time scales. However, data assimilation can still improve upon statistical reconstruction when climate covariances are a strong function of mean state, such as the Last Glacial Maximum. Paleoclimate data assimilation may also improve upon statistical reconstruction when the model has skill on longer time scales than the observed proxy.

## Acknowledgments

This research was sponsored by the National Science Foundation through Grant 0902500, awarded to the University of Washington. The first author was also by supported by fellowships from the National Defense Science and Engineering Graduate fellowship program and the Achievement Rewards for College Scientists. The first author thanks Nils Napp for insightful conversations on stochastic calculus. We thank two anonymous reviewers for constructive and thoughtful comments that greatly improved the manuscript.

### APPENDIX

#### Development and Evaluation of Error Covariance

Here we develop analytic expressions for error covariance evolution in the BB model, for both instantaneous, *C*_{total}, and time-averaged, *C*_{average}, cases. Note that the adjustment to the error covariance during the assimilation step does not require observation values, which is a convenient property of the Kalman filter. The calculation relies on stochastic calculus, an introduction to which is available in textbooks, Øksendal (2003), and Gardiner (2004), as well as an article directed toward atmospheric scientists by Penland (2003).

##### a. Instantaneous error covariance

Rearranging the system of equations in the typical style of stochastic calculus and omitting vector symbols for clarity, (6) becomes

where and *B _{t}* is a vector of white noise processes. Using the integrating factor yields the solution

We define the error of an ensemble member (or the ensemble mean) as the Euclidean distance from the true state, *T _{T}*(

*t*):

*ɛ*(

*t*) =

*T*(

*t*) −

*T*(

_{T}*t*). Defining the expectation over an infinite ensemble 〈·〉 (where ensemble members are denoted by the random variable

*B*), the ensemble error covariance is

_{s}Deviations from the ensemble mean are denoted by (·)*, so .

Noting that the expectation over a random process vanishes , the error of an ensemble member at time *t* is

Applying the expectation operators over ensembles to the squared error deviation *ɛ***ɛ**^{T}, all cross terms with only one term integrated over *dB _{s}* vanish. The error covariance is

Distributing the last transpose, using Itô’s isometry, and rearranging, the expected error covariance is

The two terms determine the error covariance evolution. The first is growth or decay of the initial error covariance, depending on the eigenvalues of . The second is the accumulation of noise, which increases the error covariance even when has only negative eigenvalues.

##### b. Error covariance of time-averaged states

The error of the time-averaged state is equal to the time-averaged error. Therefore, the error of a time-averaged state is given by

where the overbar denotes a time average.

Squaring the quantity above and applying the expectation operator, the cross terms cancel as they did for the instantaneous error covariance, leaving

Switching the order of integration and using the Itô isometry again yields

The error covariance of time-averaged states evolves in an analogous fashion to the instantaneous states, depending on initial error and the accumulation of stochastic noise.

##### c. Evaluating the integrals

Since is linear and constant for the BB model, the evaluation of the covariance integrals is straightforward. The calculation will be shown only once since all other terms can be solved in a similar manner.

The first term in the instantaneous error covariance Eq. (A6) can be evaluated directly. We define so that the second term can be written as

where the elements of are *m _{ii}* and depend on

*t*and

*s*. We calculate the elements of by doing an eigenvalue decomposition of since, if , then . Then distributing and rearranging the elements of and yields

The integral of each element of (A11) has two exponential terms: and , with all other quantities constant. This leaves the sum of many integrals of the form for *a* = *d*_{1} and *d*_{2}, which have solutions

Evaluating the error covariance of time-averaged states is similar to evaluating the instantaneous error covariances. The first term of (A9) is rewritten by distributing transposes and integrating [note that is the identity matrix and (·)^{−T} indicates an inverse and a transpose]

For the second term, we return to and its elements:

As in the instantaneous case, each of the elements of this matrix is the product of integrals over exponential terms. Each of the exponential integral parts can be broken down into

The error covariance averaged over time *L* is calculated from an initial error covariance matrix. Integration over many assimilation cycles brings the system to convergence. An example integration of the system is shown in Fig. 4.