The effect of a localized geothermal heat source on deep water formation

In a simplified two-dimensional model of a buoyancy-driven overturning circulation, we numerically study the response of the flow to a small localized heat source at the bottom. The flow is driven by differential thermal forcing applied along the top surface boundary. We evaluate the steady state solutions versus the temperature difference between the two ends of the water surface in terms of different characteristic parameters that properly describe the transition from a weak upper-layer convection state to a robust fulldepth deep convection. We conclude that a small additional bottom heat flux underneath the “cold” end of the basin is able to initiate full-depth convection even when the surface heat forcing alone is not sufficient to maintain this state.


Introduction
The engine of the Great Ocean Conveyor is the sinking of cold and saline surface waters to the bottom of the oceanic basin, referred to as Deep Water Formation (DWF).Under the present climatic conditions, DWF can occur only at a few compact polar downwelling regions, all of which are located in the Atlantic basin (van Aken, 2007).This polar sinking together with the distributed upwelling of water masses drives the Atlantic Meridional Overturning Circulation (AMOC) (Stommel, 1958), which is an important contributor to the European climate, because of its large northward heat transport of an order of magnitude of 10 15 W (Ganachaud and Wunsch, 2000).
The question arises of what makes the Atlantic the only basin where DWF can occur among the present climatic conditions.The average surface salinity of the Atlantic is significantly higher than that of all other basins, due to its larger Correspondence to: M. Vincze (vincze.m@lecso.elte.hu)evaporation-precipitation difference.The combined effect of high salinity and the heat release, which occurs in the polar regions results in an unstable stratification that makes the dense surface water parcels descend to the bottom.However, it has been shown that the amount of precipitation is quite similar over the Atlantic and the Pacific, and the evaporation difference is due to the temperature difference between the basins, which is largely a consequence of the ocean circulation itself (Huisman et al., 2009).This implies that the correct explanation of the Atlantic DWF must involve additional effects besides the widely studied role of surface buoyancy (i.e.heat or fresh water) fluxes.
Furthermore, Sandström's thermodynamic theorem (Sandström, 1908) states that a density-driven circulation can only be maintained if the spatial position of the positive heat source is below the level of the negative heat source (cooling).In the ocean this is seemingly not the case, as most of the warming and cooling occurs at the water surface.The weak geothermal heating of the seafloor is usually neglected in general ocean circulation models, as zero heat flux boundary conditions are being applied (e.g., Frankcombe et al., 2009).To fulfil Sandström's theorem, these models use large, highly anisotropic eddy diffusivity and viscosity values to parameterize turbulent mixing, which is believed to be responsible for the required upward heat and momentum "pumping".In some numerical studies (e.g., Urakawa and Hasumi, 2009), the effect of geothermal heat on a density-driven circulation was also taken into account, by applying a uniform constant heat flux at the bottom of the basins.According to in-situ measurements, however, the heat flux distribution of the oceanic crust is far from homogeneous.In Scott et al. (2001), the authors investigated the effect of an enhanced heat flow that is present in the vicinity of the mid-ocean ridges in a highly idealized ocean model, and concluded that even such an inhomogeneous thermal perturbation can cause anomalous response in the AMOC.
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
Surprisingly enough, the locations of the two most important DWF regions, namely at the Greenland Sea, and the Weddell Sea both exhibit extremely high values of geothermal heat flux (∼ 120 mW m −2 ), twice as high as the global average (Shapiro and Ritzwoller, 2004).To our best knowledge, this coincidence remained unnoticed in the literature.The scope of the present work is to investigate the response of a density-driven two-dimensional circulation to a compact, localized bottom heat source, placed underneath the "polar" end of the basin, where the surface heat flow is negative (cooling occurs).
Although the concept originates from the above discussed oceanographic problem, the parameters used here are rather unrealistic, as we intended to model a laboratory-scale setup for a better understanding of the basic underlying physics.Furthermore, we think that such parametrization supports a direct comparison with actual experiments (we compare some aspects of the ocean circulation and the numerical and experimental setups in the Appendices).Our main result is that even a weak local heat source beneath the cold end of the basin perturbs significantly the initially shallow-layer horizontal convection and markedly contributes to the DWF process.

The model setup and numerical methods
The present model is based on the two-dimensional nonhydrostatic Boussinesq equations.These were solved on a 20 × 200 equidistant array of Arakawa C cells (Arakawa and Lamb, 1977) which corresponds to a 2 m deep and 20 m long rectangular tank.The effect of salinity differences was not considered, hence the flow was only forced by the incoming heat fluxes from the surface, the bottom and the lateral sidewalls (a qualitative argument on the possible effect of salinity on the phenomena described here, is presented in Appendix B).The density ρ(T ) was assumed to obey a linear temperature dependence with reference values T ref = 25 • C, ρ ref = 997.075kg m −3 and volumetric thermal expansion coefficient α = 2.5 × 10 −4 K −1 .The kinematic viscosity ν = 10 −6 m 2 s −1 and the thermal diffusivity κ = 1.4 • 10 −7 m 2 s −1 were treated as isotropic constants, of their usual molecular values, resulting in a Prandtl number of P r = 7.14.We note, that the three orders of magnitude larger eddy viscosity and diffusivity values that are set in ocean models (e.g.ν eddy ∼ 10 −3 m 2 s −1 and κ eddy ∼ 10 −4 m 2 s −1 , as in Frankcombe et al. (2009)) yield approximately the same dimensionless ratio, therefore this part of the dynamics remains unchanged by the upscaling.
Numerical solutions were obtained using the Advanced Ocean Modeling open-source software package, written in FORTRAN95 environment (Kämpf, 2009), in which the system of PDEs are being solved with the Successive Over-Relaxation (SOR) method.According to the initial conditions, the fluid was at rest and its temperature was uniformly set to the value of T ref in the whole solution domain.At the water surface and the sidewalls slip boundary conditions were applied for the velocities, while friction was taken into account at the bottom in the form of where u denotes the horizontal velocity, with a value of u b adjacent to the bottom boundary, and the bottom-drag coefficient r was set to 0.001 as by Kämpf (2009).The flows in the domain were forced by restoring temperature boundary conditions at all boundaries, represented by a distributed heat source that is proportional to the temperature difference between the local water temperature and a spatially varying prescribed temperature field T relax (x,z), as given by: Generally, the restoring timescale τ was set to 2100 s at all boundaries, based on our laboratory measurements (details will be reported elsewhere).We note here, that for an extension of a similar model setup to oceanic scales one cannot apply the same restoring timescales at the water surface and at the other boundaries.In the case of a real ocean the heat exchange between the uppermost layer and the atmosphere is orders of magnitude more effective (τ ∼ 30 days, as in Frankcombe et al. (2009)) than it is at the abyssal regions, due to the wind-driven turbulent mixing.However, for a laboratory-scale setup that is studied here the usage of the same restoring timescales at all boundaries is a reasonable assumption.Further comparison of the thermal boundary conditions of the setup and the ocean can be found in Appendix A.
For the restoring temperature, a value of T relax (x,z) = T ref = 25 • C was defined at the vertical sidewalls and at the bottom.At the water surface, T relax (x,z) = T warm = 32 • C was prescribed in the leftmost 20 cells.A similar domain of 20 parcels was defined at the right margin of the surface as well, in which T relax (x,z) = T cold was adjusted in the range − max (white) and max (black), for actual magnitudes see Fig. 3c and f.Correspondingly, dark gray regions reflect positive (clockwise) local vorticity.In the uppermost two snapshots, labels "W" and "C" denote the warming and cooling sides.Label "HS" marks the position of the bottom "hot spot".Figures 2a,b 2f), that transfers momentum to the neighboring anticyclonic downwelling.In Fig. 2g the full-depth DWF state is present, and the counterclockwise cell already vanished.Note, that in the "hot spot" case, the marked one-cell convection is visible at T cold = 18.9 • C, unlike in the "no hot spot" case.The vertical lines in each snapshot represent the left boundary of the region over which the mean value w polar (see Fig. 3b and e) is evaluated.
of T cold = 4-24 • C for the different numerical experiments.Between these sections -that model the equatorial and polar regions of the ocean surface -the value of the restoring temperature was interpolated linearly between T warm and T cold at the water surface, as depicted in Fig. 1.
Besides T cold , our other key control parameter was the restoring temperature T relax (x,z) = T spot of a "hot spot", i.e. a 10 parcel-long section at the bottom right corner of the basin (see Fig. 1).This domain represents the aforementioned compact Atlantic seafloor regions of high geothermal flux at higher latitudes.T spot varied in an interval between 25 • C (= T ref ) and 35 • C for the different runs.
The spin-up time of the model was taken to be 12 000 s, after which a quasi-stationarity of the average surface temperature time series was achieved in every run.In order to properly describe the transition between the different con-vection states, we introduced some characteristic parameters, which exhibit jumps at the transitions between different dynamical regimes.These parameters were evaluated for the quasi-stationary part of the process.

Results
In our first series of experiments, the value of T spot was set to T ref , corresponding to the "no hot spot" scenario.In this case the flow pattern is determined by one single control parameter T cold , since T warm = 32 • C was fixed throughout all the computations.The "phase transition" from a weak upper-layer convection to a robust full-depth deep convection (i.e.DWF) is clearly visible in the left column of Fig. 2 f.The dashed and dashed-dotted vertical lines mark the critical T cold , where the onset of DWF takes place in the "no hot spot" and in the "hot spot" case, respectively.
13 Fig. 3.The T cold dependence of three parameters that characterize the transition between the convection states.As in Fig. 2, the left column corresponds to the "no hot spot" scenario, and the right column shows the T spot = 30 • C case.Figures 3a and d represent the quasi-equilibrium temperatures at the right (cold) margin of the basin at the surface and at the bottom.Figures 3b and e depict the vertical velocity ( w polar ) averaged over the region that is indicated in Fig. 2. The maxima of the streamfunctions ( max ) are shown in Figs.3c and f.The dashed and dashed-dotted vertical lines mark the critical T cold , where the onset of DWF takes place in the "no hot spot" and in the "hot spot" case, respectively.
Naturally, the onset of DWF is related to the presence of unstable density stratification at the "polar" region of the tank.Firstly, we determined the time averaged temperatures T top asympt and T bottom asympt , measured at the rightmost upper and bottom corners of the basin (solid and dashed lines with symbols in Fig. 3a).It is clearly visible in the higher T cold -range that the asymptotic bottom temperature practically reaches the prescribed value of T ref = 25 • C.This implies that the bottom layer of the fluid stays at rest in this regime.The transition -where T asympt reaches the same value at the surface and at the bottom -occurs at T * cold = 18.3 • C (denoted with dotted vertical line in Fig. 3a, b and c).
Next, we measured the mean vertical sinking velocity w polar , that is averaged in the 5 m long "polar" section identified by the vertical lines in Fig. 2. (Negative values of w polar represent downward flow.)The onset of DWF is clearly indicated as a break point at T * cold = 18.3 • C (Fig. 3b).At temperature values T cold < T * cold , the mean vertical velocities w polar have definite negative values.
Thirdly, we computed the maximal value of the timeaveraged streamfunction ( max ) for each run.The curve in Fig. 3c demonstrates the same transition: in the range T cold < T * cold , the large values of max are related to vigorous deep-water convection.
For the following numerical experiments we introduced a restoring temperature T spot > T ref at the bottom hot spot, in order to evaluate whether this extra heat flux affects the state transition in terms of the critical T cold .The right column of Fig. 3 depicts the behavior of the above mentioned order parameters in the case of T spot = 30 • C. The corresponding streamfunction patterns for the same series of experiments are shown in Fig. 2e-h    the presence of the additional "geothermal heat", the critical value of T * cold has shifted considerably upward, to the value of 20.2 • C (dashed-dotted in the right column of Fig. 3).
The next question of importance is whether there exists a certain T spot for a given T cold that maximizes the flux of the polar downwelling.Intuitively, one can argue that over a certain temperature threshold the bottom hot spot launches a rising plume -as discussed in e.g.: Stommel (1982) that would tend to reverse the direction of the circulation, thus hindering DWF rather than enhancing it.Fig. 4a shows w polar as a function of T cold for T spot = 25.0,27.2, 30.0, and 32.5 • C, the first and third curve being the same as in Figs.3b and e.There are two main observations to be emphasized: firstly, the critical crossover temperature T * cold depends weakly on the restoring hot spot temperature T spot ; secondly, the largest downward flux values in the range T cold < T * cold do not belong to the highest value of T spot .According to the expectations, high enough T spot temperatures seem to hinder DWF, see the dotted line in Fig. 4a, in the range of T cold < 15 • C (for T spot = 32.5 • C).
In order to move along an orthogonal axis of the parameter space, we evaluated the same quantity w polar as a function of the bottom heating in a broader T spot range for three fixed T cold values (Fig. 4b).When the driving horizontal temperature difference T warm − T cold is too small, DWF is hindered even when some hot spot is present (see the topmost curve of Fig. 4b, T warm = 32.00• C as in each case, T cold = 20.78• C).At large enough horizontal temperature gradients, a localized bottom heat source enhances downward DWF fluxes, exhibiting an optimal value at T opt spot = 28 • C for T cold = 19.73• C (middle curve in Fig. 4b), and T opt spot = 30 • C for T cold = 18.90 • C (lowermost curve in Fig. 4b).Both series of experiments point out that -in agreement with the qualitative reasoning -there exists a small value of localized heat flux for every given "equator-to-pole" temperature difference that is capable to maximize the flux rate of DWF.

Conclusions
We performed numerical experiments in a simplified, twodimensional laboratory-scale setup in order to capture the basic effects of a localized bottom heat source on Deep Water Formation in a convective system driven basically by top heat fluxes.As far as we know, this is the first study in a similar arrangement.Previous studies, such as Mullarney et al. (2006) incorporated uniform bottom heating and pointed out strong perturbations in the convection, however they did not study an isolated heat source.We hope that our choice of parameters promotes laboratory-scale control experiments.
Numerical simulations performed in oceanic-scale threedimensional setups (Scott et al., 2001) demonstrated that a realistically small uniform bottom heat flux (∼ 50mW m −2 ) can have a significant effect on DWF.As a consequence of the extremely weak stratification of abyssal ocean waters, even this flux can lead to an extra temperature difference of T ∼ 0.5 • C between the bottom and the surface (Scott et al., 2001).Our results support the idea that a localized hot spot with relatively weak extra heat flux is able to initiate DWF under such conditions when the surface heat forcing alone is not sufficient to maintain the deep-convection state.Therefore we believe that taking bottom heat sources into account in ocean models is clearly not an unrealistic idea.
Our results might imply that in the present climate, the equator-to-pole temperature difference is not high enough for www.nonlin-processes-geophys.net/18/841/2011/ Nonlin.Processes Geophys., 18, 841-847, 2011 the onset of a distributed polar DWF.It might happen that the current temperature gradient is sufficient to maintain only localized DWF regions in the vicinity of larger than average geothermal heat sources, where the bottom heating can enhance downwelling.This conclusion might provide a new argument in understanding the sensitivity of AMOC to the climatic conditions that has been observed in paleoclimatic data (Thornalley et al., 2011), and might explain the lack of DWF regions in the Pacific.

Heat fluxes
The incoming solar radiation plays an essential role in creating the meridional temperature and salinity differences that drive the oceanic overturning circulation.The annual average of the net irradiance received on a given surface area raises approximately from 50 W m −2 to 300 W m −2 , mainly depending on the latitude (Shapiro and Ritzwoller, 2004).
Compared to this, the geothermal heatflow at the seafloor is roughly three orders of magnitude smaller, with an average value of 50 mW m −2 .For the regions that exhibit elevated geothermal flux (e.g. the vicinity of hot spots, or mid-ocean ridges) this value can be as high as 120 mW m −2 .The locations of our particular interest, where Deep Water Formation actually occurs, are such that -being in subpolar regionsthey exhibit relatively low average irradiance at the surface and higher-than-average geothermal heating at the seafloor.We compare these realistic heatflow values to those that are present in our experimental setup.Since restoring boundary conditions (2) have been applied for the temperature in our study, once a quasi-equilibrium state is reached, the heatflow at the boundaries (in W m −2 units) can be approximated as follows, see e.g., (Frankcombe et al., 2009) or (Kämpf, 2009): where ρ 0 = 1000 kg/m 3 is the reference density, C p = 4000 J kg −1 K −1 is the heat capacity, and τ = 2100 s is the characteristic restoring timescale, as estimated by measuring the surface temperature response of an actual laboratoryscale experimental setup, after a jumpwise change in the surface heat forcing.This also sets the thermal diffusional lengthscale to δz = √ κτ ∼ 0.01 m, as taken with the molecular thermal diffusivity κ = 1.4 • 10 −7 m 2 s −1 .T relax (x,z) denotes the prescribed values of restoring temperatures at the boundaries, and T asympt (x,z) stands for the actual quasistationary temperature value which the system reaches following a transient phase.
If no differential heating, diffusion or advection took place in the setup, the restoring boundary condition alone would drive the system towards an equilibrium state of T ≡ T relax , i.e.Q = 0 (no-heatflow state).However, because of the dynamics that is present here, the actual equilibrium state exhibits non-zero fluxes at all boundaries.These can be evaluated using Eq.(A1), by detecting the values of T asympt (x,z) at grid locations in the vicinity of the different boundaries.
Substituting the values of the equilibrium temperatures of our setup, the surface heat flow is found to be on the order of 100 W m −2 for all the runs (the precise values, of course, depend on the horizontal position, and on the actual T cold and T spot boundary conditions for the given run).The magnitude of the bottom heat flow, as averaged over the whole basin length, varied between 0.01 W m −2 and 0.1 W m −2 .If averaged only over the vicinity of the "hot spot", the heat fluxes raised from 0.01 W m −2 to 10 W m −2 , depending on the value of the restoring temperature T spot .
This implies, that -as for the heat fluxes at the boundaries -the actual values of the real ocean lie within the range studied in our numerical experiments.For the better understanding of the basic dynamics though, we investigated a widerthan-realistic parameter range.
It is to be emphasized again, that this setup is far not a realistic ocean circulation model, rather it may be thought of as a "toy model" that is meant to drive the attention to a previously uncovered aspect of the dynamics in a densitydriven circulation.This effect (namely, the enhancement of downwelling by a weak localized bottom heat flow) certainly exists, but to reveal the importance of its contribution to the actual oceanic Deep Water Formation would definitely require more advanced simulations.

Salinity effects on oceanic scale
In linear approximation, the density ρ of a water parcel is determined by its temperature T and salinity S as: where ρ 0 , T 0 , and S 0 denote the reference values of density, temperature and salinity, respectively, α marks the thermal expansion coefficient and β represents the haline contraction coefficient.In our numerical setup the salinity term was neglected in order to reduce the investigated parameter space.In the case of the real ocean, however, the meridional gradient of surface salinity is an important driving force of the overturning circulation.This gradient arises because of the differential evaporation/precipitation over the basin.Intense evaporation increases the salinity of a water parcel at the Equatorial regions.Once this parcel reaches the subpolar Deep Water Formation region, this elevated salinitytogether with cooling -helps to build up a vertical density instability, that eventually forces the parcel to descend.So, taking salinity differences into account would, in general, enhance Deep Water Formation.
On the other hand, a stable vertical salinity gradient above the hot spot could theoratically counteract the destabilizing effect of this abyssal heat source, and thus suppress downwelling.However, according to field data (van Aken, 2007), the typical oceanic salinity profiles are such, that marked gradients are present only at the uppermost ∼ 100 m thick mixing layer, while in the deep ocean the salinity is basically homogeneous, therefore the buoyancy differences in the vicinity of the seafloor are determined dominantly by the temperature field.
This means, that adding realistic surface evaporationprecipitation differences and realistic density profiles to the model, we would expect the observed phenomena to be enhanced, rather than suppressed.We note however, that the evaporation-precipitation dynamics are not expected to be successfully resolved by any real laboratory-scale experiment, as the timescale of the evaporation cannot be scaled down to be comparable to the characteristic timescale set by the thermally driven part of the circulation, as in the case of the real ocean.

Fig. 1 .
Fig. 1.The schematic drawing of the setup, with the domains where the different values of T relax were applied along the boundaries.

Fig. 2 .
Fig. 2. Averaged streamfunction ( ) patterns for the quasi-stationary parts of different runs.The gray-scale values cover the interval between ,c and d correspond to the "no hot spot" scenario, i.e.T spot = T ref = 25 • C. The transition from a multi-cellular weak convection state to a DWF state is clearly visible as a function of T cold .In Fig. 2e,f,g and h the same transition is shown in the case of T spot = 30 • C. In Fig. 2e the presence of the hot spot initiates a marked two-cell convection along the whole basin, which -at T cold = 19.73• C -develops to a full-depth counterclockwise Benard cell above the hot spot (Fig.

Fig. 3 .
Fig. 3.The T cold dependence of three parameters that characterize the transition between the convection states.As in Fig.2, the left column corresponds to the "no hot spot" scenario, and the right column shows the Tspot = 30 • C case.Figures3a and drepresent the quasi-equilibrium temperatures at the right (cold) margin of the basin at the surface and at the bottom.Figures3b and edepict the vertical velocity ( w polar ) averaged over the region that is indicated in Fig.2.The maxima of the streamfunctions (Ψmax) are shown in Figures3c and

Fig. 4 .
Fig. 4. (a) w polar as a function of T cold , for four different values of Tspot, see legends.(b) Dependence of the same quantity w polar on Tspot, for three fixed values of T cold , see legends. 14

Fig. 4 .
Fig. 4. (a) w polar as a function of T cold , for four different values of T spot , see legends.(b) Dependence of the same quantity w polar on T spot , for three fixed values of T cold , see legends.