Model and observational studies have concluded that geothermal heating significantly alters the global overturning circulation and the properties of the widely distributed Antarctic Bottom Water. Here two distinct geothermal heat flux datasets are tested under different experimental designs in a fully coupled model that mimics the control run of a typical Coupled Model Intercomparison Project (CMIP) climate model. The regional analysis herein reveals that bottom temperature and transport changes, due to the inclusion of geothermal heating, are propagated throughout the water column, most prominently in the Southern Ocean, with the background density structure and major circulation pathways acting as drivers of these changes. While geothermal heating enhances Southern Ocean abyssal overturning circulation by 20%–50%, upwelling of warmer deep waters and cooling of upper ocean waters within the Antarctic Circumpolar Current (ACC) region decrease its transport by 3–5 Sv (1 Sv = 106 m3 s−1). The transient responses in regional bottom temperature increases exceed 0.1°C. The large-scale features that are shown to transport anomalies far from their geothermal source all exist in the Southern Ocean. Such features include steeply sloping isopycnals, weak abyssal stratification, voluminous southward flowing deep waters and exported bottom waters, the ACC, and the polar gyres. Recently the Southern Ocean has been identified as a prime region for deep ocean warming; geothermal heating should be included in climate models to ensure accurate representation of these abyssal temperature changes.
Geothermal heating has a nonnegligible influence on the large-scale abyssal circulation by weakening abyssal stratification and increasing the circulation of Antarctic Bottom Water (AABW) and North Atlantic Deep Water (NADW) on the order of 10%–30% (e.g., Adcroft et al. 2001; Hofmann and Morales Maqueda 2009; Emile-Geay and Madec 2009; Mashayek et al. 2013; de Lavergne et al. 2016). On a regional scale, geothermal heating can change ocean bottom temperatures by an order of magnitude more than error estimates associated with decadal abyssal temperature trends (e.g., Emile-Geay and Madec 2009; Purkey and Johnson 2010; Kouketsu et al. 2011; Wunsch and Heimbach 2014) and increase thermosteric sea level by 0.1–1 mm yr−1 (Piecuch et al. 2015). Yet geothermal heating is represented inconsistently by ocean and fully coupled models. The aim of this study is to assess how the inclusion of geothermal heat fluxes at the ocean floor influences regional ocean circulation and temperature in a coarse-resolution climate model.
Previous studies have demonstrated that geothermal heating has widespread impacts on the abyssal ocean and rivals diapycnal mixing for AABW transformation. Lighter AABW, in regions of weak stratification, is more susceptible to heat gain from geothermal input, which acts to expand the bottom water incrop area at the ocean floor (Emile-Geay and Madec 2009; de Lavergne et al. 2016). In addition, the impacts of geothermal heating vary regionally, which comes as no surprise considering that the ocean circulation patterns, water mass properties, and stratification differ between basins. For example, while the deep western boundary current contributing to the formation of NADW is warmed by geothermal heating, its temperature properties are subject to continuous buoyancy loss via surface heat fluxes (Adcroft et al. 2001). Conversely, Pacific and Indian Ocean sourced deep waters are older and have minimal contact with the ocean surface, allowing heat anomalies to accumulate over time. The conversion of geothermal heat–based gravitational potential energy to kinetic energy influences large-scale circulation differently in the Atlantic and Indo-Pacific regions (Urakawa and Hasumi 2009; Adkins et al. 2005). Here we focus on both global and regional Southern Ocean responses to geothermal heating, from the ocean surface to floor.
Few studies have focused on the Southern Ocean, where geothermal heating increases poleward heat transport by ~10% (Emile-Geay and Madec 2009) and further steepens isopycnals associated with Southern Ocean upwelling (Mashayek et al. 2013). Several model studies also lack realistic topography (essential for representation of stratification along the ocean floor and midocean ridges outputting high geothermal heat), and none to date uses a fully coupled model configuration where atmospheric–oceanic feedbacks are included. Our Southern Ocean–focused study is the first, to our knowledge, that uses a fully coupled global climate model, inclusive of realistic topography and parameterized diapycnal mixing (particularly in the bottom boundary layer), to diagnose the impacts of spatially varying geothermal heating. We assess two substantially different geothermal heat flux datasets to conclusively diagnose the transient responses in the global ocean.
The inclusion of geothermal heating in an ocean climate model is straightforward operationally. However, some models exclude this heat flux for one of two reasons. First, it has been conjectured that geothermal heating has a weak impact on the globally averaged circulation and temperature trends (Mullarney et al. 2006; Purkey and Johnson 2012). While the global mean geothermal heat flux is two to three orders of magnitude smaller than the surface heat flux, it is important to note that the geothermal heat flux is a destabilizing buoyancy source at the ocean floor, and is applied to slow moving deep ocean water masses that have limited surface outcrop regions (Emile-Geay and Madec 2009). Numerous geothermal heat sources along midocean ridges have been identified since the late 1970s (e.g., Lupton 1998; Veirs et al. 1999). The resultant hot plumes rise hundreds of meters above the ridge crest, undergo diverse chemical and physical processes, and are carried across ocean basins by the large-scale circulation (Speer and Helfrich 1995; Lupton 1998; Downes et al. 2012). As will be shown in this manuscript, geothermal heating influences the ocean from surface to floor in regions far from its midocean ridge maxima.
The second reason for the exclusion of the ocean geothermal heat flux from models is that there is no conclusive and globally complete conductive heat flow dataset. In contrast to continental regions, the oceanic heat flux estimates are not well observed. Estimates in regions of sparse data are based on experimental data and theoretical half-space cooling models that estimate the heat transfer within the oceanic lithosphere via conduction using the age of the oceanic crust (cf. Pollack et al. 1993; Hofmeister and Criss 2005). Recently, the global mean geothermal heat flux has been estimated in the range of 56–95 mW m−2 with considerable spatial variation (Hofmeister and Criss 2005; Hamza et al. 2008; Davies and Davies 2010; Goutorbe et al. 2011). Here we assess two geothermal heat flux datasets, namely that of Hamza et al. (2008) and Davies and Davies (2010), that differ primarily along midocean ridges and thus at hydrothermal vent sites.
However, model studies have shown that regardless of the magnitude of the geothermal heat flux added to the system, from 50 mW m−2 (e.g., Adcroft et al. 2001) to 100 mW m−2 (Urakawa and Hasumi 2009), whether it is applied as spatially varying or uniform (e.g., Emile-Geay and Madec 2009; Mashayek et al. 2013) or irrespective of the model configuration used (idealized to ocean-ice general circulation model), it is evident that “geothermal heating is an important factor of abyssal dynamics, and should no longer be neglected in oceanographic studies” (Emile-Geay and Madec 2009, p. 203).
Here we quantify the impact of geothermal heating across the Southern Ocean in a global coupled climate model with realistic topography and two distinct spatially varying geothermal heat flux datasets. We demonstrate that temperature changes occur far from the prime geothermal heat sources, that they propagate throughout the water column, and that they vary spatially depending on local stratification and major deep ocean circulation pathways. Furthermore, temperature differences induced by geothermal heat are similar in magnitude to (if not higher than) the deep ocean changes observed over the past few decades, implying that models without geothermal heating will have biases that exceed the changes they are designed to detect. After describing the model configuration (section 2), we quantify the global and Southern Ocean changes due to geothermal heating over centennial time scales using the Hamza et al. (2008) geothermal heat flux (sections 3a and 3b) and the Davies and Davies (2010) geothermal heat flux (section 3c), then discuss the implications of our results in section 4. In section 5, we conclude our study and reiterate the required inclusion of geothermal heating in climate models.
2. The CM2M
a. Model features
We use the National Oceanic and Atmospheric Administration (NOAA)/Geophysical Fluid Dynamics Laboratory (GFDL) CM2M—a fully coupled climate model configuration with realistic surface boundary conditions. The Earth system version of CM2M (ESM2M) contributed to the phase 5 of the Coupled Model Intercomparison Project (CMIP5) suite of model hindcasts and projections (Dunne et al. 2012). The atmospheric component (AM2; Anderson et al. 2004) is based on that used in previous GFDL coupled models, with minor bug fixes. The latitude by longitude horizontal grid resolution of AM2 is 2° × 2.5° with 24 levels in the vertical. The sea ice component is that of Winton (2000), and the land component (LM3.0) is described in Milly et al. (2014).
Previous geothermal heating studies that were run for millennial time scales were not fully coupled (e.g., Adcroft et al. 2001; Emile-Geay and Madec 2009), and thus lacked atmospheric feedbacks on the ocean that arise when geothermal heating impacts the ocean surface (e.g., Adcroft et al. 2001; Mashayek et al. 2013; Piecuch et al. 2015). Here we show that geothermal heating–induced anomalies do extend toward the surface (thus impacting surface water mass transformation), particularly in polar regions where approximated atmospheric forcing errors can be substantial (e.g., Nygard et al. 2016; Hobbs et al. 2016). In an ocean-ice model, the approximated forcing components have their own errors distinct from the model that can enhance or compete with the model errors, with no clear way to separate the two error types (cf. Griffies et al. 2009).
The ocean component is the Modular Ocean Model version 4.1 (MOM4p1; Griffies 2012) with a 1° horizontal grid resolution that tapers to ⅓° at the equator, and a bipolar grid north of 65°N. The model has 50 vertical levels using the z* geopotential coordinate, ranging from 10 m in thickness in the upper 220–250-m thick at the ocean floor. The CM2M agrees well with the observed density structure of the deep Southern Ocean. The model Antarctic Circumpolar Current (ACC) transport through Drake Passage averaged over the full control simulation (1300 yr) is 153 Sv (1 Sv = 106 m3 s−1), which is at the higher end of observationally based estimates of 135–155 Sv (cf. Cunningham et al. 2003; Griesel et al. 2012; Chidichimo et al. 2014).
Interactions between the ocean circulation and topographical features, in particular along midocean ridges, result in enhanced diapycnal mixing (Waterhouse et al. 2014). Mixing in the ocean interior and bottom is primarily governed by the breaking of internal gravity waves, lee waves, mixing associated with overflows, and mixing arising from the geothermal heating buoyancy input. Given that mixing is intricately linked to the overturning circulation (Nikurashin and Ferrari 2013; Mashayek et al. 2015), we detail how ocean mixing is parameterized herein, and discuss our results in terms of model mixing schemes in the appendix.
Turbulence arising from internal wave breaking (sourced from barotropic tidal energy) is parameterized in CM2M using the Simmons et al. (2004) scheme. In the Simmons et al. (2004) parameterization, a third of the energy dissipated from internal wave breaking is applied locally, with the remaining two-thirds of the energy dissipated nonlocally via a background diffusivity. The constant background diffusivity is 0.1–0.2 × 10−4 m2 s−1, and is stronger poleward of 30° (Dunne et al. 2012). Frictional bottom drag dissipation is parameterized following Lee et al. (2006), with the maximum associated dianeutral diffusivity set at 5 × 10−3 m2 s−1. The transport of fluid within the bottom water layer down the steep polar continental shelves is parameterized using a downslope mixing scheme detailed in Snow et al. (2015).
b. Model simulations
To assess how geothermal heating impacts Southern Ocean dynamics we run four experiments, all using the same CM2M fully coupled configuration, with or without constant geothermal heat fluxes:
A 1300-yr control simulation without geothermal heat (CTL);
A simulation including geothermal heat as defined by the spatially varying dataset of Hamza et al. (2008), starting from year 701 of the CTL, and run for 600 years (GH2); and
In all four experiments, the CM2M is run using 1860 preindustrial radiative forcing. The 1860 forcing is typically used in the spinup phase of climate models for international efforts such as CMIP5. Additionally, a constant preindustrial forcing allows us to diagnose the influence of geothermal heating in the absence of climate change. The first two experiments (CTL and GH1) are initialized from a 1-yr run of the ocean–sea ice system using Conkright et al. (2002) temperature and salinity fields. Bottom waters circulate the global ocean on centennial to millennial time scales. Thus, while the 500- (GH1) and 600-yr (GH2 and GHDD) perturbations are perturbed from a nonequilibrated ocean control state, it is an adequate length to identify the transient responses of the abyssal circulation.
Given the large range of geothermal heating applied to previous model studies (global mean of 50–100 mW m−2), we conducted an additional experiment, termed GHDD, where the GH2 case is rerun using the significantly higher geothermal heat fluxes of Davies and Davies (2010) (Fig. 1c). The results of the GHDD experiment are presented in section 3c. The Davies and Davies (2010) dataset uses a half space cooling model along midocean ridges where data are sparse, the oceanic crust is young, and hydrothermal systems are abundant; the Southern Ocean is one such region. The Davies and Davies (2010) method estimates a geothermal heat flux along midocean ridges that is at least double to that found by Hamza et al. (2008) using a spherical harmonic analysis globally (Fig. 1c). Globally, the mean geothermal heat flux in the Hamza et al. (2008) dataset is 64 mW m−2, with a mean of 95 mW m−2 in the Davies and Davies (2010) dataset.
In the results sections that follow, we compare the first 500 years of the CTL simulation (CTL1) with the GH1 case, and the CTL simulation during years 701–1300 (CTL2) with the GH2 and GHDD cases. The GH2 and GHDD cases are run for longer than the GH1 case because the CTL ocean is in a different state and we chose the centuries where the global mean temperature changes between GH1 and GH2 were similar in the respective final centuries (Fig. 2). The GH2 and GHDD cases are initialized from year 701 in the CTL simulation. While the preceding 700 CTL yr is insufficient to fully spin up the deep water mass properties, it is sufficient to allow us to examine how geothermal heating modifies deep water properties in the GH2 and GHDD cases that follow.
3. Influence of geothermal heating
Our results are divided into three main sections. Section 3a describes the global influence of geothermal heating on bottom temperature and meridional overturning circulation. We then narrow our focus to the Southern Ocean changes over the ocean floor and throughout the water column (section 3b). These first two sections compare the GH1 and GH2 perturbation experiments. Section 3c describes the GHDD perturbation results.
a. Global influence of geothermal heating
We begin by discussing the global impact of adding geothermal heating over the 500- (GH1) and 600-yr (GH2) integrations, with a primary focus on the final century of these perturbations. The global impacts of geothermal heating can be largely predicted based on the global mean rate of geothermal heating [e.g., see section D.4 of Griffies et al. (2009)]. On a global mean scale, the difference in ocean temperature at the end of the GH1 and GH2 experiments is merely ~0.06°C (Fig. 2)—an order of magnitude less than the 500-yr global mean ocean temperature change in the CTL1 case alone.
However, the regional impacts of geothermal heating are heterogeneously distributed so that temperature changes in some locations are significantly larger than the global average (Figs. 3 and 4). We find changes to the bottom temperature (and salinity) begin within the first century of GH1 and GH2. During the first three centuries, the ocean warming (excluding in the Pacific) is strongest in regions of large geothermal heating, namely along midocean ridges. However, bottom waters take at least 100 years to fill the Southern Ocean, with deep waters of northern origin not reaching the Southern Ocean for a few hundred years. Thus, it is not until the fourth century of the GH1 and GH2 runs that we find interbasin exchange of temperature and salinity anomalies via deep and bottom waters. It is important to reiterate that the transient results presented here are a function of model state and may change as the model equilibrates (discussed further in section 4).
A striking result is that regions of strongest geothermal heat flux input (along midocean ridges) do not coincide with regions where the mean ocean floor temperature difference is greatest between the geothermal heating (Figs. 3 and 4). Rather, excess heat accumulates in the deep ocean basins, primarily in the Pacific and Indian Oceans. In particular, the bottom of the ocean in the southeast Indian Ocean, the northeast Pacific, the North Atlantic, and the Pacific sector of the Southern Ocean is more than 0.15°C warmer as a result of the inclusion of geothermal heating in the final century (years 401–500 for GH1 and years 1201–1300 for GH2). The GH2 anomalies are larger than those found for the GH1 case, particularly in the Atlantic basin. The acute warming within the Antarctic Circumpolar Current region in the Pacific basin in both the GH1 and GH2 cases is likely associated with geothermal heating originating along the Pacific–Antarctic Ridge and East Pacific Rise. Hydrothermal plumes from both these midocean ridges have been observed to transport ocean properties along the ACC and south into the Ross Sea (Downes et al. 2012), and eastward into the Atlantic sector (Naveira Garabato et al. 2007). Huang (1999) suggests that the generation rate for potential energy associated with geothermal heating is greater at depth and away from shallower midocean ridges from where the greatest geothermal heat flux is found. Similarly, here we find the greatest bottom temperature anomalies are located away from the ridge axes.
Figure 5 provides a zonally integrated view of the influence of geothermal heat on the meridional overturning circulation (MOC) in density space. The MOC is defined as the combined time-mean Eulerian and parameterized eddy velocities integrated along the longitude and density coordinates [in units of Sverdrups (Sv)]. In the first century of the CTL case (Fig. 5a) the upper cell in the ACC region (maximum at ~50°S; σ2 = 1036.3 kg m−3) is 18 Sv, with 16 Sv for the subpolar cell (south of 45°S; σ2 > 1036.9 kg m−3), and 7 Sv for the lower cell where southward flowing circumpolar deep waters are returned northward via AABW (centered near 30°S; 1036.8 < σ2 < 1037 kg m−3). The overturning associated with North Atlantic Deep Water production (NADW cell) is 26 Sv. The upper and NADW ocean MOC structures are in agreement with several previous model and observationally based estimates (e.g., Lumpkin and Speer 2007; Marshall and Speer 2012); however, the model lower cell is smaller, as is often the case in coarse-resolution models (e.g., Downes et al. 2011; Downes and Hogg 2013; Downes et al. 2015; Farneti et al. 2015).
Given that the ocean state is different at year 1 versus year 701, the spatial pattern of the anomalous bottom ocean temperature in the GH1 and GH2 cases evolves differently. To understand the differences in the global MOC response between the GH1 and GH2 cases, it is important to evaluate the difference between the CTL cases in the first century of each perturbation (Fig. 5a). Specifically, by the seventh century (green contours in Fig. 5a), the upper and NADW cells have strengthened by more than 4 Sv, with the former meridionally shifting position. The counterclockwise circulating subpolar and lower cells have doubled in strength.
Applying geothermal heat at the ocean floor has minimal impact on the overturning circulation for σ2 < 1036.6 kg m−3, north of 50°S for both the GH1 (Fig. 5b) and GH2 (Fig. 5c) cases. There are large circulation changes isolated to the equatorial latitudes for σ2 > 1036.3 kg m−3, associated with a southward shift in the overturning cells. In the GH1 perturbation, the upper cell in the ACC region increases by ~10% (around 2 Sv), which is in agreement with other model studies assessing geothermal heating impacts (e.g., Adcroft et al. 2001; Mullarney et al. 2006). The lower cell, north of 50°S, strengthens largely on its northern and southern sides (i.e., the cell expands in size), attributable to local weaker stratification.
However, there is an overall large decrease of at least 3 Sv in the subpolar cell (σ2 > 1036.9 kg m−3), attributable to the upwelling of warmer waters on the Antarctic continental shelf. Approximately one-quarter of this weakening is associated with a decrease in the eddy-induced overturning circulation. The increase in the subpolar cell is double the size in the GH2 case, possibly associated with the stronger CTL2 overturning circulation state (i.e., the subpolar cell increases by similar amounts relative to its CTL years; Fig. 5a). The modeled changes in the Southern Ocean overturning circulation are of a similar magnitude to those found in multimodel studies assessing historical decadal trends and twenty-first-century projections (Downes and Hogg 2013; Farneti et al. 2015). However, in these studies, both the upper and subpolar cells strengthen in a warmer climate due to changes in surface wind and buoyancy forcing (whereas here we use a constant preindustrial atmospheric forcing).
Turning from a global to regional view, we focus on three vertical cross sections (identified in the final panels of Figs. 3 and 4) that intersect large bottom temperature anomalies, as illustrated in Fig. 6 (focusing here on the left and center panels and the right GHDD panels in section 3c). Two results are common across the three ocean cross sections in the GH1 and GH2 experiments. First, the bottom temperature anomalies upwell toward the ocean surface in the Southern Ocean. Deep ocean overturning circulation is directly linked with upper ocean cells, and thus increased overturning at depth results in upper ocean property changes. Second, the latitudes of greatest geothermal input along each section do not necessarily correlate with regions of largest temperature anomalies. We note that in the surface layers the temperature and salinity anomalies are unlikely associated with geothermal heating, but rather model drift that spatially varies in sign.
Across the 100°E Indian Ocean section (Figs. 6a,b), bottom temperature anomalies originating along the Southeast Indian, Central Indian, and Ninety East Ridges spread eastward into the region directly south of Australia, where slow circulation and the local bathymetry trap the large temperature anomaly (Figs. 3 and 4; final century panels). The GH1 and GH2 experiments result in very similar anomalies in the Indo-Pacific (Figs. 6a,b,d,e), despite differences of ~0.5°C in deep ocean temperatures in the respective CTL1 and CTL2 states resulting in a different isopycnal structure and differently upwelling pathways for temperature anomalies within and south of the ACC region. As mentioned previously, the Atlantic bottom temperature anomalies are significantly larger in the GH2 case (Figs. 6g,h), associated with a stronger NADW overturning cell (Fig. 5); however, the temperature trends south of 40°S are similar between GH1 and GH2.
Geothermal heating results in widespread warming of the ocean floor (Figs. 3 and 4). However, the local salinity response (Fig. 7) is nonuniformly distributed across the global ocean and is primarily dependent upon the background salinity field, the distribution of deep and bottom waters, and local upwelling. While the spatial pattern of the CTL case salinity (CTL1 and CTL2) and the salinity response to geothermal heating in GH1 and GH2 are similar, the magnitude of the response differs greatly. For example, in the GH1 case bottom salinity increases within the Indo-Pacific ACC region are larger (by 0.01–0.02) than those of the GH2 case (Figs. 7a–d); however, salinity changes in the central Atlantic in the GH2 case are significantly larger than and different in sign from those found in the GH1 case (Figs. 7e,f).
In the CTL case, the ocean floor salinity is freshest in the Atlantic (associated with Weddell Sea bottom water) and saltiest in the Indian and Pacific Ocean basins north of the ACC. Under geothermal heating, younger salty Indian Ocean circumpolar deep waters increase bottom salinity, and the saltier waters are upwelled toward the ocean surface. In the Atlantic, younger Weddell Sea bottom waters (representing increased AABW input) freshen the Argentine basin, and increased lower cell overturning in the Southern Hemisphere drives salty waters throughout the abyssal Atlantic basin (particularly in the GH2 case). It is particularly in the Atlantic basin where we find bottom salinity anomalies are present from upper ocean to the ocean floor within the ACC region, where the steeply sloping isopycnal structure creates favorable conditions for upwelling of anomalies (Figs. 7g,h).
Widespread freshening in the middle of the Atlantic basin is associated with the background vertical salinity gradient. In the Atlantic, salty deep waters flowing southward at middepth overlie fresher bottom waters. Thus, the increased abyssal overturning circulation resulting from geothermal heating destabilizing the already weak stratification along the ocean floor transports the fresh bottom water signatures upward (cf. Emile-Geay and Madec 2009; Mashayek et al. 2013). Conversely, in the Pacific and Indian Oceans (Figs. 6a,b,d,e and 7a,b,d,e), the southward flowing deep waters flow along the bottom of the basins until the ACC region, thus transporting a salty bottom signature upward as overturning is increased via geothermal heat input.
Overall, the global changes associated with the inclusion of geothermal heating are similar to those previously shown by similar model studies. Specifically, geothermal heating warms the global ocean in both past and present climates, and acts to strengthen the lower overturning circulation by 10%–30% (e.g., Adcroft et al. 2001; Emile-Geay and Madec 2009; Mashayek et al. 2013; Ballarotta et al. 2015). Here, as the ocean overturning cells strengthen throughout the CTL simulation, the impact of geothermal heating increases, specifically for the Southern Hemisphere lower and subpolar cells that increase by well over 30% in the GH2 case. Evidently, it is within the Southern Ocean that temperature and salinity anomalies are transported to the surface via the steeply sloping isopycnal structure inherent to the ACC region. The following section investigates these Southern Ocean temperature and circulation responses in greater detail.
b. Geothermal heating influence over the Southern Ocean
The magnitude of the changes due to geothermal heating are based on both where geothermal heat is applied and (possibly more importantly) the ambient circulation and stratification (e.g., Speer and Rona 1989; Downes et al. 2012). In terms of the background circulation, several major deep ocean pathways in the Southern Ocean advect geothermal heat signals from the local sources to the large scale, as illustrated in the previous section. The Southern Ocean deep ocean circulation is depicted in the horizontal and vertical transports illustrated in Figs. 8a–c. The major ocean pathways are the eastward flowing Antarctic Circumpolar Current (Fig. 8a), the southward flowing circumpolar deep waters and northward flowing Antarctic Bottom Waters (Fig. 8b), and the polar Ross Gyre and Weddell Gyre.
Given the strong eastward transport and upwelling (Figs. 8a,c) within the ACC region, the position of the ACC with respect to the Antarctic continent, as well as the locations of midocean ridges (and thus maximum geothermal heat input; Fig. 1a), is important for the distribution of temperature anomalies. The ACC coincides with regions of large geothermal heat input across the Southern Ocean (Fig. 1a). The ACC transport, as has been shown in observations (e.g., Sun et al. 2011), varies along its circumpolar route, and is particularly strong when constricted by topography (e.g., fracture zones) or is forced to steer around topographical features (e.g., plateaus) (see Fig. 1d for locations of highlighted topographical features). Key locations where the ACC velocity peaks have been observed near the Falkland and Crozet Plateaus, and west of the Southeast Indian and Pacific–Antarctic Ridges (Zhang et al. 2012).
In the CTL experiment, within the ACC region and to the south, the bottom meridional flow is strongest where bottom waters sourced from the Ross Sea and Adelie Land regions flow equatorward in the southwest Pacific and southeast Indian Ocean basins, and in the western half of the Atlantic Ocean where Weddell Sea AABW is exported (Fig. 8b). Deep water of northern origin flowing into the Southern Ocean is also of importance in transporting abyssal heat fluxes, with key inflow regions being the eastern Atlantic, western Indian, and the southwest Pacific Ocean basins.
Decreases in the ACC transport of 3–5 Sv with the addition of geothermal heat fluxes in the GH1 and GH2 cases are present within the first century onward (Fig. 9). These decreases are associated with a 5%–11% decrease in the meridional density gradient across the ACC region in Figs. 6. As previously mentioned, the steep isopycnal structure within the ACC region is a prime and vast upwelling zone in the global ocean (Fig. 8c). In the GH1 and GH2 perturbations, the ACC region is heated from below, upwelling warmer waters in the southern half of the ACC (Figs. 6 and 8d,f,g,i) and increasing sea ice melt (figure not shown). Increased upwelling at the depth of the midocean ridges is evident across the Southern Ocean in the GH1 and GH2 cases (Figs. 8f,i) and corresponds with the increased overturning circulation in the lower cell (Figs. 5b,c) and the upwelling of temperature and salinity anomalies south of, and within, the ACC region (Figs. 6 and 7). More downwelling in the northern half of the ACC (see upper overturning cell changes near 50°S in Figs. 5b,c) and more upwelling of salty circumpolar deep waters (Fig. 7) in the southern half result in an equatorward shift of the now denser (colder and saltier) northern boundary of the ACC.
In much of the Southern Ocean, warm waters sit beneath cool waters, with a downgradient isopycnal temperature gradient pointing upward. As first noted by Gregory (2000), and later by Hieronymus and Nycander (2013) and Morrison et al. (2013), these deep and warm temperature anomalies are transported vertically upward and southward through the action of mesoscale eddy stirring along the steeply sloped isopycnals. This process plays a sizeable role in the ventilation of heat in the Southern Ocean [see also the study from Griffies et al. (2015)]. The eddy stirring and subsequent downgradient mixing are parameterized in our coarse-resolution simulations via neutral diffusion. The action of neutral diffusion thus acts in concert with upwelling to transport warm geothermal anomalies vertically upward within the Southern Ocean.
Geothermal heating results in increased meridional deep and bottom water transport in both the GH1 (Fig. 8e) and GH2 (Fig. 8h) cases. Interestingly, it is primarily in the southwest Pacific where we find a large region of large meridional transport changes when geothermal heat is introduced in the GH1 case (Fig. 8e), where geothermal heating locally changes the direction of bottom and deep water flow. Conversely, meridional transport increases of similar magnitude in the GH2 case are distributed more broadly across the Southern Ocean (Fig. 8h).
In summary, we find that the addition of geothermal heating into the ocean floor has a nonuniform distribution of bottom temperature and salinity anomalies, dependent on the background large-scale circulation, local stratification, and vertical temperature and salinity gradients. Bottom anomalies are upwelled to the surface within the ACC region, and the ACC itself is able to propagate temperature and salinity anomalies downstream across basins. The upwelling of warmer waters in the southern region of the ACC, and cooler, fresher waters in its northern region, results in a weakening of the net ACC transport. While the ACC transport changes are relatively small compared to the strength of the ACC current, they are of a similar magnitude to that projected by climate models in the twenty-first century under extreme climate forcing (e.g., Sen Gupta et al. 2009; Meijers et al. 2013; Downes and Hogg 2013). Geothermal heating source regions (primarily midocean ridges) do not coincide with the regions of largest temperature and salinity anomalies; furthermore, these anomalies are not the same sign or magnitude throughout the water column.
c. Evaluation of the geothermal heat flux dataset
Thus far we have noted results using a relatively weak geothermal heat flux along midocean ridges (GH1 and GH2; Hamza et al. 2008), and therefore our ocean properties and circulation changes can be regarded as lower bound estimates. Previous studies have employed higher geothermal heat fluxes, both spatially varying and spatially uniform (e.g., Adcroft et al. 2001; Emile-Geay and Madec 2009; Mashayek et al. 2013), and have noted that higher geothermal heat fluxes produce stronger anomalies. To test this conclusion and to provide a quantitative range of oceanic responses to geothermal heating, we compare our year 701–1300 results from GH2 with those of the GHDD case run over the same period using the larger Davies and Davies (2010) fluxes.
Within the ocean interior, a stronger geothermal heat flux dataset (GHDD) produces a stronger temperature (Figs. 6c,f,i) and salinity (Figs. 7c,f,i) anomalies across all basins. However, the enhanced anomalies of different sign (Figs. 8j–l) amount to a very weak overall change in the Southern Ocean for the 1201–1300-averaged period. We associate this small overall change with centennial fluctuations in circulation anomalies in the perturbation experiments. For example, the Southern Ocean and southeast Pacific anomalies in Fig. 10 are stronger for the 1101–1200 period, and the ACC transport changes fluctuate on decadal to centennial time scales in the GH2 and GHDD experiments (blue and red curves in Fig. 9). Further, the global meridional overturning circulation for the GH2 and GHDD experiments changes are very similar for the 1201–1300 period (Fig. 11).
It is evident that the stronger Davies and Davies (2010) geothermal heating will change the regional circulation over the 600-yr perturbation at a different rate to the GH2 case. Using the model age tracer (initialized at 1 at the start of the CTL experiment; Fig. 12) we find that the two GH1 and GH2 experiments including the Hamza et al. (2008) dataset produce a qualitatively similar spatial pattern for the change in bottom age. Younger waters flood the Atlantic and southeast Indian Ocean basins north of the ACC, with slower southward flowing deep waters shown in the Pacific. Conversely, there is an acute decrease in the bottom age in the GHDD experiment throughout the Pacific basin.
In summary, the bottom temperature anomalies are larger with the larger geothermal heating of the Davies and Davies (2010) dataset compared with using the Hamza et al. (2008) dataset. However, we hypothesize that the centennial variation in the spatial pattern of the ocean overturning circulation varies according to the chosen geothermal heat flux constraint, particularly in the Southern Ocean.
We have provided estimates of the transient response to geothermal heating, but our methods come with notable caveats. In particular, uncertainties in our diagnosed responses to geothermal heating stem from the model configuration and parameterizations, model drift, simulation length, and (as discussed in section 3c) the chosen geothermal heat flux constraint. For a detailed discussion of the how the mixing parameterizations influence our results, please refer to the appendix. We note that the coarse model resolution, and thus coarse model topography, does not resolve individual convective plumes, such as those identified in previous hydrothermal vent and plume observational studies (e.g., Lupton 1998; Veirs et al. 1999). Testing how different geothermal heat flux datasets invoke dynamical ocean responses via plume convection in an ultra-high-resolution global configuration is recommended, but is beyond the scope (and resources) of our study.
The bottom of the ocean includes thermal and viscous boundary layers, with a wave-induced turbulence zone above (~10–100 m; cf. Armi and Millard 1976; Garrett et al. 1993). The deepest layer in the CM2M is over 250-m thick, meaning that these boundary layers are not represented. However, the model is able to capture the strengthening of stratification along the ocean floor in the polar regions due to geothermal heating, associated with the decreases in the subpolar cell (Figs. 5 and 11). A weakening of stratification due to geothermal heating in the mid-to-high latitudes is associated with a strengthening of the lower cell and upper cell.
The 500-yr model drift is not negligible; however, the difference in drift between CTL1 and GH1 is less than the diagnosed temperature anomalies, particularly below 500 m (figure not shown). The CM2M drifts warm in the deep ocean, and this warm bias is increased when geothermal heating is added. The drift is spread unevenly in the model, with the Indian and Atlantic Ocean basins showing the greatest temporal warming associated with drift (figure not shown). Excluding the warm drift from the CTL and perturbed runs does not change the regional distribution of bottom temperature changes associated with geothermal heating, however including drift in Figs. 3 and 4 would reduce the illustrated anomalies by 10%–30%.
In the upper 500 m, the model overall drifts cool relative to the initial state, thus opposing the warming anomalies associated with geothermal heating. This cool drifting is important for the surface transformation of water masses, which is the flux across the density outcrop arising from heat and freshwater input at the ocean–atmosphere interface (cf. Downes et al. 2011). For example, the subpolar cell between the CTL evolution (Fig. 5a) and the geothermal heating perturbations (Figs. 5b,c) are of opposite sign, meaning that there is a strengthening of the subpolar cell over the CTL simulation, but a weakening of the subpolar cell for the perturbations. The upwelling of warm anomalies associated with geothermal heating and consequent weakening the subpolar cell (Figs. 5b,c) is linked to a ~5% decrease in the surface transformation of deep waters into bottom waters for densities σ2 greater than 1037.0 kg m−3 (figure not shown). For lighter densities, where the geothermal heating induced temperature anomalies are not directly transported upward in the water column, the model drift and the large zonal variation in the sign and magnitude of heat fluxes dominate.
Our control simulation ends at year 1300, and the deep ocean has not equilibrated at this point. However, we have demonstrated that regional temperature anomalies throughout the water column associated with geothermal heat input are similar across year 1–500 (GH1) and year 701–1300 (GH2 and GHDD) experiments. In addition, our global meridional overturning circulation results are similar to those reported in studies using coarse-resolution ocean-only or ocean–sea ice models that have been run for more than 1500 years (e.g., Adcroft et al. 2001; Emile-Geay and Madec 2009; Urakawa and Hasumi 2009). Many of these studies also report regional temperature changes of the same order of magnitude to those found here (i.e., around 0.2°C). Thus, while we have presented transient responses to geothermal heating, the magnitude of our responses is not unlike those found in studies with an equilibrated ocean. Here we purposely have chosen to assess a transient response in an experimental design that mimics that of the preindustrial spinup simulations of global coupled climate models contributing to international efforts, such as CMIP5 (Taylor et al. 2012). The majority of these types of fully coupled coarse-resolution models currently exclude geothermal heating. Hence, while the quantitative assessment here has identified key regions of significant circulation and temperature change, it is important to note that these anomalies would likely differ in strength and region if the simulations were extended to millennial time scales.
5. Summary and conclusions
The goal of this study was to evaluate the regional transient response of ocean temperature and circulation when geothermal heating is implemented at the ocean floor. We evaluated the responses when geothermal heating is included at the start of the model spinup phase, and 700 years into this phase. Previous model studies have illustrated that geothermal heating significantly impacts ocean properties and circulation in a zonal mean sense; our study has focused on regional changes, especially in the Southern Ocean. We have explored the influence of geothermal heating for two distinct heat flux datasets with perturbations beginning at two different points in our control simulation.
Our experiments have allowed us to conclusively show that geothermal heating into the ocean floor weakens stratification, increases abyssal circulation, and influences temperature and salinity throughout the entire water column (Figs. 5–7), unlike surface ocean warming, which strengthens the local stratification. Across the three geothermal heating experiments considered here, the large and widespread bottom temperature anomalies span the range of 0.1°–0.5°C and bottom salinity anomalies span 0.005–0.02, with meridional overturning circulation of abyssal cells increasing by 20%–50%. We have illustrated that the temperature changes induced by geothermal heating are located away from local geothermal heat sources along midocean ridges (Figs. 3 and 4). Major regions of change along the ocean floor are throughout the Atlantic, across the eastern Pacific and Indian Ocean basins, and in the Southern Ocean.
While geothermal heating influences the tracer properties and circulation of the global ocean, it is primarily within the Southern Ocean where the larger anomalies are upwelled and diffused along steeply sloping isopycnals directly toward the ocean surface. The unique dynamics associated with floor-to-surface transport of heat anomalies in the Southern Ocean include the following:
The continuous source of geothermal heat (at least 60 mW m−2) along the Southeast Indian Ridge, Pacific–Antarctic Ridge, and southern tip of the East Pacific Rise.
The intersection of these geothermal heat sources with the strong Antarctic Circumpolar Current (ACC), particularly in regions where the flow is restricted by topographic features and where its zonal velocity thus peaks.
The steeply sloping isopycnals, located within the ACC that upwell and isopycnally diffuse temperature and salinity anomalies throughout the water column.
The large export of bottom waters from the Adelie region, Ross Sea, and Weddell Sea, as well as the large southward flow of salty circumpolar deep waters, along the ocean floor, and the connection of these abyssal circulation pathways with the ACC.
Deep observations are sparse, and the community thus depends upon models to dynamically interpolate temperature trends. Because of the impacts of geothermal heating on patterns of deep and abyssal temperature, we conclude that it is critical to include such heating in numerical models aiming to estimate sensitivity to increases in atmospheric greenhouse gases. We find that the evolution of geothermal heating in the GH1, GH2, and GHDD experiments differs in spatial distribution, and recommend geothermal heating constraints be implemented from the start of climate model simulations rather than post spinup after multiple centuries. Geothermal heating is not the source of abyssal warming in recent decades, with dominant terms being vertical mixing of warmer waters and advection via upwelling (e.g., Purkey and Johnson 2012). However, both mixing and geothermal heating vary spatially, and there are likely regions where the latter (despite acting on centennial time scales) is significant in the global heat budget (Purkey and Johnson 2012; de Lavergne et al. 2016; Wunsch and Heimbach 2014).
Given that the CM2M is able to adequately capture the observed Southern Ocean large-scale circulation pathways, bottom topography, and steep isopycnal structure inherent to the ACC region, there is a strong possibility that the regional geothermal heating influences described in this study are present in the real ocean. Here the model uses preindustrial forcing, whereas the real ocean is experiencing a warming climate state, which further influences the large-scale circulation, and possibly enhances (or alters) the regional impact of geothermal heating. We suggest that the influence of geothermal heating be quantified for the current climate and for future climate scenarios where the circulation is radically altered by surface warming.
We are grateful for constructive comments by Carolina Dufour and Adele Morrison on earlier versions of the manuscript, and by Ali Mashayek and two anonymous reviewers. We thank Stephen Rintoul for discussion relating to the observed Southern Ocean dynamics and circulation pathways, and Rhodri Davies for providing the Davies and Davies (2010) dataset and associated discussion. SMD was primarily supported by the ARC Centre of Excellence for Climate System Science (Grant CE110001028), and also supported by the Australian Government’s Business Cooperative Research Centres Programme through the Antarctic Climate and Ecosystems Cooperative Research Centre (ACE CRC). AMH was supported by an Australian Research Council Future Fellowship FT120100842. This research was undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government.
Evaluation of the Model Mixing Parameterization
The CM2M has a coarse-resolution horizontal and vertical grid, thus requiring parameterizations of small-scale mixing. The higher-resolution (eddy-permitting) version of CM2M indeed results in a stronger lower overturning circulation cell in the control case (Farneti and Gent 2011); however, computation time for such a run is significantly increased, and is suggested for a future study. The parameterizations of eddy-induced overturning and bottom boundary processes in the CM2M likely impact how the model responds to geothermal heating via changes in the overturning circulation. The model vertical and horizontal mixing parameterizations are a function of vertical stratification, and thus are affected by geothermal heating. Correspondingly, the associated residual overturning circulation Ψres is tightly linked with stratification.
The residual overturning circulation is the sum of the Eulerian, wind-driven circulation and the eddy-induced circulation. Following derivations by Mashayek et al. (2013, 2015), in a zonally averaged sense, the vertical advection of buoyancy Ψresbz is balanced by its downward diffusion δz(κbz). Thus, the total overturning streamfunction in the bottom boundary layer can be approximated by
where κ is the diapycnal diffusion coefficient, and buoyancy , where ρ is the local density and is the ocean mean density.
Using the relation given in Eq. (A1), and in agreement with Mashayek et al. (2013, top row of their Fig. 3), we find that the addition of geothermal heating acts to decrease both κbzz and bz, with the latter term decreasing more than the former. This results in a larger meridional streamfunction in the GH1 and GH2–GHDD cases (Figs. 5b,c). Hence, we can demonstrate that the vertical diffusivity and buoyancy combine to influence the overturning streamfunction.
Separating the influence of the vertical diffusion coefficient and buoyancy is beyond the scope of this study. We note that previous studies have indicated that using a diffusion coefficient that exponentially decays away from the bottom enhances the responses of the overturning circulation to geothermal heating (e.g., Dutay et al. 2010; Mashayek et al. 2013, 2015). Here, we made use of the Simmons et al. (2004) internal tide mixing scheme, which prescribes an exponentially decaying structure function for the diffusivities, so that diffusivities decay away from the bottom.
However, it is also feasible that our diagnosed transient responses are only influenced in a minor way by the model mixing parameterization. Dutay et al. (2010) suggest that increasing κ from 0.2 to 1.2 × 10−4 m2 s−1 at depth (in an ocean model with realistic bathymetry) increases the subpolar cell response to geothermal heating from 50% (with a constant κ) to 68%. Here, the abyssal overturning cells in the CTL2 case are markedly larger (by at least 50%) than those of the CTL1 case (Fig. 5a), and we find a larger response in the abyssal overturning cells associated with geothermal heating (Figs. 5b,c). In effect, comparing the responses in the GH1 and GH2 cases provides a test of how the ocean responds to geothermal heating with a warmer ocean (with a smaller vertical buoyancy gradient). In the CTL2(701–1300) run, the global averaged vertical buoyancy gradient for the midocean ridge depths (2000–4500 m) are larger, associated with stronger overturning circulation for the less dense half of the lower and subpolar cells. In short, while we did not explicitly test how the influence of geothermal heating is affected by a higher vertical diffusivity, we can confirm that the overturning circulation response to geothermal heating is larger for a background ocean where the vertical buoyancy gradient is smaller.