Abstract

We use the method of least squares with Lagrange multipliers to fit an ocean general circulation model to the Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) estimate of near sea surface temperature (NSST) at the Last Glacial Maximum (LGM; circa 23–19 thousand years ago). Compared to a modern simulation, the resulting global, last-glacial ocean state estimate, which fits the MARGO data within uncertainties in a free-running coupled ocean–sea ice simulation, has global-mean NSSTs that are 2°C lower and greater sea ice extent in all seasons in both the Northern and Southern Hemispheres. Increased brine rejection by sea ice formation in the Southern Ocean contributes to a stronger abyssal stratification set principally by salinity, qualitatively consistent with pore fluid measurements. The upper cell of the glacial Atlantic overturning circulation is deeper and stronger. Dye release experiments show similar distributions of Southern Ocean source waters in the glacial and modern western Atlantic, suggesting that LGM NSST data do not require a major reorganization of abyssal water masses. Outstanding challenges in reconstructing LGM ocean conditions include reducing effects from model biases and finding computationally efficient ways to incorporate abyssal tracers in global circulation inversions. Progress will be aided by the development of coupled ocean–atmosphere–ice inverse models, by improving high-latitude model processes that connect the upper and abyssal oceans, and by the collection of additional paleoclimate observations.

1. Introduction

Disagreements among general circulation model (GCM) representations of the Last Glacial Maximum [LGM; circa 23–19 thousand years ago (ka); Mix et al. (2001)] and between models and LGM paleoceanographic data (Braconnot et al. 2007; Otto-Bliesner et al. 2009; Tao et al. 2013; Dail and Wunsch 2014, hereafter DW14) illustrate a gap in our knowledge of Earth’s climate during that time period. Here we present a global ocean state estimate at the LGM, a dynamically consistent fit of an ocean general circulation model (OGCM) to surface ocean temperature proxies achieved by adjusting model initial conditions, boundary conditions, and turbulent transport parameters. This work builds on a growing body of literature combining dynamical models with proxy observations in order to interpolate between LGM observations (e.g., Winguth et al. 2013; Kurahashi-Nakamura et al. 2014; DW14; Kurahashi-Nakamura et al. 2017, hereafter KN17).

Several factors motivate studying the climate of the LGM. First, geologic evidence suggests that LGM conditions were a persistent and dramatic excursion from the present-day climate, with large ice sheets in the Northern Hemisphere, lower sea levels, and a global-mean surface air temperature reduction of several degrees Celsius (Clark et al. 2012). Second, radiocarbon dating allows measurements to be reliably placed within the LGM time frame. Finally, the LGM is a useful period to study the ocean’s role in regulating atmospheric carbon dioxide concentrations, with implications for understanding modern climate change (Sarmiento and Toggweiler 1984; Siegenthaler and Wenk 1984; Brovkin et al. 2007; Shakun et al. 2012), including the sensitivity of climate to atmospheric greenhouse gas concentrations (Schmittner et al. 2011; Hargreaves et al. 2012).

The Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) compilation of LGM surface ocean temperature estimates (Waelbroeck et al. 2009) extends the previous work of Glacial Atlantic Ocean Mapping (GLAMAP; Pflaumann et al. 2003) and Climate: Long-range Investigation, Mapping, and Prediction (CLIMAP; McIntyre et al. 1976) by including more observations from a wider range of temperature proxies. We refer to these data as representing “near” sea surface temperature (NSST) in recognition of the various depth ranges inhabited by organisms used for temperature reconstructions.

Numerous studies have used the MARGO database as a basis for comparison with numerical models, often showing qualitative disagreements on regional scales. Simulations from the Paleoclimate Modelling Intercomparison Projects (PMIP1, PMIP2, and PMIP3) used LGM boundary conditions, including global sea level, orography, greenhouse gases, and Earth’s orbital parameters (Braconnot et al. 2007), in climate models of varying complexity. Hargreaves et al. (2011) found that the intermodel spread of simulated NSSTs in PMIP1 and PMIP2 did not disagree with MARGO data within its uncertainty. However, DW14 found that, when considered individually, five PMIP2 simulations fit MARGO data poorly in the North Atlantic. In the tropical oceans, Otto-Bliesner et al. (2009) found that PMIP2 models had a similar range of global-mean NSST decrease to that estimated by MARGO and larger cooling in the Atlantic than in the Pacific, also in agreement with the observations, but that zonal gradients of LGM cooling in tropical Pacific near-surface waters were less pronounced than in MARGO. Model ensemble averages reported by Braconnot et al. (2007) and individual model results from Tao et al. (2013) show North Atlantic cooling patterns with a zonal gradient opposite that seen in the data. Data errors contributing to these disagreements could arise from chronological errors, seasonal biases, and biological effects, to name a few. Model errors include incorrectly specified initial and boundary conditions, errors in numerical solution methods, missing physics, and inaccurate parameterizations of unresolved phenomena (e.g., ocean eddies and clouds).

The Atlantic abyssal circulation may have played an important role in maintaining a climate at the LGM that was different from the modern through its role in transporting and storing heat, biological nutrients, and carbon. For instance, one interpretation of paleoceanographic data from the Atlantic is that during the LGM, deep water originating from the North Atlantic shoaled and bottom water from the Southern Ocean filled more of the abyss (e.g., Curry et al. 1988; Duplessy et al. 1988; Marchitto et al. 2002; Curry and Oppo 2005; Marchitto and Broecker 2006; Lynch-Stieglitz et al. 2007), possibly coincident with a weakening and shoaling of the upper cell of the Atlantic meridional overturning circulation (AMOC). This scenario is simulated in some climate models, but not all. While PMIP2 LGM experiments showed a broad range of strengths and depths of the upper and lower cells of the AMOC (Otto-Bliesner et al. 2007), nearly all PMIP3 simulations show deeper and stronger upper-cell AMOC transport at the LGM relative to modern simulations (Muglia and Schmittner 2015). By contrast, simplified ocean models considered by Ferrari et al. (2014) and Jansen and Nadeau (2016) point to a shallower, weaker LGM upper cell. Differences among models may arise from different model architectures, spatial resolution, bathymetry, physical parameterizations, or incomplete equilibration with surface conditions (Zhang et al. 2013; Marzocchi and Jansen 2017). Finally, estimates of LGM salinity derived from the pore fluids of sediment cores suggest that the deep ocean was not only saltier, because of the storage of freshwater in ice sheets, but also more salt stratified (Adkins et al. 2002; Insua et al. 2014). However, Miller et al. (2015) and Wunsch (2016) argue that pore fluid measurements are too few to be uniquely interpretable.

Fitting models to paleoceanographic data can improve our knowledge of model and data shortcomings. Ultimately, this approach can improve our knowledge of the ocean circulation and climate at time intervals like the LGM. Previous efforts include Dail (2012) and DW14, who obtained a state estimate of the LGM Atlantic Ocean by fitting an OGCM to Atlantic MARGO data, and KN17, who fit an OGCM to the global annual-mean MARGO data as well as oxygen and carbon isotope data in the Atlantic Ocean. Other efforts to constrain the abyssal circulation during the LGM by combining models and proxy data include LeGrand and Wunsch (1995), Gebbie and Huybers (2006), Marchal and Curry (2008), Burke et al. (2011), Gebbie (2014), and Gebbie et al. (2016). A common conclusion of these studies is the difficulty in determining past circulations uniquely because of the sparsity and noisiness of paleoceanographic measurements.

Here we present a new fit of an OGCM to the MARGO dataset of global, seasonal gridded NSST observations. This work expands upon DW14 by using (i) a global domain, (ii) a longer model integration, and (iii) atmospheric forcing derived in part from a coupled ocean–atmosphere model LGM simulation. Differences from KN17 include (i) the use of higher spatial resolution both horizontally and vertically (2° vs 3° and 50 vs 15 vertical levels, respectively) and (ii) the inclusion of seasonal MARGO data. Our state estimate is a primitive equation ocean model simulation that agrees with seasonal MARGO data within their estimated errors and allows us to analyze approximately equilibrated properties of the ocean circulation, including in the abyss. Unlike KN17, we do not use oxygen and carbon isotope data in the deep ocean to constrain the state estimate, as simulation durations required to equilibrate abyssal tracer distributions (thousands of model years) proved too computationally expensive for our state estimation framework, and fitting incompletely equilibrated model tracers to observations can lead to biased solutions (Dail 2012; Amrhein 2016). A disadvantage of this approach is that our inferences of abyssal circulation and structure are not informed by in situ measurements. A comparison of LGM state estimates in the discussion (section 4) provides insights into their uncertainties and sensitivities to different state estimation approaches.

2. Materials and methods

a. LGM NSST data

NSST data and uncertainties used in this study are from the 5° × 5° MARGO gridded products (Waelbroeck et al. 2009) constructed from microfossil and chemical measurements in ocean sediment cores representing the time interval 23–19 ka BP. The MARGO compilation includes transfer function approaches—which match past abundances of planktonic foraminifera, diatoms, dinoflagellate cysts, or radiolarians to modern analogs—and chemical thermometers based on alkenone indices and planktonic foraminiferal Mg/Ca. Gridded values are weighted means of proxy values, with weights based on data type, numbers of observations available during the time period, and calibration and instrumental errors. Three separate gridded MARGO products represent annual, January–March (JFM), and July–September (JAS) mean conditions. The spatial density of the gridded data is highest in tropical regions and at high northern latitudes, especially in the northern North Atlantic and Arctic Oceans. Data from the Southern Ocean are restricted to austral summer because of the limited seasonal representativeness of diatom assemblages, which make up most available observations in that region.

b. The MITgcm

The OGCM we fit to the MARGO data is the MITgcm, an evolved form of that described by Marshall et al. (1997) and Adcroft et al. (2004) that simulates the ocean circulation under hydrostatic and Boussinesq approximations. The model is a lower-resolution configuration of the ECCO, version 4, release 2, modern state estimation setup (ECCO; Forget et al. 2015a), with 2° horizontal resolution telescoping to higher resolution at the equator and the poles and 50 vertical levels with thicknesses ranging from a minimum of 10 m at the surface to a maximum of 456 m at depth. The MITgcm is coupled to a viscous plastic dynamic–thermodynamic sea ice model (Losch et al. 2010). Air–sea fluxes of heat, freshwater, and momentum are computed using the bulk formulae of Large and Yeager (2004). Global-mean freshwater fluxes through the sea surface are compensated at every time step by adding or subtracting a uniform freshwater flux correction that prevents drifts in global mean ocean salinity. Ocean vertical mixing is parameterized using the turbulent closure scheme of Gaspar et al. (1990). Isopycnal diffusivity is treated using the Redi (1982) scheme, and unresolved eddy advection is parameterized using the method of Gent and McWilliams (1990). Following Bugnion and Hill (2006) and Dail (2012) we use accelerated time stepping (Bryan 1984), with a tracer time step of 12 hours and a momentum time step of 20 min.

Model bathymetry for the LGM was constructed by smoothing and subsampling modern water depth estimates (Smith and Sandwell 1997) and adding the LGM minus modern bathymetry anomaly reconstructed by Peltier (2004), which has a median LGM sea level of approximately 130 m below present. A seasonal cycle of runoff is derived from Fekete et al. (2002), with runoff on the European continent between 50° and 72°N rerouted to the latitude of the English Channel, reflecting the reconstruction of Alkama et al. (2006). Sea ice and snow albedos were reduced by roughly 30% from ECCO values to prevent unrealistic sea ice growth in the LGM state estimate.

c. State estimation procedure

Procedures for obtaining data-constrained ocean state estimates used in this paper are illustrated in Fig. 1. We use the method of least squares with Lagrange multipliers (also known as the adjoint method; e.g., Wunsch 2006) to fit the MITgcm to seasonal- and annual-mean MARGO NSST data. In modern oceanography, the relative wealth of observations permits estimating the time-varying ocean state (Stammer et al. 2002; Wunsch and Heimbach 2007; Forget et al. 2015a). At the LGM, the sparsity of data motivates treating them as samples of a “seasonally steady” state—a single seasonal cycle that repeats over the interval 23–19 ka. Our goal is to generate an MITgcm simulation under annually repeating atmospheric boundary conditions that both fits the data within their uncertainties and is consistent with a quasi-steady circulation, as defined below. We will denote vectors and matrices by lower- and uppercase bold letters, respectively.

Fig. 1.

Flowchart describing the construction of simulated ocean states (MODERN and PRIOR) and data-constrained glacial state estimates (GLACIAL and GLACIAL_s) described in this study. We use 5000-yr-long simulations to equilibrate the MITgcm to different sets of model control variables (or “controls”: atmospheric conditions, turbulent transport parameters, and initial conditions). MODERN (section 3a) gives a model representation of modern oceanographic conditions; PRIOR is a first-guess glacial state that is the starting point for LGM state estimation (section 2d); GLACIAL is the LGM state estimate that is the main result of this paper; and GLACIAL_s is used to diagnose sensitivity of the state estimate to a uniform adjustment in surface air temperature (see discussion in section 4). Modern initial conditions used for MODERN and PRIOR are taken from the ECCO state estimate (Forget et al. 2015a) from the year 2007; an additional 1.1 salinity is added to every grid box in PRIOR based on the global-mean salinity change estimated at the LGM by Adkins et al. (2002). The CCSM4 refers to differences between LGM and preindustrial coupled CCSM4 simulations.

Fig. 1.

Flowchart describing the construction of simulated ocean states (MODERN and PRIOR) and data-constrained glacial state estimates (GLACIAL and GLACIAL_s) described in this study. We use 5000-yr-long simulations to equilibrate the MITgcm to different sets of model control variables (or “controls”: atmospheric conditions, turbulent transport parameters, and initial conditions). MODERN (section 3a) gives a model representation of modern oceanographic conditions; PRIOR is a first-guess glacial state that is the starting point for LGM state estimation (section 2d); GLACIAL is the LGM state estimate that is the main result of this paper; and GLACIAL_s is used to diagnose sensitivity of the state estimate to a uniform adjustment in surface air temperature (see discussion in section 4). Modern initial conditions used for MODERN and PRIOR are taken from the ECCO state estimate (Forget et al. 2015a) from the year 2007; an additional 1.1 salinity is added to every grid box in PRIOR based on the global-mean salinity change estimated at the LGM by Adkins et al. (2002). The CCSM4 refers to differences between LGM and preindustrial coupled CCSM4 simulations.

The ocean state vector at a time t, , is a complete list of the variables required to take one model time step—temperature, salinity, velocity, etc.—at all locations of the model grid. An underbar denotes a vector containing a seasonal cycle of values; for example, is a list of model variable values concatenated in time over a year. The evolution of the MITgcm under seasonally steady forcing can be written as

 
formula

where L is a nonlinear operator, is the model time step, M is a positive integer, is a vector of model parameters that are not changed in the optimization (e.g., model bathymetry), and is a vector of adjustable “control” variables (or “controls”) including fields of initial temperature and salinity, turbulent transport parameters, and monthly average atmospheric forcing (Table 1).

Table 1.

Control variables, control uncertainty standard deviations σ, sources of first-guess control values, and time periods for control variables used to derive the GLACIAL and GLACIAL_s state estimates. ECCO refers to the ECCO modern ocean state estimate in the year 2007. The CCSM4 refers to differences between LGM and preindustrial coupled CCSM4 simulations (Fig. 2). The variables , , and refer to coefficients of isopycnal diffusivity (Redi 1982), diapycnal diffusivity, and eddy diffusivity associated with the bolus velocity (Gent and McWilliams 1990). PRIOR refers to the forward simulation under ECCO + CCSM4 forcing described in section 2d.

Control variables, control uncertainty standard deviations σ, sources of first-guess control values, and time periods for control variables used to derive the GLACIAL and GLACIAL_s state estimates. ECCO refers to the ECCO modern ocean state estimate in the year 2007. The CCSM4 refers to differences between LGM and preindustrial coupled CCSM4 simulations (Fig. 2). The variables , , and  refer to coefficients of isopycnal diffusivity (Redi 1982), diapycnal diffusivity, and eddy diffusivity associated with the bolus velocity (Gent and McWilliams 1990). PRIOR refers to the forward simulation under ECCO + CCSM4 forcing described in section 2d.
Control variables, control uncertainty standard deviations σ, sources of first-guess control values, and time periods for control variables used to derive the GLACIAL and GLACIAL_s state estimates. ECCO refers to the ECCO modern ocean state estimate in the year 2007. The CCSM4 refers to differences between LGM and preindustrial coupled CCSM4 simulations (Fig. 2). The variables , , and  refer to coefficients of isopycnal diffusivity (Redi 1982), diapycnal diffusivity, and eddy diffusivity associated with the bolus velocity (Gent and McWilliams 1990). PRIOR refers to the forward simulation under ECCO + CCSM4 forcing described in section 2d.

We fit the MITgcm to the MARGO data by iteratively adjusting control variables to minimize a cost function with three terms. The first term penalizes misfits between the model and data, the second penalizes large changes to the controls, and the last imposes the dynamical constraints of the model using the Lagrange multipliers. At each iteration, the model is run forward for 100 years, cost terms are computed, and the model adjoint is run backward in time to estimate the linear sensitivity of the cost function to the controls (see the  appendix). Then control adjustments are made, the model is run again, costs are recomputed, and the cycle is repeated until an acceptable fit is found. The model–data misfit cost function term is computed over the last 20 years of the 100-yr-long forward simulations to permit the model to adjust to control changes. A period of 100 years is long enough to bring much of the surface ocean into near equilibrium with changes in atmospheric conditions but too short to equilibrate the deep ocean (Wunsch and Heimbach 2008). Thus, our results are biased against dynamical mechanisms that could reduce model–data misfits on time scales longer than a century.

After deriving controls using 100-yr-long adjoint simulations, we integrate the model for 5000 years under derived control adjustments and take the last year of the integration to be our state estimate. We perform this additional integration in order to (i) allow the abyssal ocean approximately to equilibrate to changes derived to fit NSST data and (ii) evaluate the consistency of our solution with a quasi-steady circulation. The simulation is said to be adequately steady if it fits the MARGO data within uncertainty near the beginning and end of the simulation at model years 80–100 and 4980–5000, that is, if it satisfies the sets of equations

 
formula
 
formula

Here and are the simulated seasonal cycles averaged over model years 80–100 and 4980–5000, respectively; is a matrix relating MARGO NSSTs to and ; and and are residuals to the model fit that must be consistent with magnitudes of observational errors.

d. Control variables, error covariances, and a first-guess solution

State estimation requires specifying first-guess control values that are subsequently adjusted to fit data. First guesses of atmospheric controls (Table 1) are the sums of modern ECCO fields (Forget et al. 2015a) and LGM minus preindustrial anomalies computed in the Community Climate System Model, version 4 (CCSM4; these anomalies are referred to below as CCSM4). We choose to add CCSM4 to modern ECCO fields rather than simply using CCSM4 LGM fields in an effort to mitigate potential biases from CCSM4. The CCSM4 consists of coupled ocean, atmosphere, land, and sea ice models with nominal 1° horizontal resolution. The preindustrial (PI) CCSM4 simulation (Gent et al. 2011) follows protocols for phase 5 of the Climate Model Intercomparison Project (CMIP5), while the LGM CCSM4 simulation (Brady et al. 2013) follows PMIP3 protocols, using LGM orbital parameters, greenhouse gas concentrations estimated from ice cores, modified orography due to Northern Hemisphere ice sheets, and reduced global sea level.

The CCSM4 wind stress anomalies reflect orographic changes due to the presence of Northern Hemisphere ice sheets (Brady et al. 2013; Figs. 2a,b). Surface air temperatures are everywhere reduced in the CCSM4 LGM simulation relative to the preindustrial, with especially pronounced cooling in the subpolar North Atlantic, Southern, and North Pacific Oceans (Fig. 2e). Downwelling longwave radiation (Fig. 2k) and humidity (Fig. 2i) are also lower everywhere at the LGM, likely reflecting changes in atmospheric heat content and the reduced capacity of colder air to hold moisture. Anomalies of precipitation (Fig. 2g) and shortwave downwelling radiation (Fig. 2m) show more complex patterns, possibly reflecting differences in simulated atmospheric circulation and cloud distributions as well as changes in Earth’s orbital configuration. In many regions these anomalies have the same order of magnitude as modern time-mean values.

Fig. 2.

(left) Last Glacial Maximum minus preindustrial annual-mean anomalies (ΔCCSM4) of atmospheric variables in CCSM4 (Brady et al. 2013). (right) Annual-mean adjustments to atmospheric control variables derived to fit the MITgcm to MARGO data. Panel (f) includes a uniform change in global mean surface air temperature made to fit the data.

Fig. 2.

(left) Last Glacial Maximum minus preindustrial annual-mean anomalies (ΔCCSM4) of atmospheric variables in CCSM4 (Brady et al. 2013). (right) Annual-mean adjustments to atmospheric control variables derived to fit the MITgcm to MARGO data. Panel (f) includes a uniform change in global mean surface air temperature made to fit the data.

First guesses of glacial distributions of ocean temperature and salinity are taken from a 5000-yr-long simulation of the MITgcm (PRIOR) forced by first-guess atmospheric conditions (Table 1). PRIOR is initialized with temperatures and salinities from the ECCO modern ocean state in the year 2007 plus an additional 1.1 salinity at every model grid box, based on the global-mean salinity change estimated at the LGM from changes in global-mean sea level (Fairbanks 1989; Adkins et al. 2002). The year 2007 was chosen based on the availability of modern observations; the fact that 2007 was an El Niño year may contribute to zonal Pacific temperature gradients observed in patterns of model drift. More generally, though we do not attempt to estimate it here, sensitivity to choices of first-guess conditions is an important contributor to solution uncertainty and should be prioritized in future uncertainty quantification studies.

Finally, we must assume values for the standard deviations σ of uncertainties in our choices of first-guess control variables. Following DW14, σ for shortwave and longwave downwelling radiation, humidity, and precipitation are twice those used in ECCO, and σ for surface atmospheric temperature is 4 times that in ECCO. Wind stress σ is set to 0.1 Pa, reflecting the amplitudes of CCSM4 LGM-PI wind stress changes. For initial salinity, σ is 1 on the practical salinity scale, comparable to the estimated change in ocean mean salinity over the last deglaciation. Errors for turbulent transport parameters are taken from ECCO. We assume that control variable uncertainties do not covary in space or between variables.

3. Results

This section reports results from fitting the MITgcm to MARGO LGM NSST estimates and describes properties of our best-estimate LGM ocean state (GLACIAL). We also describe the modern simulation (MODERN) used to compare to GLACIAL. Importantly, we do not claim that our state estimate is a unique fit to the data; other ocean states may exist that are qualitatively different but fit the data equally well. In particular, the abyssal ocean appears at best to be weakly constrained by the MARGO data (Kurahashi-Nakamura et al. 2014).

a. Construction of the MODERN simulation and comparison to the modern ocean

The MODERN simulation is a 5000-yr integration of the MITgcm using modern bathymetry and atmospheric conditions (Fig. 1). This simulation is used to compute LGM-modern anomalies and to identify model biases. In particular, after 5000 years of integration, annual-mean surface values (centered on 5-m water depth) of temperature and salinity in MODERN show regional deviations from modern ECCO values, which are constrained by modern observations, of over 4°C and 2, respectively (Figs. 3a,b). We attribute these differences, which accumulate on time scales of centuries to millennia (see Fig. S7 in the online supplemental material), primarily to model errors, including model “drift,” which is a common phenomenon in ocean-only models lacking atmosphere–ocean feedbacks (e.g., Griffies et al. 2009). The absence of these patterns in PRIOR–MODERN anomalies (Figs. 3c,d) indicates that similar biases are also present in PRIOR and as such will appear in the GLACIAL state estimate unless they are eliminated by fitting the model to MARGO NSST data. In addition to differences in surface water properties, the MODERN AMOC has a weaker and shallower upper cell than in modern observationally based reconstructions (Lumpkin and Speer 2007). In this study, we use MODERN as the basis for comparison with GLACIAL—rather than a modern state estimate or modern observations—because taking the difference between the two time intervals is likely to reduce the impacts of systematic model errors on our conclusions. A caveat is that wherever fitting the data adjusts GLACIAL closer to the true LGM state, common biases in GLACIAL and MODERN may no longer cancel when their anomaly is computed, leading to errors in inferred anomalies that may be as large as the model bias.

Fig. 3.

Comparison of annual-mean surface temperature and salinity at the end of the 5000-yr-long simulations MODERN and PRIOR, in the data-constrained state estimate GLACIAL, and computed as the 20-yr time average of the modern state estimate ECCO (Forget et al. 2015a). “Surface” is here taken to mean the uppermost model grid box, which is centered on 5-m water depth. (a),(b) MODERN surface temperature and salinity show basin-scale deviations from ECCO attributed to model bias due to coarse resolution and a long integration period (model drift). Because similar model biases are also present in PRIOR, the patterns evident in (a) and (b) largely cancel in (c),(d) PRIOR − MODERN. PRIOR − MODERN anomalies on the practical salinity scale in (d) have been corrected by subtracting 1.1 to account for the mean salinity increase imposed in PRIOR. (e),(f) GLACIAL − PRIOR differences show changes that result from fitting the model to the MARGO data.

Fig. 3.

Comparison of annual-mean surface temperature and salinity at the end of the 5000-yr-long simulations MODERN and PRIOR, in the data-constrained state estimate GLACIAL, and computed as the 20-yr time average of the modern state estimate ECCO (Forget et al. 2015a). “Surface” is here taken to mean the uppermost model grid box, which is centered on 5-m water depth. (a),(b) MODERN surface temperature and salinity show basin-scale deviations from ECCO attributed to model bias due to coarse resolution and a long integration period (model drift). Because similar model biases are also present in PRIOR, the patterns evident in (a) and (b) largely cancel in (c),(d) PRIOR − MODERN. PRIOR − MODERN anomalies on the practical salinity scale in (d) have been corrected by subtracting 1.1 to account for the mean salinity increase imposed in PRIOR. (e),(f) GLACIAL − PRIOR differences show changes that result from fitting the model to the MARGO data.

b. Fitting the model to data

A state estimate is considered to fit data adequately when model–data misfits normalized by observational errors have an approximately Gaussian distribution with mean 0 and standard deviation 1. By this criterion, the first-guess PRIOR simulation does not fit the MARGO data: in the annual, JAS, and JFM means, standard deviations of normalized misfits are greater than 1 (Fig. 4a). Moreover, the average value of normalized misfits is less than 0, indicating a model cold bias relative to the data. Misfits exceeding observational uncertainties are found in several regions. In both JAS and JFM, the model is warm relative to the data in the equatorial Atlantic, the northeast Atlantic, and the western Pacific, while it is too cold in the Indian, Arctic, and east Pacific Oceans (Figs. 4d,f). In JFM, the model Southern Ocean is cold relative to the data. Similarities between spatial patterns of model–data misfit and MODERN − ECCO temperature anomalies suggest that model bias is a major contributor to model–data misfit.

Fig. 4.

Model–data misfits of NSSTs are improved by the state estimation procedure. (a) Histograms of model–data misfits for annual-, JFM-, and JAS-mean values in PRIOR and GLACIAL normalized by data uncertainties. GLACIAL misfits are similar to a standard normal Gaussian distribution (gray), indicating an acceptable fit to the data. Model–data misfits normalized by data uncertainties for (left) PRIOR and (right) GLACIAL (b),(c) annual-, (d),(e) JFM-, and (f),(g) JAS-mean data. Blue (red) values indicate that the model is cold (warm) relative to the data. GLACIAL shows similar patterns of misfits as PRIOR, but with reduced amplitudes of both signs.

Fig. 4.

Model–data misfits of NSSTs are improved by the state estimation procedure. (a) Histograms of model–data misfits for annual-, JFM-, and JAS-mean values in PRIOR and GLACIAL normalized by data uncertainties. GLACIAL misfits are similar to a standard normal Gaussian distribution (gray), indicating an acceptable fit to the data. Model–data misfits normalized by data uncertainties for (left) PRIOR and (right) GLACIAL (b),(c) annual-, (d),(e) JFM-, and (f),(g) JAS-mean data. Blue (red) values indicate that the model is cold (warm) relative to the data. GLACIAL shows similar patterns of misfits as PRIOR, but with reduced amplitudes of both signs.

To reduce model–data misfits, we adjust glacial atmospheric conditions and other control variables using the method of Lagrange multipliers (see the  appendix). We found that while this approach reduced misfits of both signs, it was less effective at reducing the model cold bias. To reduce remaining biases after 10 iterations, we added a globally uniform increase of 2°C in all months to the first guess of surface air temperatures.1 After including these changes we ran 19 additional iterations for a total of 29. An additional temperature increase of 1°C was added to the control adjustments derived in January, February, and March to offset a further cold bias in that season. As a reference, a separate state estimate was produced without uniform temperature adjustments; the two solutions are compared in section 4.

Changes to atmospheric control variables are typically strongest at locations coinciding with MARGO gridded data, although large-scale changes show the ability of the data to influence the model state in regions remote from data locations (Fig. 2, right panels). Global temperature increases used to reduce the model cold bias are visible in Fig. 2f. Inferred changes to isopycnal diffusivities , diapycnal diffusivities , and eddy bolus velocity coefficients are small relative to their uncertainties σ, with changes on the order of σ at few locations (see Figs. S1–S3). Several authors have suggested that decreased sea level at the LGM may have led to increased diapycnal mixing rates in the ocean interior, as the area of shallow continental shelves where the bulk of tidal dissipation occurs in the modern ocean was reduced (Wunsch 2003; Schmittner et al. 2015). While we cannot rule out this possibility, we note that a distribution of mixing parameters similar to a modern estimate suffices to fit the MARGO data, as also pointed out by KN17. Changes to initial temperature and salinity (Figs. S4 and S5) are on the order of , as we might expect for a quasi-steady solution in which adjustments to initial conditions are not important to fit the data. In contrast, changes to air–sea fluxes of heat and freshwater play a dominant role in fitting the observations, consistent with the primary role of surface fluxes in the seasonal variability of heat and salt budgets in the upper ocean (Gill and Niller 1973). We do not claim that derived control variable changes are necessary to fit the data, only that they are sufficient and reasonable within their specified uncertainties.

Our best estimate of the glacial ocean state (GLACIAL) is the seasonal cycle of ocean variables in the MITgcm when it is run under control changes derived to fit the MARGO data, at the end of a 5000-yr-long model spinup period. Properties plotted and discussed are decadal means of the seasonal cycle in the last 10 years of the spinup period; because of the slow evolution of the model after the long integration, results are not sensitive to the choice of averaging interval. Spatial patterns of MARGO minus GLACIAL model–data misfits are similar to those for MARGO − PRIOR but with reduced amplitudes in most regions (Fig. 4, right). Average misfits in years 80–100 (not shown) and 4980–5000 (Fig. 4a) are reduced relative to PRIOR, and their normalized distribution lies close to the expected Gaussian. The ability to fit observations near the beginning and end of a 5000-yr-long integration reflects small model drifts (Fig. S7) relative to observational uncertainties. The result satisfies the data-based criteria of Eqs. (2) and (3) and supports the conclusions of DW14 and KN17 that it is possible to fit a primitive equation ocean model to the MARGO data. The fact that even our data-constrained model solution does not exactly fit the data reflects a combination of model and data errors; misfits to the state estimate are deemed acceptable in light of observational uncertainties. Subsequent adjoint iterations could further reduce model–data misfit, but at the risk of overfitting the data.

c. Analysis of the state estimate

We now describe properties of our best estimate of the LGM ocean state. When describing abyssal properties we focus on the Atlantic Ocean, where the number of paleoceanographic data is greatest.

1) The upper ocean

Differences in annual- and seasonal-mean NSSTs between GLACIAL and MODERN indicate global cooling at the LGM except for small-amplitude warming in parts of the Arctic and Southern Oceans and the equatorial Pacific (Fig. 5). The global-mean NSST difference is 2°C, similar to preceding estimates of 1.9 ± 1.8°C (MARGO), 2.2°C (KN17), and 2.4°C (in CCSM4; Brady et al. 2013). The strongest negative anomalies are found in the subpolar regions, particularly in the Northern Hemisphere. In addition to their data compilation, Waelbroeck et al. (2009) report a map of LGM minus modern surface temperature anomalies based on a nearest-neighbor interpolation algorithm. By comparison with their map, GLACIAL − MODERN anomalies resulting from our dynamical interpolation do not show pronounced zonal gradients in the equatorial Pacific and Atlantic Oceans, while in the northern North Atlantic we find that the sign of zonal gradients is reversed relative to Waelbroeck et al. (2009). Moreover, we find surface cooling, rather than warming, in both the North Pacific and the central Arctic. These disagreements arise because in addition to fitting the data, our anomaly estimates are constrained by model physics.

Fig. 5.

GLACIAL minus MODERN temperature anomalies in the uppermost grid box centered on 5-m water depth: (a) the annual mean, (b) JFM mean, and (c) JAS mean.

Fig. 5.

GLACIAL minus MODERN temperature anomalies in the uppermost grid box centered on 5-m water depth: (a) the annual mean, (b) JFM mean, and (c) JAS mean.

GLACIAL − MODERN anomalies arise from (i) control variable changes made to fit the MARGO data and (ii) the choice of a first-guess solution. Changes in ocean surface temperature resulting from fitting the MARGO data, which are illustrated by GLACIAL − PRIOR anomalies (Fig. 3e), have shorter length scales than changes arising from our choice of first guess, which are illustrated by PRIOR − MODERN anomalies (Figs. 3c,d). These differences are commensurate with mostly small-scale changes made to the control variables to fit the data compared to the larger-scale patterns of first-guess glacial atmospheric conditions arising from coupled atmosphere–ocean–ice dynamics in CCSM4 (Fig. 2). Changes in surface temperature are distinct from patterns of model bias (Figs. 3a,b), suggesting that, as used, the MARGO data are insufficient to eliminate model biases completely.

Sea ice extent in GLACIAL is greater than in MODERN (Fig. 6) in all seasons in both the Northern and Southern Hemispheres. The Arctic Ocean is filled with sea ice year-round, and winter sea ice covers much of the Nordic seas and the northwest Pacific. Winter ice thicknesses in the central Arctic are 3–5 m, with lower values in regions where ice coverage is seasonal. GLACIAL sea ice is seasonal in the northern North Atlantic, including in the Nordic seas, similar to results reported by de Vernal et al. (2006) and Waelbroeck et al. (2009). In the Southern Hemisphere, the 15% winter sea ice concentration isopleth, where concentration refers to the fractional area occupied by sea ice, is consistent with the maximum northward extent of sea ice reconstructed by Gersonde et al. (2005), whose Southern Ocean data are included in MARGO. It also falls within the range of northernmost sea ice extents simulated in PMIP3 models (Sime et al. 2016). In GLACIAL, regions where brine rejection occurs because of sea ice formation coincide with winter sea ice extent in the Southern Hemisphere (Fig. S6), and annual-mean salt fluxes due to brine rejection are increased ( kg s−1) relative to MODERN ( kg s−1).

Fig. 6.

Sea ice thickness in GLACIAL (colors) and 15% concentration isopleth in (a),(c) September and (b),(d) March for MODERN (gray contour) and GLACIAL (black contour).

Fig. 6.

Sea ice thickness in GLACIAL (colors) and 15% concentration isopleth in (a),(c) September and (b),(d) March for MODERN (gray contour) and GLACIAL (black contour).

The barotropic (vertically integrated) circulation in GLACIAL is intensified relative to MODERN (Fig. 7), especially in the Antarctic Circumpolar Current (ACC) and subpolar gyres. Annual-mean volume transport through the Drake Passage is 174 Sv (1 Sv ≡ 106 m3 s−1) in GLACIAL compared to 117 Sv in MODERN, presumably associated with differences in winds and increased production of Antarctic Bottom Water (AABW) in GLACIAL, which can steepen isopycnal slopes in the ACC (Gent et al. 2001; Hogg 2010). Like DW14, we find an increased southward return flow in the eastern interior of the North Atlantic subtropical gyre in GLACIAL relative to MODERN, though the eastward shift of the Atlantic subpolar gyre that DW14 describe is not evident. Increases in barotropic gyre circulation in GLACIAL are consistent with increased wind stress and wind stress curl.

Fig. 7.

(a) Barotropic (vertically integrated) streamfunction in GLACIAL. (b) GLACIAL − MODERN barotropic streamfunction anomaly. Barotropic transport in the ACC in GLACIAL exceeds that in MODERN by more than 100 Sv in some places.

Fig. 7.

(a) Barotropic (vertically integrated) streamfunction in GLACIAL. (b) GLACIAL − MODERN barotropic streamfunction anomaly. Barotropic transport in the ACC in GLACIAL exceeds that in MODERN by more than 100 Sv in some places.

Locations of large winter mixed layer depths (MLDs) are thought to be important for setting distributions of abyssal tracers because they determine where surface water properties are communicated to the abyssal interior (Gebbie and Huybers 2011; Amrhein et al. 2015) with possible implications for AMOC strength (Oka et al. 2012). Comparison of maximum winter MLDs in GLACIAL and MODERN reveals differences in regions of both subduction (e.g., in the model North Atlantic Current) and high-latitude convection (Fig. 8). In GLACIAL, reduced convection in the northeast North Atlantic and Arctic Oceans is due in part to (i) reduced areas of marginal seas from lower sea levels and (ii) fresher surface waters in those regions (Figs. 3d,f). Deeper winter mixed layers in the Labrador Sea are consistent with surface buoyancy losses from ocean cooling downwind of the Laurentide Ice Sheet (Fig. 2c). While MLDs are likely affected by the model drifts discussed in section 3a, differences between MODERN and GLACIAL, which are affected by similar drifts, motivate speculation that a shift of winter maximum MLDs from the eastern to western North Atlantic may contribute to differences observed in distributions of abyssal ocean tracers between the LGM and today (e.g., Keigwin 2004; Curry and Oppo 2005; Marchitto and Broecker 2006) because of a change in deep water source regions.

Fig. 8.

Mixed layer depths computed using the MLD criterion of Kara et al. (2000) for (top) MODERN and (bottom) GLACIAL in (a),(c) March and (b),(d) September.

Fig. 8.

Mixed layer depths computed using the MLD criterion of Kara et al. (2000) for (top) MODERN and (bottom) GLACIAL in (a),(c) March and (b),(d) September.

2) The abyssal Atlantic Ocean

Abyssal waters in GLACIAL are everywhere colder than in MODERN in the Atlantic, where zonal mean potential temperatures are reduced by roughly between 0.5° and 1.0°C (Figs. 9a,c). Increased salinity stratification (Figs. 9b,d) is primarily responsible for greater density stratification (Fig. 10, contours). Higher vertical salinity stratification in the GLACIAL Atlantic is consistent with larger rates of Southern Ocean brine rejection (Fig. S6), though decreased high-latitude precipitation (Figs. 2g,h) may also play a role. GLACIAL − MODERN abyssal salinity anomalies lie within the 2σ uncertainty ranges of LGM minus modern anomalies estimated from pore fluids recovered in the Pacific Ocean (Table 2; Insua et al. 2014), and are qualitatively consistent with the inference of a more salinity-stratified LGM abyssal ocean. However, we reproduce neither the relatively low pore fluid salinity anomaly estimated at Bermuda Rise (33.7°N, 57.6°W), nor the large anomaly at Shona Rise in the Southern Ocean (50.0°S, 5.9°E) that was a focus of Adkins et al. (2002). Misfits between the state estimate and pore fluid reconstructions could be due to model biases—for example, inaccurate parameterization of brine rejection—and/or to misinterpretation of the observations (Miller et al. 2015; Wunsch 2016).

Fig. 9.

Atlantic zonal mean (a),(c) potential temperature and (b),(d) salinity (on the practical salinity scale) in (top) MODERN and (bottom) GLACIAL.

Fig. 9.

Atlantic zonal mean (a),(c) potential temperature and (b),(d) salinity (on the practical salinity scale) in (top) MODERN and (bottom) GLACIAL.

Fig. 10.

(a),(c) Atlantic and (b),(d) global zonal mean meridional overturning streamfunctions in (top) MODERN and (bottom) GLACIAL. Contours denote potential density (kg m−3) minus a reference value of 1000 kg m−3. Note differences in color bars between global and Atlantic overturning and nonconstant potential density contour intervals. Potential density differences between LGM and GLACIAL reflect in part a global-mean salinity increase of 1.1 in GLACIAL.

Fig. 10.

(a),(c) Atlantic and (b),(d) global zonal mean meridional overturning streamfunctions in (top) MODERN and (bottom) GLACIAL. Contours denote potential density (kg m−3) minus a reference value of 1000 kg m−3. Note differences in color bars between global and Atlantic overturning and nonconstant potential density contour intervals. Potential density differences between LGM and GLACIAL reflect in part a global-mean salinity increase of 1.1 in GLACIAL.

Table 2.

Comparison of LGM–Holocene bottom-water salinity anomalies from pore fluid measurements and this study. Pore fluid measurements were not included in the cost function and thus provide an independent assessment of the LGM state estimate. Salinity differences are from the deepest grid box at the model grid location nearest core sites. All values are on the practical salinity scale.

Comparison of LGM–Holocene bottom-water salinity anomalies from pore fluid measurements and this study. Pore fluid measurements were not included in the cost function and thus provide an independent assessment of the LGM state estimate. Salinity differences  are from the deepest grid box at the model grid location nearest core sites. All values are on the practical salinity scale.
Comparison of LGM–Holocene bottom-water salinity anomalies from pore fluid measurements and this study. Pore fluid measurements were not included in the cost function and thus provide an independent assessment of the LGM state estimate. Salinity differences  are from the deepest grid box at the model grid location nearest core sites. All values are on the practical salinity scale.

The upper cell of the AMOC (Figs. 10a,c) is deeper and stronger in GLACIAL than MODERN, qualitatively similar to results from most PMIP3 models (Muglia and Schmittner 2015). Comparing these results to those from other studies is complicated by biases toward shoaled and weakened upper cells present in both MODERN and GLACIAL. Thus while our result of a relatively stronger, deeper GLACIAL cell contrasts with that of KN17, who found a stronger, shallower upper LGM AMOC cell, in absolute terms the LGM AMOC circulations in the two studies are similar, despite differences in state estimation procedures. Our comparison of GLACIAL and MODERN also contrasts with the idealized model of Ferrari et al. (2014), who suggested that greater sea ice extent at the LGM would shift outcropping isopycnals in the ACC equatorward and shoal the isopycnal surface separating upper and lower AMOC cells. In MODERN and GLACIAL, the upper and lower cells (separated by the zero contour in Figs. 10a,c) are separated by the 28 and 29 kg m−3 potential density anomaly isopleths, respectively. The deeper position of the dividing isopycnal in GLACIAL relative to MODERN is accompanied by steeper ACC isopycnal slopes (Figs. 10b,d), suggesting that the deeper, stronger GLACIAL upper AMOC cell is associated with stronger ACC baroclinicity. Because low-resolution models may poorly represent the role of eddies in modulating wind-driven changes to the ACC (Abernathey et al. 2011), future work should investigate LGM ACC isopycnal steepening in an eddy-resolving ocean model.

To test whether the GLACIAL circulation supports inference of a greater volume of southern-source water in the Atlantic Ocean, we perform a dye release experiment by fixing passive tracer boundary conditions in surface grid boxes to a concentration of 1 south of 60°S and to 0 elsewhere in the 5000-yr-long simulations of GLACIAL and MODERN (Fig. 11). After 5000 years, the distribution of this tracer in the Atlantic is very similar in the two simulations. From this result, we conclude that fitting the MITgcm as configured to the MARGO data does not require southern-source waters to shoal in the abyssal Atlantic, as has been suggested by interpretations of LGM abyssal tracers (e.g., Curry and Oppo 2005). This result further demonstrates the importance of including abyssal tracer data to constrain the glacial abyssal state.

Fig. 11.

Passive tracer concentrations at sections in the western Atlantic in (a) MODERN and (b) GLACIAL after 5000 years of integration. Tracer surface values are held at 1 south of 60° and 0 elsewhere.

Fig. 11.

Passive tracer concentrations at sections in the western Atlantic in (a) MODERN and (b) GLACIAL after 5000 years of integration. Tracer surface values are held at 1 south of 60° and 0 elsewhere.

4. Discussion

This paper presents a dynamical interpolation of LGM NSST observations to a seasonal cycle of gridded ocean temperature, salinity, velocities, and sea ice variables that is fully consistent with the physics of the MITgcm. While we do not claim that our glacial state estimate is a unique fit of ocean and climate variables to the data, it is a dynamically plausible hypothesis for LGM conditions. In agreement with simulations from climate models subject to glacial climate boundary conditions and with previous glacial ocean state estimates, the upper ocean at the LGM is inferred to be colder than today by 2°C in the global mean. The barotropic ocean circulation is inferred to be stronger, consistent with greater wind stress and wind stress curl. Gyre circulations, while stronger, are structurally similar to the modern circulation. Both perennial and seasonal sea ice extents are larger, and the central Arctic is filled with sea ice year-round. Regions of deep winter mixed layers differ from the modern. The abyssal ocean is more strongly salinity stratified, with an upper AMOC cell that is stronger and deeper.

Our state estimate has both similarities and differences with the state estimates of DW14 and KN17. For example, NSST fields reconstructed in KN17 are smoother than ours, which presumably reflects in part the isotropic, roughly 9° smoothing that those authors imposed on control variable adjustments. In contrast, DW14 report strong small-scale gradients in reconstructed surface ocean temperature between locations with and without LGM data. This is likely because DW14 used modern oceanographic conditions as a first guess; in contrast, our PRIOR simulation introduces large-scale patterns of cooling set by atmospheric conditions derived from coupled ocean–atmosphere–ice CCSM4 simulations. Similar to KN17, we find the largest GLACIAL minus MODERN negative anomalies in the subtropics, but our estimated surface cooling is more uniform. Like KN17, we observe a stronger salinity stratification at the LGM, which we attribute in part to greater sea ice extent; however, while KN17 find a stronger, shallower AMOC, ours is stronger and deeper. These differences may stem from a variety of factors, including the use of abyssal tracer observations in the KN17 solution, the use of seasonal NSST observations in our solution, differences in model equilibration and spatial resolution, and differences in turbulent transport parameters.

Because none of the solutions in DW14, KN17, and this work include error estimates, it is difficult to determine whether the solutions are truly in disagreement. There is currently no straightforward means to determine comprehensive uncertainties in ocean state estimates derived using the method of Lagrange multipliers. Developing tools for uncertainty quantification is an important and ongoing effort in ocean state estimation (Kalmikov and Heimbach 2014). Contrasting results in DW14, KN17, and this study suggest a sensitivity to prior choices of model controls and covariances and point to difficulties in constraining the deep ocean circulation at the LGM from available observations (e.g., LeGrand and Wunsch 1995; Huybers et al. 2007; Marchal and Curry 2008; Burke et al. 2011; Kurahashi-Nakamura et al. 2014; Gebbie et al. 2016).

We find that global-mean temperature changes are necessary to reduce an overall model cold bias. KN17 did not find such an adjustment necessary, possibly because of different choices of the first-guess ocean state and atmospheric forcing. To assess the sensitivity of our inferences to global-mean temperature changes, a separate state estimate (GLACIAL_s) was computed over six iterations without imposing such changes. Relative to our reference solution (GLACIAL), GLACIAL_s shows greater summer sea ice extent and thickness in both hemispheres (Fig. S9), a stronger reduction in NSSTs (Fig. S12), colder and saltier Atlantic bottom waters (Fig. S10), greater salinity and density stratification (Figs. S10, S11), and a marginally stronger and shallower AMOC upper cell (Fig. S11). These differences are not so large as to change our overall conclusions. To evaluate whether the mean model–data misfit arises from the first guess constructed by adding CCSM4 LGM-PI anomalies to modern ECCO atmospheric conditions, we ran an additional state estimate (not shown) using CCSM4 LGM conditions as a first guess. We find a similar model–data bias, suggesting that the first-guess choice was not a major factor. A similar result (not shown) was obtained for a first guess derived from a different coupled model LGM simulation (MIROC; Sueyoshi et al. 2013). Ultimately, the mean model–data misfit may be due to biases in the data, the MITgcm, our choice of first-guess boundary conditions, the coupled models used to generate first guesses, and/or the choice of boundary conditions used to force coupled models. Resolving the origin of this bias is important given the use of LGM climate to infer climate sensitivity (Schmittner et al. 2011; Hargreaves et al. 2012) and the use of large-scale atmospheric cooling to simulate LGM conditions in idealized models (Jansen 2017).

The GLACIAL solution is consistent with a seasonally steady state (i.e., a single repeating annual cycle) insofar as it is taken from a 5000-yr-long model integration that fits LGM observations near its beginning and end. The fact that control adjustments derived over a relatively short period (100 years) can still fit observations after a longer integration (5000 years) is not surprising given that the data are fit largely by local changes in surface heat fluxes, to which the upper ocean adjusts on time scales shorter than 100 years; Forget et al. (2015b) found a similar result for the ECCO state estimate. Adjoint integration times longer than afforded here could reveal other, longer-time-scale mechanisms that also permit the model to fit the data.

5. Perspectives

This work points to several ways forward to improve paleoceanographic state estimation. First, there is the issue of how changes are made to atmospheric conditions in order to fit paleoceanographic observations. In ECCO, first-guess atmospheric conditions come from reanalysis products, which are constrained both by satellite observations and by coupled models. The assumption in ECCO is that the reanalysis products are sufficiently accurate that changes to accommodate ocean observations will have small amplitudes that are uncorrelated over large spatial scales. In contrast, first guesses of atmospheric conditions for preinstrumental periods are poorly constrained—here, for instance, we use the quasi equilibrium of a coupled climate model (CCSM4)—and we should expect that they differ from true atmospheric states on all spatial scales, reflecting the full range of coupled ocean–atmosphere–ice dynamics. Instead, in our state estimate, we infer “patchy” adjustments to atmospheric controls (Fig. 2, right) whose length scales reflect data availability and ocean dynamics and are not informed by atmospheric or coupled dynamics. While KN17 mitigated this patchiness by smoothing control variables in space, it is not obvious that this approach yields more accurate atmospheric fields. A separate issue is that different choices of atmospheric controls can have similar effects on the ocean state, leading to degeneracies; for instance, similarities between changes in shortwave and longwave radiation inferred to fit the data (Figs. 2l,n) reflect the inability of the data and model to differentiate between different sources of ocean heating. Finally, the absence of feedbacks between the ocean and atmosphere in the presence of large changes to atmospheric forcing can lead to unphysical patterns of heating and cooling that contribute to model drift. These caveats urge caution in attempting to rationalize inferred atmospheric conditions physically. They point to a need for more accurate estimates of how uncertainties in control variables covary in space and time and, ultimately, to a need for coupled ocean–atmosphere–ice state estimation.

Second, assuming a steady or seasonally steady LGM ocean circulation at once provides a strong constraint on the state estimate and poses technical challenges for reaching model equilibrium. In this work, we found that long simulations intended to equilibrate the deep ocean to reconstructed surface conditions led to model biases, due in part to model drifts. While drifts can be reduced by relaxing ocean surface values of temperature and salinity to fixed climatological values (Danabasoglu et al. 2014), relaxation generates undesirable sources and sinks of temperature and salinity in the state estimate. One approach to mitigate this problem could be to adjust patterns and time scales of relaxation in the control vector to fit the data. More broadly, the extent to which the ocean circulation is ever in equilibrium (including at the LGM) is unclear. Paleoceanographic data provide an important arena for challenging assumptions about climate stationarity, and steadiness should only be assumed when absolutely necessary. Satisfying a version of Eqs. (2) and (3) yields a solution that is only as steady as the data require and provides a less restrictive modeling criterion for the steadiness of the LGM and other geologic intervals.

Third, this work raises the question of how well suited the current generation of ocean models is to paleoceanographic state estimation, particularly for abyssal properties. Unlike in the modern state estimation problem, there are no direct measurements of ocean hydrography at the LGM. While the MARGO data can constrain some features in the surface ocean, the impact of surface temperature data on inferences of abyssal properties is mediated by deep-water formation processes that are typically parameterized and that occur in poorly sampled regions such as Antarctic shelves and the Labrador Sea. Locations and rates of high-latitude deep water formation are important for setting abyssal values of temperature, salinity, and passive tracers (Amrhein et al. 2015). We expect that improving model representations of high-latitude processes will be effective at increasing the accuracy of reconstructed abyssal ocean conditions at the LGM.

Finally, LGM state estimation will benefit from an increased number, spatial coverage, and diversity of proxy observations, as well as greater understanding of how to represent those observations in numerical models. Of particular utility is the inclusion of abyssal tracer measurements, which have inspired many hypotheses about LGM watermass reorganizations (e.g., Curry and Oppo 2005). While KN17 take the important step of including carbon and oxygen stable isotope measurements in their state estimation, realizing the full potential of these measurements is challenging because of the long time scales of tracer equilibration (Wunsch and Heimbach 2008), which necessitates running long and computationally expensive adjoint simulations. Dail (2012), Amrhein (2016), and KN17 describe technical improvements on this front that should be explored in future work. Ultimately, the goal is to derive a state estimate using all possible observations from the LGM and to include new observations as they become available.

Acknowledgments

The authors acknowledge Patrick Heimbach, Raffaele Ferrari, Mick Follows, David McGee, and LuAnne Thompson for helpful comments on earlier versions of the work, as well as thoughtful comments from three reviewers that improved the manuscript. An Nguyen, Jeff Scott, and Jean-Michel Campin provided important assistance with the MITgcm adjoint. Computing time was provided by NASA Advanced Supercomputing and Ames Research Center. DEA was supported by a NSF Graduate Research Fellowship and NSF Grant OCE-1060735. OM acknowledges support from the NSF. GF was supported by NASA Award 1553749 and Simons Foundation Award 549931. The adjoint model was generated using TAF (Giering and Kaminski 1998). Output from the state estimate is available by contacting the corresponding author.

APPENDIX

Seasonal State Estimation by the Method of Least Squares with Lagrange Multipliers

Seasonal state estimation seeks a set of annually repeating control variables that minimizes a scalar cost function with three contributions. The first contribution, Jdata, is the squared, weighted model–data misfit, which has three terms,

 
formula

where the L are numbers of observations available in each time period (annual, JFM, and JAS); , , and are observations; the matrices , , and have the form of , where angle brackets denote the expected value, and are observational noise covariances constructed from MARGO uncertainty estimates; and the matrices , , and relate model variables across space and time to the data. Multiplication by for JFM and JAS divides by the number of model monthly means included in the cost function; for annual observations, this factor is 1. The second contribution is

 
formula

where is a vector of Lagrange multipliers that serves to impose the MITgcm model equations upon the solution. The vector can more generally represent model errors as well. The last contribution,

 
formula

penalizes control adjustments, where is the error covariance of the control variables. Here is assumed to be zero except for diagonal values that are equal to the squared standard deviations σ assumed for control variable uncertainties (Table 1).

Minimization of the total cost function, , is a problem of constrained nonlinear optimization. The dimension of the state vector and the complexity of L preclude an analytical solution. Instead, automatic differentiation of the MITgcm code (Giering and Kaminski 1998) is used to adjust the control variables iteratively in the direction of locally steepest descent using a quasi-Newton algorithm (Gilbert and Lemaréchal 1989). After each iteration, the cost function and local sensitivities are recomputed and the procedure is repeated until the distribution of model–data misfits, normalized by observational uncertainty, approximates a normal (Gaussian) distribution with zero mean and unity variance. At this point, the state estimate is considered acceptable, so long as the control adjustments are also acceptable. For a more detailed discussion see Wunsch and Heimbach (2007).

REFERENCES

REFERENCES
Abernathey
,
R.
,
J.
Marshall
, and
D.
Ferreira
,
2011
:
The dependence of Southern Ocean meridional overturning on wind stress
.
J. Phys. Oceanogr.
,
41
,
2261
2278
, https://doi.org/10.1175/JPO-D-11-023.1.
Adcroft
,
A.
,
C.
Hill
,
J.-M.
Campin
,
J.
Marshall
, and
P.
Heimbach
,
2004
: Overview of the formulation and numerics of the MIT GCM. Proc. ECMWF Seminar on Recent Developments in Numerical Methods for Atmospheric and Ocean Modelling, Reading, United Kingdom, ECMWF,
139
149
, https://www.ecmwf.int/sites/default/files/elibrary/2004/7642-overview-formulation-and-numerics-mit-gcm.pdf.
Adkins
,
J. F.
,
K.
McIntyre
, and
D. P.
Schrag
,
2002
:
The salinity, temperature, and O of the glacial deep ocean
.
Science
,
298
,
1769
1773
, https://doi.org/10.1126/science.1076252.
Alkama
,
R.
,
M.
Kageyama
, and
G.
Ramstein
,
2006
:
Freshwater discharges in a simulation of the Last Glacial Maximum climate using improved river routing
.
Geophys. Res. Lett.
,
33
,
L21709
, https://doi.org/10.1029/2006GL027746.
Amrhein
,
D.
,
2016
: Inferring ocean circulation during the Last Glacial Maximum and last deglaciation using data and models. Ph.D. thesis, Massachusetts Institute of Technology–WHOI, 192 pp.
Amrhein
,
D.
,
G.
Gebbie
,
O.
Marchal
, and
C.
Wunsch
,
2015
:
Inferring surface water equilibrium calcite δ18O during the last deglacial period from benthic foraminiferal records: Implications for ocean circulation
.
Paleoceanography
,
30
,
1470
1489
, https://doi.org/10.1002/2014PA002743.
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
, https://doi.org/10.5194/cp-3-261-2007.
Brady
,
E. C.
,
B. L.
Otto-Bliesner
,
J. E.
Kay
, and
N.
Rosenbloom
,
2013
:
Sensitivity to glacial forcing in the CCSM4
.
J. Climate
,
26
,
1901
1925
, https://doi.org/10.1175/JCLI-D-11-00416.1.
Brovkin
,
V.
,
A.
Ganopolski
,
D.
Archer
, and
S.
Rahmstorf
,
2007
:
Lowering of glacial atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry
.
Paleoceanography
,
22
,
PA4202
, https://doi.org/10.1029/2006PA001380.
Bryan
,
K.
,
1984
:
Accelerating the convergence to equilibrium of ocean–climate models
.
J. Phys. Oceanogr.
,
14
,
666
673
, https://doi.org/10.1175/1520-0485(1984)014<0666:ATCTEO>2.0.CO;2.
Bugnion
,
V.
, and
C.
Hill
,
2006
:
Equilibration mechanisms in an adjoint ocean general circulation model
.
Ocean Dyn.
,
56
,
51
61
, https://doi.org/10.1007/s10236-005-0052-z.
Burke
,
A.
,
O.
Marchal
,
L. I.
Bradtmiller
,
J. F.
McManus
, and
R.
François
,
2011
:
Application of an inverse method to interpret 231Pa/230Th observations from marine sediments
.
Paleoceanography
,
26
,
PA1212
, https://doi.org/10.1029/2010PA002022.
Clark
,
P. U.
, and Coauthors
,
2012
:
Global climate evolution during the last deglaciation
.
Proc. Natl. Acad. Sci. USA
,
109
,
E1134
E1142
, https://doi.org/10.1073/pnas.1116619109.
Curry
,
W. B.
, and
D. W.
Oppo
,
2005
:
Glacial water mass geometry and the distribution of C of CO2 in the western Atlantic Ocean
.
Paleoceanography
,
20
,
PA1017
, https://doi.org/10.1029/2004PA001021.
Curry
,
W. B.
,
J. C.
Duplessy
,
L. D.
Labeyrie
, and
N. J.
Shackleton
,
1988
:
Changes in the distribution of δ13C of deep water ΣCO2 between the last glaciation and the Holocene
.
Paleoceanography
,
3
,
317
341
, https://doi.org/10.1029/PA003i003p00317.
Dail
,
H. J.
,
2012
: Atlantic Ocean circulation at the Last Glacial Maximum: Inferences from data and models. Ph.D. thesis, Massachusetts Institute of Technology–WHOI, 236 pp.
Dail
,
H. J.
, and
C.
Wunsch
,
2014
:
Dynamical reconstruction of upper-ocean conditions in the Last Glacial Maximum Atlantic
.
J. Climate
,
27
,
807
823
, https://doi.org/10.1175/JCLI-D-13-00211.1.
Danabasoglu
,
G.
, and Coauthors
,
2014
:
North Atlantic simulations in Coordinated Ocean-Ice Reference Experiments phase II (CORE-II). Part I: Mean states
.
Ocean Modell.
,
73
,
76
107
, https://doi.org/10.1016/j.ocemod.2013.10.005.
de Vernal
,
A.
,
A.
Rosell-Melé
,
M.
Kucera
,
C.
Hillaire-Marcel
,
F.
Eynaud
,
M.
Weinelt
,
T.
Dokken
, and
M.
Kageyama
,
2006
:
Comparing proxies for the reconstruction of LGM sea-surface conditions in the northern North Atlantic
.
Quat. Sci. Rev.
,
25
,
2820
2834
, https://doi.org/10.1016/j.quascirev.2006.06.006.
Duplessy
,
J. C.
,
N. J.
Shackleton
,
R. G.
Fairbanks
,
L.
Labeyrie
,
D.
Oppo
, and
N.
Kallel
,
1988
:
Deepwater source variations during the last climatic cycle and their impact on the global deepwater circulation
.
Paleoceanography
,
3
,
343
360
, https://doi.org/10.1029/PA003i003p00343.
Fairbanks
,
R. G.
,
1989
:
A 17 000-year glacio-eustatic sea level record: Influence of glacial melting rates on the Younger Dryas event and deep-ocean circulation
.
Nature
,
342
,
637
642
, https://doi.org/10.1038/342637a0.
Fekete
,
B. M.
,
C. J.
Vörösmarty
, and
W.
Grabs
,
2002
:
High-resolution fields of global runoff combining observed river discharge and simulated water balances
.
Global Biogeochem. Cycles
,
16
,
1042
, https://doi.org/10.1029/1999GB001254.
Ferrari
,
R.
,
M. F.
Jansen
,
J. F.
Adkins
,
A.
Burke
,
A. L.
Stewart
, and
A. F.
Thompson
,
2014
:
Antarctic sea ice control on ocean circulation in present and glacial climates
.
Proc. Natl. Acad. Sci. USA
,
111
,
8753
8758
, https://doi.org/10.1073/pnas.1323922111.
Forget
,
G.
,
J.-M.
Campin
,
P.
Heimbach
,
C.
Hill
,
R. M.
Ponte
, and
C.
Wunsch
,
2015a
:
ECCO version 4: An integrated framework for non-linear inverse modeling and global ocean state estimation
.
Geosci. Model Dev.
,
8
,
3071
3104
, https://doi.org/10.5194/gmd-8-3071-2015.
Forget
,
G.
,
D.
Ferreira
, and
X.
Liang
,
2015b
:
On the observability of turbulent transport rates by Argo: Supporting evidence from an inversion experiment
.
Ocean Sci.
,
11
,
839
853
, https://doi.org/10.5194/os-11-839-2015.
Gaspar
,
P.
,
Y.
Grégoris
, and
J.-M.
Lefevre
,
1990
:
A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and long-term upper ocean study site
.
J. Geophys. Res.
,
95
,
16 179
16 193
, https://doi.org/10.1029/JC095iC09p16179.
Gebbie
,
G.
,
2014
:
How much did glacial North Atlantic water shoal?
Paleoceanography
,
29
,
190
209
, https://doi.org/10.1002/2013PA002557.
Gebbie
,
G.
, and
P.
Huybers
,
2006
:
Meridional circulation during the Last Glacial Maximum explored through a combination of South Atlantic O observations and a geostrophic inverse model
.
Geochem. Geophys. Geosyst.
,
7
,
Q11N07
, https://doi.org/10.1029/2006GC001383.
Gebbie
,
G.
, and
P.
Huybers
,
2011
:
How is the ocean filled?
Geophys. Res. Lett.
,
38
,
L06604
, https://doi.org/10.1029/2011GL046769.
Gebbie
,
G.
,
G. J.
Streletz
, and
H. J.
Spero
,
2016
:
How well would modern-day oceanic property distributions be known with paleoceanographic-like observations?
Paleoceanography
,
31
,
472
490
, https://doi.org/10.1002/2015PA002917.
Gent
,
P. R.
, and
J. C.
McWilliams
,
1990
:
Isopycnal mixing in ocean circulation models
.
J. Phys. Oceanogr.
,
20
,
150
155
, https://doi.org/10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2.
Gent
,
P. R.
,
W. G.
Large
, and
F. O.
Bryan
,
2001
:
What sets the mean transport through Drake Passage?
J. Geophys. Res.
,
106
,
2693
2712
, https://doi.org/10.1029/2000JC900036.
Gent
,
P. R.
, and Coauthors
,
2011
:
The Community Climate System Model version 4
.
J. Climate
,
24
,
4973
4991
, https://doi.org/10.1175/2011JCLI4083.1.
Gersonde
,
R.
,
X.
Crosta
,
A.
Abelmann
, and
L.
Armand
,
2005
:
Sea-surface temperature and sea ice distribution of the Southern Ocean at the EPILOG Last Glacial Maximum—A circum-Antarctic view based on siliceous microfossil records
.
Quat. Sci. Rev.
,
24
,
869
896
, https://doi.org/10.1016/j.quascirev.2004.07.015.
Giering
,
R.
, and
T.
Kaminski
,
1998
:
Recipes for adjoint code construction
.
ACM Trans. Math. Software
,
24
,
437
474
, https://doi.org/10.1145/293686.293695.
Gilbert
,
J. C.
, and
C.
Lemaréchal
,
1989
:
Some numerical experiments with variable-storage quasi-Newton algorithms
.
Math. Program.
,
45
,
407
435
, https://doi.org/10.1007/BF01589113.
Gill
,
A. E.
, and
P. P.
Niller
,
1973
:
The theory of the seasonal variability in the ocean
.
Deep-Sea Res. Oceanogr. Abstr.
,
20
,
141
177
, https://doi.org/10.1016/0011-7471(73)90049-1.
Griffies
,
S. M.
, and Coauthors
,
2009
:
Coordinated Ocean-Ice Reference Experiments (COREs)
.
Ocean Modell.
,
26
,
1
46
, https://doi.org/10.1016/j.ocemod.2008.08.007.
Hargreaves
,
J. C.
,
A.
Paul
,
R.
Ohgaito
,
A.
Abe-Ouchi
, and
J. D.
Annan
,
2011
:
Are paleoclimate model ensembles consistent with the MARGO data synthesis?
Climate Past
,
7
,
917
933
, https://doi.org/10.5194/cp-7-917-2011.
Hargreaves
,
J. C.
,
J. D.
Annan
,
M.
Yoshimori
, and
A.
Abe-Ouchi
,
2012
:
Can the Last Glacial Maximum constrain climate sensitivity?
Geophys. Res. Lett.
,
39
,
L24702
, https://doi.org/10.1029/2012GL053872.
Hogg
,
A. M.
,
2010
:
An Antarctic Circumpolar Current driven by surface buoyancy forcing
.
Geophys. Res. Lett.
,
37
,
L23601
, https://doi.org/10.1029/2010GL044777.
Huybers
,
P.
,
G.
Gebbie
, and
O.
Marchal
,
2007
:
Can paleoceanographic tracers constrain meridional circulation rates?
J. Phys. Oceanogr.
,
37
,
394
407
, https://doi.org/10.1175/JPO3018.1.
Insua
,
T. L.
,
A. J.
Spivack
,
D.
Graham
,
S.
D’Hondt
, and
K.
Moran
,
2014
:
Reconstruction of Pacific Ocean bottom water salinity during the Last Glacial Maximum
.
Geophys. Res. Lett.
,
41
,
2914
2920
, https://doi.org/10.1002/2014GL059575.
Jansen
,
M. F.
,
2017
:
Glacial ocean circulation and stratification explained by reduced atmospheric temperature
.
Proc. Natl. Acad. Sci. USA
,
114
,
45
50
, https://doi.org/10.1073/pnas.1610438113.
Jansen
,
M. F.
, and
L.-P.
Nadeau
,
2016
:
The effect of Southern Ocean surface buoyancy loss on the deep-ocean circulation and stratification
.
J. Phys. Oceanogr.
,
46
,
3455
3470
, https://doi.org/10.1175/JPO-D-16-0084.1.
Kalmikov
,
A. G.
, and
P.
Heimbach
,
2014
:
A Hessian-based method for uncertainty quantification in global ocean state estimation
.
SIAM J. Sci. Comput.
,
36
,
S267
S295
, https://doi.org/10.1137/130925311.
Kara
,
A. B.
,
P. A.
Rochford
, and
H. E.
Hurlburt
,
2000
:
An optimal definition for ocean mixed layer depth
.
J. Geophys. Res.
,
105
,
16 803
16 821
, https://doi.org/10.1029/2000JC900072.
Keigwin
,
L. D.
,
2004
:
Radiocarbon and stable isotope constraints on Last Glacial Maximum and Younger Dryas ventilation in the western North Atlantic
.
Paleoceanography
,
19
,
PA4012
, https://doi.org/10.1029/2004PA001029.
Kurahashi-Nakamura
,
T.
,
M.
Losch
, and
A.
Paul
,
2014
:
Can sparse proxy data constrain the strength of the Atlantic meridional overturning circulation?
Geosci. Model Dev.
,
7
,
419
432
, https://doi.org/10.5194/gmd-7-419-2014.
Kurahashi-Nakamura
,
T.
,
A.
Paul
, and
M.
Losch
,
2017
:
Dynamical reconstruction of the global ocean state during the Last Glacial Maximum
.
Paleoceanography
,
32
,
326
350
, https://doi.org/10.1002/2016PA003001.
Large
,
W. G.
, and
S. G.
Yeager
,
2004
: Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies. NCAR Tech. Note NCAR/TN-460+STR, 105 pp., https://doi.org/10.5065/D6KK98Q6.
LeGrand
,
P.
, and
C.
Wunsch
,
1995
:
Constraints from paleotracer data on the North Atlantic circulation during the Last Glacial Maximum
.
Paleoceanography
,
10
,
1011
1045
, https://doi.org/10.1029/95PA01455.
Losch
,
M.
,
D.
Menemenlis
,
J.-M.
Campin
,
P.
Heimbach
, and
C.
Hill
,
2010
:
On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations
.
Ocean Modell.
,
33
,
129
144
, https://doi.org/10.1016/j.ocemod.2009.12.008.
Lumpkin
,
R.
, and
K.
Speer
,
2007
:
Global ocean meridional overturning
.
J. Phys. Oceanogr.
,
37
,
2550
2562
, https://doi.org/10.1175/JPO3130.1.
Lynch-Stieglitz
,
J.
, and Coauthors
,
2007
:
Atlantic meridional overturning circulation during the Last Glacial Maximum
.
Science
,
316
,
66
69
, https://doi.org/10.1126/science.1137127.
Marchal
,
O.
, and
W. B.
Curry
,
2008
:
On the abyssal circulation in the glacial Atlantic
.
J. Phys. Oceanogr.
,
38
,
2014
2037
, https://doi.org/10.1175/2008JPO3895.1.
Marchitto
,
T. M.
, and
W. S.
Broecker
,
2006
:
Deep water mass geometry in the glacial Atlantic Ocean: A review of constraints from the paleonutrient proxy Cd/Ca
.
Geochem. Geophys. Geosyst.
,
7
,
Q12003
, https://doi.org/10.1029/2006GC001323.
Marchitto
,
T. M.
,
D. W.
Oppo
, and
W. B.
Curry
,
2002
:
Paired benthic foraminiferal Cd/Ca and Zn/Ca evidence for a greatly increased presence of Southern Ocean water in the glacial North Atlantic
.
Paleoceanography
,
17
,
1038
, https://doi.org/10.1029/2000PA000598.
Marshall
,
J.
,
A.
Adcroft
,
C.
Hill
,
L.
Perelman
, and
C.
Heisey
,
1997
:
A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers
.
J. Geophys. Res.
,
102
,
5753
5766
, https://doi.org/10.1029/96JC02775.
Marzocchi
,
A.
, and
M. F.
Jansen
,
2017
:
Connecting Antarctic sea ice to deep-ocean circulation in modern and glacial climate simulations
.
Geophys. Res. Lett.
,
44
,
6286
6295
, https://doi.org/10.1002/2017GL073936.
McIntyre
,
A.
,
N. G.
Kipp
,
A. W.
,
T.
Crowley
,
T.
Kellogg
,
J. V.
Gardner
,
W.
Prell
, and
W. F.
Ruddiman
,
1976
: Glacial North Atlantic 18 000 years ago: A CLIMAP reconstruction. Investigation of Late Quaternary Paleoceanography and Paleoclimatology, GSA Memoirs, Vol. 145, Geological Society of America, 43–76.
Miller
,
M. D.
,
M.
Simons
,
J. F.
Adkins
, and
S. E.
Minson
,
2015
:
The information content of pore fluid δ18O and [Cl−]
.
J. Phys. Oceanogr.
,
45
,
2070
2094
, https://doi.org/10.1175/JPO-D-14-0203.1.
Mix
,
A. C.
,
E.
Bard
, and
R.
Schneider
,
2001
:
Environmental processes of the ice age: Land, oceans, glaciers (EPILOG)
.
Quat. Sci. Rev.
,
20
,
627
657
, https://doi.org/10.1016/S0277-3791(00)00145-1.
Muglia
,
J.
, and
A.
Schmittner
,
2015
:
Glacial Atlantic overturning increased by wind stress in climate models
.
Geophys. Res. Lett.
,
42
,
9862
9868
, https://doi.org/10.1002/2015GL064583.
Oka
,
A.
,
H.
Hasumi
, and
A.
Abe-Ouchi
,
2012
:
The thermal threshold of the Atlantic meridional overturning circulation and its control by wind stress forcing during glacial climate
.
Geophys. Res. Lett.
,
39
,
L09709
, https://doi.org/10.1029/2012GL051421.
Otto-Bliesner
,
B. L.
,
C. D.
Hewitt
,
T. M.
Marchitto
,
E.
Brady
,
A.
Abe-Ouchi
,
M.
Crucifix
,
S.
Murakami
, and
S. L.
Weber
,
2007
:
Last Glacial Maximum ocean thermohaline circulation: PMIP2 model intercomparisons and data constraints
.
Geophys. Res. Lett.
,
34
,
L12706
, https://doi.org/10.1029/2007GL029475.
Otto-Bliesner
,
B. L.
, and Coauthors
,
2009
:
A comparison of PMIP2 model simulations and the MARGO proxy reconstruction for tropical sea surface temperatures at Last Glacial Maximum
.
Climate Dyn.
,
32
,
799
815
, https://doi.org/10.1007/s00382-008-0509-0.
Peltier
,
W. R.
,
2004
:
Global glacial isostasy and the surface of the ice-age Earth: The ICE-5G (VM2) model and GRACE
.
Annu. Rev. Earth Planet. Sci.
,
32
,
111
149
, https://doi.org/10.1146/annurev.earth.32.082503.144359.
Pflaumann
,
U.
, and Coauthors
,
2003
:
Glacial North Atlantic: Sea-surface conditions reconstructed by GLAMAP 2000
.
Paleoceanography
,
18
,
1065
, https://doi.org/10.1029/2002PA000774.
Redi
,
M. H.
,
1982
:
Oceanic isopycnal mixing by coordinate rotation
.
J. Phys. Oceanogr.
,
12
,
1154
1158
, https://doi.org/10.1175/1520-0485(1982)012<1154:OIMBCR>2.0.CO;2.
Sarmiento
,
J. L.
, and
J. R.
Toggweiler
,
1984
:
A new model for the role of the oceans in determining atmospheric PCO2
.
Nature
,
308
,
621
624
, https://doi.org/10.1038/308621a0.
Schmittner
,
A.
,
N. M.
Urban
,
J. D.
Shakun
,
N. M.
Mahowald
,
P. U.
Clark
,
P. J.
Bartlein
,
A. C.
Mix
, and
A.
Rosell-Melé
,
2011
:
Climate sensitivity estimated from temperature reconstructions of the Last Glacial Maximum
.
Science
,
334
,
1385
1388
, https://doi.org/10.1126/science.1203513.
Schmittner
,
A.
,
J. A. M.
Green
, and
S.-B.
Wilmes
,
2015
:
Glacial ocean overturning intensified by tidal mixing in a global circulation model
.
Geophys. Res. Lett.
,
42
,
4014
4022
, https://doi.org/10.1002/2015GL063561.
Shakun
,
J. D.
, and Coauthors
,
2012
:
Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation
.
Nature
,
484
,
49
54
, https://doi.org/10.1038/nature10915.
Siegenthaler
,
U.
, and
T.
Wenk
,
1984
:
Rapid atmospheric CO2 variations and ocean circulation
.
Nature
,
308
,
624
626
, https://doi.org/10.1038/308624a0.
Sime
,
L. C.
,
D.
Hodgson
,
T. J.
Bracegirdle
,
C.
Allen
,
B.
Perren
,
S.
Roberts
, and
A. M.
de Boer
,
2016
:
Sea ice led to poleward-shifted winds at the Last Glacial Maximum: The influence of state dependency on CMIP5 and PMIP3 models
.
Climate Past
,
12
,
2241
2253
, https://doi.org/10.5194/cp-12-2241-2016.
Smith
,
W. H. F.
, and
D. T.
Sandwell
,
1997
:
Global sea floor topography from satellite altimetry and ship depth soundings
.
Science
,
277
,
1956
1962
, https://doi.org/10.1126/science.277.5334.1956.
Stammer
,
D.
, and Coauthors
,
2002
:
Global ocean circulation during 1992–1997, estimated from ocean observations and a general circulation model
.
J. Geophys. Res.
,
107
,
3118
, https://doi.org/10.1029/2001JC000888.
Sueyoshi
,
T.
, and Coauthors
,
2013
:
Set-up of the PMIP3 paleoclimate experiments conducted using an Earth system model, MIROC-ESM
.
Geosci. Model Dev.
,
6
,
819
836
, https://doi.org/10.5194/gmd-6-819-2013.
Tao
,
W.
,
L.
Yi
, and
H.
Wei
,
2013
:
Last Glacial Maximum sea surface temperatures: A model-data comparison
.
Atmos. Ocean. Sci. Lett.
,
6
,
233
239
, https://doi.org/10.3878/j.issn.1674-2834.13.0019.
Waelbroeck
,
C.
, and Coauthors
,
2009
:
Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum
.
Nat. Geosci.
,
2
,
127
132
, https://doi.org/10.1038/ngeo411.
Winguth
,
A. M. E.
,
D.
Archer
,
E.
Maier-Reimer
, and
U.
Mikolajewicz
,
2013
: Paleonutrient data analysis of the glacial Atlantic using an adjoint ocean general circulation model. Inverse Methods in Global Biogeochemical Cycles, Geophys. Monogr., Vol. 114, Amer. Geophys. Union,
171
183
.
Wunsch
,
C.
,
2003
:
Determining paleoceanographic circulations, with emphasis on the Last Glacial Maximum
.
Quat. Sci. Rev.
,
22
,
371
385
, https://doi.org/10.1016/S0277-3791(02)00177-4.
Wunsch
,
C.
,
2006
: Discrete Inverse and State Estimation Problems: With Geophysical Fluid Applications. Cambridge University Press, 384 pp.
Wunsch
,
C.
,
2016
:
Last Glacial Maximum and deglacial abyssal seawater oxygen isotopic ratios
.
Climate Past
,
12
,
1281
1296
, https://doi.org/10.5194/cp-12-1281-2016.
Wunsch
,
C.
, and
P.
Heimbach
,
2007
:
Practical global oceanic state estimation
.
Physica D
,
230
,
197
208
, https://doi.org/10.1016/j.physd.2006.09.040.
Wunsch
,
C.
, and
P.
Heimbach
,
2008
:
How long to oceanic tracer and proxy equilibrium?
Quat. Sci. Rev.
,
27
,
637
651
, https://doi.org/10.1016/j.quascirev.2008.01.006.
Zhang
,
X.
,
G.
Lohmann
,
G.
Knorr
, and
X.
Xu
,
2013
:
Different ocean states and transient characteristics in Last Glacial Maximum simulations and implications for deglaciation
.
Climate Past
,
9
,
2319
2333
, https://doi.org/10.5194/cp-9-2319-2013.

Footnotes

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

a

Current affiliation: School of Oceanography and Department of Atmospheric Sciences, University of Washington, Seattle, Washington.

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

1

The method reduces a cost function by search methods. At any stage of the search, estimates of the optimized state can and should be introduced to speed convergence.

Supplemental Material