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.
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.
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
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).
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
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.
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.
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.
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.
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.
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).
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.
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.
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).
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.
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.
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.
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.
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,
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
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,
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).
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-17-0769.1.s1.
Current affiliation: School of Oceanography and Department of Atmospheric Sciences, University of Washington, Seattle, Washington.
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.