Despite slow rates of ocean mixing, observational and modeling studies suggest that buoyancy is redistributed to all depths of the ocean on surprisingly short interannual to decadal time scales. The mechanisms responsible for this redistribution remain poorly understood. This work uses an Earth system model to evaluate the global steady-state ocean buoyancy (and related steric sea level) budget, its interannual variability, and its transient response to a doubling of CO2 over 70 years, with a focus on the deep ocean. At steady state, the simple view of vertical advective–diffusive balance for the deep ocean holds at low to midlatitudes. At higher latitudes, the balance depends on a myriad of additional terms, namely mesoscale and submesoscale advection, convection and overflows from marginal seas, and terms related to the nonlinear equation of state. These high-latitude processes rapidly communicate anomalies in surface buoyancy forcing to the deep ocean locally; the deep, high-latitude changes then influence the large-scale advection of buoyancy to create transient deep buoyancy anomalies at lower latitudes. Following a doubling of atmospheric carbon dioxide concentrations, the high-latitude buoyancy sinks are suppressed by a slowdown in convection and reduced dense water formation. This change is accompanied by a slowing of both upper and lower cells of the global meridional overturning circulation, reducing the supply of dense water to low latitudes beneath the pycnocline and the commensurate flow of light waters to high latitudes above the pycnocline. By this mechanism, changes in high-latitude buoyancy are communicated to the global deep ocean on relatively fast advective time scales.
Given measured and modeled rates of turbulent mixing in the open ocean, the time scale for a property anomaly to mix across the thermocline is many centuries.1 Yet, recent observational analyses suggest there is significant variability in global average deep ocean temperature and buoyancy in the last 20 years, a phenomenon reproduced even in low-mixing models (Palmer et al. 2011). For example, an analysis of repeat hydrographic sections shows basin-scale temperature gains below 2000 m during the 1990–2000s (Purkey and Johnson 2010; Kouketsu et al. 2011). Likewise, ocean heat uptake beneath 700 m is required to reconcile satellite-observed sea level rise with estimates of change in ocean mass and near-surface steric height (Domingues et al. 2008; Song and Colberg 2011). Deep ocean warming is important in the planetary heat budget, which cannot be balanced without accounting for the ocean beneath 700 m (Katsman and van Oldenborgh 2011; Meehl et al. 2011). In agreement with these observational analyses, a recent modeling study suggests that heat can be redistributed to all depths of the modeled ocean on decadal time scales, although the physical mechanisms for this redistribution are unclear (Palmer et al. 2011). Thus, the appearance of deep buoyancy anomalies raises the same question that has been asked for decades about the steady-state buoyancy budget [see Kuhlbrodt et al. (2007) for a review]: What mechanisms transport buoyancy beneath the ocean’s thermocline?
The primary hypothesis offered to explain the variability of deep ocean buoyancy gives a leading role to deep and bottom water formation (Purkey and Johnson 2010; Meehl et al. 2011). It is assumed that a slowdown in the formation of cold water masses at high latitudes would impact the ocean’s deep buoyancy by reducing the convective mixing that removes buoyancy from the ocean interior. However, this explanation leaves intact a key question: What mechanisms warm the vestigial deep water mass, allowing the slowdown in high-latitude convection to rapidly influence the global deep ocean? Recent studies have suggested that circulation changes are of paramount importance in determining change to the ocean’s heat budget under global warming (Banks and Gregory 2006; Lowe and Gregory 2006; Pardaens et al. 2011; Xie and Vallis 2012; Winton et al. 2013).
Here we add to this body of work by evaluating the physical processes that set the ocean’s steady-state buoyancy (and, thus, the steric sea level) budget, its variability, and its transient response to a doubling of CO2 in two simulations of an earth system model [the Geophysical Fluid Dynamics Laboratory’s (GFDL’s) Earth System Model with the Modular Ocean Model (ESM2M)]. In the control simulation, atmospheric CO2 is held at preindustrial concentrations of 286 ppm. In the CO2-doubling simulation, atmospheric CO2 increases at a compounded rate of 1% yr−1 until doubling at 70 years, after which the concentration is held steady. During the doubling of CO2, the warming of the ocean is surface-intensified, but its signature extends to the ocean bottom (Fig. 1). Over 20% of the globally integrated steric sea level change is due to the warming of the water column beneath the pycnocline (here taken as 1270 m), similar to observation-based estimates (Church et al. 2011; Ponte 2012).
One novelty of this study is our access to the full diagnosis of terms in the buoyancy budget, including those due to nonlinearities in the equation of state, which cannot be accurately calculated offline from temperature, salinity, and transport fields. While related process-based budgets have been calculated for temperature (e.g., Gregory 2000; Banks and Gregory 2006; Wolfe et al. 2008), these buoyancy diagnostics are uniquely available for the GFDL ocean model, so our results are limited to the study of this model. As such, some results are likely to be model dependent, with choices of tuned parameterizations influencing ocean circulation and the buoyancy budget. A desirable follow-up would be to compare our results with other simulations, especially those that resolve the mesoscale motions entirely parameterized here.
In section 2, we describe the model and these diagnostic tools. In section 3, we address three related questions. First, what processes set the global ocean’s steady-state buoyancy budget (section 3a)? Second, what controls variability of deep ocean buoyancy (section 3b)? Finally, which processes cause the deep ocean to gain buoyancy under a doubling of CO2 (section 3c)? Conclusions are presented in section 4.
a. The model and the buoyancy budget diagnostics
The global model used in this study is ESM2M, which was developed at GFDL for the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) (Dunne et al. 2012). The ocean component in ESM2M evolved from the ocean model used in CM2.1, as documented and compared against observations in several studies (Griffies et al. 2005; Delworth et al. 2006; Gnanadesikan et al. 2006; Dunne et al. 2012). The ocean model’s horizontal grid is nominally 1° square and is initialized from the Levitus climatology (Levitus et al. 1994). Results shown here are taken from 100-yr simulations following a spinup of more than 2000 yr (see Dunne et al. 2012), after it has reached an approximate steady state (as further described in section 3a). We analyze the 100-yr means of the control experiment to understand the steady-state buoyancy budget, and the interannual variability of annual mean quantities to characterize the variability on shorter time scales. In our analysis of the changes due to a doubling of CO2, we compare the average of years 61–80 of the CO2-doubling experiment (where year 70 is the final year of rising atmospheric CO2) to the same years from the control experiment. The model configuration is briefly described in the appendix; for additional details on the ocean model and its buoyancy budget, see Griffies and Greatbatch (2012).
This study takes advantage of a suite of model diagnostics that expose the processes affecting the time derivative of potential temperature Θ, salinity S, and thus locally referenced potential density ρ. The tracer equations for potential temperature and salinity take the material form
where J consists of subgrid-scale fluxes arising from lateral and diasurface mixing fluxes, with neutral and vertical diffusion providing the most important examples. Here, SΘ and SS are the sum of sources and sinks of potential temperature and salinity, with the nonlocal transport term from the K-profile parameterization (KPP) of Large et al. (1994) being an important example. The material derivative is defined as
where v† is the three-dimensional velocity field including both resolved motions (v) as well as a parameterization of unresolved motion, often called eddy-induced or quasi-Stokes velocity (v*). Consequently, it is the effective velocity, also called the residual mean velocity,
that transports tracers in our ocean model.
The model evaluates this expression at every time step and for every process that influences potential temperature and salinity by evaluating the derivatives, and , at the grid cell center and then multiplying by the finite-volume integrated expressions for the flux convergence and sources (Griffies 2012). Thus, the diagnostics give the tendency of the locally referenced potential density for every process that contributes to this budget.
Buoyancy, b, has long been used to study ocean circulation and energetics [for a useful review, see Vallis (2006, chapter 16)]; it is related to density via
where g is the gravitational acceleration and ρo is a reference density, and where the density tendency arises just from changes in temperature and salinity. To facilitate comparisons with Gnanadesikan et al. (2005), we omit division by the reference density in our buoyancy budgets, so our depictions of the buoyancy budget differ from the density budget only by a factor of −g.
Terms that influence the ocean’s density field also change the steric sea level, ηsteric, according to the following equation [see Eq. (13) in Griffies and Greatbatch (2012) for a derivation]:
Equation (8) describes the evolution of local steric sea level as the full water column integral of the density tendency from the ocean’s bottom (−H) to its free surface (η). This evolution is affected by changes in Θ, S, and pressure, with our focus on the influence of Θ and S alone, as we are interested in the buoyancy budget and its relation to sea level changes. When discussing the deep ocean contribution to the steric sea level tendency, the right-hand side of Eq. (8) is broken into the sea level tendency due to the surface ocean above and below a chosen depth, −D:
We have chosen a depth D = 1270 m, as it is the top of the shallowest model grid cell found below the global pycnocline. Because sea level change is an ocean property of consequence to human activity and is directly observable by satellite, many of the following results are presented in terms of the process-driven budget for steric sea level.
b. A note about subgrid-scale processes
Equation (5) describes the temporal evolution of ocean density arising from a variety of resolved and parameterized processes, and each process is further described in the appendix. This section offers brief descriptions essential for understanding the modeled buoyancy budget.
The eddy-induced or quasi-Stokes velocity (v*) is parameterized by two schemes, both of which appear as restratifying streamfunctions that adiabatically overturn isopycnals from the vertical to the horizontal. The effect of mesoscale eddies is parameterized by Gent and McWilliams (1990, hereafter GM) as modified by Ferrari et al. (2010), and uses a spatially and temporally variable thickness diffusivity. Submesoscale mixed layer eddies are parameterized using the Fox-Kemper et al. (2011) scheme (hereafter referred to as submesoscale).
The model’s diagnostic for the vertical diffusion term [included in J in Eq. (5)] includes contributions from background vertical diffusion with a time-invariant diffusivity, parameterizations of tidal mixing, and enhanced vertical diffusivity in regions of shear and gravitational instability. Therefore, convection influences the vertical diffusion term, as specified in the KPP boundary layer scheme of Large et al. (1994), as does tidally enhanced mixing in regions of rough topography (Simmons et al. 2004; Lee et al. 2006). To assess the relative importance of the background diffusion relative to the parameterizations dealing with convection and tidal mixing, an offline estimate of this term is calculated as
where kυ is the background vertical diffusivity, which is equal to 10−5 m2 s−1 equatorward of 35° latitude, and makes a tanh transition to 1.5 × 10−5 m2 s−1 at higher latitudes (Dunne et al. 2012). Because kυ is constant in time, the vertical diffusion of heat and salt is linear and can be accurately diagnosed from temperature and salinity output. However, the partial derivatives of density with respect to temperature and salinity are temporally variable, so the resultant density tendency from this offline diagnostic neglects this source of nonlinearity.
When surface buoyancy forcing drives convection above a stratified water column, buoyant water is entrained from the pycnocline into the mixed layer (Large et al. 1994), an effect not captured by enhanced diffusivity in regions of unstable stratification. This characteristic of convection motivates the introduction of a nonlocal transport, where “nonlocal” refers to the region of the water column where turbulent fluxes do not depend on the local gradients in tracers. The KPP scheme of Large et al. (1994) includes this term over a depth of the ocean determined by a bulk Richardson number criterion. Hence, convective mixing directly impacts two separate terms: the familiar, downgradient, vertical diffusion of buoyancy and the nonlocal term associated with convection (hereafter the nonlocal convective term).
Neutral diffusion [also included in J in Eq. (5)] is often referred to as isopycnal mixing; it represents turbulent tracer mixing in the direction of constant locally referenced potential density (McDougall 1987a). In these simulations, a constant neutral diffusivity of 600 m2 s−1 was chosen because it reduced biases in the climate model (Griffies et al. 2005). Given that neutral diffusion mixes tracers in the direction of constant locally referenced potential density, it is surprising that such mixing can play a role in the temporal evolution of density. However, seawater’s nonlinear equation of state gives rise to density source and sink terms in the ocean interior when mixing crosses temperature, salinity, and pressure gradients on a neutral surface. These processes are referred to as cabbeling and thermobaricity. Cabbeling describes the fact that mixing two parcels of equal density but different temperature and salinity creates a parcel that is denser than either one. Thermobaricity arises due to the pressure dependence of the coefficient of thermal expansion and can create a source or sink of density when heat is mixed along sloping isopycnals, but it integrates globally to a large buoyancy sink (Klocker and McDougall 2010). Both processes are defined as occurring only due to mixing in neutral directions (McDougall 1987b), according to
where An is the neutral diffusion coefficient, ∇nθ and ∇np are the neutral gradients of temperature and pressure, and C and T are the cabbeling and thermobaric coefficients. In addition to the source/sink terms from cabbeling and thermobaricity, our implementation of neutral diffusion transitions exponentially to horizontal diffusion in regions of steep neutral slopes, thus including a dianeutral component in these regions.
3. Results and discussion
a. The global ocean buoyancy budget near steady state
1) The steric sea level budget for the global ocean
For the average over the final 100 years of the control simulation, the global steric sea level is very nearly at steady state: it increases at a global average rate of 0.04 mm yr−1, associated with an upward temperature drift of approximately 5 × 10−5 K yr−1. However, the sea level tendency diagnosed online according to Eq. (5) and averaged over the 100-yr simulation would suggest that sea level is falling by 2 mm yr−1 (Table 1), an error arising from the approximation of the partial derivatives of density with respect to temperature and salinity [see Eq. (5)]. These derivatives are evaluated at the fixed temperature, salinity, and pressure of the grid cell center, an approximation that neglects higher-order terms of a Taylor series expansion. Although these errors place an important caveat on our study, the budget remains useful to understand which processes are quantitatively important to the ocean’s sea level budget and their variability, as the leading individual terms are generally much larger than the diagnostic error, and the individual processes influencing buoyancy [i.e., the right-hand side of Eq. (5)] sum with high precision to the total time tendency of buoyancy (i.e., the left-hand side). Furthermore, because the diagnosed density tendency drifts to a similar degree in both the control and CO2-doubling simulations, their difference allows a meaningful assessment of the processes governing the ocean’s response to warming.
The trend in the diagnosed sea level tendency in Table 1 results from the residual of several large, nearly compensating processes, namely a steric sea level rise due to surface heat fluxes and its reduction due to transport processes. Freshwater crossing the ocean’s surface and salt fluxes from melting sea ice change salinity in the upper ocean and thus impact density, but generally to a much lesser degree than the heat fluxes.
The tendency for surface heating to increase steric sea level by 9 mm yr−1 (Table 1) may appear surprising at first, given that the control simulation’s ocean heat budget is very near steady state. The steric sea level increase caused by surface heating is a consequence of the strong temperature dependence of the thermal expansion coefficient, α, where
Because of the meridional gradient in sea surface temperature, α is more than 10 times larger at low latitudes than high latitudes (Fig. 2a). Thus, surface heat fluxes that are nearly balanced via ocean warming at low latitudes and cooling at high latitudes provide a large source of buoyancy to the global ocean. This process is equivalent to adding volume to the ocean at the impressive rate of roughly 0.1 Sv (1 Sv = 106 m3 s−1), about equal to half the inflow of the Amazon River (similar to the value of 0.07 Sv estimated from observational reanalysis products of air–sea buoyancy fluxes; Schanze and Schmitt 2013).
In the case of a steady state for ocean volume, the source of buoyancy by surface forcing is balanced by densification due to processes that transport heat across gradients of α (McDougall and Garrett 1992; Griffies and Greatbatch 2012). Cabbeling is not the principal buoyancy sink when defined due only to mixing along isopycnals [Eq. (10)], as in McDougall (1987b). Rather, in the control simulation, vertical diffusion provides the biggest buoyancy sink by mixing heat out of the warm surface ocean (where α is large) and into the cold ocean interior (where α tends to be smaller). When the heat is mixed beneath the thermocline into the unstratified interior, this buoyancy sink becomes less effective, since α increases with pressure beneath the thermocline (see Fig. 2b).
2) The globally integrated buoyancy budget as a function of depth
The ocean’s buoyancy budget offers a useful framework for considering the energetics of ocean circulation. Analyzing the buoyancy budget, Gnanadesikan et al. (2005) revealed that the transport of buoyancy by the large-scale, long-term mean flow supplies buoyancy to the deep ocean by tilting isopycnals into the vertical, as was suggested by Gregory (2000). The time-varying anomalies from this mean flow (i.e., mesoscale eddies) remove this buoyancy from the deep ocean by laying those isopycnals flat, as schematized in Fig. 3. The competition between the mean and eddying components of deep ocean advective heat transport was later confirmed in models that explicitly resolve, rather than parameterize, relevant mesoscale motions (Wolfe et al. 2008; F. Bryan 2013, unpublished manuscript).
The terms of the globally integrated buoyancy budget as a function of depth and separated by temperature and salinity is shown in Fig. 4. Positive values reflect a supply of buoyancy below each depth. In agreement with previous studies (Gregory 2000; Gnanadesikan et al. 2005; Wolfe et al. 2008), the large-scale mean advection of warm waters supplies buoyancy to the ocean interior (thin blue line in Fig. 4a). These warm waters are apparently saltier than the waters they displace, such that the large-scale advection of salinity makes the ocean interior denser (thin blue line in Fig. 4b), albeit at a slower rate than the buoyancy increase due to the temperature transport. The total buoyancy supply by the resolved advection (Fig. 4c) is larger than the vertical diffusive supply between 300 and 1000 m and nearly the same size beneath 1000 m. The sum of the parameterized mesoscale and submesoscale advection (v*) removes buoyancy at a rate that exceeds the supply by the resolved motions (Fig. 4c). Cabbeling and thermobaricity, which cause the majority of the buoyancy sink due to neutral diffusion (see the end of section 2b), destroy buoyancy in the ocean interior. The component of the convective removal of buoyancy that is represented by the nonlocal KPP term is small beneath 400 m. Geothermal heating (not shown) adds buoyancy at a slower rate than any of the terms depicted in Fig. 4, although it is thought to play an indirect role in the buoyancy budget by reducing the strength of abyssal stratification (Mashayek et al. 2013).
Our estimate of the buoyancy sink due to cabbeling is about an order of magnitude smaller than that estimated by Gnanadesikan et al. (2005). In contrast with our approach, this earlier study associated large convergence of vertical density fluxes in the interior with cabbeling. In so doing, their analysis separated the impact of lateral neutral fluxes and vertical neutral fluxes, which in our analysis cancel each other significantly. A recent study provides a valuable comparison between an observationally based estimate of cabbeling and thermobaricity (Klocker and McDougall 2010) and our model diagnostics. The global integral of the velocities from the observational estimate is approximately 5 Sv of dense water formed at a neutral density of 1027.4 kg m−3 (with this neutral surface having a mean depth of 1000 m) from these nonlinear processes. For comparison, we estimate dianeutral velocities from cabbeling and thermobaricity, w, from the density tendency diagnostics, that is, dividing the sum of Eqs. (11) and (12) by the local stratification:
The dianeutral velocity due to cabbeling and thermobaricity integrates globally to roughly a fifth of the Klocker and McDougall (2010) estimate: 3 Sv at 270 m, decreasing to 1 Sv at 800 m and 0.1 Sv at 2600 m. Because the observational estimates of cabbeling and thermobaric velocities are linearly dependent on the choice of diffusion coefficient [assumed to be a constant 1000 m2 s−1 in Klocker and McDougall (2010)], they would be 40% lower if calculated using the diffusivity in our simulations (600 m2 s−1), and somewhat closer to our results.
Although our simulations may suggest a smaller role for cabbeling and thermobaricity than inferred in previous analyses (Gnanadesikan et al. 2005; Klocker and McDougall 2010), their sum remains first order in the buoyancy budget for the deep Southern Ocean and can play a role in moderating local rates of local sea level rise, as discussed in the following sections.
3) The zonal-mean steric sea level budget
The ocean’s role in the global redistribution of heat—transporting heat received at low latitudes toward the poles—is apparent in a depiction of the column integrated, zonal-mean ocean steric sea level budget (Fig. 5a). This figure differs from a simple depiction of heat transport because multiplication by a highly variable thermal expansion coefficient creates pronounced meridional gradients in steric sea level anomaly. Between 15°S and 20°N the heat flux from the atmosphere expands the water column. In steady state, advection transports this light water to higher latitudes where heat exchange with the atmosphere removes ocean buoyancy. At all latitudes north of 20°N, the ocean is made denser via heat loss to the atmosphere. In contrast, between 40° and 60°S, the ocean is warmed and freshened: this region of high-latitude ocean buoyancy gain is a reflection of the wind-driven upwelling of deep water masses in the Antarctic Circumpolar Current (ACC) that are cooler than the overlying atmosphere.
The zonal-mean steric sea level budget for the deep ocean (>1270 m) reveals the physical processes that transport buoyancy below the pycnocline (Figs. 5b,c). Throughout the low latitudes (~30°S–40°N), the steady-state buoyancy budget is qualitatively as Munk (1966) hypothesized: vertical diffusion lightens the deep ocean at the same rate that resolved advection makes it denser (Fig. 5b). At high latitudes, the balance is dramatically different.
In the Southern Ocean, the resolved and parameterized components of the advective buoyancy transport strongly compensate one another, with resolved advection supplying buoyancy to the deep ocean and parameterized eddies removing it (Fig. 5c), consistent with (Wolfe et al. 2008). South of 40°S, the balance tips slightly in favor of the resolved motions such that the resultant residual mean advection (v†) is generally a supply term of buoyancy to the deep Southern Ocean, as is the vertical diffusion. As first hypothesized on the basis of observational data (Gill 1973), the cabbeling and thermobaric destruction of buoyancy is an important sink term in the buoyancy budget of the Southern Ocean (Fig. 5b).
In the simulated high-latitude North Atlantic, vertical diffusion supplies buoyancy to the deep ocean more than an order of magnitude faster than elsewhere (visible in the zonal mean of Fig. 5b). The large role for the vertical diffusion may be tied to the fact that Mediterranean Outflow Water allows the middepths of the North Atlantic to fill with warm, salty waters (Lozier et al. 1995). This warm Mediterranean-influenced layer sits on top of the cold, relatively fresh waters formed in subpolar regions. Therefore, the midlatitude North Atlantic has the world’s strongest vertical temperature and salinity gradients at intermediate depths, where parameterizations of tidal mixing raise the vertical diffusivities. Furthermore, because the thermal expansion coefficient, α, is pressure dependent (and this pressure dependence is the principal control on α outside of the tropical thermocline; see Fig. 2b), the deeper heat is mixed, the more it adds to the water column’s buoyancy. The steric sea level tendency due to background diffusion (diagnosed offline; see section 2a) provides little of the total sea level tendency due to vertical diffusion (recall Fig. 4); thus, the components of vertical diffusion due to tidal parameterizations and water column instability must be of leading importance.
The role of deep water formation is apparent in the important buoyancy sink provided by the nonlocal convective term and the mixing with dense marginal seawater. At steady state at 60°N, the buoyancy source due to the strong vertical diffusion is balanced in equal parts by the buoyancy sinks from the nonlocal convective term, parameterized mixing with the dense overflows from the marginal seas ringing the North Atlantic, and GM and submesoscale advection. North of 60°N, a large advective source of buoyancy to the deep North Atlantic is balanced in steady state by mixing with the marginal seas (Fig. 5b).
b. Causes of temporal variability in deep steric sea level
The importance of the deep ocean for the temporal evolution of steric sea level has been raised on the basis of several observational studies (Domingues et al. 2008; Purkey and Johnson 2010; Song and Colberg 2011). In agreement with these observational studies, the globally integrated deep steric sea level of our control simulation regularly swings by as much as 0.3 mm yr−1 (Fig. 6). The time series of global average deep (>1270 m) steric sea level tendency shows that resolved advection dominates deep steric sea level budget (Fig. 6): the correlation between the globally averaged steric sea level tendency due to resolved advection and the total time rate of change is r2 = 0.81 in the control simulation.
To evaluate which terms locally set the variability in the deep ocean buoyancy budget, we examine the independent linear regressions between each process term and the local time rate of change, holding the other terms constant, according to
where the left-hand side is the time series of local buoyancy tendency and the right-hand side is the time series of each buoyancy transport process or source/sink multiplied by its regression coefficient, m. A regression coefficient close to unity signals that the size of the transport or source/sink term is similar to the size of the total time rate of change, and the correlation coefficient gives an indication of how well a linear equation describes the relationship. Figure 7 synthesizes these statistics for the global deep (>1270 m) ocean. The key message from these global maps is that resolved advection sets variability in the buoyancy tendency budget almost everywhere in the deep ocean, except in convective regions where many additional terms are important. A second interesting point from Fig. 7 is that thermobaricity is as important to the deep ocean buoyancy budget and its variability as vertical diffusion.
It is not surprising that variability in the budget for deep ocean buoyancy beneath a fixed depth is dominated almost everywhere by the three-dimensional resolved advection, which causes heaving of the pycnocline (Fig. 7). Other terms—particularly parameterized mesoscale advection and vertical diffusion—are quantitatively important to the budget throughout the global deep ocean, as indicated by the slope of the regression (i.e., the size of the point), but do not contribute significantly to the temporal variability. In contrast, at high latitudes, subgrid-scale processes compete with resolved advection to set the buoyancy variability of the deep ocean. Most of this variability stems from the North Atlantic and Southern Ocean. Therefore, we focus on the balances in these key regions and their influence on lower latitudes.
1) The North Atlantic
The origins of the deep North Atlantic steric sea level anomalies can be traced to atmosphere–ocean heat exchange at high latitudes. The buoyancy sink due to surface heat forcing averaged in the North Atlantic north of 55°N is significantly correlated with the local deep steric sea level tendency at zero lag (Fig. 8). How the signal of the surface buoyancy forcing penetrates the deep ocean is a complicated question. As we saw in the zonal mean deep steric sea level budgets (Fig. 5), several buoyancy transport processes play important roles in the northern high latitudes and each of these processes varies on interannual time scales.
The buoyancy budget in convective regions such as the Labrador Sea is strongly influenced by sink terms: parameterized mesoscale and submesoscale advection, convection, and thermobaricity (Fig. 7). In the region between Iceland and the British Isles, a reduction of the parameterized overflow connecting the marginal seas and the open North Atlantic (called downslope mixing), submesoscale advection, and convection all allow for the accumulation of deep buoyancy on interannual time scales. It is the suppression of these deep buoyancy sink terms, as well as variability in the advective supply of buoyancy to the deep ocean, that causes buoyancy to accumulate in the deep, high-latitude North Atlantic on interannual time scales (Fig. 7).
Having penetrated to the deep layers of the high-latitude North Atlantic, buoyancy anomalies then propagate southward to a latitude of 20°N, evident in the high correlation between surface forcing north of 55°N and deep buoyancy tendency extending southward at increasing lag (Fig. 8). This propagation is more directly seen by tracing the Atlantic zonal mean buoyancy tendency as a function of latitude and time (Fig. 9a). The southward propagation of the buoyancy anomalies is clearly linked to the large-scale advection (Fig. 9b), which varies because of anomalies in the velocity field and the buoyancy it transports. A number of modeling and observational studies have shown similar coherent propagation on these time scales (Curry et al. 1998; Zhang 2010; van Sebille et al. 2011; Mauritzen 2012). This coherence suggests the possibility of significant lead time in predicting changes in the deep ocean contribution to regional steric sea level evolution.
2) The Southern Ocean
Resolved advection is the dominant process controlling interannual variability in deep buoyancy of the Southern Ocean (Fig. 7). Subgrid-scale processes play a first-order role only in limited regions. A related result was reported by Ito et al. (2010), who showed, using a ⅙° model simulation, that the large-scale mean Eulerian advection controlled the intra-annual variability in Southern Ocean anthropogenic CO2 budget, and on these time scales was not effectively counterbalanced by eddies, to the degree that they were resolved.
In addition to the important role of advection, interannual swings in processes associated with convection are evident (Fig. 7h) in the Weddell Sea, where deep convection is a common feature of many IPCC climate models, including GFDL's coupled climate model, similar to the one used here (de Lavergne et al. 2013, manuscript submitted to Nat. Climate Change). Here, thermobaricity also plays a leading role in setting buoyancy variability (Fig. 7f).
c. Processes governing the CO2-induced steric sea level change
The change in the buoyancy budget due to warming is diagnosed as the difference between 20-yr means in the control simulation and those in the CO2-doubling experiment during the final years of rising atmospheric CO2. As prescribed atmospheric CO2 makes its final climb to 572 ppm, global average steric sea level rises at 2 mm yr−1 faster than in the control simulation, about double the observational estimates of late-twentieth-century rates of steric sea level rise (0.88 ± 0.33 from 1993 to 2008, as observed atmospheric CO2 topped 360 ppm; Church et al. 2011). As noted in the introduction, there is a buoyancy gain at every depth of the globally integrated ocean. We can now see that the vertical penetration of buoyancy is due principally to an increase in the supply from resolved advection (Fig. 10c), and that the budget change is dominated by temperature transport (Fig. 10a) rather than salinity (Fig. 10b).
The global average vertical diffusion of buoyancy decreases slightly beneath 1000 m under a doubling of CO2 (Fig. 10), due mostly to a reduction at high latitudes (Fig. 11b). Assessing the precise cause of the slowdown is not trivial, as the vertical diffusion term is the sum of several parameterized processes. The background mixing cannot explain the slowdown in the vertical diffusive supply of buoyancy, as the constant diffusivity acting on stronger buoyancy gradients increases the supply to deep layers (Fig. 10). The buoyancy tendencies due to tidal parameterizations are generally too small to cause the slowdown, and, importantly, the spatial distribution of the changes to these terms, as deduced from changes in their associated diffusivities (not shown), does not resemble the pattern of change in the total vertical diffusion term. Thus, we are left to speculate that the most likely process causing a slowdown is a reduction of the region with enhanced diffusivity brought about by shear instability. Where the water column is stably stratified (i.e., N2 > 0) and enhanced diffusivities are triggered by a strong velocity shear, the vertical diffusion term supplies buoyancy to deeper layers of the ocean. We conjecture that the enhanced stratification under warming reduces the region where the threshold for such instability is achieved, thereby reducing diffusivities in these regions, and slowing the vertical diffusive supply of buoyancy to the deep ocean. The diagnostics needed to deduce the effects just from shear instability were not saved, so we cannot uncover for sure the distinction between each process.
There is significant meridional variability of the zonal mean steric sea level rise, with spatial patterns set by a competition between increasing atmospheric input of heat (or a decrease in the ocean loss of heat to the atmosphere) and the oceanic transport of the resultant buoyancy gain (Fig. 11a). The largest change in buoyancy forcing is in the northern high latitudes, where a reduction in ocean to atmosphere heat loss would force sea level to rise by up to 25 mm yr−1 (Fig. 11a). However, in these regions, the steric sea level rise is less than 3 mm yr−1, because the advective supply of buoyancy is reduced by a nearly equal amount (Fig. 11a). Other regions (e.g., 10°–40°S) see similar rates of steric sea level rise in spite of an unchanged regional surface buoyancy flux (Fig. 11): here, sea level rise is primarily caused by an increase in the advective supply of buoyancy. It is interesting to note that changes in the cabbeling and thermobaric terms can influence steric sea level change by as much as 3 mm yr−1 between 45° and 65°N.
Of the 2 mm yr−1 global average steric sea level rise caused by a doubling of CO2, 23% is due to a reduction in density below the pycnocline (i.e., beneath 1270 m). Observationally based estimates of the fraction of total sea level due to thermal expansion in the deep ocean are highly uncertain: a recent study suggested that the deep ocean beneath 700 m accounts for 19% of the total, but with uncertainties approximately as large as the estimate itself (Church et al. 2011), bringing it into the same range as the model estimate. Here, we are interested in what processes may cause the accumulation of buoyancy in the deep ocean. From 40°S to 50°N the buoyancy accumulation in the deep ocean is entirely due to changes in the heat and freshwater advection by the model’s resolved currents (Figs. 11b,c). At high latitudes a more complicated picture emerges.
A doubling of CO2 slows convection in the North Atlantic, as reflected by the nonlocal convective term. The end of the warmer century likewise sees a reduction in the mixing of dense water from the overflows connecting the Greenland, Norwegian, and Iceland Seas to the North Atlantic (reflected in the downslope mixing term that parameterizes the overflows; Fig. 11b). These buoyancy sink terms are strongly suppressed, but the vertical diffusive supply of buoyancy is reduced by almost the same amount. In the deep layer of the high northern latitudes, there is a nontrivial difference between the change in density tendency as diagnosed from the online approximation and the same change diagnosed from the evolution of the ocean’s density field; we believe the general evolution of each term remains relevant, despite that the compensation between the terms is more complete than the diagnostics indicate.
In the deep Southern Ocean, south of 40°S, there is a sea level trend ranging from 0 to 0.5 mm yr−1 (Fig. 11b), with contributions from changes in the advection and convection. With a doubling of CO2, the Southern Hemisphere westerly winds strengthen and move poleward (Fig. 12) with consequences for ocean circulation. Yet, for the final 20-yr mean of the CO2-doubling experiment, changes in the buoyancy supply due to mean residual advection are small south of 40°S (Fig. 12b), because changes in the GM buoyancy transport compensate changes in the resolved advection almost perfectly (Fig. 11c). This compensation on the century time scale is distinct from the analysis of interannual variability, in which the parameterized eddies do not seem to keep pace with variability in resolved advection [as discussed in section 3b(2)]. Importantly, the strength of this compensation may exert a leading control on ocean heat uptake efficiency across different climate models (Kuhlbrodt and Gregory 2012). The question of how the GM parameterization responds to increasing wind strength has been raised on the basis of model studies at various resolutions, but often with a much stronger wind perturbation than the one predicted here (Hallberg and Gnanadesikan 2006; Farneti et al. 2010; Farneti and Gent 2011). Recent results from a coupled model simulation with a ° ocean suggest that the GM parameterization is effective at mimicking the net effect of mesoscale eddies on ocean heat transport and their changes due to a doubling of atmospheric CO2, albeit not in every detail (Bryan et al. 2013).
In addition to the change in the advective terms, a reduction in convection, evident in the decrease of the nonlocal convective term (Fig. 11b), causes an increase in ideal age of the deep Southern Ocean (Fig. 12b). Ideal age is an ocean tracer that is set to zero everywhere at the surface and then ages at a rate of 1 yr yr−1 in the interior. The entire deep Southern Ocean is older at the end of the warmer century, while the surface ~1400 m are younger to the north of 50°S, consistent with a slowdown in deep convection and attendant reduction in the upward flux of old waters south of 65°S. It appears that the weakening convection is slamming the door to the deep ocean shut, overwhelming the ability of the strengthened winds to prop it open (Russell et al. 2006). Thus, we may need to adjust our expectations regarding the response of Southern Ocean heat and carbon storage to a warming climate.
Zonal sections showing the impact of advection on the buoyancy budget overlaid with overturning streamfunctions help to summarize the key impacts of CO2 doubling on buoyancy (Figs. 13a,b). The shallow overturning cells above 500 m transport large amounts of buoyancy from the equator to midlatitudes; this transport strengthens slightly at the end of the warmer century. Below this surface layer and above 3000 m, the upper cell of the meridional overturning circulation (MOC) converges dense water (i.e., reduces buoyancy) in the subsurface ocean at all latitudes between 40°S and 60°N (Fig. 13a). With a doubling of CO2, this cell of the MOC slows to the north of 20°S, effectively reducing the spread of dense water and allowing the accumulation of buoyancy (Fig. 13b), similar to conclusions of previous studies (Gregory 2000; Banks and Gregory 2006; Winton et al. 2013). To the south of 20°S, the upper cell of the MOC strengthens slightly, and advection becomes a supply term of buoyancy to the deep ocean. The bottom cell of the MOC, transporting Antarctic Bottom Water to the low latitudes, slows and reduces the supply of dense water below 3000 m. Changes in vertical diffusion are large and highly variable above 1000 m, but in the deep ocean these changes are smaller than changes in advection (Figs. 13c,d). In the control simulation, the sum of the other subgrid-scale processes provides a buoyancy sink almost everywhere (Fig. 13e); that sink is substantially reduced at high latitudes in both hemispheres after a doubling of CO2 (Fig. 13f).
An analysis of the buoyancy budget in the ocean component of an Earth system model at quasi-steady state and in response to a doubling of CO2 has yielded a number of insights into the time-mean circulation, the temporal variability of steric sea level, and its evolution under a doubling of CO2, many of which are schematized in Fig. 14. Although both observational (Domingues et al. 2008; Church et al. 2011; Katsman and van Oldenborgh 2011; Song and Colberg 2011) and modeling studies (Palmer et al. 2011) have shown that buoyancy can be redistributed throughout the full water column on time scales of years to decades, the mechanisms responsible for this redistribution have remained poorly understood.
The apparent importance of mesoscale and submesoscale advection and turbulent vertical and lateral mixing implies that we are far from having the capacity to diagnose the complete buoyancy budget on large spatial scales from observations, and models will continue playing a crucial role for advancing our understanding. Still, our results are subject to the limitations of both the model and the density tendency diagnostics, which use linear approximations that give rise to error. Results are likely sensitive to tuned parameterizations: for instance, the global buoyancy budget is strongly dependent on the choice of coefficient for the GM parameterization and vertical diffusivity (Gnanadesikan 1999; Xie and Vallis 2012).
Vertical diffusion is a leading source of buoyancy to the deep ocean in steady state, consistent with Munk (1966). Although advection is classically framed as the buoyancy sink that balances the diffusion, we find that this framework is only relevant for the global ocean only if we consider the residual mean advection (Gnanadesikan et al. 2005; Wolfe et al. 2008). The total advective sink of buoyancy in the deep ocean is due to the imbalance between buoyancy supply by the large-scale mean circulation and its removal due to mesoscale and submesoscale motions (Fig. 14a; see the sea level tendencies beneath 1270 m). These motions are parameterized in this model, but similar results were found in the eddy-resolving models of Wolfe et al. (2008) and F. Bryan (2013, unpublished manuscript).
Throughout the low and midlatitudes, interannual to decadal variability of deep ocean buoyancy (and, thus, the deep contribution to steric sea level) is controlled entirely by resolved advection. Vertical diffusion constantly supplies buoyancy to the deep low-latitude ocean. In steady state, this diffusive buoyancy supply is balanced by the advective injection of dense waters (Fig. 14a; see the deep layer between 18°S and 25.5°N). So, an anomaly need not be mixed across the pycnocline to increase the buoyancy of the deep low-latitude ocean. Rather, changes in high-latitude buoyancy poleward of the thermocline’s wintertime surface outcrop are communicated with the global deep ocean on much faster advective time scales.
In the North Atlantic’s Deep Water formation regions, variability in surface heat fluxes creates buoyancy anomalies that penetrate the deep ocean via the slowing or cessation of a number of buoyancy sink terms, namely open-ocean convection, mixing with the overflow waters from the Greenland–Norwegian Seas, mesoscale and submesoscale advection, and the thermobaric destruction of buoyancy. These perturbations then influence the circulation, such that buoyancy anomalies propagate to low latitudes due to changes in the resolved advection. This advective perturbation creates spatially coherent variability in the deep North Atlantic on decadal time scales.
With a doubling of CO2, atmospheric buoyancy forcing changes dramatically in the high-latitude North Atlantic and near the equator, but steric sea level grows at similar rates at all latitudes, due to changes in the meridional advection (Fig. 14b). A slowdown in convection and dense water formation in marginal seas suppresses the principal deep ocean buoyancy sink in the high latitudes of both hemispheres. This change is accompanied by a slowdown in the lower cell of the overturning circulation and an attendant reduction in the supply of Antarctic Bottom Water to the deep low latitudes, in spite of an increase in the Southern Hemisphere westerlies. Similarly, a slowdown in the upper cell of the overturning circulation in the Northern Hemisphere denies the deep ocean its supply of dense North Atlantic Deep Water, thereby increasing the buoyancy of the deep low latitudes. In contrast, the upper cell of the overturning circulation is slightly strengthened in the Southern Ocean, and here the residual mean advection becomes a source of buoyancy to the deep ocean, by removing dense waters adiabatically. While the vertical mixing of heat is incapable of rapidly communicating a buoyancy anomaly to the deep ocean, perturbations to the global circulation allow for the redistribution of buoyancy at all depths of the ocean on interannual to decadal time scales.
The authors thank Mike Winton, Carolina Dufour, Daniele Bianchi, Jonathan Gregory, Frank Bryan, and an anonymous reviewer for helpful suggestions. JBP gratefully acknowledges funding from NOAA’s Global Carbon Cycle program and Canada’s NSERC Discovery Program.
Ocean Model Configuration
The ocean component of the ESM2M simulations used here is the Modular Ocean Model (MOM4p1), configured using the nominal 1° horizontal resolution of GFDL Climate Model version 2.1 (CM2.1; Delworth et al. 2006) and ESM2M (Dunne et al. 2012) and generalized level coordinates (Stacey et al. 1995; Adcroft and Campin 2004; Griffies 2012). We use the volume-conserving Boussinesq configuration with a z* vertical coordinate, defined as
where z = −H(x, y) is the ocean bottom and z = η(x, y, t) is the ocean sea surface height. The latitudinal resolution is refined to ⅓° at the equator, thus leading to 200 latitudinal rows in total. The ocean has 50 vertical levels, with 22 levels in the upper 220 m.
A more detailed discussion of each buoyancy tendency term and its impact on global sea level in a quasi-steady state, Coordinated Ocean-Ice Reference Experiment (CORE)-forced (Griffies et al. 2009; Large and Yeager 2009) ocean-ice model is offered by Griffies and Greatbatch (2012). The purpose of this appendix is to briefly describe the processes in the ocean model that influence buoyancy, dominated by the following seven terms:
1) Boundary fluxes of heat and salt. The boundary flux of freshwater does not influence ocean density until it mixes with the top ocean layer, thereby changing its density. Shortwave radiation is allowed to penetrate through the water column using a climatological chlorophyll field as updated from the methods detailed in Sweeney et al. (2005). The shortwave attenuation with ocean depth is based on the seawater optics of Manizza et al. (2005). Here we present the sum of the boundary fluxes of temperature and the penetrating component of shortwave radiation as surface heating (Table 1; Figs. 5a, 11a, and 14).
2) The mean advective transport of heat and salt. The model assumes that the advective transport can be decomposed into a mean component, which is resolved in the coarse-grained model, and components due to mesoscale and submesoscale motions, which are parameterized.
3) The quasi-Stokes transport from mesoscale eddies. This transport is implemented using the skew flux approach of Griffies (1998). The quasi-Stokes streamfunction is computed via a boundary value problem extending across the full ocean column according to Ferrari et al. (2010), which contrasts with the local approach of Gent and McWilliams (1990) and Gent et al. (1995). The horizontal variation of the eddy diffusivity is based on vertically averaged flow properties (Griffies et al. 2005) with an allowable range of values specified between 100 and 800 m2 s−1.
4) The parameterization of submesoscale eddy-induced mixed layer restratification (Fox-Kemper et al. 2011).
5) Neutral diffusion is based on Griffies et al. (1998) with a constant diffusivity of 600 m2 s−1 and the neutral slope tapering scheme of Danabasoglu and McWilliams (Danabasoglu and McWilliams 1995), using for the maximum slope. Neutral diffusion contributes to the evolution of neutral density due to the effects of nonlinearities in the equation of state, namely the processes of cabbeling and thermobaricity. In addition, our implementation of neutral diffusion provides for an exponential transition to horizontal diffusion in regions of steep neutral slopes, thus including a dianeutral component in such regions. This transition to horizontal diffusion is made in regions where the surface boundary layer is encountered, following the recommendations from Treguier (1992), Ferrari et al. (2008), and Ferrari et al. (2010) and next to solid walls, following from the recommendations of Gerdes et al. (1991).
6) Ocean vertical mixing. Vertical mixing is parameterized in various schemes: the K-profile parameterization scheme for the ocean boundary layer from Large et al. (1994), also described in Griffies (Griffies 2012); parameterizations of tide mixing from Simmons et al. (2004) and Lee et al. (2006); enhanced vertical diffusivity in regions of gravitational or shear instability; and a specified background diffusivity of 10−5 m2 s−1 in low latitudes and transitioning to 1.5 × 10−5 m2 s−1 poleward of 35°. The KPP scheme also provides a nonlocal mixing term that redistributes heat and salt in the boundary layer in the presence of negative surface buoyancy forcing. The vertical tracer and momentum diffusivity are enhanced in regions of gravitational or shear instability. There is no further convective adjustment scheme.
7) Tracer mixing due to parameterizations of overflow processes and inland sea exchange, with details given in Griffies et al. (2005) and Griffies (2012). These processes contribute a negligible effect to global mean sea level (see Table 1), but can be important in the deep ocean (Fig. 14), particularly in the North Atlantic, which is surrounded by marginal seas.
Additional terms influence the simulated density field, but are far smaller than the dominant terms listed above when averaged globally. They are as follows:
8) As river water enters the ocean, it is mixed into the upper four grid cells with the ambient water. The river water is assumed to have zero salinity and have the same temperature as the sea surface. This mixing generally increases the seawater buoyancy, and so provides a positive but small contribution to sea level.
9) Because of the computational null mode present on a B grid (Mesinger 1973; Killworth et al. 1991), the sea surface is smoothed using a Laplacian filter. This smoothing necessitates a corresponding flux of buoyancy in order to conserve tracer [see section 12.9 of Griffies (2004)].
10) The sea ice model assumes a uniform salinity of five parts per thousand [see the appendix to Delworth et al. (2006)]. Thus exchange with the sea ice model can provide a source or sink of salt and freshwater with a corresponding impact on density.
This scale analysis is based on the vertical diffusion equation, where a time rate of change in temperature, T, arises as the result of turbulent diffusion alone with a coefficient kυ = 10−5 m2 s−1: . The time, Δt, for a temperature anomaly to mix across a distance, Δz, scales as , requiring on the order of 1000 years for an anomaly to spread 1000 m.