Abstract

The Bern3D coupled three-dimensional dynamical ocean–energy balance atmosphere model is introduced and the atmospheric component is discussed in detail. The model is of reduced complexity, developed to perform extensive sensitivity studies and ensemble simulations extending over several glacial–interglacial cycles. On large space scales, the modern steady state of the model compares well with observations. In a first application, several 800 000-yr simulations with prescribed orbital, greenhouse gas, and ice sheet forcings are performed. The model shows an increase of Atlantic meridional overturning circulation strength at glacial inceptions followed by a decrease throughout the glaciation and ending in a circulation at glacial maxima that is weaker than at present. The sensitivity of ocean temperature to atmospheric temperature, Atlantic meridional overturning circulation (AMOC), and Antarctic bottom water (AABW) strength is analyzed at 23 locations. In a second application the climate sensitivities of the modern and of the Last Glacial Maximum (LGM) state are compared. The temperature rise for a doubling of the CO2 concentration from LGM conditions is 4.3°C and thus notably larger than in the modern case (3°C). The relaxation time scale is strongly dependent on the response of AABW to the CO2 change, since it determines the ventilation of the deep Pacific and Indian Ocean.

1. Introduction

With increasing computer power, the realism of climate models could be increased and thus models have become more complex. Spatial resolution has been refined to better resolve small-scale phenomena, and more processes have been included. Nonetheless, computer power still limits the most complex type of models, the coupled atmosphere–ocean general circulation models (AOGCMs), from permitting simulations exceeding a few hundred years [e.g., the National Center for Atmospheric Research (NCAR) Community Climate System Model (CCSM; Collins et al. 2006) and third climate configuration of the Met Office Unified Model (HadCM3; Gordon et al. 2000)]. At the time of writing, Liu et al. (2009) have presented the longest simulations with an AOGCM of almost 8 kyr (1 kyr = 1000 yr) using CCSM3. Otherwise, for simulations on a multimillenial time scale, so-called earth system models of intermediate complexity (EMICs) have been built. They are typically based on simplified physics and parameterizations of a larger number of processes [e.g., the Bern2.5D global model (Stocker et al. 1992), the University of Victoria (UVic) earth system climate model (Weaver et al. 2001), ECBilt-CLIO (Goosse et al. 2005), the Climate and Bisphere Model version 3α (Climber3α) (Montoya et al. 2005), and the Grid Enabled Integrated Earth system model (GENIE) (Edwards and Marsh 2005)].

To the present day, only few model simulations have been done spanning more than one glacial cycle of approximately 100 000 years. Those studies used models with zonally averaged ocean basins (Tuenter et al. 2005) or with a three-dimensional atmosphere but a slab ocean component and accelerated variations of the orbital configuration (Jackson and Broccoli 2003). Here we present a very cost-efficient coupled three-dimensional dynamical ocean–energy balance atmosphere intermediate complexity model, which permits simulations on glacial-to-interglacial time scales. Currently 50 000 model years per day can be run on a single personal computer CPU. With this, the model is considerably more efficient than most three-dimensional EMICs. Thus, extensive long time-scale parameter sensitivity studies or ensemble simulations become feasible. Because the model also includes a prognostic formulation of the carbon cycle and a palette of other tracers, it is a powerful tool for comprehensive paleoceanographic model studies.

In this paper, we introduce the atmospheric component of the coupled model in detail, present the modern steady state of the coupled ocean–atmosphere model, and perform several coupled 800 000-yr simulations with prescribed orbital, greenhouse gas, and ice sheet forcings.

2. Model description

a. The ocean model component

The Bern3D ocean component is a seasonally forced three-dimensional frictional geostrophic global ocean model with coarse spatial resolution. It consists of 36 × 36 grid boxes in the horizontal direction and 32 vertical layers. The year is divided into 48 time steps, corresponding to about a week per time step. It is based on the ocean model of Edwards et al. (1998) and described in detail by Müller et al. (2006). A new feature is the possibility of barotropic flow around the American continent and Australia. In the modern control state, there is 0.5 Sv (1 Sv = 106 m3 s−1) northward flow through the Bering Strait and 23 Sv Indonesian Throughflow from the Pacific to the Indian Ocean. In ocean-only simulations, the model is run under restoring surface boundary conditions for temperature and salinity. Temperature fields are taken from Levitus and Boyer (1994) and salinity fields from Levitus et al. (1994). The model also contains a prognostic carbon cycle (Parekh et al. 2008; Tschumi et al. 2008). The Bern3D ocean component has been used for a range of applications (Gerber and Joos 2010; Gerber et al. 2009; Ritz et al. 2008; Parekh et al. 2008; Tschumi et al. 2008; Müller et al. 2008; Siddall et al. 2007; Muscheler et al. 2007).

b. The energy balance model of the atmosphere

The single-layer energy balance is described in spherical coordinates using ϕ ∈ {0; 2π} for the longitude and ϑ ∈ {−π/2; π/2} for the latitude and is similar to the model described by Weaver et al. (2001). The spatial and temporal resolutions are equal to the resolution of the ocean model. Notation and values of the model parameters are given in Table 1. Depth-integrated horizontal heat fluxes are parameterized in terms of eddy-diffusive fluxes with uniform zonal, Kϕ, and meridional, Kϑ, diffusivities. The vertical energy fluxes consist of shortwave (sw) and longwave (lw) fluxes at the top of the atmosphere (TOA) and across the atmosphere–ocean (AO), atmosphere–sea ice (AI), and atmosphere–land (AL) boundaries, respectively:

 
formula

where Ta is the atmospheric temperature. For model stability reasons, horizontal advective transport is not taken into account. The fluxes at the bottom of the atmosphere (BOA) are the heat gains of the ice-free ocean, sea ice, and land surface, respectively. They are separated into a shortwave and longwave radiation term, a sensible heat flux, and a latent heat flux term:

 
formula
 
formula

(the emissivities of water and ice are very similar for infrared wavelengths), and

 
formula

where To is the surface ocean temperature, Ti the surface sea ice temperature, and Tl the surface temperature over land. Note that since evaporation is included in and in , it needs to be subtracted in Eq. (1). In this version, water is not stored on land, and therefore evaporation is zero on land boxes. Land temperatures are calculated by solving

 
formula

The land surface scale height hl is chosen to be 2 m. This corresponds to the depth to which temperature is affected by seasonality (Hartmann 1994). For reasons of simplicity, except for Antarctica (AA), the parameters for the land surface are chosen to be global and correspond to sandy, saturated soil (Martin 2002). For Antarctica, . These choices help to decrease atmospheric temperature seasonality in this region.

Table 1.

Parameter values for energy and moisture transport in the atmosphere and for the sea ice model.

Parameter values for energy and moisture transport in the atmosphere and for the sea ice model.
Parameter values for energy and moisture transport in the atmosphere and for the sea ice model.

Incoming solar radiation is calculated following the algorithm of Berger (1978). A parameterization of Bintanja (1996) determines how much of the incoming radiation is transmitted through the atmosphere to the bottom of the atmosphere , how much of it is reflected back into space , and how much is absorbed by the atmosphere:

 
formula

where are the radiation fluxes transmitted through the atmosphere for clear-sky and overcast conditions, respectively; ξ(ϑ, n) denotes the zonally averaged fractional cloud amount climatology taken from the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40), with n being the model time step within a year; and α(ϕ, ϑ, n) denotes the surface albedo. Analogously,

 
formula

Here include approximations for the absorptive and reflective properties of the atmospheric constituents, the solar zenith angle and the surface elevation [taken from 5-minute gridded elevations/bathymetry for the world (ETOPO5); see http://www.ngdc.noaa.gov/mgg/global/etopo5.html] and are calculated as described by Bintanja (1996). However, the following change has been made: In contrast to Bintanja (1996), cloud optical depth τ, a measure of cloud transparency, is not set to a constant value but rather depends on the liquid water path, W = qaρahq, the integrated amount of water in the atmospheric column:

 
formula

with W0 = 1 kg m−2 (Stephens 1978); qa is the surface specific humidity. It satisfies a balance equation (see below).

Following Weaver et al. (2001), the parameterization for outgoing planetary infrared irradiance for clear-sky conditions at TOA of Thompson and Warren (1982) is used and extended by a parameterization for the radiative forcing owing to deviations of atmospheric CO2 concentrations from a reference value. Additionally, a simplified term for CH4 greenhouse gas forcing is added, as well as a term representing the water vapor feedback:

 
formula

where pCO2,0 = 278 ppm is the preindustrial atmospheric carbon dioxide concentration, ΔF2×CO2 = 5.35 W m−2 (Myrhe et al. 1998), pCH4,0 = 700 ppb is the preindustrial atmospheric methane concentration, and υ = 0.036 W m−2 (Shi 1992). The coefficients ai depend on the relative humidity rh according to

 
formula

The coefficients bj,i are taken from Thompson and Warren (1982); rh = qa/qs(Ta), where qs(Ta) is the saturation specific humidity (explained below). The last term of Eq. (9) is a very simple approximation for the water vapor feedback. We define the seasonally dependent temperature deviation from the modern control state ΔTa = Ta(n, y) − TaCTRL(n), where Ta(n, y) is the global mean atmospheric temperature at model time step n of year y; TaCTRL is the temperature average of a 5-kyr control run. The feedback parameter λ is tuned (λ = 1 W m−2 K−1) to produce an equilibrium climate sensitivity (global temperature rise for a doubling of the atmospheric CO2 content) of 3°C for a modern steady state. It is found that the effect of clouds on the longwave radiation needs to be accounted for, especially at high latitudes, where cloud cover is high compared to the global mean value. Thus, the following simple parameterization is implemented into the model: Since for the presence of clouds the location of emitted longwave radiation is the cloud-top level instead of the earth’s surface, the graybody radiation temperature is lower and thus also the outgoing radiation (Hartmann 1994). Therefore, we reduce the graybody radiation temperature of Eq. (9) to

 
formula

depending on fractional cloud cover ξ. Here Tact expresses the temperature at cloud top and the temperature reduction when ξ = 1; is chosen such that a reasonable global atmospheric temperature is obtained.

Following Weaver et al. (2001), evaporation at the ocean surface EAO is calculated according to

 
formula

where |u| is the surface wind speed at 10-m height from ERA-40 reanalysis, and CE = CE(ϕ, ϑ, n) is the Dalton number. It is diagnosed during a 50-yr initialization run by solving Eq. (12) for CE and using monthly evaporation fields from ERA-40 reanalysis as well as a Ta climatology from ERA-40 and ocean temperatures from the ocean-only simulation. As proposed by Isemer et al. (1989), 6.0 × 10−5CE ≤ 2.19 × 10−3. Also, qs(To) is the saturation specific humidity at the ocean surface (g water per kg air). It is calculated using the parameterization of Bolton (1980):

 
formula

where c1 = 3.80 g kg−1, c2 = 17.67, and c3 = 243.5 K; Ts is either ocean temperature To or the atmospheric temperature Ta. Sublimation over sea ice EAI is parameterized as evaporation [Eq. (12)], except for replacing surface ocean saturation specific humidity qs(To) by sea ice surface saturation specific humidity

 
formula

where c1 = 3.80 g kg−1, c4 = 21.87, and c5 = 265.5 K.

Precipitation occurs when relative humidity rises above a maximum value rh,max. Precipitation stops when rh is equal to rh,precip:

 
formula

A fraction of the precipitation over land is instantly transported to the ocean as runoff. Its pathway is defined by a pseudoelevation map: Each land box contains information about the direction of runoff (north, south, west, or east). This fraction χ is taken from observations (Hartmann 1994) and is different for every continent. It ranges from 0.33 to 0.43, except for Africa, where reevaporation is high (χ = 0.16), and Antarctica, where the ice sheet is in broad equilibrium and the mass increase by snowfall is compensated by snowmelt and iceberg calving at the ocean margin (χ = 0.83). The rest, which in the real world would be stored on land or reevaporated to the atmosphere, is distributed uniformly into every surface ocean box. Note that it would physically make more sense to redirect (1 − χ) · P back to the atmosphere, since this is the fraction that is taken up by the soils and eventually reevaporated. However, because of the coarse resolution and the parameterization of the model, a large part of the reevaporated moisture would be precipitated again in the following time step. This would strongly increase the atmosphere–ocean moisture turnover and thus the amount of runoff. Evaporation and precipitation in m s−1 are converted into a latent heat flux by multiplying the reference density of water ρo and the latent heat of evaporation Lυ.

At the atmosphere–ocean and atmosphere–sea ice interface the parameterization of Weaver et al. (2001) for the sensible heat flux is used:

 
formula

where CH = 0.94CE is the Stanton number and Ts is either surface ocean or sea ice temperature. Sensible heat fluxes over land surface are parameterized as

 
formula

using a constant bulk coefficient Dl (Martin 2002).

In contrast to the energy balance Eq. (1), moisture is transported by diffusion and advection. Meridional advection is important in the tropical regions where moisture is transported equator ward by intertropical convergence (ITC), where precipitation occurs. In a diffusive-only scheme, moisture would not converge but diverge in this region. Zonally averaged monthly wind velocity fields from ERA-40 reanalysis are applied (note that zonally resolved winds would lead to convergence in various boxes in the mid and high latitudes and thus, as a direct effect, to an unrealistically high amount of precipitation. These boxes negatively affect the state of the model. In the three-dimensional real world, convergence only leads to precipitation when the rising air masses cool sufficiently). The wind fields are vertically density weighted and averaged up to the moisture scale height. The vertical fluxes are given by evaporation and precipitation. The moisture balance equation is formulated as follows:

 
formula

Eddy diffusivities Kϕq and Kϑq are assumed to be globally constant.

The sea ice model component is based on work by Semtner (1976) and Hibler (1979) and is similar to the sea ice model of Edwards and Marsh (2005). The model determines three variables: fractional sea ice area Ai, ice thickness Hi, and surface sea ice temperature Ti. Note that Hi is averaged over the ice and the open-ocean fractions. Ice dynamics are kept very simple: Ice flows with the surface ocean currents, and processes induced by horizontal gradients are parameterized by diffusion with a diffusivity Ki.

Vertical heat fluxes are separated into the atmosphere–ocean heat flux across ice-free areas, [Eq. (2)], the atmosphere–ice heat flux across ice-covered areas, [Eq. (3)], and the ice-ocean heat flux

 
formula

which brings To back to the freezing temperature Tf by either melting or growing ice. Note that Tf is salinity dependent and is parameterized as

 
formula

(Doronin and Kheisin 1977). Also, uτ is the skin friction velocity at the ice–ocean boundary, and ch is an empirical constant after McPhee (1992). The parameter values are given in Table 1. The total heat flux from the atmosphere is , where Ai is the ice-cover fraction.

Fractional ice area satisfies a balance equation that consists of a horizontal flow term, an ice area production term, and an ice area destruction term. The ice produced in the open-ocean area is uniformly spread using a minimal thickness H0. In the ice area destruction term, the ice is assumed to be distributed uniformly between height 0 and 2Hi/Ai over the ice-covered fraction Ai (Fig. 1). Following simple intercept theorems, the formed open-ocean fraction dAi can be derived from

 
formula

where dHi is the thickness of the melted ice layer. Hence,

 
formula

Thus, fractional sea ice area is calculated by solving the following equation:

 
formula

Similarly,

 
formula

If the new ice thickness is smaller than a minimal ice thickness H0, then Hi and Ai are set to zero. The surface temperature of the ice is calculated by equating the radiative incoming heat flux with the conductive heat flux through the ice, thus assuming a linear temperature profile within the ice:

 
formula

where Icond is the thermal conductivity of ice. In the model, Ti = min(Ti, Tf).

Fig. 1.

Outline of how sea ice is distributed for the ice-melting term in Eq. (23); Ai is the ice-covered fraction of the box, Hi the height of the ice, dHi the melted ice layer, and dAi the ice fraction that is melted away.

Fig. 1.

Outline of how sea ice is distributed for the ice-melting term in Eq. (23); Ai is the ice-covered fraction of the box, Hi the height of the ice, dHi the melted ice layer, and dAi the ice fraction that is melted away.

A very simple parameterization for ocean and sea ice albedo is used:

 
formula

The coefficients are chosen to match values of ice-free and fully ice-covered ocean areas of the Kukla and Robinson (1980) ocean–sea ice–albedo climatology. Over land, the zonally averaged land–albedo climatology of Kukla and Robinson (1980) is used.

Note that the presence of sea ice requires the sensible heat flux Fsh and evaporation E to be separated into an ice-covered and an open-ocean fraction. Finally, the heat flux into the ocean is calculated as

 
formula

As in Eq. (23), ice is formed over the open-ocean fraction when . In this case the released heat of fusion is considered in the first term of Eq. (27). Also, Qm is the heat of fusion of the additional amount of meltwater that is added to the ocean when the ice thickness falls below the minimal thickness (Hi < H0 but in the previous time step Hi,t−Δt > 0). The freshwater flux is calculated as

 
formula

where R is runoff and Qm/(Lfρo) is the additional amount of freshwater added to the ocean when Hi < H0 but Hi,t−Δt > 0. In the model, not freshwater but salt is added and taken out of the ocean, respectively: , with Sref = 34.78 psu being a reference salinity for the surface ocean.

Because of the absence of dynamics in the atmosphere, an Atlantic-to-Pacific freshwater flux (Zaucker et al. 1994) must be prescribed. We apply 0.17 Sv to increase the strength of the Atlantic meridional overturning circulation (AMOC). The freshwater is taken from the North Atlantic from the 30° to 71°N basin and distributed into the North Pacific from 46° to 71°N.

An additional freshwater flux out of the ocean is applied to the Ross Sea and the Weddell Sea (two boxes each) in order to stimulate deep water formation in these regions. This correction is due to the fact that the small-scale processes in these regions are not resolved by the model. We apply a flux of 0.2 Sv divided into these four boxes. This flux correction is compensated by adding the same freshwater amount to the remaining ocean boxes around Antarctica (63°–71°S).

The horizontal transport term in the energy and moisture balance Eqs. (1) and (18) is solved implicitly. Thus, the time step for the energy balance model (EBM) can be chosen equal to the ocean time step. Since the velocities and diffusivities in Eqs. (1) and (18) do not vary interannually, the linear systems of equations need to be inverted only during the first year of every simulation. This makes the EBM very efficient. Since the sea ice has low flow speeds and diffusion, it is not solved implicitly. We halved the time step for the land temperature [Eq. (5)] for numerical stability reasons. Equation (25) cannot be solved analytically. Note that Ti4 is linearized and discretized using a first-order Taylor approximation so that , where Ti,t−Δt is the temperature at the previous time step.

Atmospheric temperature values of the boxes closest to the poles (poleward of 71°) are averaged longitudinally after every time step to increase numerical stability. The discretization schemes Euler forward, centered differences, and variable upwind (for zonal advection) are used.

3. Present-day simulations

a. Ocean

The parameters described in the model section—particularly the relative humidity after precipitation rh,precip, zonal and meridional eddy diffusivities Kϕ and Kϑ, the temperature reduction at cloud top for overcast conditions , the Dalton number CE, and the freshwater correction fluxes from the North Atlantic to the Pacific and in the Ross and Weddell Seas—have been tuned such that a good representation of the modern climate is achieved. The model tuning was done on the basis of observational fields and Taylor diagrams. Special emphasis was placed on atmospheric and surface ocean temperature, sea surface salinity, sea ice cover, Atlantic and Pacific zonal mean temperature, salinity, and radiocarbon concentration. Global relative standard deviations and correlations between the mentioned quantities and observations were calculated and optimized.

A steady state of the coupled atmosphere–ocean model is obtained by spinning up the ocean-only model for 10 kyr, followed by a 50-yr EBM initialization run, where evaporation patterns of ERA-40 reanalysis are approached by diagnosing the Dalton number CE seasonally at every grid point. Finally, ocean and atmosphere are coupled and a follow-on 10-kyr spinup is performed. The result is a stable and steady model state: Global mean net TOA radiation fluxes converge to zero. The 100-yr average of globally integrated radiation fluxes at TOA is 0.005 PW, which corresponds to an average flux of 0.010 W m−2.

The modern steady-state annual mean overturning circulation of the Atlantic and Pacific basin and of the global ocean are shown in Fig. 2. North Atlantic Deep Water (NADW) reaches down to 3–4-km depth before flowing southward. Because of the coarse resolution of the model, NADW is formed south of Greenland, one box row south of the Greenland–Iceland–Norwegian (GIN) Seas, where deep water should be formed. The AMOC strength with a maximum of approximately 14 Sv is low compared to other models (Randall et al. 2007). Observations of the Atlantic radiocarbon content (Key et al. 2004) are, however, consistent with the AMOC strength (Fig. 3). The global overturning shows the Southern Ocean overturning cell with a strength of approximately 18 Sv. Again, this strength leads to a deep Pacific radiocarbon concentration that compares well with the observations (Fig. 4).

Fig. 2.

Modern annual mean Atlantic, Pacific, and global overturning circulation (in Sv). The Indonesian Throughflow at 10°S in the Pacific (gray bar) produces a discontinuity in the contour lines.

Fig. 2.

Modern annual mean Atlantic, Pacific, and global overturning circulation (in Sv). The Indonesian Throughflow at 10°S in the Pacific (gray bar) produces a discontinuity in the contour lines.

Fig. 3.

Atlantic zonally and annually averaged latitude–depth fields of (top) ocean temperature, (second row) salinity, (third row) natural radiocarbon (bomb-produced radiocarbon is filtered out in the observations), and (bottom) phosphate concentrations. (left) Model results are compared to (right) observations from Levitus and Boyer (1994) for temperature, Levitus et al. (1994) for salinity, GLODAP (Key et al. 2004) for radiocarbon, and the World Ocean Atlas 2001 (WOA01; Conkright et al. 2002) for phosphate.

Fig. 3.

Atlantic zonally and annually averaged latitude–depth fields of (top) ocean temperature, (second row) salinity, (third row) natural radiocarbon (bomb-produced radiocarbon is filtered out in the observations), and (bottom) phosphate concentrations. (left) Model results are compared to (right) observations from Levitus and Boyer (1994) for temperature, Levitus et al. (1994) for salinity, GLODAP (Key et al. 2004) for radiocarbon, and the World Ocean Atlas 2001 (WOA01; Conkright et al. 2002) for phosphate.

Fig. 4.

As in Fig. 3, but for the Pacific.

Fig. 4.

As in Fig. 3, but for the Pacific.

Zonally and annually averaged latitude–depth plots of ocean temperature, salinity, radiocarbon, and phosphate distributions are shown for the Atlantic (Fig. 3) and for the Pacific (Fig. 4). Atlantic temperatures agree well with observations of Levitus and Boyer (1994). The Pacific below 1-km depth is about 1°C too cold. The salinity fields in both Atlantic and Pacific are too fresh at the surface and too salty at depth. This deficiency is possibly linked to the distribution and thus to the parameterization of evaporation and precipitation in the atmosphere.

Besides radiocarbon, constraints on the state of the overturning circulation can be inferred from distributions of nutrients and other biogeochemical tracers. Therefore, the model is run with the prognostic carbon cycle, allowing us to compare the model output to observations from the World Ocean Atlas 2001 (WOA01; Conkright et al. 2002) for phosphate and silicate, and the Global Data Analysis Project (GLODAP) data (Key et al. 2004) for dissolved inorganic carbon, alkalinity, and chlorofluorocarbon (CFC-11). The modeled phosphate distribution is in fair agreement with the observations. The largest deficiencies are found in the Pacific, where surface concentrations are too high and the surface-to-deep gradient is too small. In the Atlantic, the highest observational values are found in the equatorial upwelling region of the African Margin (Gulf of Angola). In the model, the highest values are also found at the African Margin, but the region extends to the Gulf of Guinea. Therefore, a sharp change to lower values is found in the zonally averaged field (Fig. 3). The phosphate accumulation in the deep ocean of the African Margin is a consequence of underestimated ventilation of the deepest boxes. For a quantitative comparison between the model and observations, a Taylor diagram is constructed providing information about the relative standard deviation and the correlation between annually averaged model and data fields (Fig. 5a). For CFC-11, a transient simulation was performed prescribing atmospheric concentrations starting at year 1931 A.D. and ending at year 2000. Because the observational data from GLODAP (Key et al. 2004) were collected during multiple years of the 1990s, the data of every cruise station were interpolated and regridded on to the model depth grid and then compared to the corresponding year of the model output. Additionally, the model results for temperature, salinity, and radiocarbon are compared to the results of the ocean-only simulation described by Müller et al. (2006). It is expected that the ocean-only simulation performs better than the coupled simulation because the surface–ocean boundary is restored to observations, while in the coupled run here the surface–ocean conditions follow from the energy balance model. Thus, it is remarkable that the temperatures of the coupled run agree nearly as well with the observations as those of the ocean-only simulation. Correlation and relative standard deviation of the salinity, however, is quite poor. As described earlier, this is primarily due to the poorly represented surface-to-deep gradient of salinity. On the other hand, radiocarbon of the coupled run compares better with observations, indicating that the time scales of surface-to-deep transport are realistic in the model.

Fig. 5.

(a) Taylor diagram showing correlation and relative standard deviation of modeled annual mean ocean temperature, salinity, and various tracer fields compared to observations. The filled symbols and the plus sign (+) represent the coupled model; the outlined symbols indicate the ocean-only model, as in Fig. 8 of Müller et al. (2006). Observations are taken from Levitus and Boyer (1994) and Levitus et al. (1994) for temperature and salinity, respectively; GLODAP (Key et al. 2004) for radiocarbon, chlorofluorocarbon CFC-11, dissolved inorganic carbon (DIC) and alkalinity; and WOA01 (Conkright et al. 2002) for phosphate and silicic acid. Model results for CFC-11 were determined from a transient simulation, in which atmospheric concentrations were prescribed from years 1931 to 2000 A.D. Perfect agreement with the data is at point (1, 0) in the Taylor diagram. (b) Correlation and relative standard deviation of modeled atmospheric temperature, evaporation and precipitation with ERA-40 reanalysis data. Black symbols indicate the annual mean; gray symbols, January fields; outlined symbols, July fields. Only evaporation over the ocean is compared. Symbols for annual mean and January precipitation are on top of each other.

Fig. 5.

(a) Taylor diagram showing correlation and relative standard deviation of modeled annual mean ocean temperature, salinity, and various tracer fields compared to observations. The filled symbols and the plus sign (+) represent the coupled model; the outlined symbols indicate the ocean-only model, as in Fig. 8 of Müller et al. (2006). Observations are taken from Levitus and Boyer (1994) and Levitus et al. (1994) for temperature and salinity, respectively; GLODAP (Key et al. 2004) for radiocarbon, chlorofluorocarbon CFC-11, dissolved inorganic carbon (DIC) and alkalinity; and WOA01 (Conkright et al. 2002) for phosphate and silicic acid. Model results for CFC-11 were determined from a transient simulation, in which atmospheric concentrations were prescribed from years 1931 to 2000 A.D. Perfect agreement with the data is at point (1, 0) in the Taylor diagram. (b) Correlation and relative standard deviation of modeled atmospheric temperature, evaporation and precipitation with ERA-40 reanalysis data. Black symbols indicate the annual mean; gray symbols, January fields; outlined symbols, July fields. Only evaporation over the ocean is compared. Symbols for annual mean and January precipitation are on top of each other.

b. Surface ocean and atmosphere

Annual mean atmospheric and sea surface temperature (SST), sea surface salinity (SSS), evaporation, and precipitation are shown in Fig. 6 and compared to observations from Levitus and Boyer (1994) for SST and Levitus et al. (1994) for SSS and to ERA-40 reanalysis data for atmospheric temperature, evaporation, and precipitation [for atmospheric temperature, values at standard sea level pressure (1013.25 hPa) are used for the comparison]. Atmospheric temperatures in the tropics and the latitudinal gradient are well represented by the model. The largest deficiency is found in the GIN Seas, where temperatures are too low. Because deep water formation occurs too far south in the model, the warm ocean currents do not contribute to the atmospheric temperature in this region, leading to too cold temperatures. Because in zonal direction advection is not considered in the model and diffusion is constant, zonal gradients are smaller than in ERA-40 data. Sea surface temperatures are in good agreement with the observations. In the tropics, temperatures are too low. This deficit is associated with the precipitation pattern. As described before, sea surface salinity is generally too low compared to the data of Levitus et al. (1994). In the initialization phase of the energy balance model, the Dalton number CE is diagnosed while prescribing ERA-40 evaporation. Thus, the simulated evaporation pattern matches the data well. The total amount of evaporation (and precipitation) is 20% lower than in ERA-40. Because winds are zonally averaged, model precipitation has large differences to the data of ERA-40. When using zonally resolved winds, the global precipitation pattern improves, but unrealistically high precipitation occurs in various boxes (as a direct consequence of convergence), which negatively affects the state of the model.

Fig. 6.

Model results of annual mean atmospheric temperature, sea surface temperature, sea surface salinity, evaporation, and precipitation compared to ERA-40 reanalysis data [for atmospheric temperature, values at standard sea level pressure (1013.25 hPa) are used for the comparison]. In the third column the model – observation differences are displayed.

Fig. 6.

Model results of annual mean atmospheric temperature, sea surface temperature, sea surface salinity, evaporation, and precipitation compared to ERA-40 reanalysis data [for atmospheric temperature, values at standard sea level pressure (1013.25 hPa) are used for the comparison]. In the third column the model – observation differences are displayed.

Because NADW is formed too far south, modeled sea ice extent is too large in the Arctic Ocean east of Greenland as compared to the dataset of Rayner et al. (2003). On the other hand, Southern Ocean sea ice extent is slightly too low throughout the year (Fig. 7).

Fig. 7.

Seasonal variation of modeled sea ice cover compared to the dataset of Rayner et al. (2003). The dataset is averaged over years 1958–2001, the time span of the ERA-40 data. (a),(b) Time series of Northern and Southern Hemisphere sea ice cover. The ice cover in (b) was calculated from the original dataset (1° × 1° resolution, solid line) and after regridding to the Bern3D grid (dashed line). (c)–(f) February and September fractional ice area, respectively.

Fig. 7.

Seasonal variation of modeled sea ice cover compared to the dataset of Rayner et al. (2003). The dataset is averaged over years 1958–2001, the time span of the ERA-40 data. (a),(b) Time series of Northern and Southern Hemisphere sea ice cover. The ice cover in (b) was calculated from the original dataset (1° × 1° resolution, solid line) and after regridding to the Bern3D grid (dashed line). (c)–(f) February and September fractional ice area, respectively.

The seasonal variability of temperature in the atmosphere reproduces the larger July-to-January temperature difference over continents compared to the ocean (Fig. 8). The amplitude is too large over the GIN Seas because deep-water formation, which “attracts” the warm surface currents and which occurs too far south in the model, is a winter-only phenomenon. The amplitude is also too large over South America and the Southern Ocean. Evaporation and precipitation vary strongly during the season. The model reproduces the seasonal shift of precipitation in the equatorial zone (Fig. 9).

Fig. 8.

Seasonal variations of atmospheric temperature. Model results of (left) July and (middle) January temperatures, and (right) July–January differences, are compared to ERA-40 reanalysis data. In the bottom row the model – observation differences are displayed.

Fig. 8.

Seasonal variations of atmospheric temperature. Model results of (left) July and (middle) January temperatures, and (right) July–January differences, are compared to ERA-40 reanalysis data. In the bottom row the model – observation differences are displayed.

Fig. 9.

Seasonal variations of evaporation and precipitation. Model results of January and July (top two rows) evaporation and (bottom two rows) precipitation are compared to ERA-40 reanalysis data. ERA-40 evaporation over land surface is not shown in this figure. In the third column the model – observation differences are displayed.

Fig. 9.

Seasonal variations of evaporation and precipitation. Model results of January and July (top two rows) evaporation and (bottom two rows) precipitation are compared to ERA-40 reanalysis data. ERA-40 evaporation over land surface is not shown in this figure. In the third column the model – observation differences are displayed.

Correlation and relative standard deviation of modeled annual mean, January, and July fields of atmospheric temperature, evaporation, and precipitation with ERA-40 data are displayed in Fig. 5b.

Meridional heat transport in the model is shown in Fig. 10. The ocean heat fluxes are compared to observationally based estimates from Lumpkin and Speer (2007), Fasullo and Trenberth (2008), and Trenberth and Caron (2001). Modeled Atlantic meridional transport is 0.2 to 0.6 PW lower than found in the observations. In the Indo-Pacific, they agree well. Atmospheric heat flux is divided into sensible and latent heat flux. Compared to the data-based estimates of Fasullo and Trenberth (2008), atmospheric heat transport is by up to 1.5 PW too low in the Northern Hemisphere and by up to 2 PW too low in the Southern Hemisphere. The total heat flux, the sum of atmospheric and ocean heat flux, is also shown.

Fig. 10.

(a) Northward transport of heat in the Atlantic, Pacific, and global ocean. The model (thick lines) is compared to data from Lumpkin and Speer (2007) (squares) and to data of Trenberth and Caron (2001) derived from ECMWF atmospheric fields (thin lines). The total ocean heat flux is also compared to newer estimates inferred from satellite retrievals from the Earth Radiation Budget Experiment (ERBE) and Clouds and Earth’s Radiant Energy System (CERES) fields (Fasullo and Trenberth 2008) (dashed line). (b) Modeled northward heat transport of the atmosphere, the ocean, and in total (thick lines). The atmospheric transport is separated into sensible and latent heat flux. The model is compared to the satellite-based data of Fasullo and Trenberth (2008) (thin lines).

Fig. 10.

(a) Northward transport of heat in the Atlantic, Pacific, and global ocean. The model (thick lines) is compared to data from Lumpkin and Speer (2007) (squares) and to data of Trenberth and Caron (2001) derived from ECMWF atmospheric fields (thin lines). The total ocean heat flux is also compared to newer estimates inferred from satellite retrievals from the Earth Radiation Budget Experiment (ERBE) and Clouds and Earth’s Radiant Energy System (CERES) fields (Fasullo and Trenberth 2008) (dashed line). (b) Modeled northward heat transport of the atmosphere, the ocean, and in total (thick lines). The atmospheric transport is separated into sensible and latent heat flux. The model is compared to the satellite-based data of Fasullo and Trenberth (2008) (thin lines).

c. Sensitivities

In a first sensitivity test, modeled atmospheric temperatures are compared to ERA-40 data at 2 m above ground by taking into account land topography based on ETOPO5 (data available online at http://www.ngdc.noaa.gov/mgg/global/etopo5.html). Atmospheric temperature at altitude z is calculated by applying a constant lapse rate Γ = 5 K km−1. We replace Ta by Taalt = Ta − Γz in the atmosphere–land flux Eq. (4), the equation for sensible heat flux over land (17) and in Eq. (9) for the longwave radiation flux at TOA, thus updating Eq. (11) to . As in the standard case, the temperature at sea level is used for the horizontal transport. The value of Γ was chosen such that the root-mean-square deviation between modeled temperature and ERA-40 is minimized.

Global and annual mean atmospheric temperature at sea level is 16.3°C when topography is taken into account (run henceforth called ALTI), which compares to 15.4°C in the standard case (called CTRL); the mean ERA-40 temperature at sea level is 15.7°C. Root-mean-square deviations between atmospheric temperature of ALTI and ERA-40 at 2 m above ground are 2.8°C for the annual mean, 5.5°C in January, and 4.3°C in July. In comparison, root-mean-square deviations between atmospheric temperature of CTRL and ERA-40 at sea level are 3.0°C for the annual mean, 5.4°C in January, and 4.5°C in July. So generally, the model compares better to ERA-40 when topography is taken into account. Also, temperatures over the Eurasian continent are closer to ERA-40. On the other hand, temperatures over Antarctica and the Southern Ocean are too warm. As a consequence, AABW is too weak.

Modeled Northern Hemisphere summer temperatures deviate most from ERA-40 in the vicinity of the Hudson Bay (Fig. 8). Thus, in a second sensitivity test, the Hudson Bay is taken into account in the model by replacing two land boxes by ocean boxes in North America such that the Hudson Bay is connected to the Arctic Ocean. The presence of the Hudson Bay indeed leads to significantly cooler July air temperatures in this area by approximately 3°C compared to the standard case. However, it cannot explain the deviation from ERA-40 of approximately 13.5°C. Also, the Hudson Bay hardly affects the adjacent boxes and it does not affect the overall state of the model.

4. Application of the coupled Bern3D model: Paleosimulations

a. Last Glacial Maximum

To analyze the circulation state of the model under Last Glacial Maximum (LGM) conditions, a 10-kyr steady-state simulation is performed, setting orbital parameters to values at 20 kyr B.P. (before present; i.e., before 1950 A.D.), atmospheric CO2 to 180 ppm, atmospheric CH4 to 350 ppb, and increasing surface albedos to account for the presence of Northern Hemispheric ice sheets during the LGM. Since ice sheets are not simulated explicitly by the model, their extent is taken from Peltier (1994). During the first 2 kyr of the simulation, the ice sheet is linearly scaled from the modern to the LGM state, taking into account the relocation of the freshwater from the ocean onto the land and the additional latent heat flux (see appendix for details). As described in section 2, the freshwater relocation is done by adding salt to the ocean instead of lowering the sea level. As a consequence, ocean gateways such as the Bering Strait remain open.

The glacial state of the overturning circulation has been qualitatively reconstructed by analyzing estimates of LGM nutrient distributions inferred from Cd/Ca ratios of benthic foraminifera found in marine sediment cores (Lynch-Stieglitz et al. 2007; Marchitto and Broecker 2006). As in the present-day simulation, the model is run with the prognostic carbon cycle. Additionally, LGM aeolian iron fluxes have been applied (Mahowald et al. 2006). Modeled phosphate concentrations are converted to cadmium concentrations following the linear relationship described by Elderfield and Rickaby (2000). In the modeled glacial state, AMOC is shallower and somewhat weaker compared to the modern state, forming the so-called Glacial North Atlantic Intermediate Water (GNAIW) in agreement with paleoclimatic reconstructions (Labeyrie et al. 1992; Gherardi et al. 2009). On the other hand, Antarctic Bottom Water (AABW) becomes stronger and penetrates farther to the north (Fig. 11). This is also seen in the reconstruction of the cadmium distribution. Compared to the paleoceanographic estimates, model concentrations are by approximately 0.1 to 0.2 nmol kg−1 too high. Also, too strong nutrient trapping in the equatorial surface ocean creates too high cadmium values in this region. However, as discussed for the modern phosphate distribution, these high concentrations occur at the African Margin, where no paleoceanographic data are available (Marchitto and Broecker 2006). Modeled equatorial concentrations in the central surface and deep Atlantic are by approximately 0.1 nmol kg−1 lower. Note that paleoceanographic data used for the reconstruction are very sparse.

Fig. 11.

(a) Last Glacial Maximum (LGM) Atlantic overturning circulation. (b) Zonally averaged modeled Atlantic cadmium concentrations for the LGM. Modeled cadmium concentrations are obtained by applying the following equation to the PO4 model output: [Cd] = 1.2 nmol kg−1/(2(3300 nmol kg−1/[PO4] − 1) + 1) (Elderfield and Rickaby 2000). (c) LGM cadmium concentration estimates from the ratio of Cd/Ca in shells of benthic foraminifera (Marchitto and Broecker 2006). The data locations are indicated in the figure. (d) Modeled cadmium field, where model data are only taken at the data locations of Marchitto and Broecker (2006). The interpolation is done as in (c).

Fig. 11.

(a) Last Glacial Maximum (LGM) Atlantic overturning circulation. (b) Zonally averaged modeled Atlantic cadmium concentrations for the LGM. Modeled cadmium concentrations are obtained by applying the following equation to the PO4 model output: [Cd] = 1.2 nmol kg−1/(2(3300 nmol kg−1/[PO4] − 1) + 1) (Elderfield and Rickaby 2000). (c) LGM cadmium concentration estimates from the ratio of Cd/Ca in shells of benthic foraminifera (Marchitto and Broecker 2006). The data locations are indicated in the figure. (d) Modeled cadmium field, where model data are only taken at the data locations of Marchitto and Broecker (2006). The interpolation is done as in (c).

LGM sea ice cover is larger in the North Pacific and the Southern Ocean compared to the modern model state, whereas the North Atlantic shows little difference (Fig. 12). The Southern Ocean sea ice margin extends farthest in the Atlantic sector to about 45°S in winter and 60°S in summer. This is in rough agreement with reconstructions from Gersonde et al. (2005), who also find the largest sea ice extent in the Atlantic sector to be about 48°S in winter and 52°S in summer.

Fig. 12.

Seasonal variation of LGM sea ice cover. (a) Time series of Northern and Southern Hemisphere sea ice cover; (b) February and (c) September fractional ice area.

Fig. 12.

Seasonal variation of LGM sea ice cover. (a) Time series of Northern and Southern Hemisphere sea ice cover; (b) February and (c) September fractional ice area.

In a sensitivity test the Bering Strait was closed. Compared to the standard setup, this results in a (2 Sv) weaker and shallower AMOC and stronger AABW, and the Arctic Ocean becomes fresher. When additionally closing the passage between North America and Greenland, the Arctic freshens more. Global mean atmospheric temperature is not significantly affected.

Comprehensive climate models such as NCAR CCSM and HadCM3 also show a shallower and weaker AMOC in the LGM state as compared to the modern state. Others, such as the Model for Interdisciplinary Research on Climate (MIROC; Hasumi and Emori 2004) and ECBilt-CLIO, show a stronger and deeper AMOC (Otto-Bliesner et al. 2007). Weber et al. (2007) developed a metric to analyze the reasons for the different behaviors of the models and applied it to various models of the Paleoclimate Modeling Intercomparison Project (PMIP2). Here, we apply the same metric to the Bern3D model. First, the components of the Atlantic freshwater budget according to the balance equation Msurf = Mov + Maz + Mdiff + MBS (in equilibrium) for the modern state and their difference to the LGM state are calculated. Here Msurf is the basin integral of evaporation, precipitation, continental runoff, ice melt, brine rejection, and flux corrections; Mov is the meridional overturning component and Maz the azonal component of the oceanic freshwater transport through the southern boundary. The Bering Strait and diffusive contributions are determined as a residual term R = Mdiff + MBS. Other quantities to determine Southern Ocean controls versus Atlantic processes are ρAtl, the density difference between the northern and southern ends of the Atlantic basin (taken at 55°N and 30°S, respectively and averaged over the top 1500 m), and ρNS, the density contrast between NADW and AABW (taken at 55°N and 55°S, respectively and averaged over all depth levels). The results of the Bern3D model are summarized and compared to other models in Tables 2 and 3.

Table 2.

Values for the AMOC strength ΨAMOC, the Atlantic basin-integrated result of evaporation, precipitation, and other freshwater fluxes Msurf, the meridional overturning component Mov and the zonal component Maz of the oceanic freshwater transport through the South Atlantic boundary, and a rest term R that includes Bering Strait and diffusive contributions to close the Atlantic freshwater budget (in Sv). The values shown are for the modern model state (mod) and for the LGM to modern difference (Δ = LGM − mod). The Bern3D values are compared to observations (Pardaens et al. 2003; Weijer et al. 1999) and to various PMIP2 models (Weber et al. 2007).

Values for the AMOC strength ΨAMOC, the Atlantic basin-integrated result of evaporation, precipitation, and other freshwater fluxes Msurf, the meridional overturning component Mov and the zonal component Maz of the oceanic freshwater transport through the South Atlantic boundary, and a rest term R that includes Bering Strait and diffusive contributions to close the Atlantic freshwater budget (in Sv). The values shown are for the modern model state (mod) and for the LGM to modern difference (Δ = LGM − mod). The Bern3D values are compared to observations (Pardaens et al. 2003; Weijer et al. 1999) and to various PMIP2 models (Weber et al. 2007).
Values for the AMOC strength ΨAMOC, the Atlantic basin-integrated result of evaporation, precipitation, and other freshwater fluxes Msurf, the meridional overturning component Mov and the zonal component Maz of the oceanic freshwater transport through the South Atlantic boundary, and a rest term R that includes Bering Strait and diffusive contributions to close the Atlantic freshwater budget (in Sv). The values shown are for the modern model state (mod) and for the LGM to modern difference (Δ = LGM − mod). The Bern3D values are compared to observations (Pardaens et al. 2003; Weijer et al. 1999) and to various PMIP2 models (Weber et al. 2007).
Table 3.

Modern to LGM differences (LGM − modern) of ΨAMOC, Msurf, the density difference between the northern and southern ends of the Atlantic ρatl, and the density difference between NADW and AABW ρNS in the Bern3D model and in various PMIP2 models (Weber et al. 2007). The potential controlling process for the change in AMOC strength in each model is indicated in the last column. None of these processes could be attributed as a controlling factor in the Bern3D model as well as in ECBilt and UVic, since they are either anticorrelated to ΨAMOC or insignificant.

Modern to LGM differences (LGM − modern) of ΨAMOC, Msurf, the density difference between the northern and southern ends of the Atlantic ρatl, and the density difference between NADW and AABW ρNS in the Bern3D model and in various PMIP2 models (Weber et al. 2007). The potential controlling process for the change in AMOC strength in each model is indicated in the last column. None of these processes could be attributed as a controlling factor in the Bern3D model as well as in ECBilt and UVic, since they are either anticorrelated to ΨAMOC or insignificant.
Modern to LGM differences (LGM − modern) of ΨAMOC, Msurf, the density difference between the northern and southern ends of the Atlantic ρatl, and the density difference between NADW and AABW ρNS in the Bern3D model and in various PMIP2 models (Weber et al. 2007). The potential controlling process for the change in AMOC strength in each model is indicated in the last column. None of these processes could be attributed as a controlling factor in the Bern3D model as well as in ECBilt and UVic, since they are either anticorrelated to ΨAMOC or insignificant.

As in most PMIP2 models, compared to observations (Pardaens et al. 2003; Weijer et al. 1999) Msurf is overestimated and Maz underestimated in the modern state of the Bern3D model. The overturning component Mov is negative but too small compared to the observations. The sign of Mov simulated by NCAR CCSM, HadCM3, and the UVic model is inconsistent with the observations. The diffusive transport is overestimated, as can be expected from a coarse-resolution model. The differences between modern and LGM values for these quantities are small in the Bern3D model compared to the PMIP2 models, except for UVic. Although Msurf remains the same, Mov switches to a positive value. This change is compensated by the azonal component. The differences in Mov and Maz between glacial and interglacial state are consistent with CCSM, HadCM3, MIROC, and ECBilt. Weber et al. (2007) conclude that none of these quantities are good indicators for the changes in AMOC strength, basically because their LGM-to-modern differences are all in line and independent of the differences in AMOC strength. This is also true for the Bern3D model, especially because Msurf exhibits negligible changes and ΔMov is positive but ΔΨAMOC negative.

In a second step, Weber et al. (2007) analyze the correlation of ΔΨAMOC with Δρatl and ΔρNS. The processes that correlate positively with ΔΨAMOC indicate potential drivers for the AMOC strength. In the Bern3D model, ΔMsurf ≈ 0, ΔρNS ≈ 0, and ΔΨMOC and Δρatl are anticorrelated. Therefore, as for the other EMICs (ECBILT-CLIO and UVic) analyzed by Weber et al. (2007), the changes in the analyzed quantities are too small to attribute them as controlling processes for the change in AMOC strength (Table 3).

b. Transient simulations over the past 800 000 years

1) Variations in greenhouse gas concentrations and ice sheet extent

Simulations over ice age cycles are possible with the Bern3D model. Here we present seven transient simulations over the past 800 000 years (Table 4), where orbital forcing, atmospheric greenhouse gases CO2 and CH4 and ice sheets are prescribed. In simulation ORBI, only orbital forcing is applied. In addition to the orbital forcing, simulation PCO2 varies CO2, simulation PCH4 varies CH4, and simulation ICE1 varies ice sheets. In ALL1, all forcings are taken into account. For CO2 and CH4, data of Lüthi et al. (2008) and Loulergue et al. (2008), respectively, are used (Fig. 13c). As before, ice sheets for the modern and LGM states are taken from Peltier (1994) and scaled using the benthic δ18O stack of Lisiecki and Raymo (2005), which is a proxy of global ice volume (Fig. 13d; see appendix for details). The simulations that involve ice sheet changes are computed twice: once taking into account albedo changes, the freshwater relocation from the ocean onto the land and vice versa, and the additional latent heat flux when ice sheets are formed and wasted (simulations ICE1 and ALL1) and once only accounting for albedo changes (ICE2 and ALL2). Snow–albedo feedback is not included in the simulations.

Table 4.

Summary of the 800-kyr model simulations. Quantities not mentioned in the “prescribed forcing” column are kept constant at preindustrial values. These are land mask, winds, atmospheric and ocean diffusivities, terrestrial surface albedo (unless ice sheets are present), cloud cover, atmospheric CO2 and CH4, and ice sheets.

Summary of the 800-kyr model simulations. Quantities not mentioned in the “prescribed forcing” column are kept constant at preindustrial values. These are land mask, winds, atmospheric and ocean diffusivities, terrestrial surface albedo (unless ice sheets are present), cloud cover, atmospheric CO2 and CH4, and ice sheets.
Summary of the 800-kyr model simulations. Quantities not mentioned in the “prescribed forcing” column are kept constant at preindustrial values. These are land mask, winds, atmospheric and ocean diffusivities, terrestrial surface albedo (unless ice sheets are present), cloud cover, atmospheric CO2 and CH4, and ice sheets.
Fig. 13.

A series of 800-kyr simulations. (a) Global and annual mean atmospheric temperatures for a simulation with varying solar insolation (green), insolation and CH4 (orange), insolation and CO2 (red), insolation and ice sheets (dark blue when freshwater relocation from the ocean onto the land and the additional latent heat flux is taken into account, light blue when it is not and only the surface albedo changes), and the combination of the three (black and gray, respectively, when the freshwater relocation is not taken into account). For the parameters held constant, preindustrial values are used. (b) Changes of AMOC strength for the runs described above. For ORBI, PCH4, PCO2, and ICE2, the 500-yr running average is plotted to disambiguate the panel. (c)–(e) Atmospheric CO2 (Lüthi et al. 2008), the benthic δ18O stack of Lisiecki and Raymo (2005) (which is a proxy for global ice volume), and the δD Antarctic deuterium record (a proxy for atmospheric temperature in Antarctica) (Jouzel et al. 2007), respectively. The shaded vertical bars indicate interglacial periods (Augustin et al. 2004).

Fig. 13.

A series of 800-kyr simulations. (a) Global and annual mean atmospheric temperatures for a simulation with varying solar insolation (green), insolation and CH4 (orange), insolation and CO2 (red), insolation and ice sheets (dark blue when freshwater relocation from the ocean onto the land and the additional latent heat flux is taken into account, light blue when it is not and only the surface albedo changes), and the combination of the three (black and gray, respectively, when the freshwater relocation is not taken into account). For the parameters held constant, preindustrial values are used. (b) Changes of AMOC strength for the runs described above. For ORBI, PCH4, PCO2, and ICE2, the 500-yr running average is plotted to disambiguate the panel. (c)–(e) Atmospheric CO2 (Lüthi et al. 2008), the benthic δ18O stack of Lisiecki and Raymo (2005) (which is a proxy for global ice volume), and the δD Antarctic deuterium record (a proxy for atmospheric temperature in Antarctica) (Jouzel et al. 2007), respectively. The shaded vertical bars indicate interglacial periods (Augustin et al. 2004).

The changes in atmospheric global and annual mean temperature and AMOC strength are shown in Fig. 13. In the model, the effects of orbital forcing only or in combination with CH4 forcing are small. The CO2 and orbital forcing change global mean atmospheric temperatures by 1.5° to 2°C, where the higher value stands for large glacial-to-interglacial changes such as the LGM-to-modern change. The lower value represents smaller glacial-to-interglacial changes (e.g., at 500 kyr B.P.). The AMOC is not significantly affected by this forcing. Ice sheet forcing provokes glacial-to-interglacial temperature changes of 1° to 1.5°C. In ICE1, when freshwater relocation from the ocean to the continent is taken into account, the AMOC increases rapidly at glacial inceptions and decreases again at the terminations. The reason for this behavior is the densification of the surface ocean due to the ice sheet buildup. In the North Atlantic, this leads to stronger convection and hence to stronger AMOC. Analogously, the AMOC weakens during glacial terminations. This result is in agreement with an earlier study by Meissner and Gerdes (2002), who also find an AMOC intensification during glaciation.

In ALL1, when all forcings are combined, glacial-to-interglacial temperatures vary by 3° to 4°C. This temperature difference is at the lower end of the 4° to 7°C cooling estimated for the LGM (Jansen et al. 2007). Modeled glacial-to-interglacial temperatures vary by 2° to 3°C in the tropics, by 6° to 8°C in the southern high latitudes, and by 3.5° to 6°C in the northern high latitudes. Pollen records also indicate to a 2° to 3°C cooler climate in the tropics during the LGM (Farrera et al. 1999). For the high latitudes, temperature reconstructions from polar ice cores indicate a LGM-to-modern temperature increase of about 9°C in Antarctica (Stenni et al. 2001) and about 22°C in Greenland (Dahl-Jensen et al. 1998). We tentatively attribute the large difference between reconstructed and modeled northern high-latitude temperatures to the too low sea ice sensitivity and to the absence of topographic effects of the North American and Eurasian ice sheets, which changes the course of the jet stream. Besides the different magnitudes of the change, the evolution of temperature at the poles and in the tropics and of the global mean are similar, suggesting a fast glacial-to-interglacial transition and a more gradual interglacial-to-glacial transition, as also seen in the atmospheric CO2 (Fig. 13c) and the δ18O (Fig. 13d) records. The relatively fast interglacial-to-glacial transitions found in the δD Antarctic ice core record (Jouzel et al. 2007) (Fig. 13e), a proxy for atmospheric temperature in Antarctica, are not represented by the model.

As in ICE1, AMOC in ALL1 increases at the glacial inception, but its strength decreases again with the cooling. At the glacial maxima, the AMOC is weaker than during the interglacial periods. The weakening of the AMOC is also observed when the ice sheet freshwater relocation is not taken into account. Global atmospheric temperatures are similar in ALL1 and ALL2. Note that the parameterization of the ice sheets does not allow ablation at the ice sheet margin in an ice sheet buildup phase. Also, the switch from an ice sheet melting to an ice sheet build-up phase and vice versa occurs globally at the same time. This might lead to an overestimation of the abrupt AMOC changes that are induced by the switch from ice sheet melting to ice sheet buildup.

2) Local ocean temperature sensitivity to changes in global mean atmospheric temperature, AMOC, and AABW

We focus on ALL1, the most realistic simulation of this study, and determine the sensitivity of ocean temperature in various regions to changes in global mean atmospheric temperature, AMOC, and AABW strength. Therefore, the assumption is made that the ocean temperature time series at any location in the ocean can be reconstructed by a linear combination of global mean atmospheric temperature Tatm, AMOC, and AABW strength (ΨAMOC, Fig. 14c, and ΨAABW, Fig. 14d) time series. The purpose is to test whether the inversion of this approach would be feasible, that is, to determine important global-scale atmosphere (Tatm) and ocean (Ψ) variables using a specific combination of reconstructions of local temperature from paleoceanographic records.

Fig. 14.

(a),(b) Atlantic and Pacific equatorial deep ocean temperature of simulation ALL1 compared to reconstructions of this region (Martin et al. 2002). (c) AMOC strength. (d) AABW strength. (e) Normalized atmospheric temperature (gray) and deep equatorial Atlantic temperature (black). The normalization was done as follows: = (TT)/ max(TT), where T is the temperature average. The red curve is a fit of the function to the normalized local temperature; ΨAMOC and ΨAABW are the AMOC and AABW strengths, respectively, as in (c) and (d). The coefficients as well as the correlations between atmospheric temperature and local temperature RTatm, and between the fit and the local temperature Rfit, are listed. (f) As in (e), but for the deep equatorial Pacific.

Fig. 14.

(a),(b) Atlantic and Pacific equatorial deep ocean temperature of simulation ALL1 compared to reconstructions of this region (Martin et al. 2002). (c) AMOC strength. (d) AABW strength. (e) Normalized atmospheric temperature (gray) and deep equatorial Atlantic temperature (black). The normalization was done as follows: = (TT)/ max(TT), where T is the temperature average. The red curve is a fit of the function to the normalized local temperature; ΨAMOC and ΨAABW are the AMOC and AABW strengths, respectively, as in (c) and (d). The coefficients as well as the correlations between atmospheric temperature and local temperature RTatm, and between the fit and the local temperature Rfit, are listed. (f) As in (e), but for the deep equatorial Pacific.

First the quantities are normalized as follows: = (TT)/ max(TT), where T is the average over 800 kyr, and likewise for Ψ. Then the linear combination

 
formula

is fitted to the local ocean temperature time series oc by minimizing . The calculation returns values for the coefficients aTatm, aAMOC, and aAABW, which give information on the sensitivity of the local temperature to these quantities. Note that Tatm, ΨAMOC, and ΨAABW are not orthogonal, but their correlations are small: 0.46 between Tatm and ΨAMOC, −0.44 between Tatm and ΨAABW, and 0.07 between ΨAMOC and ΨAABW.

Modeled deep ocean temperatures can be validated by comparing the results to reconstructions from marine sediment cores. However, because of major difficulties in reconstructing deep-sea temperatures, only few time series exist currently. Martin et al. (2002) provide reconstructions from benthic foraminiferal Mg/Ca records for the last 330 kyr in the eastern tropical Atlantic and for the past 230 kyr in the eastern tropical Pacific (Figs. 14a,b). Because of the large uncertainties of the reconstructions of about 1.5°C, we cannot go into detail when comparing the data with the model results. However, as already discussed for the modern ocean, the modeled deep Pacific is too cold by about 1°C. One major feature of the reconstruction, the early warming of the Pacific during glacial times, is not reflected in the model. Model and reconstruction compare somewhat better in the Atlantic, where both show a gradual temperature decrease during glaciation.

The fit for the deep equatorial Atlantic temperature time series correlates surprisingly well with the modeled temperature (Rfit = 0.99; Fig. 14e). For a comparison, the correlation between atmospheric and ocean temperature is substantially lower, with RTatm = 0.90. Thus, the circulation plays an important role in this region. The calculated coefficients are aTatm = 0.75, aNADW = 0.53, and aAABW = −0.13. In other words, a 100% change in atmospheric temperature explains a 75% change in the deep equatorial Atlantic temperature. A 100% increase in AABW strength explains a local temperature decrease by 13%. At this location, all three components contribute to changes in the local temperature, although the role of AABW is minor.

The same exercise for deep equatorial Pacific temperature results in a correlation of Rfit = 0.97, and the coefficients read aTatm = 0.70, aNADW = 0.15, and aAABW = −0.13 (Fig. 14f). As expected, AABW is more important here than in the Atlantic, and correspondingly the influence of AMOC is less important. However, it is surprising that the AMOC is not negligible in this region: when AMOC is not taken into account in the fitting, the correlation Rfit is reduced to 0.96.

The described procedure has been applied to 23 locations in the ocean. The chosen regions are the northern, equatorial, and southern Atlantic and Pacific and the equatorial and southern Indian Ocean. In every region an average of several boxes is taken (light gray area in Fig. 15) and separated into upper ocean (400–650-m depth), intermediate ocean (1750–2300-m depth), and deep ocean (>2300 m, depending on the topography).

Fig. 15.

Sensitivity of ocean temperature at various regions to changes in atmospheric temperature, AMOC, and AABW strength in the ALL1 simulation. Every region is an average of several boxes (light gray area) and is separated into upper ocean (400–650-m depth), intermediate ocean (1750–2300-m depth), and deep ocean (>2300-m depth depending on the topography). For reasons of topography, the deep North Atlantic location does not represent bottom waters in this region. The bars depict the values of the coefficients in Eq. (29), given in %. They therefore give the percentage of a local temperature change per 100% change in atmospheric temperature, AMOC strength, and AABW strength, respectively. In brackets the local temperature change is given in °C per °C change in atmospheric temperature (and °C per Sv change in AMOC and AABW, respectively). The R value is the correlation between modeled local temperature time series and the linear combination of atmospheric temperature, AMOC, and AABW time series according to the given coefficients.

Fig. 15.

Sensitivity of ocean temperature at various regions to changes in atmospheric temperature, AMOC, and AABW strength in the ALL1 simulation. Every region is an average of several boxes (light gray area) and is separated into upper ocean (400–650-m depth), intermediate ocean (1750–2300-m depth), and deep ocean (>2300-m depth depending on the topography). For reasons of topography, the deep North Atlantic location does not represent bottom waters in this region. The bars depict the values of the coefficients in Eq. (29), given in %. They therefore give the percentage of a local temperature change per 100% change in atmospheric temperature, AMOC strength, and AABW strength, respectively. In brackets the local temperature change is given in °C per °C change in atmospheric temperature (and °C per Sv change in AMOC and AABW, respectively). The R value is the correlation between modeled local temperature time series and the linear combination of atmospheric temperature, AMOC, and AABW time series according to the given coefficients.

In all locations, the sensitivity is highest to the atmospheric temperature. AMOC has a strong influence in the Southern Hemisphere intermediate and deep oceans and in the deep Atlantic and Indian Oceans. In most cases, an increase in AMOC leads to a warming. Exceptions are only the upper equatorial and southern Atlantic, where a stronger AMOC exports more heat to the North and thus leads to a local cooling. In the Pacific and Indian Ocean, stronger AABW leads to cooling. In the intermediate Atlantic, where the AABW cell returns south, a stronger AABW induces warming.

Using the local sensitivity coefficients calculated above, the model can potentially be used to extract information of the ocean circulation from paleoceanographic temperature reconstructions. This requires ocean temperature time series at multiple locations.

5. Application of the coupled Bern3D model: Climate sensitivity

Here we investigate the possibility of model state-dependent equilibrium climate sensitivity by comparing the temperature response of the model to a doubling and quadrupling of atmospheric CO2 concentration for the modern and the LGM states. Atmospheric pCO2 is doubled from 278 to 556 ppm and quadrupled to 1112 ppm in the modern case and from 180 to 360 ppm and to 720 ppm, respectively, in the LGM case. The model is run for 10 000 years to a new equilibrium. We compare the change in global mean atmospheric temperature and the associated relaxation time scales.

As described earlier, the modern state equilibrium climate sensitivity to a doubling of atmospheric CO2 is tuned to 3°C. Modeled climate sensitivity of the LGM state is 4.3°C. This higher sensitivity is probably due to the presence of more sea ice in the LGM. When sea ice melts, surface albedo decreases strongly and more radiation is absorbed. Note that we assume that the water vapor feedback parameter λ remains constant. Also, land ice remains constant. A quadrupling of atmospheric CO2 from the modern state results in a temperature increase of 5.5°C and from the LGM state in an increase of 7.3°C.

To quantify the relaxation time scales, a least squares exponential fit of the form with the constraint is performed, where T0 is the equilibrium climate sensitivity and τi are characteristic time scales of the transient global mean response. The results for the four simulations are summarized in Table 5 and Fig. 16. In the 2 × CO2 modern case, τ1 = 3.0 yr, τ2 = 40 yr, and τ3 = 600 yr with the coefficients a1 = 0.48, a2 = 0.30, and a3 = 0.22. How long it takes for the climate to equilibrate depends on the efficiency of ocean heat uptake. The empirical time scales cannot be assigned to individual processes, but in general τ1 covers short-term processes such as wind-driven mixing at the surface ocean and convection into the intermediate and deep ocean; τ2 and τ3 cover intermediate and long-term processes that transport the temperature signal to the deep ocean. To a large extent, the equilibration time scale of the deep ocean is determined by the AMOC and AABW strength. When comparing the modern 2 × CO2 and 4 × CO2 simulations, it is notable that a3 and τ3 of the 4 × CO2 run are substantially reduced. While the ocean circulation is hardly affected by the CO2 doubling, the AMOC collapses in the 4 × CO2 model world and AABW strength increases from 18 to 28 Sv. The relaxation time scales of the LGM 2 × CO2 and 4 × CO2 simulations are similar. However, the coefficient a3 is reduced in the 4 × CO2 run, indicating a faster overall response in the 4 × CO2 run compared to the 2 × CO2 run. Compared to the modern 2 × CO2 simulation the time scales are longer, though. In particular, τ3 is twice as large. The steady-state LGM ocean circulation after the perturbation shows a strong AMOC of 17 Sv (2 × CO2) and 18 Sv (4 × CO2), respectively, filling the entire Atlantic. AABW strength decreases to 14 Sv in both runs.

Table 5.

Equilibrium climate sensitivity and the relaxation time scales for a doubling and quadrupling of atmospheric CO2 in the modern and LGM climate state. The 4 × CO2 case is compared to results from the HadCM3 model (Li and Jarvis 2009) (2σ confidence interval in brackets).

Equilibrium climate sensitivity and the relaxation time scales for a doubling and quadrupling of atmospheric CO2 in the modern and LGM climate state. The 4 × CO2 case is compared to results from the HadCM3 model (Li and Jarvis 2009) (2σ confidence interval in brackets).
Equilibrium climate sensitivity and the relaxation time scales for a doubling and quadrupling of atmospheric CO2 in the modern and LGM climate state. The 4 × CO2 case is compared to results from the HadCM3 model (Li and Jarvis 2009) (2σ confidence interval in brackets).
Fig. 16.

The 2 × CO2 and 4 × CO2 climate sensitivity simulations for the modern and LGM case as summarized in Table 5. The temperature responses are normalized to compare the rates of change. Also shown is a 4 × CO2 simulation from the HadCM3 model for the modern case (Li and Jarvis 2009). Because the model was only run for 1000 years, the uncertainties are large.

Fig. 16.

The 2 × CO2 and 4 × CO2 climate sensitivity simulations for the modern and LGM case as summarized in Table 5. The temperature responses are normalized to compare the rates of change. Also shown is a 4 × CO2 simulation from the HadCM3 model for the modern case (Li and Jarvis 2009). Because the model was only run for 1000 years, the uncertainties are large.

It appears that the response time of the climate system to a strong CO2 perturbation strongly depends on the reaction of its components. Since in the Bern3D model winds are not affected by the warming, τ1 and a1 are similar in all simulations; however, τ2 and τ3 anticorrelate with AABW strength, which determines how fast the temperature signal is transported to the remotest regions of the ocean. The role of the AMOC is minor.

Li and Jarvis (2009) have also calculated the relaxation time scales to a 4 × CO2 perturbation of the modern climate using the HadCM3 AOGCM (Table 5). While HadCM3 shows a similar response time on the short time scale (τ1 = 4.5 yr), τ2 = 140 yr and τ3 = 1476 yr are approximately 4 times larger than in the Bern3D model. We attribute this slower response in atmospheric temperature of HadCM3 to a slower surface-to-deep transport in the ocean than in the Bern3D model. Note that the uncertainties are quite large in HadCM3 because the model was run for only 1000 years.

6. Conclusions

We have developed an efficient coupled three-dimensional dynamical ocean–energy balance atmosphere model capable of performing sensitivity studies and ensemble simulations on glacial–interglacial time scales. In this work, an energy balance component has been added to the previously developed physical ocean circulation model (Edwards et al. 1998; Müller et al. 2006) and the modules describing the penetration of transient tracers (CFCs, anthropogenic carbon, bomb-produced and natural radiocarbon; Müller et al. 2006, 2008); the cycling of carbon, carbon isotopes, alkalinity, phosphate, iron, silica, calcite, and oxygen (Parekh et al. 2008; Tschumi et al. 2008); 231Pa/230Th (Siddall et al. 2007); the ecosystem model Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) and the cycle of aragonite (Aumont and Bopp 2006); a sediment module (Tschumi 2009); and an ensemble Kalman filter framework for inverse analyses (Gerber et al. 2009; Gerber and Joos 2010). Although the resolution of the model is coarse, zonal gradients within ocean basins are resolved. Thorough comparisons of the modern model state to observations and reanalysis data have been made to validate the model.

Several 800 000-yr simulations were performed where the model was forced with orbital parameters, greenhouse gases, and ice sheets. We find that changes in these parameters have direct effects on the global overturning circulation. Atlantic meridional overturning varies between approximately 8 and 19 Sv in these transient simulations, with the lowest values during glacial maxima and terminations and highest values during glacial inceptions. When ice sheets, CO2, and CH4 are prescribed, the model shows glacial–interglacial global temperature variations of 3° to 4°C. This is on the lower end of the 4° to 7°C cooling estimated for the LGM (Jansen et al. 2007). A low bias is expected as the cooling influence of forcing by dust and by vegetation changes is not considered here. The CO2-only forcing accounts for 1.5° to 2°C and the ice sheet-only forcing for 1° to 1.5°C. The model shows an increase of AMOC strength at the glacial inception due to the relocation of freshwater from the surface ocean onto the land for the ice sheet buildup. In the glacial state, NADW weakens and shallows as indicated by paleoceanographic reconstructions.

We analyzed the sensitivity of the ocean temperature at 23 locations to changes in atmospheric temperature, AMOC, and AABW strength. It is surprising how well the local temperature correlates with a linear combination of these three components. Although atmospheric temperature has the greatest influence at all analyzed locations, depending on the region, AMOC and AABW also have considerable influence.

Modeled LGM climate is, with a global mean temperature change of 4.3°C, more sensitive to a doubling of CO2 than the modern climate (3°C climate sensitivity). How long it takes for the ocean to equilibrate strongly depends on the reaction of AABW to the CO2 change. In the LGM case, the CO2 doubling weakens AABW strength. Thus, the relaxation time scale is approximately twice as large as in the modern case. On the other hand, the AMOC collapses in a modern 4 × CO2 world, AABW strength increases, and equilibration is faster.

With the recently completed coupling of a sediment model (Tschumi 2009) and the incorporation of the PISCES ecosystem model, the Bern3D model is approaching a comprehensive EMIC. The next step will be to couple the Lund–Potsdam–Jena dynamic global vegetation model (Sitch et al. 2003; Strassmann et al. 2008) to the Bern3D model.

Future work could include fine tuning of the model using an ensemble Kalman filter (Gerber et al. 2009). Applications of the coupled model could include comparing past long-time scale local climate reconstructions such as SST from marine sediment cores (Martrat et al. 2007) with the model. This could further validate the model and provide a quantitative understanding of the proxy behavior and of the past global climate by bringing the variations found in the reconstruction into a global context. Also, the step can be made to estimate past AMOC and AABW changes using the model and bottom water temperature time series. However, the result will have large uncertainties due not only to the model but to the difficulties in reconstructing bottom water temperatures. We are now also in a position to tackle transient glacial–interglacial changes in the atmospheric and oceanic carbon isotopes Δ14C and δ13C. Again, the model results will be compared to observations (Reimer et al. 2004; Marchitto et al. 2007; Elsig et al. 2009). One of the main goals of paleoclimate research is to quantitatively explain glacial–interglacial variations of atmospheric CO2. Here, the addition of the EBM to the ocean component will allow us to progress beyond earlier work with the Bern3D model (Parekh et al. 2008; Tschumi et al. 2008). The Bern3D model will serve as a tool with which ensemble simulations of the Earth System over many 100 000-yr periods are possible.

Acknowledgments

This study was funded by GRACCIE (CONSOLIDER-INGENIO 2010) and by the Swiss National Science Foundation. F. Joos acknowledges support by the European Community’s Seventh Framework Programme (FP7/2007–2013) through the European Project on Ocean Acidification (EPOCA) and the Project Past4Future and by the Swiss Staatssekretariat für Bildung und Forschung (Cost Action 735). We thank J. Lynch-Stieglitz for providing the LGM cadmium data, M. Gerber for sharing the regridded CFC data, and N. Edwards for valuable comments concerning the EBM. We also thank N. Edwards, S. A. Müller, T. Tschumi, P. Parekh, R.Gangstø, M. Gerber, and A. Bass for their contributions in the development of the Bern3D model and three anonymous reviewers, whose comments have considerably improved the paper. ERA-40 reanalysis fields used in the EBM are provided by the European Centre for Medium-Range Weather Forecasts (ECMWF).

REFERENCES

REFERENCES
Augustin
,
L.
, and
Coauthors
,
2004
:
Eight glacial cycles from an Antarctic ice core.
Nature
,
429
,
623
628
.
Aumont
,
O.
, and
L.
Bopp
,
2006
:
Globalizing results from ocean in situ iron fertilization studies.
Global Biogeochem. Cycles
,
20
,
GB2017
.
doi:10.1029/2005GB002591
.
Berger
,
A. L.
,
1978
:
Long-term variations of daily insolation and quaternary climatic changes.
J. Atmos. Sci.
,
35
,
2362
2367
.
Bintanja
,
R.
,
1996
:
The parameterization of shortwave and longwave radiative fluxes for use in zonally averaged climate models.
J. Climate
,
9
,
439
454
.
Bolton
,
D.
,
1980
:
The computation of equivalent potential temperature.
Mon. Wea. Rev.
,
108
,
1046
1053
.
Collins
,
W. D.
, and
Coauthors
,
2006
:
The Community Climate System Model version 3 (CCSM3).
J. Climate
,
19
,
2122
2143
.
Conkright
,
M. E.
,
R. A.
Locarnini
,
H. E.
Garcia
,
T. D.
O’Brien
,
T. P.
Boyer
,
C.
Stephens
, and
J. I.
Antonov
,
2002
:
World Ocean Atlas 2001: Objective analyses, data statistics, and figures—CD-ROM Documentation.
National Oceanographic Data Center, 17 pp. [Available online at ftp://ftp.nodc.noaa.gov/pub/data.nodc/woa/PUBLICATIONS/readme.pdf]
.
Dahl-Jensen
,
D.
,
K.
Mosegaard
,
N.
Gundestrup
,
G. D.
Clow
,
S. J.
Johnsen
,
A. W.
Hansen
, and
N.
Balling
,
1998
:
Past temperatures directly from the Greenland ice sheet.
Science
,
282
,
268
271
.
Doronin
,
Y. P.
, and
D.
Kheisin
,
1977
:
Sea Ice.
Amerind, 323 pp
.
Edwards
,
N. R.
, and
R.
Marsh
,
2005
:
Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model.
Climate Dyn.
,
24
,
415
433
.
Edwards
,
N. R.
,
A. J.
Willmott
, and
P. D.
Killworth
,
1998
:
On the role of topography and wind stress on the stability of the thermohaline circulation.
J. Phys. Oceanogr.
,
28
,
756
778
.
Elderfield
,
H.
, and
R. E. M.
Rickaby
,
2000
:
Oceanic Cd/P ratio and nutrient utilization in the glacial Southern Ocean.
Nature
,
405
,
305
310
.
Elsig
,
J.
, and
Coauthors
,
2009
:
Stable isotope constraints on Holocene carbon cycle changes from an Antarctic ice core.
Nature
,
461
,
507
510
.
Enting
,
I. G.
,
1987
:
On the use of smoothing splines to filter CO2 data.
J. Geophys. Res.
,
92
,
10977
10984
.
Farrera
,
I.
, and
Coauthors
,
1999
:
Tropical climates at the Last Glacial Maximum: A new synthesis of terrestrial palaeoclimate data. I. Vegetation, lake levels and geochemistry.
Climate Dyn.
,
15
,
823
856
.
Fasullo
,
J. T.
, and
K. E.
Trenberth
,
2008
:
The annual cycle of the energy budget. Part II: Meridional structures and poleward transports.
J. Climate
,
21
,
2313
2325
.
Gerber
,
M.
, and
F.
Joos
,
2010
:
Carbon sources and sinks from an ensemble Kalman filter ocean data assimilation.
Global Biogeochem. Cycles
,
24
,
GB3004
.
doi:10.1029/2009GB003531
.
Gerber
,
M.
,
F.
Joos
,
M.
Vazquez-Rodriguez
,
F.
Touratier
, and
C.
Goyet
,
2009
:
Regional air-sea fluxes of anthropogenic carbon inferred with an ensemble Kalman filter.
Global Biogeochem. Cycles
,
23
,
GB1013
.
doi:10.1029/2008GB003247
.
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
.
Gherardi
,
J. M.
,
L.
Labeyrie
,
S.
Nave
,
R.
Francois
,
J. F.
McManus
, and
E.
Cortijo
,
2009
:
Glacial–interglacial circulation changes inferred from 231Pa/230Th sedimentary record in the North Atlantic region.
Paleoceanography
,
24
,
PA2204
.
doi:10.1029/2008PA001696
.
Goosse
,
H.
,
T. J.
Crowley
,
E.
Zorita
,
C. M.
Ammann
,
H.
Renssen
, and
E.
Driesschaert
,
2005
:
Modelling the climate of the last millennium: What causes the differences between simulations?
Geophys. Res. Lett.
,
32
,
L06710
.
doi:10.1029/2005GL022368
.
Gordon
,
C.
,
C.
Cooper
,
C. A.
Senior
,
H.
Banks
,
J. M.
Gregory
,
T. C.
Johns
,
J. F. B.
Mitchell
, and
R. A.
Wood
,
2000
:
The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments.
Climate Dyn.
,
16
,
147
168
.
Hartmann
,
D. L.
,
1994
:
Global Physical Climatology.
Academic Press, 411 pp
.
Hasumi
,
H.
, and
S.
Emori
,
2004
:
K-1 coupled GCM (MIROC) description.
K-1Tech. Rep. 1, Center for Climate System Research, University of Tokyo, 34 pp
.
Hibler
,
W. D.
,
1979
:
A dynamic thermodynamic sea ice model.
J. Phys. Oceanogr.
,
9
,
815
846
.
Isemer
,
H.
,
J.
Willebrand
, and
L.
Hasse
,
1989
:
Fine adjustment of large-scale air–sea energy flux parameterizations by direct estimates of ocean heat transport.
J. Climate
,
2
,
1173
1184
.
Jackson
,
C. S.
, and
A. J.
Broccoli
,
2003
:
Orbital forcing of Arctic climate: Mechanisms of climate response and implications for continental glaciation.
Climate Dyn.
,
21
,
539
557
.
Jansen
,
E.
, and
Coauthors
,
2007
:
Palaeoclimate.
Climate Change 2007: The Physical Science Basis, S. Solomon et al., Eds., Cambridge University Press, 433–497
.
Jouzel
,
J.
, and
Coauthors
,
2007
:
Orbital and millennial Antarctic climate variability over the past 800,000 years.
Science
,
317
,
793
796
.
Key
,
R. M.
, and
Coauthors
,
2004
:
A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP).
Global Biogeochem. Cycles
,
18
,
GB4031
.
doi:10.1029/2004GB002247
.
Kukla
,
G.
, and
D.
Robinson
,
1980
:
Annual cycle of surfae albedo.
Mon. Wea. Rev.
,
108
,
56
68
.
Labeyrie
,
L. D.
,
J-C.
Duplessy
,
J.
Duprat
,
A.
Juillet-Leclerc
,
J.
Moyes
,
E.
Michel
,
N.
Kallel
, and
N. J.
Shackleton
,
1992
:
Changes in the vertical structure of the North Atlantic Ocean between glacial and modern times.
Quat. Sci. Rev.
,
11
,
401
413
.
Levitus
,
S.
, and
T.
Boyer
,
1994
:
Temperature.
Vol. 4, World Ocean Atlas 1994, NOAA Atlas NESDIS 4, 117 pp
.
Levitus
,
S.
,
R.
Burgett
, and
T.
Boyer
,
1994
:
Salinity.
Vol. 3, World Ocean Atlas 1994, NOAA Atlas NESDIS 3, 99 pp
.
Li
,
S.
, and
A.
Jarvis
,
2009
:
Long run surface temperature dynamics of an A-OGCM: The HadCM3 4×CO2 forcing experiment revisited.
Climate Dyn.
,
33
,
817
825
.
Lisiecki
,
L. E.
, and
M. E.
Raymo
,
2005
:
A Pliocene–Pleistocene stack of 57 globally distributed benthic δ18O records.
Paleoceanography
,
20
,
PA1003
.
doi:10.1029/2004PA001071
.
Liu
,
Z.
, and
Coauthors
,
2009
:
Transient simulation of last deglaciation with a new mechanism for Bølling-Allerød warming.
Science
,
325
,
310
314
.
Loulergue
,
L.
, and
Coauthors
,
2008
:
Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years.
Nature
,
453
,
383
386
.
Lumpkin
,
R.
, and
K.
Speer
,
2007
:
Global ocean meridional overturning.
J. Phys. Oceanogr.
,
37
,
2550
2562
.
Lüthi
,
D.
, and
Coauthors
,
2008
:
High-resolution carbon dioxide concentration record 650,000–800,000 years before present.
Nature
,
453
,
379
382
.
Lynch-Stieglitz
,
J.
, and
Coauthors
,
2007
:
Atlantic meridional overturning circulation during the Last Glacial Maximum.
Science
,
316
,
66
69
.
Mahowald
,
N. M.
,
D. R.
Muhs
,
S.
Levis
,
P. J.
Rasch
,
M.
Yoshioka
,
C. S.
Zender
, and
C.
Luo
,
2006
:
Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates.
J. Geophys. Res.
,
111
,
D10202
.
doi:10.1029/2005JD006653
.
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
.
doi:10.1029/2006GC001323
.
Marchitto
,
T. M.
,
S. J.
Lehman
,
J. D.
Ortiz
,
J.
Flückiger
, and
A.
van Geen
,
2007
:
Marine radiocarbon evidence for the mechanism of deglacial atmospheric CO2 rise.
Science
,
316
,
1456
1459
.
Martin
,
F.
,
2002
:
Klimamodellierung mit einem vereinfachten Modell höherer Auflösung. M.S. thesis, Abteilung für Klima- und Umweltphysik, Physikalisches Institut der Universität Bern, 89 pp
.
Martin
,
P. A.
,
D. W.
Lea
,
Y.
Rosenthal
,
N. J.
Shackleton
,
M.
Sarnthein
, and
T.
Papenfuss
,
2002
:
Quaternary deep sea temperature histories derived from benthic foraminiferal Mg/Ca.
Earth Planet. Sci. Lett.
,
198
,
193
209
.
Martrat
,
B.
,
J. O.
Grimalt
,
N. J.
Shackleton
,
L.
de Abreu
,
M. A.
Hutterli
, and
T. F.
Stocker
,
2007
:
Four climate cycles of recurring deep and surface water destabilizations on the Iberian Margin.
Science
,
317
,
502
507
.
McPhee
,
M.
,
1992
:
Turbulent heat flux in the upper ocean under sea ice.
J. Geophys. Res.
,
97
,
5365
5379
.
Meissner
,
K. J.
, and
R.
Gerdes
,
2002
:
Coupled climate modelling of ocean circulation changes during ice age inception.
Climate Dyn.
,
18
,
455
473
.
Montoya
,
M.
,
A.
Griesel
,
A.
Levermann
,
J.
Mignot
,
M.
Hofmann
,
A.
Ganopolski
, and
S.
Rahmstorf
,
2005
:
The earth system model of intermediate complexity CLIMBER-3α: Part I: Description and performance for present-day conditions.
Climate Dyn.
,
25
,
237
263
.
Müller
,
S. A.
,
F.
Joos
,
N. R.
Edwards
, and
T. F.
Stocker
,
2006
:
Water mass distribution and ventilation time scales in a cost-efficient, three-dimensional ocean model.
J. Climate
,
19
,
5479
5499
.
Müller
,
S. A.
,
F.
Joos
,
G-K.
Plattner
,
N. R.
Edwards
, and
T. F.
Stocker
,
2008
:
Modeled natural and excess radiocarbon: Sensitivies to the gas exchange formulation and ocean transport strength.
Global Biogeochem. Cycles
,
22
,
GB3011
.
doi:10.1029/2007GB003065
.
Muscheler
,
R.
,
F.
Joos
,
J.
Beer
,
S. A.
Müller
,
M.
Vonmoos
, and
I.
Snowball
,
2007
:
Solar activity during the last 1000 yr inferred from radionuclide records.
Quat. Sci. Rev.
,
26
,
82
97
.
Myrhe
,
G.
,
E. J.
Highwood
,
K. P.
Shine
, and
F.
Stordal
,
1998
:
New estimates of radiative forcing due to well mixed greenhouse gases.
Geophys. Res. Lett.
,
25
,
2715
2718
.
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
.
doi:10.1029/2007GL029475
.
Pardaens
,
A. K.
,
H. T.
Banks
,
J. M.
Gregory
, and
P. R.
Rowntree
,
2003
:
Freshwater transports in HadCM3.
Climate Dyn.
,
21
,
177
195
.
Parekh
,
P.
,
F.
Joos
, and
S. A.
Müller
,
2008
:
A modeling assessment of the interplay between aeolian iron fluxes and iron-binding ligands in controlling carbon dioxide fluctuations during Antarctic warm events.
Paleoceanography
,
23
,
PA4202
.
doi:10.1029/2007PA001531
.
Peltier
,
W. R.
,
1994
:
Ice age paleotopography.
Science
,
265
,
195
201
.
Randall
,
D. A.
, and
Coauthors
,
2007
:
Climate models and their evaluation.
Climate Change 2007: The Physical Science Basis, S. Solomon et al., Eds., Cambridge University Press, 589–662
.
Rayner
,
N. A.
,
D. E.
Parker
,
E. B.
Horton
,
C. K.
Folland
,
L. V.
Alexander
,
D. P.
Rowell
,
E. C.
Kent
, and
A.
Kaplan
,
2003
:
Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century.
J. Geophys. Res.
,
108
,
4407
.
doi:10.1029/2002JD002670
.
Reimer
,
P. J.
, and
Coauthors
,
2004
:
IntCal04 terrestrial radiocarbon age calibration, 0–26 cal kyr BP.
Radiocarbon
,
46
,
1029
1058
.
Ritz
,
S. P.
,
T. F.
Stocker
, and
S. A.
Müller
,
2008
:
Modeling the effect of abrupt ocean circulation change on marine reservoir age.
Earth Planet. Sci. Lett.
,
268
,
202
211
.
Semtner
,
A. J.
,
1976
:
A model for the thermodynamic growth of sea ice in numerical investigations of climate.
J. Phys. Oceanogr.
,
6
,
379
389
.
Shi
,
G.
,
1992
:
Radiative forcing and greenhouse effect due to the atmospheric trace gases.
Sci. China
,
35B
,
217
229
.
Siddall
,
M.
,
T. F.
Stocker
,
G. M.
Henderson
,
F.
Joos
,
M.
Frank
,
N. R.
Edwards
,
S. P.
Ritz
, and
S. A.
Müller
,
2007
:
Modeling the relationship between 231Pa/230Th distribution in North Atlantic sediment and Atlantic meridional overturning circulation.
Paleoceanography
,
22
,
PA2214
.
doi:10.1029/2006PA001358
.
Sitch
,
S.
, and
Coauthors
,
2003
:
Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model.
Global Change Biol.
,
9
,
161
185
.
Stenni
,
B.
,
V.
Masson-Delmotte
,
S.
Johnsen
,
J.
Jouzel
,
A.
Longinelli
,
E.
Monnin
,
R.
Rothlisberger
, and
E.
Selmo
,
2001
:
An oceanic cold reversal during the last deglaciation.
Science
,
293
,
2074
2077
.
Stephens
,
G. L.
,
1978
:
Radiation profiles in extended water clouds. II: Parametrization schemes.
J. Atmos. Sci.
,
35
,
2123
2132
.
Stocker
,
T. F.
,
D. G.
Wright
, and
L. A.
Mysak
,
1992
:
A zonally averaged, coupled ocean–atmosphere model for paleoclimate studies.
J. Climate
,
5
,
773
797
.
Strassmann
,
K. M.
,
F.
Joos
, and
G.
Fischer
,
2008
:
Simulating effects of land use changes on carbon fluxes: Past contributions to atmospheric CO2 increases and future commitments due to losses of terrestrial sink capacity.
Tellus
,
60B
,
583
603
.
Thompson
,
S. L.
, and
S. G.
Warren
,
1982
:
Parameterization of outgoing infrared radiation derived from detailed radiative calculations.
J. Atmos. Sci.
,
39
,
2667
2680
.
Trenberth
,
K. E.
, and
J. M.
Caron
,
2001
:
Estimates of meridional atmosphere and ocean heat transports.
J. Climate
,
14
,
3433
3443
.
Tschumi
,
T.
,
2009
:
Modeling the ocean’s contribution to past and future changes in global carbon cycling. Ph.D. thesis, University of Bern, 150 pp
.
Tschumi
,
T.
,
F.
Joos
, and
P.
Parekh
,
2008
:
How important are Southern Hemisphere wind changes for low glacial carbon dioxide? A model study.
Paleoceanography
,
23
,
PA4208
.
doi:10.1029/2008PA001592
.
Tuenter
,
E.
,
S.
Weber
,
F. J.
Hilgen
,
L. J.
Lourens
, and
A.
Ganopolski
,
2005
:
Simulation of climate phase lags in response to precession and obliquity forcing and the role of vegetation.
Climate Dyn.
,
24
,
279
295
.
Weaver
,
A. J.
, and
Coauthors
,
2001
:
The UVic Earth System Climate Model: Model description, climatology, and application to past, present, and future climates.
Atmos.–Ocean
,
39
,
361
428
.
Weber
,
S. L.
, and
Coauthors
,
2007
:
The modern and glacial overturning circulation in the Atlantic Ocean in PMIP coupled model simulations.
Climate Past
,
3
,
51
64
.
Weijer
,
W.
,
W. P. M.
De Ruijter
,
H. A.
Dijkstra
, and
P. J.
van Leeuwen
,
1999
:
Impact of interbasin exchange on the Atlantic overturning circulation.
J. Phys. Oceanogr.
,
29
,
2266
2284
.
Zaucker
,
F.
,
T. F.
Stocker
, and
W. S.
Broecker
,
1994
:
Atmospheric freshwater fluxes and their effect on the global thermohaline circulation.
J. Geophys. Res.
,
99
,
12443
12457
.

APPENDIX

Ice Sheet Scaling

Ice sheets for the modern and LGM states are taken from Peltier (1994) and linearly scaled using the benthic δ18O stack of Lisiecki and Raymo (2005). High-frequency variability of the δ18O dataset was removed by applying a spline fit according to Enting (1987) with a cutoff period of 10 kyr. The ice sheet is assumed to have a constant height hice. At the ice sheet margin the height goes abruptly from hice to 0. Here hice is calculated by solving

 
formula

where Aice is the fractional ice sheet–covered area per box for the LGM and modern case, respectively, noc is the number of surface ocean boxes, and Δhoc = 120 m is the sea level rise since the LGM.

For the ice sheet interpolation we distinguish between the Northern and Southern Hemisphere. In the Northern Hemisphere, the sum of the fractional ice sheet amount at longitude ϕ and time t is calculated as follows:

 
formula

where η = (δ18O(t) − δ18Omod)/(δ18OLGMδ18Omod) and δ18Omod and δ18OLGM are the modern and LGM δ18O values. When η is between 0 and 1 then the ice sheet is interpolated between LGM and modern conditions; when η < 0 or η > 1 the ice sheet is extrapolated. The ice is latitudinally distributed starting at the northernmost land box. If, for example, ∑ϑ>0 Aice (ϕeg, teg) = 2.3, then Aice of the two northernmost land boxes at ϕ = ϕeg is set to 1 and Aice of the box to the south is set to 0.3. The method is analogously used in the Southern Hemisphere.

In the simulations where freshwater relocation from the ocean onto the land (and vice versa) is taken into account, the freshwater required for the ice sheet buildup is taken out of the ocean globally, whereas melted ice flows along the precipitation runoff path.

Footnotes

Corresponding author address: Stefan P. Ritz, University of Bern, Physics Institute, Climate and Environmental Physics, Sidlerstr. 5, Bern 3012, Switzerland. Email: ritz@climate.unibe.ch