On the Movements of the North Atlantic Subpolar Front in the Preinstrumental PastAbstract
Three sediment records of sea surface temperature (SST) are analyzed that originate from distant locations in the North Atlantic, have centennial-to-multicentennial resolution, are based on the same reconstruction method and chronological assumptions, and span the past 15 000 yr. Using recursive least squares techniques, an estimate of the time-dependent North Atlantic SST field over the last 15 kyr is sought that is consistent with both the SST records and a surface ocean circulation model, given estimates of their respective error (co)variances. Under the authors’ assumptions about data and model errors, it is found that the 10°C mixed layer isotherm, which approximately traces the modern Subpolar Front, would have moved by ~15° of latitude southward (northward) in the eastern North Atlantic at the onset (termination) of the Younger Dryas cold interval (YD), a result significant at the level of two standard deviations in the isotherm position. In contrast, meridional movements of the isotherm in the Newfoundland basin are estimated to be small and not significant. Thus, the isotherm would have pivoted twice around a region southeast of the Grand Banks, with a southwest–northeast orientation during the warm intervals of the Bølling–Allerød and the Holocene and a more zonal orientation and southerly position during the cold interval of the YD. This study provides an assessment of the significance of similar previous inferences and illustrates the potential of recursive least squares in paleoceanography.
The surface circulation in the North Atlantic Ocean appears to be dominated by three connected currents that transport tropical waters to subpolar latitudes (e.g., Schmitz and McCartney 1993; Rossby 1996; Fratantoni 2001; Reverdin et al. 2003) (Fig. 1). The first is the Gulf Stream, which flows to the north and east after it separates from the North American coast at Cape Hatteras. The second is the North Atlantic Current (NAC), which constitutes a branch of the Gulf Stream after it splits near 50°W and which flows northward along the continental rise east of the Grand Banks. The third is the Subpolar Front, which is the northern and eastern extension of the NAC after it turns anticyclonically in the “northwest corner” (notice that many authors refer to the continuing flow to the east as the NAC or North Atlantic Drift as well). The Subpolar Front is not a sharply defined front in the sense of the Gulf Stream after it detaches from the coast, but a relatively wide region that separates the subtropical gyre from the subpolar gyre and where the main thermocline shoals to the surface (Rossby 1996).
|
The North Atlantic Current–Subpolar Front (NAC–SPF) system transports mass and heat far north into the North Atlantic (e.g., Pérez-Brunius et al. 2004). These waters become colder and denser as they flow northward and eastward, both by mixing with colder surrounding waters and by heat release to the atmosphere. The heat release to the atmosphere is thought to be largely responsible for the mild winters of western and northern Europe (e.g., Krauss 1986), although this notion has not gone unchallenged (Seager et al. 2002). It has been speculated that, should the poleward ocean heat transport decrease or turn zonal at midlatitudes, the mean temperature of the North Atlantic and European landmasses would drop precipitously (Rossby 1996).
Analyses of hydrographic, current meter, float, and satellite observations showed that the NAC–SPF system is characterized by significant variability on monthly and interannual time scales (e.g., Colin de Verdière et al. 1989; Sy et al. 1992; White and Heywood 1995; Belkin and Levitus 1996; Bower and von Appen 2008; Read et al. 2010; Walter and Mertens 2013). For example, Sy et al. (1992) relied on CTD data obtained in 1981–84, XBT profiles, and long-term current meter moorings to study the NAC above and east of the Mid-Atlantic Ridge (MAR). The NAC was found to be composed of clearly defined branches (jets). Whereas the northernmost branch appeared to be topographically fixed at the Charlie–Gibbs Fracture Zone (cutting the MAR at 52°–53°N), the number and intensity of the remaining branches, as well as their path over the MAR, were subject to intense variability. Bottom topography, wind stress curl, and mesoscale eddies appear to all be important for the NAC–SPF system and for the meridional shifts of the SPF (e.g., White and Heywood 1995; Bower and von Appen 2008; Read et al. 2010).
Inferences about the position of large-scale fronts in the North Atlantic have also been drawn for times in the preinstrumental past, in particular for the Last Glacial Maximum (LGM), a time interval around 20 kyr before present (BP) (Mix et al. 2001; Clark et al. 2009). The CLIMAP (1976) project published a global map of summer (August) sea surface temperature (SST) at the LGM, which was prepared from planktonic foraminiferal counts in a number of sediment cores. The authors estimated that the Gulf Stream shifted slightly southward off the northeastern North American coast and swept almost directly across the basin to the Iberian Peninsula. The reconstructed steep thermal gradient at 42°N in the Atlantic was interpreted as marking the southern edge of “polar” waters (the “polar front”). The SST distribution, the position of the Gulf Stream, and the position of the polar front in the glacial North Atlantic have been discussed in subsequent studies (e.g., Pflaumann et al. 2003; Matsumoto and Lynch-Stieglitz 2003; MARGO 2009; Eynaud et al. 2009). In a recent work, Barker et al. (2015, their Fig. 1d) located the front to south of about 50°N, based on the proportion of the polar species Neogloboquadrina pachyderma (left coiling) in the planktonic foraminiferal assemblage measured in several sediment cores.
The movements of the North Atlantic polar front during the last deglaciation have been estimated from a synthesis of sediment records (Ruddiman and McIntyre 1981). These authors inferred that the front retreated from its glacial position to the northwest at ~13 14C kyr BP (roughly coincident with the warm episode of the Bølling–Allerød) but then readvanced back to near that position at ~11 14C kyr BP [coincident with the onset of the cold episode of the Younger Dryas (YD)]. At ~10 14C kyr BP (approximately, the end of the YD), they reported that the front would have moved again to the northwest to a position along the mouth of the Labrador Sea and close to the east Greenland coast. The idea of a pivoting polar front during abrupt climate changes has been echoed in subsequent publications (e.g., Zahn 1994; Barker et al. 2015).
The time scales characterizing the deglacial movements of the polar front have been discussed by some authors (e.g., Bard et al. 1987; Lehman and Keigwin 1992). For example, Bard et al. (1987) reported deglacial reconstructions of winter (February) SST from two cores in the eastern North Atlantic: core CH73–139C from the Rockall Plateau (54°38′N, 16°21′W) and core SU81–18 from the Iberian margin (37°46′N, 10°11′W) (Fig. 1). They concluded that the front moved between these two locations in less than 1000 yr for the earliest retreat (dated at 12.0–12.5 14C kyr BP) and in less than 400 yr at the onset and termination of the YD.
Studies on the deglacial movements of North Atlantic fronts, however, suffer from several limitations. First, the reconstructed SSTs used to infer these movements have estimated errors of O(1°C) on average (e.g., Waelbroeck et al. 1998). This value is smaller but not negligible compared to the estimated changes, from LGM to today, in North Atlantic SSTs (annual mean, winter, or summer values) derived from multiple proxies (MARGO 2009). Second, different sediment indicators can result in noticeably different SST estimates (e.g., Chapman et al. 1996; Marchal et al. 2002; MARGO 2009), particularly at high latitudes where planktonic foraminiferal counts may suggest low temperatures relative to other proxy estimates (e.g., Bard 2001; Waelbroeck et al. 2014). Third, although the last deglacial period is within the range of 14C and other dating techniques (e.g., Mix et al. 2001), the chronology of deep-sea sediment records is difficult to establish, especially for cores originating from high latitudes (e.g., Waelbroeck et al. 2001). Fourth, despite the generally high deposition rates in the North Atlantic compared to other oceanic basins, many sediment records do not have the temporal resolution that would be needed to properly document the front migrations (e.g., Bard et al. 1987; Lehman and Keigwin 1992). Finally, inferences about past frontal movements are generally made in the absence of an ocean circulation model, implying that the dynamical consistency of the reported movements has not been demonstrated.
In this paper, inverse methods are applied in an effort to infer the movements of the North Atlantic SPF over the past 14.5 kyr (i.e., from the beginning of the Bølling to the present). Frontal movements are estimated by fitting a surface ocean circulation model to three sediment SST records using sequential methods of optimal estimation theory: a Kalman filter and a related smoother. These methods allow us to address some, but not all, of the limitations listed above. First, they allow us to account for the uncertainties in the reconstructed SSTs (and in the model) used for their interpretation. As a result, the significance of the frontal movements estimated from the combination of the data and the model is evaluated. Second, they allow us to infer a time-evolving SST field that is consistent (in the least squares sense) with both the data records and the physics of surface circulation as represented in the model. Our reconstruction of frontal movements is thus aimed at extending previous studies where these movements were estimated solely on the basis of sediment records.
This paper is organized as follows. The SST records and the surface circulation model are described in sections 2 and 3, respectively. The sequential methods used to combine the records and the model are presented in section 4. In section 5, these methods are applied in order to estimate the time-dependent SST field in a region in the North Atlantic over the past 14.5 kyr. In section 6, the significance and paleoceanographic implications of our results are discussed. Conclusions and possible extensions of this work follow in section 7.
The records of SST considered in this study1 have been generated from sediment cores raised from three locations in the North Atlantic (Fig. 1): core SU81–18 from the western Iberian margin (37°46′N, 10°11′W, water depth 3135 m), core CH69-K09 from the southern Newfoundland basin (41°45′N, 47°21′W, 4100 m), and core NA87–22 from the Rockall Trough (55°29′N, 14°41′W, 2161 m).
The chronologies of the three records have been established from 14C dates measured on monospecific planktonic foraminifera (Bard et al. 1987; Duplessy et al. 1992; Labeyrie et al. 1999). The number of planktonic foraminiferal 14C dates amounts to 30 for core SU81–18, 13 for core CH69-K09, and 25 for core NA87–22. Radiocarbon ages were converted into calendar age estimates using (i) the CALIB 4.1 software (Stuiver and Reimer 1993), (ii) the 310-yr moving average marine calibration curve of Stuiver et al. (1998), and (iii) reservoir ages constrained from the assumption that the two deglacial warming events apparent in each record were synchronous with the onset of the Bølling (
kyr BP) and the end of the YD (
kyr BP) in the GRIP ice core, central Greenland [dates from Waelbroeck et al. (2001)]. Further details about the sediment core chronologies can be found in that paper.
SST estimates for the three sediment cores were derived from planktonic foraminiferal counts. The faunal counts were converted into SST estimates for the “cold” and “warm” seasons using the revised analog method (RAM) developed by Waelbroeck et al. (1998). The SST values used in this study are obtained by averaging the cold and warm season SSTs reported in Waelbroeck et al. (2001) (Fig. 2), with the acceptance that the resulting averages might be biased estimates of the annual mean values. The uncertainties in the averaged SSTs are calculated by propagating the errors in the seasonal values (Bevington and Robinson 1992), assuming no error correlation. In the averaged SST records, the most recent value amounts to 
C for core SU81–18 (estimated calendar age of 0 yr BP), 
C for core CH69-K09 (550 yr BP), and 
C for core NA87–22 (530 yr BP), where the errors represent one standard deviation [note that an error of 0.02°C is likely a strong underestimate, reflecting the standard deviation of the core top values used to compute reconstructed SSTs (Waelbroeck et al. 1998)]. These values are not significantly different, at the level of two standard deviations, from the annual mean SSTs at the corresponding closest points of the World Ocean Atlas [cf. values defined as “statistical means” and “averaged decades” in the World Ocean Atlas 2013 (Locarnini et al. 2013)] (Fig. 2). The number of SST estimates for the past 14.5 kyr BP is equal to 24 for core SU81–18, 103 for core CH69-K09, and 96 for core NA87–22. These values imply a mean temporal resolution of, respectively, 604, 141, and 151 yr for these cores.
|
The SST records derived from the three sediment cores show noticeable temporal variations during the last deglaciation (Fig. 2). Common to all three records are the occurrences of relatively (i) high SSTs during the Bølling–Allerød (ca. 14.7–13.0 kyr BP), (ii) low SSTs during the YD (ca. 13.0–11.6 kyr BP), and (iii) high SSTs during the Holocene (ca. 11.6–0 kyr BP), which at least partly reflects the common chronological assumptions for the records. The average error in the reconstructed SSTs for the last 14.5 kyr amounts to 0.65°C for core SU81–18, 1.54°C for core CH69-K09, and 0.56°C for core NA87–22.
The model considered in this study is an advective model of the mixed layer (ML), where temperature and salinity are vertically uniform in the layer but where the effect of horizontal advection is retained. The motive for using an advective ML model is twofold. First, SST variability on interannual and longer time scales observed in instrumental records has been interpreted in terms of processes represented in such models (e.g., Deser et al. 2010). Second, a more complete model that also describes the dynamics of deeper layers would render the application of sequential methods impractical given the available resources. None of these justifications is fully satisfying. In particular, the path of the NAC–SPF system appears to be strongly constrained by topography (e.g., Sy et al. 1992; Bower and von Appen 2008), and it is clear that an ML model without any interior dynamics cannot provide a complete description of this system (e.g., Rossby 1996). In an effort to account for the limitations of an ML model, the model equations are not imposed exactly, but only in the mean square, in the estimation of the SST fields (section 4).
The ML extends from the sea surface at 
to the base of the mixed layer at 
(Fig. 3). Below the ML is a transition zone with a thickness δ and characterized by a relatively strong vertical density gradient. The transition zone separates the strongly turbulent ML above from the weakly turbulent stratified interior below. It is assumed to be thin compared to the ML (i.e., 
). Moreover, the horizontal motion at the base of the transition zone is taken to be small relative to that in the ML, as in previous models (e.g., de Szoeke 1980; Welander 1981; Cushman-Roisin 1981).
Our ML model is similar to the model of Welander (1981) with two modifications: the tendency term is retained in the ML temperature equation, and the geostrophic flow in the ML is considered in addition to the Ekman flow.2 The effects of both temperature T and salinity S on density ρ are taken into account through a linear equation of state:
are reference values, α is a thermal expansion coefficient, and β is a haline contraction coefficient.The equation for ML temperature is the governing equation of the model. It is derived from the vertical average of the temperature equation from 
to 
:
is the ML average of horizontal velocity, 
is the interior vertical velocity at the depth 
, 
is the interior temperature at that depth, 
is a parameter that describes the strength of heat exchange between the atmosphere and the ML, and 
is an apparent atmospheric temperature (Welander 1981; de Ruijter 1983). Also, t is time and 
is the horizontal gradient operator. According to (2), warming or cooling in the ML is caused by horizontal heat advection, heat exchange with water from below the ML, and (or) heat exchange with the atmosphere.
The horizontal velocity in the ML, 
, is the sum of a geostrophic velocity, 
, and an Ekman velocity, 
(i.e., 
), where
is the Coriolis parameter, where 
is the Earth angular velocity and ϕ is latitude, 
is a unit vector pointing upward, 
is a horizontal stress, and p is pressure. Clearly, 
. Consider first 
, the vertical average of the Ekman velocity in the ML. This is easily obtained by integrating the Ekman momentum equation from 
to 
, which leads to
is the surface wind stress, and the stress at the bottom of the transition zone has been assumed to be negligible.Consider then 
, the vertical average of the geostrophic velocity in the ML. The second equation in (3) is differentiated with respect to depth. This yields, using the hydrostatic approximation, 
, where g is the acceleration due to gravity, and the equation of state (1),
. Averaging this equation from 
to 
yields, using (5),
, has been assumed to be small compared to the average geostrophic velocity in the ML, 
, as in previous models (e.g., Cushman-Roisin 1981). This assumption implies that the effects of geostrophic motion below the ML are neglected and is perhaps the strongest assumption of our model.Note that, under the approximations above, the thermal contribution to the geostrophic flow,
vanishes identically. As a result, heat is carried horizontally in the model only by the Ekman flow, 
, and by the saline contribution to the geostrophic flow,
and 
The equation for ML temperature (2) includes, in addition to the effect of horizontal advection, the heat exchange with water from below the ML and the heat exchange with the atmosphere. The first term depends on the interior vertical velocity 
, whereas the second depends on 
. Consider first 
. This is obtained by integrating the statement of volume conservation, 
, from 
to 
:
[an assumption implicit in (2)]. The vertical velocity near the base of the ML, 
, is thus equal to the divergence of the total horizontal transport in the ML.Consider then 
. Haney (1971) proposed, as a surface thermal boundary condition for ocean circulation models, that the surface heat flux be set equal to 
, where 
is a parameter with units of W m−2 °C−1. Based on zonal- and time-averaged quantities, this author reported observational estimates of 
in the range from about 25 W m−2 °C−1 to about 50 W m−2 °C−1. The parameter 
would therefore vary between 
and 
m s−1, considering a density 
kg m−3 and a heat capacity at constant pressure 
J kg−1 °C−1 for seawater. Unless stipulated otherwise, a constant value 
m s−1 is used in this study. The parameters of the ML model are listed in Table 1.
|
Table 1. Parameters of the ocean mixed layer model.
The domain of the ML model is a region between 36°–62°N and 11°–47°W in the North Atlantic (Fig. 1). The temperature values at the boundaries of the domain are assumed to be governed solely by the heat exchange with the atmosphere,
In this section, we illustrate the ML temperature distribution and surface circulation in the modern North Atlantic, which are calculated by the model and compare briefly the model solution with observational estimates. The ML temperature equation (2) is integrated to steady state using the numerical method described in appendix A and with values of 
derived from modern climatologies (appendix B). The initial conditions are provided by modern annual mean SSTs (Locarnini et al. 2013) averaged in model grid cells.
The simulated ML temperature distribution (Fig. 4) shows a modest adjustment compared to the initial conditions (Fig. 1), which reflects the fact that 
is based on the modern SSTs (appendix B) and the importance of surface flux in the ML heat balance (2). The importance of surface heat flux suggests that the boundary condition (11) provides a reasonable approximation of the ML heat balance in the interior of the domain.
|
The simulated distribution of ML total velocity, 
, shows a generally southeastward flow (Fig. 4a). This flow reveals the Ekman contribution to 
in our North Atlantic domain, where the zonal component of surface wind stress is predominantly eastward (Risien and Chelton 2008). Temperature and salinity tend to offset their effects on horizontal density gradients and hence on the geostrophic contribution to 
(not shown). Indeed, observational studies have shown that large-scale fronts in the North Atlantic are partially compensated, with waters tending to be warm and salty on one side of the fronts and cool and fresh on the other side (Stommel 1993; Chen 1995).
Observational estimates of North Atlantic surface circulation have been derived from drifters which have drogues attached at depth (15 or 100 m) to minimize the effects of the winds, whether direct or indirect through the Ekman currents (Rossby 1996; Fratantoni 2001). To compare consistently these estimates with our model, we consider the simulated distribution of 
, the geostrophic part of the total ML velocity (Fig. 4b). The simulated field of 
reveals a generally northeastward flow, which agrees qualitatively with the drifter observations.
Notice that agreement is not found throughout the domain, however. In the region between 40°–50°W and 40°–50°N, the simulated velocity 
(Fig. 4b) does not reproduce the northward flow along the continental rise and the anticyclonic eddy [the Mann eddy (Mann 1967)] southeast of the Grand Banks, as inferred from drifters (e.g., Fratantoni 2001). In this region, the mean speed of surface current estimated from drifters [38 cm s−1 (Fratantoni 2001)] exceeds the simulated values of 
by one order of magnitude. Among the possible sources of disagreement are the model limitations, in particular the lack of deep ocean dynamics and the coarse horizontal resolution.
Sequential methods are applied in order to combine the sediment SST records (section 2) with the surface circulation model (section 3). They will allow us to produce an estimate of the time-dependent state of the surface North Atlantic over the past 14.5 kyr, which is consistent both observationally and dynamically. For conciseness, the problem of estimating past ocean states is referred below to as the inverse problem.
The elements of the state vector are the variables that define the state of the surface North Atlantic in terms of the model. They are the unknowns of the inverse problem, including the distributions of the ML temperature, the apparent atmospheric temperature, the interior temperature, the ML depth, and the Ekman and saline contributions to 
. The thermal contribution to 
is not considered as a separate variable, since it can be determined diagnostically from the ML temperature field using (8).
The state vector (
below) is more specifically defined as follows. Let 
and 
be, respectively, the zonal and meridional component of the vector sum 
. The apparent atmospheric temperature 
, the interior temperature 
, the ML depth 
, and the velocity components 
are each a function of longitude λ, latitude ϕ, and time t. The inclusion in the state vector of the gridded field of all these variables, however, would lead to a vector 
with a dimension of O(
), which would make the application of sequential methods impractical given the available resources. To lower the dimension of 
, a “state reduction” approach is adopted, whereby some of the variables in 
are approximated by analytic functions of spatial coordinates (e.g., Gaspar and Wunsch 1989). Specifically, the variables 
are each represented as the sum of a polynomial function and a residual term:
is a time-dependent coefficient, 
and 
, where 
W and 
N (the point 
is near the “center” of the domain), 
are integral exponents, and 
is a residual term that varies with 
. The coefficients 
determine the time-dependent part of the distributions of 
. Their evolution is described in probabilistic terms (random walk):
is a purely random process (noise) with zero mean and constant variance. The effect of this description is to impose some temporal covariance and nonstationarity on the fields of 
. The state is then defined by the gridded field of ML temperature (including the T values at the boundaries) and by the coefficients 
, 
:
by their respective polynomial coefficients can reduce the dimension of the state vector very significantly, depending on K. Unless stated otherwise, these fields are approximated with 
terms (Table 2), which reduces the dimension of the state vector to 306 elements and entails some loss of horizontal resolution.
|
Table 2. Integral exponents of the polynomial approximations.
The state vector 
satisfies two equations, both of which are written below in discrete form [e.g., 
is written as 
where i is a time index]. The first equation is the observation equation
is a vector of observations, 
is a matrix that relates the state elements to the observations, and 
is a vector of observational errors. In this study, the observations are the sediment records of SST (Fig. 2) and, for the modern time (0 yr BP), observational estimates of the elements of 
derived from modern data. These estimates come from (i) the annual mean SSTs reported in Locarnini et al. (2013) and (ii) the coefficients 
, 
, determined from modern observations (appendix C).The second equation satisfied by the state 
is the dynamic equation,
is a vector of functions of 
, and 
is another vector of errors. The functions in 
are obtained from (i) the finite-difference approximations of the ML temperature equation (2) and of the boundary condition (11), and (ii) the probability model for the polynomial coefficients (13) (see appendix A). Note that some of the functions in 
contain products of state elements, which makes our inverse problem nonlinear.Estimates of SST derived from the sediment (section 2) and the ocean model considered for their interpretation (section 3) have sizeable uncertainties and limitations. Accordingly, neither the SST estimates nor the model equations should be imposed perfectly when analyzing the SST records. The following assumptions are made about the data errors 
and the model errors 
. Both 
and 
are assumed to have zero mean and to show no temporal correlation. Moreover, they are taken to be mutually uncorrelated at any time. Collectively,
is the expected value (the mean) and 
is the Kronecker delta (i.e., 
if 
and 
otherwise). The quantities 
and 
are the covariance matrices for the observational errors and the model errors, respectively.The matrix 
is constructed as follows. Its diagonal elements are the variances of the errors in the observations defined in 
: that is, in the sediment SSTs and in the modern observational estimates of the state elements. The error variances for the sediment SSTs are set equal to the square of the errors displayed in Fig. 2. The error variances for the modern SSTs are calculated by propagating the standard errors in the annual mean SSTs (Locarnini et al. 2013) contributing to the gridcell averages (appendix B), assuming no error correlation. The error variances for the modern estimates of 
are determined as described in appendix C. The off-diagonal elements of 
are the covariances of the errors in the observations in 
. They are set to zero, except for the modern estimates of 
(appendix C).
The matrix 
is established as follows. Its diagonal elements are the variances of the errors in the equations defined in 
, and its off-diagonal elements are the covariances between these errors. The error variances are assumed to be the same at every time, and the error covariances are taken as zero, so that 
is time invariant (
) and diagonal. The error variances for the discrete forms of the ML temperature equation (2) and of the boundary condition (11) are arbitrarily taken as proportional to the spatial average of annual mean SSTs in the model domain according to modern observations (Locarnini et al. 2013). This average is noted 
, where 
refers to the modern time (0 yr BP). In the same vein, the variance of the noise in the probability model (13) is set proportional to the magnitude of 
as estimated from modern observations. Thus,
The conventional Kalman filter assumes that the equations in 
and the observations in 
are linear functions of the state elements in 
(e.g., Anderson and Moore 1979). Consequently, it cannot be applied here given the nonlinearities in some of the equations. Different versions of the Kalman filter have been proposed to account for the presence of nonlinearities in the system dynamics and (or) in the observations [for a short review, see Wunsch (2006b)]. Two of these versions are considered in this paper: the extended Kalman filter (EKF) and the linearized Kalman filter (LKF). Both versions rely on the linearization of the equations in 
, but around a different state. In the EKF, the equations in 
are linearized about the most recent state estimate at time i. In the LKF, they are linearized about a reference state, which in this study is a constant state constrained from modern observations. The ML temperatures estimated from the EKF and LKF can be considered as weighted least squares estimates, with the weighting provided by 
and 
. Details about the two filters are provided in appendix D.
The following notation is adopted in this paper. The Kalman filter estimate of the state at time i is denoted as 
(+), where the plus sign indicates that the estimate considers the data prior to and at time i. This notation is used to contrast 
(+) from 
(−), which is the estimate at time i that only considers the data prior to this time. The error covariance matrix, or uncertainty, of 
(
) is 
, where 
is the true state at time i.
For both the EKF and LKF, the initial state estimate at 14 500 yr BP, 
(+), is determined from modern observations, and its uncertainty, 
, includes relatively large error (co)variances to account for the fact that the state of the surface North Atlantic at 14.5 kyr BP may have been different from the modern state. Consider first 
(+). The initial ML temperatures are obtained by averaging in model grid cells the modern annual mean SSTs at 1° × 1° resolution of Locarnini et al. (2013). The initial values of the polynomial coefficients 
, 
, are determined from modern observations, as described in appendix C. Consider then 
. The error variances for the initial ML temperatures are set equal to the spatial variance of these temperature values, which amounts to 
. The error covariances for the initial ML temperatures are taken as zero. The error (co)variances for the initial coefficients 
, 
, are taken as 4 times the error (co)variances for the modern estimates of these coefficients (appendix C). Thus, the errors at 14.5 kyr BP in the spatial patterns of 
, which are described by (12), are set equal to twice the corresponding errors at 0 kyr BP.
We find that the EKF and LKF lead to very similar estimates of ML temperature and of their errors near the sediment core locations (appendix D). Thus, the estimation of the North Atlantic ML temperature field and of its errors does not seem to depend appreciably on whether the equations in 
are linearized about the most recent state estimate, 
(+), or about the initial state estimate, 
(+). Solutions obtained from the LKF are only considered in the remainder of the paper.
The application of a Kalman filter allows us to derive estimates of the time-evolving state of the surface North Atlantic that are constrained by both the SST records and the ML model. A Kalman filter, however, constitutes an incomplete analysis, in the sense that the filter estimates are based solely on past and present observations: the filter estimate 
(+) is determined from data prior to and at time i but is not constrained by observations posterior to time i. To account for posterior observations, a smoother should be used (e.g., Anderson and Moore 1979).
A linearized smoother is applied in order to estimate states of the surface North Atlantic that are consistent with the entire SST records and the modern observations, as well as with the model physics. Since the state of the surface North Atlantic is to be estimated at times within a fixed interval, a so-called fixed-interval smoother is applied (e.g., Anderson and Moore 1979). As for the filter estimates, the ML temperatures estimated from the smoother can be viewed as weighted least squares estimates, but with the weighting provided by 
, 
, and 
. Details about the smoother can be found in appendix E.
In this section, the LKF is applied to estimate the time-dependent state of the surface North Atlantic from 14 500 yr BP (
) to 0 yr BP (
). Our goal is to illustrate the effect of model errors on the state estimates and to produce a solution that will constitute the first part of the smoothing solution to be discussed in section 5b.
The application of a Kalman filter such as LKF requires knowledge of the covariance matrices for the data and model errors (appendix D). Whereas the data errors are relatively well understood, the model errors are more difficult to constrain. However, the sequence 
(appendix D) contains information about the ability of the model to replicate observations and hence about model errors. Consider a linear and time-invariant system: that is, a system with 
and 
, where 
are constant matrices and 
is a deterministic forcing. A filter of such a system is optimal provided that 
is a normal white noise sequence with covariance 
(Mehra 1970). The properties of 
, known as the innovation sequence, can therefore be consulted in order to assess the optimality of the filtering solution. Although the system (15)–(16) is not linear and not time invariant (e.g., 
is variable), the innovation properties are nonetheless useful for evaluating the ability of the model to explain the observations and for constraining a plausible choice of model errors in 
, as shown below.
We consider three solutions of the LKF obtained from different assumptions about the model errors in 
: 
, 
, and 
. For each solution, we illustrate both the elements of 
(Fig. 5) and the distribution function of the innovation elements normalized to their standard deviations (Fig. 6). Note that the elements of 
, which includes modern observations, are removed from both figures in order to isolate the ability of the model to replicate the sediment SST records.
|
|
We find that, for the solution with 
(relatively small model errors), the elements of 
average −0.8°C, with a standard error of 0.1°C, and reconstructed SSTs between ca. 11 and 13 kyr BP are strongly overestimated (Fig. 5a). This solution is apparently biased: the SST records are poorly replicated due to the prescription of too small model errors (Fig. 6a). In contrast, for the solution with 
(higher model errors), the innovation elements average −0.1°C with a standard error of 0.1°C, and the reconstructed SST values between ca. 11 and 13 kyr BP are better reproduced (Fig. 5c), revealing no apparent bias. However, the close fit of the filter temperatures to the sediment SSTs that is obtained in this case does not seem to be warranted given the errors in the sediment SSTs (Fig. 6c), suggesting that the choice 
is giving too much confidence in the data. The intermediate choice 
leads to a solution that shows both no clear bias (Fig. 5b) and no clear sign of under- or overfit to the SST records (Fig. 6b).
The time series of ML temperature near the sediment core locations as estimated by the LKF for 
are shown in Fig. 7. As expected, the temperatures estimated from filtering are generally consistent with those reconstructed from the faunal counts for each core: the differences between the filtering and reconstructed temperatures are generally less than one standard deviation in the filtering estimates.
|
Notice the discontinuities in the time series of the filter temperatures and of their errors, which occur when SST observations are available (Fig. 7; see also Figs. D1 and D2). When observations are available, the filter combines the model forecast δ
(−) with the observations 
to produce the filter estimate δ
(+) according to (D9) (appendix D). The state then deviates from its evolution driven by the model to a degree that depends on the Kalman gain (i.e., on the relative influence of 
and 
). Furthermore, the errors in the filter temperature estimates are reduced when observations are available (Fig. 7; see also Fig. D2). This behavior can be understood from the following equation:
(appendix D) using the matrix inversion lemma (Liebelt 1967). As shown by (20), the uncertainty of the state after measurement, 
, is never larger than 
, since 
is at least positive definite (Bryson and Ho 1975). In contrast, when observations are not available, the evolution of the state is entirely governed by the model (D8a), and the state uncertainty evolves according to the error propagation equation (D8b).In this section, the linearized smoother is applied to estimate the time-dependent state of the surface North Atlantic over the past 14 500 yr (with 
). The time series of ML temperature near the sediment core locations as estimated by the smoother are shown in Fig. 8. In contrast to the filtering estimates (Fig. 7), each smoothing estimate of temperature is constrained by the entire dataset, including by the data posterior to the time of the estimate. The consideration of posterior data has two noticeable consequences. First, compared to the filtering estimates, the smoothing estimates of temperature present reduced discontinuities at times when data are available (hence the name “smoother”). Second, the errors in the smoothing estimates are lower than the errors in the filtering estimates, as anticipated from (E4) (cf. Fig. 8 with Fig. 7).
|
We consider the position of the 10°C ML isotherm in the North Atlantic domain during the Younger Dryas, which is estimated in the smoothing solution (Fig. 9). In the modern North Atlantic, the 10°C surface represents reasonably well the 
kg m−3 surface in the main thermocline, and its rapid shoaling near 50°N indicates the path of the NAC–SPF system (Rossby 1996). As shown in Fig. 9, the 10°C ML isotherm is estimated to have been more zonal and located more to the south at 12 000 yr BP compared to today. Surface waters warmer than 10°C would have been confined to south of 48°N throughout the basin in the annual mean sense. The errors in the smoothing estimates of ML temperature at 12 000 yr BP vary from 0.54° to 0.94°C, with minima in the latitude band between the sediment core locations and maxima outside the band (dashed lines in Fig. 9). This pattern results from a combination of model dynamics 
, model errors 
, data distribution 
, and data errors 
, as shown by the error covariance equations of the filter (D8b) and the smoother (E4).
|
The estimated meridional shift of the 10°C ML isotherm is much larger in the east than in the west of the studied domain (Fig. 9), despite the fact that the three SST records show comparable warming from the YD to the present time (Fig. 2). The smoothing solution is generally consistent with the YD data for each record (Fig. 8) but features a much greater shift of the isotherm near the longitudes of eastern cores SU81–18 and NA87–22 than near the longitude of western core CH69-K09 (Fig. 9). We interpret this result as being due to the location of the western core CH69-K09 within a region of large meridional SST gradients (Fig. 1).
In this section, we discuss in more detail the deglacial movements of the 10°C ML isotherm in the North Atlantic that are estimated in the smoothing solution. This is done with the understanding that a single isotherm may not be appropriate to characterize the path of the NAC–SPF system over the entire North Atlantic, particularly because (i) temperature may vary along the path, and (ii) the system may be composed of different branches (e.g., Sy et al. 1992; White and Heywood 1995; Read et al. 2010). We then clarify the implications of our results for North Atlantic deglacial oceanography.
We consider the estimated meridional movements of the 10°C isotherm at three different longitudes: 13°W (close to the longitudes of cores SU81–18 and NA87–22), 29°W (approximately the central longitude of the domain), and 47°W (close to the longitude of core CH69-K09). At each of these longitudes, the latitude of the 10°C isotherm, 
, is determined by linear interpolation from the encompassing temperatures 
and 
,
(
) is the latitude of the grid point where 
(
). The uncertainty in 
, 
, is calculated by error propagation with due regard for error correlation,
(
) is the error (standard deviation) of 
(
) and 
is the covariance between these errors. The partial derivatives in (22) are evaluated analytically from (21).The meridional displacements of the 10°C isotherm that are estimated by the smoother are very distinct at 13°, 29°, and 47°W (Fig. 10). In the eastern North Atlantic (13°W), the displacements at the start and end of the YD are estimated to reach about 15° (Fig. 10a). They are significant at the level of two standard deviations in the isotherm position 
, the particular significance level depending on the precise pair of dates for which 
values are compared. The meridional variations of 
associated with the YD are clearly the most pronounced and most significant ones over the past 14.5 kyr. In the central North Atlantic (29°W), the 10°C isotherm would have experienced meridional displacements of only about 5° across the YD (Fig. 10b). Note the relatively large uncertainties in the isotherm position over the past 9–10 kyr at this longitude. They stem primarily from the fact that in the smoothing solution the meridional thermal gradients near the latitude of the 10°C isotherm are small during this time interval (not shown), leading to enhanced values of the partial derivatives in (22) and hence to enhanced values of 
. Finally, in the western North Atlantic (47°W), the 10°C isotherm is estimated to have remained within a few degrees of the latitude of 45°N (Fig. 10c). Collectively, these results imply that the 10°C isotherm would have pivoted twice around a region southeast of the Grand Banks, with a SW–NE orientation during the Bølling–Allerød and the Holocene and a more zonal orientation and southerly position during the YD. This pivotal motion is consistent with previous depictions of the polar front movements during past climate changes (e.g., Ruddiman and McIntyre 1981; Zahn 1994; Barker et al. 2015).
|
We focus on the pronounced meridional variations of the 10°C isotherm at 13°W that are inferred at the onset and termination of the YD (Fig. 11). An apparent speed of migration of the isotherm during the two periods is estimated from a weighted least squares fit of 
versus time, with the weighting provided by the variable 
. We find that the isotherm would have moved southward at an apparent speed of 
km yr−1 from 13.5 to 13.0 kyr BP and northward at an apparent speed of 
km yr−1 from 11.6 to 11.0 kyr BP, where the errors are standard errors. These values are comparable to the apparent speed of advance (retreat) of the polar front at the onset (termination) of the YD between cores CH73–139C and SU81–18 (Bard et al. 1987). According to these authors, the front would have moved in less than 400 yr over the distance of about 1930 km between these two cores (Fig. 1): that is, at an apparent speed of less than 4.8 km yr−1.
|
The front migration speeds estimated here and in prior work should be interpreted with caution. The finite deposition rates, discrete sampling along the core, and bioturbation imply that the sediment records have a bounded resolution, limited here to centennial and lower-frequency variability. Besides, the front migration speeds estimated from sediment records arise from relatively small phase differences between the records and hence may be particularly sensitive to the assumptions about core chronologies.
We consider the effects on our results of (i) the parameter 
, which determines the strength of air–sea heat exchange, and (ii) the number K of terms retained in the polynomial functions for 
(12).
Consider first the effect of 
. We produce two other smoothing solutions with 
m s−1 and 
m s−1, which approximate the range of values of 
reported by Haney (1971) (section 3d). For both solutions, the estimated latitudes of the 10°C isotherm, 
, show small differences compared to those estimated in our reference solution (section 5b): for each longitude (13°, 29°, and 47°W), the differences in the estimated latitude are always less than 0.5° (notice that 
values at times when 
is estimated to be present at more than one latitude at a given longitude are not considered in the comparison). In the solution with 
m s−1, the speed of isotherm movement amounts to 
and 
km yr−1 between 13.0 and 13.5 kyr BP and 11.0 and 11.6 kyr BP, respectively. In the solution with 
m s−1, they amount to 
and 
km yr−1, respectively. Thus, our results do not seem to be sensitive to 
, provided that 
is in the range determined from modern observations.
Consider then the effect of K, the number of terms retained in (12). A smoothing solution is obtained with 
, thereby increasing by 5 the number of terms in each polynomial (Table 2). The estimated latitudes of the 10°C isotherm are comparable to those of the reference solution: for each longitude (13°, 29°, and 47°W), the absolute differences in the estimated latitude average to less than 1°, with a maximum value of 5° (again, excluding 
values at times when 
is present at more than one latitude at a given longitude). The speed of isotherm movement reaches 
and 
km yr−1 between 13.0 and 13.5 kyr BP and 11.0 and 11.6 kyr BP, respectively. As for 
, the smoothing solution does not seem to be very sensitive to K, at least within the range of values being considered.
The past states of the surface North Atlantic that are inferred in this study are intended to represent annual averages. Although estimates of cold- and warm-season SSTs are available at the core locations (section 2), no attempt is made to infer the SST field at subannual time scales. Thus, the sediment SSTs considered are averages of the seasonal values (section 2), and the observations used to establish the initial and final states are annual means (sections 3–4). In fact, the ML temperatures estimated in this work remain above −1.9°C, the freezing point of seawater for a salinity of 35 and zero hydrostatic pressure (Millero 1978) (the sole exception is for the EKF solution discussed in appendix D, for which the ML temperature at the northern- and western-most grid point of the domain is −2.1°C at 13 140 yr BP and −1.9°C thereafter). This result does not imply that sea ice was not present but that the (annual mean) SST records being considered do not require its presence in our North Atlantic domain. Below, two specific mechanisms of large-scale frontal movements in the deglacial North Atlantic are discussed that do not explicitly involve a role of sea ice.
Keffer et al. (1988) proposed that the Laurentide Ice Sheet (LIS) may have modified the North Atlantic atmospheric circulation in such a manner as to cause the line of zero wind stress curl to become more zonal and, in turn, force a more zonal boundary between the subtropical and subpolar gyres during the LGM. It is not clear, however, that changes in the North Atlantic wind field associated with LIS changes were fast enough to cause the front to move over the relatively short period of 500–1000 yr (Fig. 11). Keffer et al. (1988) stressed that this mechanism would account only for oceanic changes over the relatively long time scales characterizing the buildup and decay of continental ice sheets. It would not explain the YD event (Keffer et al. 1988), although the effects of a varying size or shape of continental ice sheets have more recently been implicated in abrupt climate changes (e.g., Schulz et al. 1999; Jackson 2000; Schulz 2002; Wunsch 2006a; Zhang et al. 2014).
Rossby and Nilsson (2003) proposed an ocean current switching mechanism that involves a rapid movement of the front between the two gyres at the termination of the YD. They suggested the operation of a salinity-driven feedback that amplifies the production of North Atlantic Deep Water. According to this alternative mechanism, a drop in steric height or sea level at the Iceland–Scotland Ridge would be transmitted southward and westward by topographic and planetary waves. The NAC near the Grand Banks (the Gulf Stream in their paper) would respond to the resulting perturbation along the western boundary by transporting warm salty water north along a topographically defined path, thereby reinforcing deep water production. The rapid switch of the NAC–SPF system from its glacial zonal to present meridional course, the authors argued, would have led to the abrupt end of the YD as inferred from Greenland ice core records (e.g., Alley et al. 1993).
We compare the time scales of signal transmission by wave propagation and salt advection with the time scale of front migration as inferred in this study. Rossby and Nilsson (2003) estimated that the initial perturbation at the Iceland–Scotland Ridge would be signaled by topographic and planetary waves to the northeastern coast of North America in about 5 months and that the transit of salty water from the tail of the Grand Banks to Iceland would take about 2 yr. Although these figures are rough estimates, as pointed out by the authors, the collective time scale of about 2.5 yr is smaller by two to three orders of magnitude than the time scale of about 500–1000 yr inferred in this study (Fig. 11). The mismatch would provide evidence against the current switching envisioned by Rossby and Nilsson (2003), with the caveat that the bounded resolution of sediment records and systematic errors in their phase relationships may bias the estimates of isotherm migration speed (section 6a).
In this paper, a coarse-resolution model of surface ocean circulation is fitted to three sediment records of SST using sequential methods. The result of the fit is an estimate of model variables, such as mixed layer temperature, and of their errors in a region of the North Atlantic over the past 14 500 yr. In contrast to previous studies, the estimated evolution of ML temperatures is consistent with a dynamical description of surface flow, and formal uncertainties are provided. Emphasis is placed on the meridional shifts of the 10°C isotherm, which roughly coincides with the modern pathway of the North Atlantic Current–Subpolar Front system.
We find that the 10°C isotherm would have pivoted twice around a region southeast of the Grand Banks, with a SW–NE orientation during the Bølling–Allerød and the Holocene and a more zonal orientation and southerly position during the Younger Dryas. This result is broadly consistent with previous inferences about the “polar front” in the paleoceanographic literature (e.g., Ruddiman and McIntyre 1981). In the eastern North Atlantic (13°W), the estimated deglacial movements of the isotherm span over 15° of latitude and are significant at the level of two standard deviations in the isotherm position. At this longitude, the isotherm would have migrated southward at a speed of 
km yr−1 between 13.0 and 13.5 kyr BP and northward at a speed of 
km yr−1 between 11.0 and 11.6 kyr BP (reference solution), where the errors reflect the error (co)variances of the estimated ML temperature fields.
The time scale of the inferred frontal movements would amount to less than 500–1000 yr, assuming negligible errors in the phase differences between the sediment records. It appears too small to support the LIS-induced mechanism envisioned by Keffer et al. (1988), although more recent ideas have invoked a role of varying topography of continental ice sheets in millennial-scale climate changes (e.g., Jackson 2000; Wunsch 2006a; Zhang et al. 2014). On the other hand, it is greater by two to three orders of magnitude than the time scale (ca. 2.5 yr) involved in the current switching mechanism to explain the abrupt termination of the YD (Rossby and Nilsson 2003). The bounded resolution of sediment records and the uncertainties in their phase relationships, however, prevent testing this mechanism with high confidence.
Extensions of this work appear endless. Needless to say, the present experiments should be repeated with a larger number of SST records, which ideally should meet the following criteria. First, the records should be based on well-constrained age models with a consistent set of chronological assumptions (radiocarbon calibration, reservoir ages, etc.). Second, they should possess sufficiently high resolution in order to document SST changes associated with relatively rapid events such as the YD. Finally, they should rely on the same reconstruction method in order to avoid offsets arising from the consideration of different temperature indicators and (or) different calibrations to SST. Alternately, different temperature proxies could be combined in an attempt to reduce the bias in SST estimates derived from sediment indicators (MARGO 2009). The chronology of the sediment records considered in this study strongly depends on the assumed synchronism with Greenland records [Waelbroeck et al. (2001); for a recent discussion of this approach, see Austin and Hibbert (2012)]. Particularly desirable would be records that are both free of such an assumption and based on reliable reservoir ages. Records from sites that would bracket or delineate the Subpolar Front at different times would be particularly useful, such as from the northern Newfoundland basin closer to the Northwest Corner and from the Irminger Sea. The analysis of a large number of records meeting the above criteria, when available, might lead to different results, especially in the western North Atlantic.
Likewise, experiments should be conducted with more complete ocean models, in particular models that represent motion both within and below the mixed layer. A more complete description of the NAC–SPF system and of its variability would require the consideration of Sverdrup dynamics and bottom topography (e.g., Sy et al. 1992; Bower and von Appen 2008; but see Wunsch 2011). The combination of complex models with (long) paleoceanographic records, however, would involve substantial resources given the computational and storage demands of sequential methods. In fact, conventional methods such as the EKF and LKF may not be practical in this case, and other approaches to filtering and smoothing may be needed (e.g., Fukumori 2002; Evensen 2003).
The present experiments could also be extended to include time intervals prior to the Bølling, in particular the LGM. A major difficulty in this regard is the prescription of a consistent initial state: whereas a number of reconstructed SSTs exist for the LGM (e.g., MARGO 2009), other key aspects of the glacial North Atlantic, such as sea surface salinity, mixed layer depth, wind stress components, etc., are still very poorly known. While use could be made of results from climate models subject to glacial boundary conditions (e.g., Braconnot et al. 2007) or dynamical reconstructions derived from the combination of LGM data with a model (e.g., Dail and Wunsch 2014), estimates of state error covariance are generally not available. The difficulty of establishing a consistent initial state motivated us to start the experiments at a time (14 500 yr BP) during the Bølling, with the assumption that the state of the surface North Atlantic was then comparable to the modern one.
Finally, assumptions in the filtering and smoothing procedures would need to be explored. The state estimates obtained in this study are only as reliable as the error covariance matrices 
and 
from which they have been derived. Future work will also need to assess the consequences of representing property fields with analytic functions of longitude and latitude, as these should alter the observability of the dynamic system (e.g., Chen 1999; Marchal 2014).
Acknowledgments
The authors express their gratitude to Daniel Amrhein, Geoffrey Gebbie, Lloyd Keigwin, Delia Oppo, and Carl Wunsch for useful and detailed comments on an earlier version of the manuscript. Discussions with Katherine Allen, Brian Farrell, Eli Tziperman, and Carl Wunsch have also been very helpful. The constructive and detailed remarks from Thomas Rossby and two anonymous reviewers have significantly improved the manuscript. Sarah Shipley worked as a guest student on several aspects of this project. OM acknowledges support from the U.S. National Science Foundation. CW acknowledges support from the European Research Council ERC Grant ACCLIMATE 339108.
| Alley, R. B., and Coauthors, 1993: Abrupt increase in Greenland snow accumulation at the end of the Younger Dryas event. Nature, 362, 527–529, doi:https://doi.org/10.1038/362527a0. Crossref, Google Scholar | |
| Anderson, B. D. O., and G. B. Moore, 1979: Optimal Filtering. Prentice Hall, 368 pp. Google Scholar | |
| Austin, W. E. N., and F. D. Hibbert, 2012: Tracing time in the ocean: A brief review of chronological constraints (60–8 kyr) on North Atlantic marine event-based stratigraphies. Quat. Sci. Rev., 36, 28–37, doi:https://doi.org/10.1016/j.quascirev.2012.01.015. Crossref, Google Scholar | |
| Bard, E., 2001: Comparison of alkenone estimates with other paleotemperature proxies. Geochem. Geophys. Geosyst., 2, 1002, doi:https://doi.org/10.1029/2000GC000050. Crossref, Google Scholar | |
| Bard, E., M. Arnold, P. Maurice, J. Duprat, J. Moyes, and J.-C. Duplessy, 1987: Retreat velocity of the North Atlantic polar front during the last deglaciation determined by 14C accelerator mass spectrometry. Nature, 328, 791–794, doi:https://doi.org/10.1038/328791a0. Crossref, Google Scholar | |
| Barker, S., J. Chen, X. Gong, L. Jonkers, G. Knorr, and D. Thornalley, 2015: Icebergs not the trigger for North Atlantic cold events. Nature, 520, 333–336, doi:https://doi.org/10.1038/nature14330. Crossref, Google Scholar | |
| Belkin, I. M., and S. Levitus, 1996: Temporal variability of the Subarctic Front near the Charlie-Gibbs Fracture Zone. J. Geophys. Res., 101, 28 317–28 324, doi:https://doi.org/10.1029/96JC02794. Crossref, Google Scholar | |
| Bevington, P. R., and D. K. Robinson, 1992: Data Reduction and Error Analysis for the Physical Sciences. 2nd ed. McGraw-Hill, 328 pp. Google Scholar | |
| Bower, A. S., and W.-J. von Appen, 2008: Interannual variability in the pathways of the North Atlantic Current over the Mid-Atlantic Ridge and the impact of topography. J. Phys. Oceanogr., 38, 104–120, doi:https://doi.org/10.1175/2007JPO3686.1. Link, Google Scholar | |
| Braconnot, P., and Coauthors, 2007: Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum—Part 1: Experiments and large-scale features. Climate Past, 3, 261–277, doi:https://doi.org/10.5194/cp-3-261-2007. Crossref, Google Scholar | |
| Bryson, A. E., and Y.-C. Ho, 1975: Applied Optimal Control. Taylor and Francis, 482 pp. Google Scholar | |
| Chapman, M., N. J. Shackleton, M. Zhao, and G. Eglinton, 1996: Faunal and alkenone reconstructions of subtropical North Atlantic surface hydrography and paleotemperature over the last 28 kyr. Paleoceanography, 11, 343–357, doi:https://doi.org/10.1029/96PA00041. Crossref, Google Scholar | |
| Chen, C. T., 1999: Linear System Theory and Design. Oxford University Press, 334 pp. Google Scholar | |
| Chen, L. G., 1995: Mixed layer density ratio from the Levitus data. J. Phys. Oceanogr., 25, 691–701, doi:https://doi.org/10.1175/1520-0485(1995)025<0691:MLDRFT>2.0.CO;2. Link, Google Scholar | |
| Clark, P. U., and Coauthors, 2009: The Last Glacial Maximum. Science, 325, 710–714, doi:https://doi.org/10.1126/science.1172873. Crossref, Google Scholar | |
| CLIMAP, 1976: The surface of the Ice-Age Earth. Science, 191, 1131–1137, doi:https://doi.org/10.1126/science.191.4232.1131. Crossref, Google Scholar | |
| Colin de Verdière, A., H. Mercier, and M. Arhan, 1989: Mesoscale variability transition from the western to the eastern Atlantic along 48°N. J. Phys. Oceanogr., 19, 1149–1170, doi:https://doi.org/10.1175/1520-0485(1989)019<1149:MVTFTW>2.0.CO;2. Link, Google Scholar | |
| Cushman-Roisin, B., 1981: Effects of horizontal advection on upper-ocean mixing: A case of frontogenesis. J. Phys. Oceanogr., 11, 1345–1356, doi:https://doi.org/10.1175/1520-0485(1981)011<1345:EOHAOU>2.0.CO;2. Link, Google Scholar | |
| Dail, H., and C. Wunsch, 2014: Dynamical reconstruction of upper-ocean conditions in the Last Glacial Maximum. J. Climate, 27, 807–823, doi:https://doi.org/10.1175/JCLI-D-13-00211.1. Link, Google Scholar | |
| de Ruijter, W. P. M., 1983: Effects of velocity shear in advective mixed-layer models. J. Phys. Oceanogr., 13, 1589–1599, doi:https://doi.org/10.1175/1520-0485(1983)013<1589:EOVSIA>2.0.CO;2. Link, Google Scholar | |
| Deser, C., M. A. Alexander, S.-P. Xie, and A. S. Phillips, 2010: Sea surface temperature variability. Annu. Rev. Mar. Sci., 2, 115–143, doi:https://doi.org/10.1146/annurev-marine-120408-151453. Crossref, Google Scholar | |
| de Szoeke, R. A., 1980: On the effects of horizontal variability of wind stress on the dynamics of the ocean mixed layer. J. Phys. Oceanogr., 10, 1439–1454, doi:https://doi.org/10.1175/1520-0485(1980)010<1439:OTEOHV>2.0.CO;2. Link, Google Scholar | |
| Duplessy, J.-C., L. Labeyrie, M. Arnold, M. Paterne, J. Duprat, and T. C. E. Van Weering, 1992: Changes in surface salinity of the North Atlantic Ocean during the last deglaciation. Nature, 358, 485–488, doi:https://doi.org/10.1038/358485a0. Crossref, Google Scholar | |
| Evensen, G., 2003: The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dyn., 53, 343–367, doi:https://doi.org/10.1007/s10236-003-0036-9. Crossref, Google Scholar | |
| Eynaud, F., and Coauthors, 2009: Position of the Polar Front along the western Iberian margin during key cold episodes of the last 45 ka. Geochem. Geophys. Geosyst., 10, Q07U05, doi:https://doi.org/10.1029/2009GC002398. Crossref, Google Scholar | |
| Fraser, D. C., 1967: A new technique for the optimal smoothing of data. Ph.D. dissertation, Massachusetts Institute of Technology, 195 pp. Google Scholar | |
| Fratantoni, D., 2001: North Atlantic surface circulation during the 1990’s observed with satellite-tracked drifters. J. Geophys. Res., 106, 22 067–22 093, doi:https://doi.org/10.1029/2000JC000730. Crossref, Google Scholar | |
| Fukumori, I., 2002: A partitioned Kalman filter and smoother. Mon. Wea. Rev., 130, 1370–1383, doi:https://doi.org/10.1175/1520-0493(2002)130<1370:APKFAS>2.0.CO;2. Link, Google Scholar | |
| Gaspar, P., and C. Wunsch, 1989: Estimates from altimetric data of barotropic waves in the northwestern Atlantic Ocean. J. Phys. Oceanogr., 19, 1821–1844, doi:https://doi.org/10.1175/1520-0485(1989)019<1821:EFADOB>2.0.CO;2. Link, Google Scholar | |
| Gelb, A., J. F. Kasper, R. A. Nash, C. F. Price, and A. A. Sutherland, 1974: Applied Optimal Estimation. M.I.T. Press, 374 pp. Google Scholar | |
| Haney, R. L., 1971: Surface thermal boundary conditions for ocean general circulation models. J. Phys. Oceanogr., 1, 241–248, doi:https://doi.org/10.1175/1520-0485(1971)001<0241:STBCFO>2.0.CO;2. Link, Google Scholar | |
| Jackson, C., 2000: Sensitivity of stationary wave amplitude to regional changes in Laurentide ice sheet topography in single-layer models of the atmosphere. J. Geophys. Res., 105, 24 443–24 454, doi:https://doi.org/10.1029/2000JD900377. Crossref, Google Scholar | |
| Jazwinski, A. H., 1970: Stochastic Processes and Filtering Theory. Academic Press, 376 pp. Google Scholar | |
| Keffer, T., D. G. Martinson, and B. H. Corliss, 1988: The position of the Gulf Stream during Quaternary glaciations. Science, 241, 440–442, doi:https://doi.org/10.1126/science.241.4864.440. Crossref, Google Scholar | |
| Krauss, W., 1986: The North Atlantic Current. J. Geophys. Res., 91, 5061–5074, doi:https://doi.org/10.1029/JC091iC04p05061. Crossref, Google Scholar | |
| Labeyrie, L., and Coauthors, 1999: Temporal variability of the surface and deep waters of the North West Atlantic Ocean at orbital and millennial scales. Mechanisms of Global Climate Change at Millennial Time Scales, Geophys. Monogr., Vol. 112, Amer. Geophys. Union, 77–98, doi:https://doi.org/10.1029/GM112p0077. Google Scholar | |
| Lehman, S. J., and L. D. Keigwin, 1992: Sudden changes in North Atlantic circulation during the last deglaciation. Nature, 356, 757–762, doi:https://doi.org/10.1038/356757a0. Crossref, Google Scholar | |
| Liebelt, P. B., 1967: An Introduction to Optimal Estimation. Addison-Wesley, 273 pp. Google Scholar | |
| Locarnini, R. A., and Coauthors, 2013: Temperature. Vol. 1, World Ocean Atlas 2013, NOAA Atlas NESDIS 73, 40 pp. Google Scholar | |
| Mann, C. R., 1967: The termination of the Gulf Stream and the beginning of the North Atlantic Current. Deep-Sea Res., 14, 337–359, doi:https://doi.org/10.1016/0011-7471(67)90077-0. Crossref, Google Scholar | |
| Marchal, O., 2014: On the observability of oceanic gyres. J. Phys. Oceanogr., 44, 2498–2523, doi:https://doi.org/10.1175/JPO-D-13-0183.1. Link, Google Scholar | |
| Marchal, O., and Coauthors, 2002: Apparent long-term cooling of the sea surface in the northeast Atlantic and Mediterranean during the Holocene. Quat. Sci. Rev., 21, 455–483, doi:https://doi.org/10.1016/S0277-3791(01)00105-6. Crossref, Google Scholar | |
| MARGO, 2009: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum. Nature, 2, 127–132, doi:https://doi.org/10.1038/ngeo411. Google Scholar | |
| Matsumoto, K., and J. Lynch-Stieglitz, 2003: Persistence of Gulf Stream separation during the Last Glacial Period: Implications for current separation theories. J. Geophys. Res., 108, 3174, doi:https://doi.org/10.1029/2001JC000861. Crossref, Google Scholar | |
| Mehra, R. K., 1970: On the identification of variances and adaptive Kalman filtering. IEEE Trans. Autom. Control, 15, 175–184, doi:https://doi.org/10.1109/TAC.1970.1099422. Crossref, Google Scholar | |
| Millero, F. J., 1978: Freezing point of seawater. UNESCO Tech. Papers in Marine Science 28, 29–31. Google Scholar | |
| Mix, A., E. Bard, and R. Schneider, 2001: Environmental Processes of the Ice Age: Land, Oceans, and Glaciers (EPILOG). Quat. Sci. Rev., 20, 627–657, doi:https://doi.org/10.1016/S0277-3791(00)00145-1. Crossref, Google Scholar | |
| Monterey, G., and S. Levitus, 1997: Seasonal variability of mixed layer depth for the World Ocean. NOAA Atlas NESDIS 14, 96 pp. Google Scholar | |
| Pérez-Brunius, P., T. Rossby, and D. R. Watts, 2004: Absolute transports of mass and temperature for the North Atlantic Current–subpolar front system. J. Phys. Oceanogr., 34, 1870–1883, doi:https://doi.org/10.1175/1520-0485(2004)034<1870:ATOMAT>2.0.CO;2. Link, Google Scholar | |
| Pflaumann, U., and Coauthors, 2003: Glacial North Atlantic: Sea-surface conditions reconstructed by GLAMAP 2000. Paleoceanography, 18, 774, doi:https://doi.org/10.1029/2002PA000774. Crossref, Google Scholar | |
| Rauch, H. E., F. Tung, and C. T. Striebel, 1965: Maximum likelihood estimates of linear dynamic systems. AIAA J., 3, 1445–1450, doi:https://doi.org/10.2514/3.3166. Crossref, Google Scholar | |
| Read, J. F., R. T. Pollard, P. I. Miller, and A. C. Dale, 2010: Circulation and variability of the North Atlantic Current in the vicinity of the Mid-Atlantic Ridge. Deep-Sea Res. I, 57, 307–318, doi:https://doi.org/10.1016/j.dsr.2009.11.010. Crossref, Google Scholar | |
| Reverdin, G., P. P. Niiler, and H. Valdimarsson, 2003: North Atlantic Ocean surface currents. J. Geophys. Res., 108, 2-1–2-21, doi:https://doi.org/10.1029/2001JC001020. Crossref, Google Scholar | |
| Risien, C. M., and D. B. Chelton, 2008: A global climatology of surface wind and wind stress fields from eight years of QuickSCAT scatterometer data. J. Phys. Oceanogr., 38, 2379–2413, doi:https://doi.org/10.1175/2008JPO3881.1. Link, Google Scholar | |
| Rossby, T., 1996: The North Atlantic Current and surrounding waters: At the crossroads. Rev. Geophys., 34, 463–481, doi:https://doi.org/10.1029/96RG02214. Crossref, Google Scholar | |
| Rossby, T., and J. Nilsson, 2003: Current switching as the cause of rapid warming at the end of the Last Glacial Maximum and Younger Dryas. Geophys. Res. Lett., 30, 423, doi:https://doi.org/10.1029/2002GL015423. Crossref, Google Scholar | |
| Ruddiman, W. F., and A. McIntyre, 1981: The mode and mechanism of the last deglaciation: Oceanic evidence. Quat. Res., 16, 125–134, doi:https://doi.org/10.1016/0033-5894(81)90040-5. Google Scholar | |
| Schmitz, W. J., and M. S. McCartney, 1993: On the North Atlantic circulation. Rev. Geophys., 31, 29–49, doi:https://doi.org/10.1029/92RG02583. Crossref, Google Scholar | |
| Schulz, M., 2002: The tempo of climate change during Dansgaard-Oeschger interstadials and its potential to affect the manifestation of the 1470-year climate cycle. Geophys. Res. Lett., 29, 2-1–2-4, doi:https://doi.org/10.1029/2001GL013277. Crossref, Google Scholar | |
| Schulz, M., W. H. Berger, M. Sarnthein, and P. M. Grootes, 1999: Amplitude variations of 1470-year climate oscillations during the last 100,000 years linked to fluctuations of continental ice mass. Geophys. Res. Lett., 26, 3385–3388, doi:https://doi.org/10.1029/1999GL006069. Crossref, Google Scholar | |
| Seager, R., D. S. Battisti, J. Yin, N. Gordon, N. Naik, A. C. Clement, and M. A. Cane, 2002: Is the Gulf Stream responsible for Europe’s mild winters? Quart. J. Roy. Meteor. Soc., 128, 2563–2586, doi:https://doi.org/10.1256/qj.01.128. Crossref, Google Scholar | |
| Stommel, H. M., 1993: A conjectural regulating mechanism for determining the thermohaline structure of the ocean mixed layer. J. Phys. Oceanogr., 23, 142–150, doi:https://doi.org/10.1175/1520-0485(1993)023<0142:ACRMFD>2.0.CO;2. Link, Google Scholar | |
| Stuiver, M., and P. J. Reimer, 1993: Extended 14C data base and revised CALIB 3.0 14C age calibration program. Radiocarbon, 35, 215–230. Crossref, Google Scholar | |
| Stuiver, M., and Coauthors, 1998: Intcal98 radiocarbon age calibration, 24,000–0 cal BP. Radiocarbon, 40, 1041–1083. Crossref, Google Scholar | |
| Sy, A., U. Schauer, and J. Meincke, 1992: The North Atlantic Current and its associated hydrographic structure over and eastwards of the mid-Atlantic Ridge. Deep-Sea Res., 39, 825–853, doi:https://doi.org/10.1016/0198-0149(92)90124-C. Crossref, Google Scholar | |
| Waelbroeck, C., L. Labeyrie, J.-C. Duplessy, J. Guiot, M. Labracherie, H. Leclaire, and J. Duprat, 1998: Improving past sea surface temperature estimates based on planktonic fossil faunas. Paleoceanography, 13, 272–283, doi:https://doi.org/10.1029/98PA00071. Crossref, Google Scholar | |
| Waelbroeck, C., J.-C. Duplessy, E. Michel, L. Labeyrie, D. Paillard, and J. Duprat, 2001: The timing of the last deglaciation in North Atlantic climate records. Nature, 412, 724–727, doi:https://doi.org/10.1038/35089060. Crossref, Google Scholar | |
| Waelbroeck, C., L. Labeyrie, D. Paillard, and J. Duprat, 2014: Constraints on surface seawater oxygen isotopic change between the Last Glacial Maximum and the Late Holocene. Quat. Sci. Rev., 105, 102–111, doi:https://doi.org/10.1016/j.quascirev.2014.09.020. Crossref, Google Scholar | |
| Walter, M., and C. Mertens, 2013: Mid-depth mixing linked to North Atlantic Current variability. Geophys. Res. Lett., 40, 4869–4875, doi:https://doi.org/10.1002/grl.50936. Crossref, Google Scholar | |
| Welander, P., 1981: Mixed layers and fronts in simple ocean circulation models. J. Phys. Oceanogr., 11, 148–152, doi:https://doi.org/10.1175/1520-0485(1981)011<0148:MLAFIS>2.0.CO;2. Link, Google Scholar | |
| White, M. A., and K. J. Heywood, 1995: Seasonal and interannual changes in the North Atlantic subpolar gyre from Geosat and TOPEX/POSEIDON altimetry. J. Geophys. Res., 100, 24 931–24 941, doi:https://doi.org/10.1029/95JC02123. Crossref, Google Scholar | |
| Wunsch, C., 2006a: Abrupt climate change: An alternative view. Quat. Res., 65, 191–203, doi:https://doi.org/10.1016/j.yqres.2005.10.006. Crossref, Google Scholar | |
| Wunsch, C., 2006b: Discrete Inverse and State Estimation Problems. Cambridge University Press, 371 pp. Google Scholar | |
| Wunsch, C., 2011: The decadal mean ocean circulation and Sverdrup balance. J. Mar. Res., 69, 417–434, doi:https://doi.org/10.1357/002224011798765303. Crossref, Google Scholar | |
| Zahn, R., 1994: Core correlations. Nature, 371, 289–290, doi:https://doi.org/10.1038/371289a0. Google Scholar | |
| Zhang, X., G. Lohmann, and C. Purcell, 2014: Abrupt glacial climate shifts controlled by ice sheet changes. Nature, 512, 290–294, doi:https://doi.org/10.1038/nature13592. Crossref, Google Scholar | |
| Zweng, M. M., and Coauthors, 2013: Salinity. Vol. 2, World Ocean Atlas 2013, NOAA Atlas NESDIS 74, 39 pp. Google Scholar |
1 Data are available upon request to the corresponding author.
2 The computer code of the model and a manual including details about model derivation are available at https://www.ncdc.noaa.gov/paleo/study/19300.
D1 The derivatives are listed in the manual at https://www.ncdc.noaa.gov/paleo/study/19300.
The model consists of the governing equation for ML temperature (2), together with the expressions for the Ekman velocity (4), the geostrophic velocity (7) with its thermal (8) and saline contributions (9), and the interior vertical velocity 
(10). The boundary condition is provided by (11).
The model equations are solved using a finite-difference method. The temporal and spatial derivatives in (2) and (8)–(11) are approximated with finite differences derived for a staggered grid with uniform longitude and latitude spacings 
(Fig. 1). A time step 
yr is used for all model integrations.
The grid is arranged as follows. To lighten notation, the bar over the velocities that is used in the text to denote vertical averages is omitted [the notation 
is used for another purpose in this appendix]. The zonal (meridional) component of 
is noted by 
, and the zonal (meridional) component of surface wind stress is designated by 
. The values of 
are carried by grid points at the center of 2° × 2° cells. The values of 
are carried by points along the meridional boundaries of the cells and at the same latitudes as the 
-carrier points. Conversely, the values of 
are carried by points along the zonal boundaries of the cells and at the same longitudes as the 
-carrier points.
With this arrangement of the grid, the finite-difference approximations of (2), (4), and (8)–(11) are as follows. Consider first (2), the equation for ML temperature. The advection term, 
, is expressed as
vanishes in the difference 
. Equation (2) is then discretized using an upstream advection scheme:
stands for 
, and r is Earth’s radius. The advective fluxes in (A4a)–(A4b) are evaluated at locations halfway between the T-carrier points:
. Note that in (A6) all the variables are implicitly defined at the time level n. The same convention is adopted below, unless specified otherwise.Consider (4), the Ekman velocity in the mixed layer. The Ekman velocity components are obtained from
Consider (8), the thermal contribution to the geostrophic velocity. The zonal and meridional components of this contribution are computed from, respectively,
in (A9b) is given by 
. The zonal and meridional components of the saline contribution to the geostrophic velocity (9) are computed from similar equations.The vertical interior velocity (10) is computed from
is governed by
The apparent atmospheric temperatures 
are obtained from the annual mean (statistical means, averaged decades) SSTs (
m) with 1° × 1° resolution available in Locarnini et al. (2013). The SST values from Locarnini et al. (2013) are averaged in 2° × 2° cells centered at the T-carrier points of the grid to produce the 
values. The interior oceanic temperatures 
are derived from the same procedure, except that these temperatures are reduced by 0.5°C, consistent with our criterion for ML depth (see below). The values of salinity S are obtained from the annual mean (statistical means, averaged decades) surface salinities (
m) with 1° × 1° resolution reported in Zweng et al. (2013). The salinity values from Zweng et al. (2013) are averaged in 2° × 2° cells centered at the S-carrier points of the model grid to produce the S values. The ML depths h are derived from a climatology of monthly mean ML depths based on a 0.5°C potential temperature criterion and with 1° × 1° resolution (Monterey and Levitus 1997). From this climatology, a distribution of annual mean ML depth with the same resolution is computed. The annual mean ML depths are then averaged in 2° × 2° cells centered at the h-carrier points to produce the h values. Finally, the surface wind stresses are obtained from the Scatterometer Climatology of Ocean Winds (Risien and Chelton 2008). The monthly averages of zonal and meridional wind stresses, available at ¼° × ¼° resolution, are averaged to produce annual means. The annual mean values of zonal (meridional) stress are then averaged in 2° × 2° cells centered at the 
- (
-) carrier points.
Consider the estimation of 
, with 
. Equations (12) lead to systems of linear albegraic equations, one system for each variable and one equation for each pair 
). For instance, the system for the variable 
can be written as
, 
is a coefficient matrix, 
is a vector of modern data, and 
is a vector that includes the errors in these data. The data in this case are observational estimates of the apparent atmospheric temperature. They are obtained by averaging in 2° × 2° model grid cells the annual mean SSTs with 1° × 1° resolution reported in Locarnini et al. (2013) (appendix B of this paper). Their errors are calculated by propagating the standard errors in the annual mean SSTs contributing to the gridcell averages, neglecting error correlation. The squares of these errors are included along the diagonal of a square matrix 
. The off-diagonal elements of 
are set to zero. An estimate of 
that minimizes the quadratic form 
is then sought. This weighted least squares estimate and its uncertainty are given by (e.g., Wunsch 2006b)
are estimated along the same lines. For 
, the data in 
are as described in appendix B, and their errors represented in 
are the same as for 
. For 
, the data in 
are as described in appendix B, and their errors represented in 
are assumed to be 10 m. Finally, for 
and 
, the data in 
are zonal and meridional velocity components, 
and 
, computed from the wind stress and salinity data described in appendix B, using the numerical scheme detailed in appendix A. Their errors represented in 
are set equal to 0.1 cm s−1. As for 
, the off-diagonal elements of 
are set to zero for the estimation of 
.In this appendix, a brief overview of the EKF and LKF is provided; textbooks should be consulted for real understanding (e.g., Jazwinski 1970; Gelb et al. 1974; Bryson and Ho 1975; Anderson and Moore 1979).
In the EKF, the equations in 
are linearized about the more recent estimate of the state, 
(+), in the calculation of state error propagation:
is a matrix of partial derivatives evaluated at 
.D1 Like other versions of the Kalman filter, the EKF consists of two sets of equations. In the first set (the time-update equations), the model is used to extrapolate the state and its uncertainty from time 
to time i:
, whereas an unstable system will experience unbounded error growth in the absence of observations (Gelb et al. 1974). The second term on the right-hand side represents the effect of model errors. For example, for 
in (19), the error in ML temperature would tend to increase by an amount equal to 
in one time step in the absence of observations. With a time step of 0.1 yr−1 [the value used for the model integrations (appendix A)], the ML temperature errors tend to grow by an amount equal to 
in 100 yr, should no observations be available over this interval.In the second set (the data-update equations), the extrapolated state and its uncertainty are modified in order to account for the presence of observations [if no observation is available, then the state and its uncertainty keep evolving according to (D2a)–(D2b)]. A state estimate is sought that minimizes the objective function:
, and the second represents the deviation from the observations 
, weighted by the error covariances 
and 
, respectively. The resulting state estimate and its uncertainty are given by the data-update equations:
and 
(e.g., Bryson and Ho 1975). Notice that (D4b) is used here rather than the simpler and mathematically equivalent form 
, as it tends to promote nonnegativity of 
(e.g., Anderson and Moore 1979).The EKF has been found to yield accurate estimates in number of applications and reduces to the conventional Kalman filter when the system dynamics and the observations are linear. For these reasons, the EKF is usually one of the first methods to be tried for any nonlinear filtering problem (e.g., Gelb et al. 1974). On the other hand, the EKF may be prone to instabilities (e.g., Wunsch 2006b). In this study, another version of the Kalman filter that can deal with nonlinearities, the LKF, is considered as an alternative method to combine the SST records with the ocean model.
In the LKF, the state vector is represented as the sum of a nominal or reference state 
and a perturbation 
:
The time-update equations of the LKF are
, which has the same form as (D4b). In (D9), the gain 
is computed from (D5), and 
is the vector difference 
.In this study, the reference state 
is taken as time independent and is obtained from modern observations. Specifically, 
(section 4d), so that 
. A more plausible reference state would require prior knowledge about the past distributions of SST, ML depth, surface wind stress, etc., which is not available. In general, the LKF is not expected to yield accurate state estimates if the state drifts very far from the reference state. In such a situation, the EKF can lead to more accurate results.
We compare the ML temperature estimates obtained from the EKF and LKF near the sediment core locations for the time interval between 10 and 14.5 kyr BP. This interval is characterized by the largest temperature changes over the past 14.5 kyr BP in the sediment records (Fig. 2). It is therefore the interval when the largest differences between the results from both filters can reasonably be expected, given the varying assumptions in both methods. To ensure comparability of the results, the equations of the EKF and LKF are solved with the same assumptions about 
(section 4c) and 
(
; section 5a).
We find that the ML temperature estimates derived from the EKF and LKF show small differences: the estimates derived from the extended filter differ from those derived from the linearized filter by less than one standard deviation of either the EKF or LKF estimates (Figs. D1, D2). Likewise, the standard deviations in the ML temperature estimates obtained from both methods are generally very close (Fig. D2).
|
|
The smoother equations used in this study are based on an algorithm originally developed for linear systems and reported in Bryson and Ho (1975, p. 392) (see also Fraser 1967). An estimate of the state perturbation 
is sought that minimizes
. The resulting smoothing estimate is given by
and with 
. The first term on the right-hand side of (E2) is the state perturbation estimated by the LKF (with 
; section 5a). The second is a correction proportional to a vector 
driven by posterior data according to the recursion equation (E3). The state estimate (E2) is thus constrained by the entire set of observations, from the initial time (
) to the terminal time (
). As for the Kalman filter estimates, the smoother estimates can be regarded as weighted least squares estimates, but with the weighting now provided by 
, 
, and 
over the entire interval (e.g., Bryson and Ho 1975).The uncertainty in the smoothing estimate (E2) is given by
. The first term on the right-hand side of (E4) is the state uncertainty obtained from the LKF. The second is always positive definite, indicating that the errors in the state estimates generally decrease as posterior data are considered in the analysis (more precisely, these errors never increase). It depends on the matrix 
, which accounts for the errors in the posterior data represented in 
and which is calculated by recursion (E5).The algorithm (E2)–(E5) is equivalent to the Rauch–Tung–Striebel (RTS) algorithm (Rauch et al. 1965), a conventional method for fixed-interval smoothing (Gelb et al. 1974; Bryson and Ho 1975; Anderson and Moore 1979). The two algorithms differ in the matrix inversions: whereas the RTS requires the calculation of 
, (E2)–(E5) rely on 
. Since in the present study 
is generally of much smaller order than 
, (E2)–(E5) offer a clear numerical advantage.
