A dynamical reconstruction of the Last Glacial Maximum ocean state constrained by global oxygen isotope data

Combining ocean general circulation models with proxy data via data assimilation is a means to obtain estimates of past ocean states that are consistent with model physics as well as with proxy data. The climate during the Last Glacial Maximum (LGM, 19–23 ka) was substantially different from today. Even though boundary conditions are comparatively well known, the large-scale patterns of the ocean circulation during this time remain uncertain. Previous efforts to combine ocean models with proxy data have shown dissimilar results regarding the state of the ocean, in particular of the Atlantic Meridional 5 Overturning Circulation. Here, we present a new LGM ocean state estimate that extents previous estimates by using global benthic as well as planktic data on the oxygen isotopic composition of calcite. It is further constrained by global seasonal and annual sea surface temperature (SST) reconstructions. The estimate shows an Atlantic Ocean that is similar to the Late Holocene Atlantic Ocean but with a reduced formation of Antarctic Bottom Water, in contrast to results of previous studies. The results indicate that SST and oxygen isotopic data alone do not require the presence of a shallower North Atlantic Deep 10 Water and a more extensive Antarctic Bottom Water, and highlight the need for more proxy data of different types to obtain reliable ocean state estimates. Additional adjoint sensitivity experiments reveal that data from the deep North Atlantic and from the global deep Southern Ocean are most important to constrain the Atlantic Meridional Overturning Circulation.

We use a general circulation model enhanced with a water isotope module (Völpel et al., 2017(Völpel et al., , 2018 that enables it to simulate the oxygen isotopic composition of seawater (δ 18 O sw ) in the complete water column. This study investigates the robustness of previous ocean state estimates and the influence of the additional data constraint placed by global surface data and deep-ocean data from outside of the Atlantic Ocean on the estimate, in particular, on the AMOC strength and the vertical extent of the NADW. We additionally perform adjoint sensitivity experiments based on the 400-year estimate that provide the sensitivity of 5 the AMOC with respect to the global ocean temperature and salinity and, therefore, indicate which areas of the global ocean are most important for constraining the AMOC.

Model
We used the coupled ocean-sea-ice Massachussets Institute of Technology general circulation model (MITgcm; Marshall et al., 1997;MITgcm Group, 2016) in the same global configuration as used by Breitkreuz et al. (2018). The configuration uses a cubed-sphere grid (Ronchi et al., 1996) with a horizontal resolution of about 2.8 • and 15 vertical levels. Outgoing radiation, wind stress, and evaporation were computed by bulk formulae (Large and Yeager, 2004) and a GM/Redi scheme (Redi, 1982;Gent and Mcwilliams, 1990) was used to parameterize subgrid-scale mixing. In this study the model was driven by monthly air temperature, meridional and zonal wind velocities, wind speed, specific humidity, precipitation, downward shortwave radiation, downward longwave flux, and river run-off, based on a fully coupled LGM simulation of the Community Climate System Model Version 3 (Merkel et al., 2010). Following Völpel et al. (2018) and Breitkreuz et al. (2019, in review), the Mediterranean Sea was excluded from the model domain, because a shallow passage trough the Strait of Gibraltar was not possible due to the coarse vertical resolution. The isotopic composition of seawater with respect to the Vienna Standard Mean Ocean Water (VSMOW, R VSMOV = 2, 005.2 · 10 −6 , Gonfiantini, 1978) can be 10 computed from the concentration of the isotopes H 16 2 O and H 18 2 O, which are included in the model by a water isotope module (Völpel et al., 2017). The climatological isotopic composition of precipitation and water vapor need to be prescribed and were obtained from a water isotope-enabled LGM simulation with the Community Atmosphere Model version 3.0 (IsoCAM3.0, Tharammal et al., 2013). To apply the adjoint method, the adjoint of the model code needs to be obtained. The MITgcm is tailored to automatic differentiation and the adjoint code can be generated with a source-to-source translator (Giering and and water masses have distinct δ 18 O sw signals (Breitkreuz et al., 2018). The oxygen isotopic composition of foraminiferal calcite preserves the δ 18 O sw and temperature signal of the past ocean. With the use of anomalies, we aim at eliminating species-specific vital effects.
We used the SST anomaly reconstructions by the MARGO Project Members (2009), which are based on a combination of multiple microfossil-based and geochemical proxies. In this study we used the annual as well as the seasonal temperature 5 estimates for July, August, September (JAS) and January, February, March (JFM). The MARGO Project Members (2009) additionally provide estimates of the respective uncertainties of the reconstructed temperature and a mean reliability index.
Following Breitkreuz et al. (2019, in review), we averaged the raw data and their respective uncertainties onto our model grid according to the reliability index as proposed by the MARGO Project Members (2009). In total, 1,120 grid cells were filled with SST data including the seasonal and the annual data. 10 We used the δ 18 O c compilation including uncertainties used by Breitkreuz et al. (2019, in review), which include global δ 18 O c data from planktic as well as from benthic species from different sources. An overview over the δ 18 O c data sets they used can be found in Table 1. Breitkreuz et al. (2019, in review) used the same model grid and averaged the data and their uncertainties onto the model grid. The planktic δ 18 O c data were assigned to the first depth level of the model covering 0-50 m. The depth levels for the benthic data were determined by subtracting 130 m from the core depth to account for the mean sea 15 level change during the LGM (Clark et al., 2009). We used their averaged data set containing 200 grid cells with benthic and 136 grid cells with planktic data. Figure 1 shows the locations of the combined LGM SST and δ 18 O c data from benthic and planktic foraminifera.  Figure 1. Locations of proxy data used in this study: δ 18 Oc data from planktic and benthic foraminifera and annual, JFM, and JAS SST.

Optimization
We used the adjoint method (Wunsch, 1996;Errico, 1997), also called 4-dimensional variational method (4D-Var), or method of Lagrange multipliers, to obtain an estimate of the LGM ocean that is consistent with the observational data in a leastsquares sense as well as with the model physics. The adjoint method minimizes a cost function measuring the model-data misfit by adjusting defined control variables. The control variables can include, for example, initial or boundary conditions, or 5 internal model parameters. Iteratively, the gradient of the cost function close to the current best guess of the control variables is computed from the adjoint model code and used by a quasi-Newton descent algorithm (Nocedal, 1980;Gilbert and Lemaréchal, 1989) to compute adjustments to the control variables.
To create initial conditions for the optimization, the model was first spun-up for 3,000-years using the the original firstguess control variables without the two isotopic tracers and subsequently for another 3,000 years with the isotopic tracers.

10
The adjoint method is typically used to optimize simulations covering 10-50 years (e.g., Köhl et al., 2007;Köhl and Stammer, 2008;Forget et al., 2015) and the longer the simulation, the more difficult it is to achieve a sufficient reduction of the cost function. This might be due to the non-linearity of the model and the increasingly non-convex shape of the cost function with increasing simulation length (Evensen, 2009). To obtain a comparatively long optimized simulation of 400 years, we used a carry-over technique (Dail, 2012;Kurahashi-Nakamura et al., 2017;Breitkreuz et al., 2018), where first, a short simulation 15 is optimized and, subsequently, the optimized control variables are used as a first guess for the optimization of a longer run, eventually reaching the desired length. Following Kurahashi-Nakamura et al. (2017) and Breitkreuz et al. (2018), we chose the atmospheric forcing fields (i.e., air temperature, specific humidity, precipitation, zonal and meridional wind velocities, downward shortwave radiation, downward longwave flux, isotopic composition in precipitation and water vapor), the initial conditions for the physical tracers (salinity and temperature), and the spatially-varying vertical diffusivity as control variables.
The optimization with the carry-over method was started with a 50-year long run followed by a 100-, a 150-, a 200-, and finally a 400-year run. We were not able to obtain a significant reduction of the SST model-data misfit in the first 50-year optimization while using the isotopic control variables. This might be due to the non-convex shape of the cost function and the possibility that the optimization came to a halt in a local minimum. We, therefore, included the isotopic control variables only 5 in the 100-year optimization. Following the 100-year optimization, we excluded the isotopic control variables once again as we already found comparatively big local changes in those control variables.
Following Breitkreuz et al. (2018), the cost function consists of three parts J = J misfit +J ctrl +J eq , which quantify the modeldata misfit (J misfit ), the deviation from the first guess of the control variables (J ctrl ), and the model drift (J eq ). The first term is given by Here, SST sim , SST obs , δ 18 O sim  Table 2). The weighting matrices W δ 18 Oc/SST are the inverse of the error-covariance matrices of the respective proxy data. The uncertainties of the proxy data are assumed to be spatially uncorrelated such that the matrices are diagonal.
The simulated δ 18 O c values were computed from simulated δ 18 O sw and temperature T according to the paleo-temperature 20 equation by Shackleton (1974) Beforehand, an offset of 1.1 ‰ was added to simulated δ 18 O sw to account for the global mean change during the LGM (Duplessy et al., 2002) and the values were transferred from the VSMOW to the Vienna Peedee belemnite (VPDB) standard by subtracting 0.27 ‰ (Hut, 1987). The simulated LGM-LH anomalies were computed using a state estimate of the modern ocean 25 (Breitkreuz et al., 2018) that was obtained by fitting the MITgcm in the same configuration to modern climatological salinity, temperature and δ 18 O sw data with the adjoint method. The use of anomalies enabled us to use one paleo-temperature equation for all species because species-specific vital effects cancel out. Additionally, systematic model errors in the simulated LH and LGM state are assumed to largely neutralize each other. The second part of the cost function has the form It measures the deviation of the adjusted (adj) from the original (orig) control variables weighted by the inverse of the respective error covariance matrices W. The control variables and their prior uncertainties are given in Table 3. The value of this part of the cost function increases during the optimization with growing adjustments to the control variables and, therefore, limits the deviation from the first-guess control variables. This term serves as a regularization of the problem as the number of control variables is substantially bigger than the number of observations and the problem would be highly under-determined without 10 this term. As proposed by Kurahashi-Nakamura et al. (2017), the factor N data · N −1 ctrl is included to balance the contributions of J misfit and J ctrl . The number of control variables is significantly higher than the number of model-data comparisons and a reduction of J misfit would be prohibited by the penalty terms without the balancing factor. Table 3. Control variables, assumed prior uncertainties (one standard deviation, σ), and normalized deviation from the first-guess values after the final phase of the carry-over process (J ctrl ).

Symbol Control variable
Unit Isotopic ratio of water vapor -1.00 · 10 −3 9.3 · 10 −6 K Vertical diffusion coefficient m 2 /s 1.00 · 10 −5 6.8 · 10 −1 Following Breitkreuz et al. (2018), the third term of the cost function is given by where T is the annual mean temperature, η the annual mean global sea surface elevation, and AMOC the annual mean strength of the southward transport at 45 • N in the Atlantic Ocean measured at each year 1, .., t end in the cost function interval. This term penalizes the drift in these three quantities and supports that the obtained state estimate is stable and equilibrated. Following  (Table 4) empirically such that a sufficient reduction of the cost function was still possible. η Annual mean sea surface elevation 9-point spatial smoothing to the adjustments of the control variables, and set minimum and maximum values for precipitation, specific humidity, downward shortwave radiation, air temperature, and the isotopic ratios in precipitation and water vapor (Table 5). The following results refer to the mean over the last 100 years of the final 400-year optimization, which will be termed LGM400. Isotopic ratio in water vapor -2.13 · 10 −3 2.215 · 10 −3

Adjoint Sensitivity Experiments
To investigate the sensitivity of the AMOC to data constraints in different regions of the global ocean, we performed the following adjoint sensitivity experiments. An adjoint simulation with the MITgcm provides the sensitivity of the cost function 15 with respect to all state variables. We implemented a cost function that measures the 10-year mean strength of the AMOC streamfunction at 45 • N to investigate the sensitivity of the AMOC strength to the spatially-varying ocean state. The sensitivity of the AMOC with respect to the ocean state is in turn an indication of how sensitive the AMOC strength is with regard to data constraints in certain regions of the global ocean. We re-ran the 400-year optimized simulation and its adjoint simulation with the AMOC cost function to obtain the sensitivities. 3 Results

Optimization
The optimization with the adjoint method greatly reduced the model-data misfit as measured by the cost function ( Table 2). The first guess simulation showed a good agreement with the SST data, but a very large misfit for the planktic and benthic δ 18 O c data. During the carry-over process the model-data misfit for the SST, planktic, and benthic data were reduced by 58.8 %, 5 95.0 %, and 51.3 %, respectively ( Table 2). The normalized model-data misfit (J misfit ) for the SST has a value below one, which indicates that the optimized simulated SST field agrees on average with the data within their uncertainties. Only single data points show a mismatch in areas where other data points indicate agreement (Fig. 2). The simulated first guess surface δ 18 O c showed a high model-data misfit to the planktic δ 18 O c data in the low latitudes and in the western North Atlantic.
The optimization greatly improves the simulated surface δ 18 O c but some small discrepancies remain, mainly in the North 10 Atlantic (Fig. 3). The simulated deep-ocean δ 18 O c shows an overall agreement with the benthic δ 18 O c data but again some discrepancies remain after the optimization, mostly at shallower depth above 2,000 m (Fig. 4). As the control variables were adjusted during the optimization, the deviation of the control variables from their first guess (J ctrl ) was zero before the optimization and increased during the carry-over process ( Table 2). The normalized value (J ctrl ) stayed well below one, which indicates an average change in the control variables within their assumed uncertainties. Except  (Table 2).

LGM Ocean State
The estimated LGM surface ocean shows an overall cooling compared to the LH (Fig. 2). The global mean cooling in our  The dye tracer experiments provide a similar picture of the water masses in the Atlantic Ocean in LGM400 (Fig. 7). In both 5 the LGM and the LH estimate the center of the northern-source NADW is at a depth of about 2,200 m, but in the LGM estimate the NADW penetrates much deeper and further south than during the LH. Accordingly, the southern-source waters reach much further north in the LH estimate.
The estimated AMOC is only slightly weaker compared to the estimated LH AMOC, and the southward transport reaches a similar depth in both estimates (Fig. 5). The deep counter-clockwise turning cell of the AMOC is stronger in the LH reaching 10 a maximum of 2.1 Sv compared to 0.9 Sv in the LGM estimate. To summarize, LGM400 provides an estimate of the LGM ocean that is similar to the modern one, except for an AABW that is much smaller in its extent.

Adjoint AMOC Sensitivities
The adjoint sensitivities in Figs

LGM Ocean State Estimate
According to a number of previous studies (e.g., Duplessy et al., 1988;Curry and Oppo, 2005;Lynch-Stieglitz et al., 2007;10 Muglia et al., 2018) the deep Atlantic Ocean was filled with a southern-source deep water mass overlaid by a shallow NADW during the LGM. With the adjoint method, we successfully found an LGM ocean state that is consistent with the employed data as well as with the physics of our model. However, in contrast to previous studies, our estimate shows a deep NADW and only a very weak extension of the southern-source AABW.
The modern Atlantic Ocean is dominated by a relatively salty NADW with waters originating from the tropics. The very 15 deep southern Atlantic Ocean is filled with the fresher but colder AABW. For the LGM Atlantic Ocean it has been proposed that the AABW was much saltier compared to today due to increased sea ice formation and brine rejection, and that the glacial  NADW was likely less salty relative to the glacial AABW (Adkins et al., 2002). Recently, however, a critical assessment of these results showed that while the proposed ocean state was possible, it was not required by the data (Miller et al., 2015;Wunsch, 2016). The global deep ocean is thought to have been overall substantially colder than the modern ocean (Adkins et al., 2002).
In our estimate these expectations are only partly met. The AABW is slightly saltier relative to the NADW (Fig. 6) and the 5 deep ocean is mostly cooler, however, not as much as previously indicated (Adkins et al., 2002). Additionally, the negative temperature anomaly of the NADW is about 2 • C larger compared to the AABW (Fig. 6). Relative higher salinity and colder temperatures make the density change of the NADW in our estimate larger compared to the density change in the AABW.
This coincides with the surface density changes in our LGM estimate, which show higher positive values for the region south of Iceland where deep water is formed in our LGM estimate, but a smaller increase in density in the Southern Ocean in the 10 Weddell Sea (Fig. 10). These density changes are, for example, in contrast to the results of Paul and Schäfer-Neth (  AABW into the Atlantic Ocean is very limited in our estimate. Additionally, the week counter-clockwise turning AMOC cell might contribute to the small northward extent of the AABW. The AMOC in our estimate shows a very similar strength and depth compared with the modern ocean (Ganachaud, 2003;Breitkreuz et al., 2018). Previous studies using models Weber et al., 2007;Muglia and Schmittner, 2015), a combination of models and proxy data via data assimilation (Kurahashi-Nakamura et al., 2017;Breitkreuz et al., 5 2019, in review), or 231 Pa / 230 Th data (McManus et al., 2004;Gherardi et al., 2005;Lippold et al., 2016) targeting the LGM AMOC strength have shown dissimilar results. Most recent studies inferred a weaker AMOC from the best fit of a physicalbiogeochemical ocean model to δ 13 C, radiocarbon, and δ 15 N data (Muglia et al., 2018) and from 231 Pa / 230 Th and neodymium isotopes (Lippold et al., 2016;Howe et al., 2016).  the AABW and colder temperatures close to the freezing point in the deep Atlantic Ocean (e.g., Adkins et al., 2002;Curry and Oppo, 2005).

Data Constraint
Above approximately 2,000 m the benthic δ 18 O c anomalies show some remaining model-data misfit (Fig. 4) The issue that δ 18 O c is influenced by temperature and δ 18 O sw and that the effects of both tracers might cancel each other out has been previously acknowledged (e.g., Paul et al., 1999) and exacerbates the constraint of the circulation through only δ 18 O c data. A likely reason for a very salty LGM AABW is increased sea ice formation (Adkins et al., 2002), but again, δ 18 O c 15 data is not a good constraint for these salinity changes because δ 18 O sw is not as much influenced by sea ice formation and brine rejection as salinity is. Salinity estimates derived from pore-water measurements would be crucial in this regard but are still sparse and have high uncertainties (Wunsch, 2016 Similarly, the AMOC strength seems to be not well constrained in our estimate. During the carry-over process the estimated AMOC strength varies between 11.0 Sv and 21.3 Sv (Table 2) 2017) additionally used benthic carbon isotopic data, whose constraint is 10 missing in our estimate. However, they did not use any planktic isotopic data, such that the simulated surface values for δ 13 C and δ 18 O c were not constrained by proxy data. Additionally, the benthic data they used were confined to the Atlantic Ocean.
Their optimized control variables were also used with a similar model configuration to obtain an equilibrium LGM simulation that showed a shallower and weaker AMOC, again highlighting the influence of the simulation length or of slight changes in the model configuration on the result (Völpel et al., 2018). 15 Carbon isotopic data has shown to place a stronger constraint on the ocean circulation than oxygen isotopic data in previous studies (e.g., Marchal and Curry, 2008). The carbon isotopic composition in the deep ocean, however, depends on the remineralization of organic carbon. The implementation of carbon isotopes into the MITgcm, including 14 C, will be an important next step and allow to assimilate δ 13 C and 14 C data, which will likely place a better constraint on the LGM ocean circulation. Breitkreuz et al. (2019, in review) obtained an LGM ocean state estimate from the same observational data and the same 20 model version, but their approach was based on a state reduction approach for the control variables and a Kalman smoother method. Their estimate shows a shallower NADW and a weaker AMOC with a maximum strength of 13.8 Sv. However, in their estimate a substantially bigger model-data misfit remained and they found a second ocean state with similar model-data misfit but without an active overturning circulation in the Atlantic Ocean. The differences in the estimates suggest that results from ocean state estimation highly depend on the assimilation approach, the choice of control variables, the remaining model-data Additionally, the adjoint method as well as the Kalman smoother method used by Breitkreuz et al. (2019, in review) are highly first-guess dependent.

30
Additional to the type of proxy data, the location of the proxy data can be essential for the success in constraining the ocean circulation. Figures 8 and 9 indicate how sensitive the AMOC is to changes in potential temperature and salinity in different regions of the global ocean and, therefore, provide information on which areas are most important to be constrained by proxy data. The adjoint sensitivities for δ 18 O sw do not hold any information because δ 18 O sw is not an active tracer, that is, a change in δ 18 O sw does not have an influence on the circulation or the AMOC. Nevertheless, the sensitivities to temperature and salinity changes indicate in which regions a constraint through different proxy data has the biggest effect.
Most important seem to be data in the North Atlantic Ocean and in the deep Southern Ocean. The North Atlantic region, where deep water is formed due to the density gradient between surface and deep water, is crucial for the strength of the AMOC and hence, high sensitivities are visible for the surface as well as for the deep ocean in this region. While there is a good coverage of SST reconstructions for the North Atlantic (Fig. 1), the coverage of planktic and benthic isotopic data is not as extensive and more isotopic data from this region might place a better constraint on the simulated density gradient and improve future ocean state estimates.
The second region with highest sensitivities is the Southern Ocean. The Southern Ocean links all three major oceans via the Antarctic Circumpolar Current (ACC). The density of the water in the Southern Ocean in the Atlantic region is, therefore,  (Fig. 1). 15 Regarding the surface ocean, a constraint of the Atlantic Ocean seems to be more important than of the Indian or Pacific Oceans, supporting the results of Breitkreuz et al. (2019, in review). Whereas they could not determine the importance of deep-ocean data from outside the Atlantic Ocean from their results, the adjoint sensitivities specifically highlight the need for more benthic data in the southern part of the global ocean.
In general, sensitivities are highest 150 years before the AMOC is computed. In this and previous LGM ocean state estimates 20 the cost function covered only the last 5 years (Dail andWunsch, 2014), 10 years (Kurahashi-Nakamura et al., 2017), 60 years (this study), or 100 years (Breitkreuz et al., 2019, in review) of the optimized simulation. The adjoint sensitivities indicate that an extension of the cost function to a longer time interval might help in obtaining better constrained ocean state estimates.
An extension of the cost function in an adjoint simulation, however, creates high recomputation and storage costs because the required state variables in each time step during the cost function interval are needed for the adjoint calculation.

Conclusions
We presented a new ocean state estimate that was obtained from an oxygen isotope-enabled general circulation model and a successful application of the adjoint method. The adjoint method provides a solution that is consistent with the employed data and the model physics. The estimate extends previous state estimates by using benthic as well as planktic data on the oxygen isotopic composition of calcite. It is additionally constrained by global annual and seasonal SST reconstructions. It shows an 30 LGM Atlantic Ocean that is dominated by a northern-source deep water mass corresponding to the modern NADW. The deep southern-source water mass is less extensive than during the LH ocean, in contrast to previous studies. The unexpected results might be due to an insufficient constraint of the available oxygen isotopic and SST data. Comparisons with previous ocean state study will be submitted to PANGAEA (https://www.pangaea.de/) by the authors of Breitkreuz et al. (2019, in review).