A conceptual model of oceanic heat transport in the Snowball Earth scenario

Geologic evidence suggests that the Earth may have been completely covered in ice in the distant past, a state known as Snowball Earth. This is still the subject of controversy, and has been the focus of modeling work from low-dimensional models up to state-of-the-art general circulation models. In our present global climate, the ocean plays a large role in redistributing heat from the equatorial regions to high latitudes, and as an important part of the global heat budget, its role in the initiation a Snowball Earth, and the subsequent climate, is of great interest. To better understand the role of oceanic heat transport in the initiation of Snowball Earth, and the resulting global ice covered climate state, the goal of this inquiry is twofold: we wish to propose the least complex model that can capture the Snowball Earth scenario as well as the present-day climate with partial ice cover, and we want to determine the relative importance of oceanic heat transport. To do this, we develop a simple model, incorporating thermohaline dynamics from traditional box ocean models, a radiative balance from energy balance models, and the more contemporary “sea glacier” model to account for viscous flow effects of extremely thick sea ice. The resulting model, consisting of dynamic ocean and ice components, is able to reproduce both Snowball Earth and present-day conditions through reasonable changes in forcing parameters. We find that including or neglecting oceanic heat transport may lead to vastly different global climate states, and also that the parameterization of under-ice heat transfer in the ice–ocean coupling plays a key role in the resulting global climate state, demonstrating the regulatory effect of dynamic ocean heat transport.


Introduction
It is well known that the albedo difference between seawater and sea ice leads to a crucial climatic feedback.In a warming climate, melting of high-albedo sea ice exposes low-albedo seawater, thus increasing the fraction of incoming solar radiation that is absorbed and thereby amplifying warming.Conversely, in a cooling climate the growth of sea ice increases the planetary albedo, which amplifies cooling.Classic energy balance models (EBMs) demonstrate how this well-known ice-albedo feedback can lead to multiple steady climate states (Budyko, 1969;Sellers, 1969).Given forc-ings which resemble present-day conditions, these models are bistable, with one possible steady state having a partial ice cover and another being completely ice-free.With a sufficient reduction in the solar input or the greenhouse effect these energy balance models yield completely ice-covered steady states, reminiscent of the "Snowball Earth" episodes of the Neoproterozoic era.
The virtue of simple models, of course, is that they make it possible to explore ranges of relevant parameters easily.First-order energy balances at climatic scales are overwhelmingly radiative.The first energy balance models were not intended to model meridional energy transport, and if included, it was incorporated as a diffusive process.Diffusiondominated transport tends to mitigate the tendency for sea ice to grow in a cooling climate: formation of extra sea ice would increase the temperature contrast between the ice-covered high latitudes and the low latitudes, which in turn would increase the rate of heat transport to the high latitudes.This then slows the growth of the ice and exports some of the excess cooling due to the ice-albedo feedback to lower latitudes.
While the small ice cap instability has been a common phenomenon in EBMs (Held and Suarez, 1974), this bistability has not been seen in more sophisticated general circulation models (GCMs) (Armour et al., 2011).In a recent study, Wagner and Eisenman (2015) showed the bistability does not appear when an EBM was combined with a single column model for ice thermodynamics, suggesting sufficient complexity, including a seasonal cycle and diffusive heat transport, eliminates this bifurcation.While the small ice cap instability is not the primary focus of this study, we will use our model to investigate this small ice cap bifurcation.
Another class of simple models, box ocean models, have been critical to understanding meridional ocean circulation, upwelling and mixing, among other processes (Stommel, 1961).In particular, they have suggested important climate questions to pursue by means of more sophisticated models, data, or both.Models and measurements established; however, the importance of the basic thermohaline circulation and the more complete meridional ocean circulation.With it, models for energy transport via diffusive processes were replaced by a more heterogeneous and dynamic mechanism.The Atlantic Meridional Overturning Circulation (AMOC), in particular, contributes a net energy transfer into the North Atlantic equivalent to several percent of the total incoming shortwave solar radiation incident to the region.Models of the meridional overturning circulation, from the simplest to the most detailed, agree that under present-day conditions the circulation is bistable (Rahmstorf, 2000;Rahmstorf et al., 2005).In addition to the thermally dominated steady state that is currently observed, a second, salinity-dominated steady state is also possible, with a much weaker circulation which flows in the direction opposite to the present flow.It is clearly important to understand the interaction between this oceanic circulation and the distribution of sea ice.
Since the present-day AMOC transports heat into the North Atlantic, it tends to reduce the extent of sea ice.A strengthening of the circulation would then reduce ice cover, while a weakening would cause it to expand.The circulation itself is driven by wind stress in and near the Southern Ocean, as well as in part by mechanical mixing from tidal action (Munk and Wunsch, 1998), and it is sustained by variations in the density of circulating water as it exchanges heat and fresh water with the atmosphere as it flows along the surface.As water flows northward along the surface through the trop-ics it is warmed, and its salinity increases as a result of excess evaporation.The increase in temperature decreases the water density, while the increase in salinity increases it.As the water passes to higher latitudes, it cools and freshens as precipitation exceeds evaporation.Again, the two effects tend to change the density in opposite directions.In the present configuration of the AMOC, the thermal effects dominate, so the water becomes denser as it moves through the subpolar latitudes and ultimately sinks, returning southward at depth.
The presence of sea ice affects the processes that change the density of circulating seawater at high latitudes.An ice cover isolates the water from the atmosphere and thus cuts off the precipitation that otherwise would reduce the salinity of the water and lower the rate at which its density increases.It also insulates the water thermally from the atmosphere.This by itself would not have much effect on the water density, since the water temperature cannot fall below the freezing temperature anyway.However, it allows the atmosphere to become colder than it would be if it were in contact with the seawater.Heat transfer to the cold atmosphere through the ice layer results in freezing of seawater at the base of the ice.Brine rejection then increases the salinity and hence the density of the remaining water.
A simple model of ocean circulation with sea ice may be useful in studying the Snowball Earth episodes in the Neoproterozoic era.Since the discovery of geologic evidence suggesting that glaciation occurred in the tropics at least twice in the Neoproterozoic (Kirschvink, 1992;Hoffman and Schrag, 2002), some 710 and 635 million years ago (Ma), there has been substantial debate about whether the Earth was ever in a completely ice-covered Snowball Earth state (Lubick, 2002).These events are thought to have lasted several million years, raising such questions as how life could have survived a long period if the Earth were in a completely ice-covered state (McKay, 2000).Thus, in trying to explain these signs of apparent tropical glaciation in the context of global climate dynamics, alternative hypotheses have been proposed that leave some portion of the ocean either free of ice or covered only in thin ice (Hyde et al., 2000;Pollard and Kasting, 2005;Abbot and Pierrehumbert, 2010;Abbot et al., 2011).Poulsen et al. (2001) counts among the first studies to suggest dynamic ocean heat transport is important to the Snowball Earth hypothesis.Using the fully coupled Fast Ocean-Atmosphere Model (FOAM), initialized to Neoproterozoic parameter values to facilitate snowball conditions, the authors found that a global ice cover was produced when using a mixed-layer ocean model that parameterized heat transport through diffusion.In fully coupled experiments with the ocean component, on the other hand, the ice margin would retreat to high latitudes.Other studies have considered the Snowball Earth problem in an EBM framework, with ocean heat transport typically parameterized by a diffusion process (Pollard and Kasting, 2005;Rose and Marshall, 2009).In particular, Pollard and Kasting (2005) examine the feasi-bility of a tropical thin-ice solution, incorporating detailed treatment of optical properties of ice and a nonlinear internal ice temperature profile, as well as a separate snow layer and an evaporation minus precipitation term to facilitate surface melt/accumulation.
GCMs that have more detailed ocean physics have also been used to study the initiation of a Snowball Earth (Yang et al., 2012a, b).A recent study by Voigt et al. (2011) uses the state-of-the-art atmosphere-ocean model ECHAM5/MPI-OM to study the Snowball Earth scenario.They implement a Marinoan (635 Ma) land mask in their coupled GCM simulations, as well as the lower insolation of a younger, weaker sun.In addition to ocean dynamics, their study also included sea ice dynamics (albeit with thin ice) and interactive clouds.All three of these had previously been found to be essential for snowball initiation (Poulsen et al., 2001;Poulsen and Jacob, 2004;Lewis et al., 2003Lewis et al., , 2007)).Voigt et al. (2011) were able to achieve snowball initiation, and also to prevent snowball initiation in the same setting by doubling carbon dioxide levels.Stability analysis of an EBM analog, based on the 0-D model of global mean ocean temperature developed in Voigt and Marotzke (2010), indicates an insolation bifurcation point for Snowball Earth in the Marinoan setting of about 95-96 % of pre-industrial levels, in agreement with their computational results.In their experiments that resulted in partial ice cover, the ice margin was around 30 to 40 • latitude, with maximum stable sea ice extent of 55 % of ocean cover observed in their experiments.
Sea ice in a global ice cover can be very thick, to the extent that flow by plastic deformation under its own weight should be considered.Thus, its non-Newtonian fluid dynamics must be considered in addition to its thermodynamics.Goodman and Pierrehumbert (2003) (henceforth GP03) first considered these flow effects in the Snowball Earth scenario.(The same framework was used in Abbot and Pierrehumbert (2010) and Li and Pierrehumbert (2011) to transport dust to low latitudes in the mudball scenario).Their model runs outside a global circulation model, using FOAM output for forcing data, and it has neither an active ocean component nor a parameterization for oceanic heat transport.They use the term "sea glacier" to describe their modeled ice, to distinguish it from present-day sea ice, which only grows to thicknesses on the order of meters, and from land ice or ice shelves.The sea glacier is formed in the ocean, yet it achieves the thickness of a land ice sheet without the land-ice interface, and its non-Newtonian rheology is taken into account in the calculation of its flow.They are able to achieve both partial glaciation and a full Snowball Earth state through changes in the atmospheric forcing (surface temperature and precipitation minus evaporation).They find that the additional viscous flow term is highly effective at allowing the ice margin to penetrate low-latitude regions of melting, thus encouraging Snowball Earth initiation.
A recent study by Ashkenazy et al. (2013) found a dynamic ocean in a Snowball Earth scenario with strong cir-culation, in contrast to a stagnant ocean typically expected due to ice cover serving as an insulation layer to atmospheric forcing.Their model was forced with geothermal heat, which was spatially varying with a peak near the Equator, averaging to 0.1 W m −2 .In their 2-D and 3-D ocean simulations coupled with a 1-D ice model extending upon that of GP03, they found that the ocean plays a larger role in determining ice thickness than the atmosphere, and that geothermal heat forcing plays a dominant role in ice-covered ocean dynamics.This was expanded upon in Ashkenazy et al. (2014), where a dynamic ocean with strong equatorial jets and a strong overturning circulation was found in simulations of a steady-state globally glaciated Earth.
Atmospheric dynamics and cloud cover undoubtedly play a large role in such climate systems, as demonstrated by Voigt et al. (2011).It is, however, difficult to isolate the role played by oceanic transport in these coupled simulations, due to the necessary inclusion of the complex dynamics of the atmosphere and cloud distribution as well as sea ice dynamics.In fact, in Voigt and Abbot (2012) it was found that ocean heat transport has no effect on the critical sea ice cover that leads to snowball initiation.This motivated us to consider a simpler model that includes oceanic heat transport coupled to ice dynamics.We aim to extend the framework laid out in GP03 to include ocean heat transport effects, including under the ice layer.Realistic oceanic transport undoubtedly leads to highly nonuniform heat distributions, likely with local consequences on the global snowball scenario.However, these local effects are beyond the scope of our study, and indeed beyond the scope of any low-dimensional model.By omitting atmospheric effects, we aim to get an assessment of the effects from oceanic heat transport alone.
The purpose of this paper is to investigate the interaction between sea ice and the meridional overturning circulation.By now there are several studies that have used complex circulation models to confirm that ocean transport is an important component of any explanation of how sea ice recedes and grows along with changes in forcing, albedo, biogeochemistry, etc.With a simple model, however, it is possible to efficiently test our understanding and propose questions critical to our being able to further understand the complexities of radiation, ice cover, and oceanic/atmospheric transport, such as the Snowball Earth hypothesis.
To this end we have combined a one-dimensional energy balance model with a box model of the meridional overturning circulation and a dynamic ice component.Our model is described in detail in Sect.2, a key element of which is the under-ice heat exchange with the ocean.In Sect. 3 we present model results in different climate regimes.The sensitivity of the model to key parameterizations is studied in Sect.4, and concluding remarks are in Sect. 5.
Our model consists of a four-box ocean model with transport, similar to that first proposed by Stommel (1961), coupled to a one-dimensional EBM similar to that of Budyko (1969) and a dynamic model for sea ice coverage (GP03), and is depicted in Fig. 1.The ocean component is a hemispheric model with thermohaline dynamics.While the ocean model uses a traditional transport equation for the salinity, it differs from traditional box models in the use of an energy conservation model to capture the temperature dynamics.This allows us to couple the ice and radiation components to the ocean dynamics.The ice layer is zonally averaged, so its thickness is taken to depend only on latitude θ and time t.The ice margin evolves dynamically, and we include non-Newtonian flow so that the model can accommodate an ice layer thick enough to be appropriate to Snowball Earth conditions.The surface absorbs incoming solar radiation, with an albedo which takes into account whether the surface is open ocean or ice, and emits longwave radiation with a specified emissivity.Where ice overlies the ocean, heat conducts through the ice layer and is exchanged with the ocean at the ocean-ice interface, where melting and freezing can occur, while also supplying heat for redistribution through ocean circulation, a departure from typical EBMs.We also account for geothermal heat forcing, found to be the dominant forcing in a Snowball Earth ocean (Ashkenazy et al., 2013).
In addition to physical properties of ice and seawater, geometrical features of the ocean basin, and the distribution of insolation, there are three parameters which are important to our studies.In the energy balance model we model the greenhouse effective by including an effective emissivity, ε, in terms of modeling outgoing longwave radiation.The box model requires a hydraulic coefficient, k, which relates the strength of circulation to the densities of the water in the various boxes.The third parameter quantifies the thermal coupling between the sea ice layer and the water beneath; we express this as an effective thermal boundary layer thickness, D, with the rate of heat transfer from the water to the ice being proportional to the temperature difference between the box water and the base of the ice, divided by the thickness D. Among the issues we will investigate is the question of how, and indeed whether, the coupling relates the bistability of the energy balance model and the independent bistability of the box model.

Ocean and energy components
The ocean component of our model is similar to the one proposed by Griffies and Tziperman (1995) and Kurtze et al. (2010).Four boxes are used to represent the ocean in one hemisphere, from pole to Equator, with each box representing a zonal average across longitude.Referring to Fig. 1, we define Box 1 as the tropical surface ocean box, Box 2 its polar counterpart, Box 3 below Box 2, and Box 4 below Ice cover

Colatitude θ
Ice margin = η(t) • simply because we are interested in investigating climate regimes ranging from global ice cover to zero ice cover, so there is no advantage in trying to confine the ice cover to a polar box only.We have found that changing the location of this boundary within midlatitudes does not qualitatively affect our results.The dynamic ice margin is at η(t), the ice cover thickness is given by h(θ, t), and its poleward meridional velocity is given by v(θ, t).Model parameter values, including geometric quantities pertaining to the box structure, are given in Table 2.For clarity, model variables will be subscripted with letters rather than numbers, using ut = upper (surface)  Each box has a (well-mixed) temperature T j (t) and salinity S j (t), where j = ut, up, lp, lt, which determine the density of each box by a linear equation of state: Here ρ 0 is a reference density corresponding to a reference temperature and salinity T 0 , S 0 .β S and β T are the expansion coefficients associated with salinity and temperature.The density-driven flow between the boxes is denoted by f , where we adopt the convention that f < 0 is surface poleward flow (from Box 1 to Box 2).As in Griffies and Tziper-man (1995), the (buoyancy-driven) transport rate is where k = k 0 = 8 × 10 4 /ρ 0 Sv is the hydraulic constant which governs the strength of the density-driven flow.In Sect.4.2 we explore the model's sensitivity to this parameter.The flux f is purely thermohaline-driven.We could modify this flux to include the effect of wind stresses in a crude manner by an additive correction to the flux, however, this has been omitted in this study.
The equations for each box's salinity and temperature will depend on the direction of the mean meridional flow.The salinity equations when f < 0, corresponding to poleward Here V j is the volume of the j box and M(θ, t) is the total production/melting rate of ice, with M > 0 corresponding to ice production and M < 0 corresponding to melting, described in Sect.2.2.The terms involving the circulation rate f correspond to fluxes across the box boundaries.We assume that ice sits atop the surface ocean boxes, and that the mass of ice is much less than the total mass of ocean water.The surface ocean box volumes V j are kept constant, and any changes to the deep ocean box volumes are negligible.An important assumption is that ice that is formed is freshwater ice, which rejects brine into the ocean.The integral term represents the change in salinity due to net freshwater added/removed through ice melting/production.The bounds of integration represent the portion of each box covered in ice, and r E is the Earth's radius.The ice component of the model serves as a saline forcing on the ocean box model component.There are similar equations for when ocean circulation is in the reverse direction.
In contrast to traditional Stommel box models, rather than using transport equations for the temperature, we opt instead for thermal balance equations.The box temperatures change as heat transfers between boxes with the flow f , and as it transfers into the box via net radiation or conduction through overlying ice.Thus, for f < 0 the temperature equations are Here c w is the ocean water heat capacity.The first term in each equation accounts for net energy accumulation due to fluxes across box boundaries.The second term in Eqs. ( 3) and ( 4) represents the radiative balance, in a similar form as appears in EBMs dating back to Budyko (1969).Here α w is the ocean water albedo, and F s (θ ) is insolation, for which we use the parameterization of McGehee and Lehman (2012): Here e is the eccentricity of the Earth's orbit (presently at 0.0167), β is the obliquity (presently at 23.5 • ), and γ is longitude.This parameterization is annually averaged (no seasonal cycle), but with time-dependent orbital parameters that allow for accounting for the Milankovitch cycles.However, in our results we found that these were not strong enough to qualitatively affect the resulting model state, regarding the location of the ice margin, or Snowball Earth initiation or deglaciation, so only present-day insolation values are used.The final terms in Eqs. ( 5) and ( 6) represent a uniform geothermal heating forcing, with default F g = 0.05 W m −2 , as in Ashkenazy et al. (2013).
In Eqs.
(3) and (4), insolation is balanced by outgoing longwave radiation, integrated over the exposed ocean portion of the box, where we assume blackbody radiation from the surface using the full Stefan-Boltzmann law, with σ the Stefan-Boltzmann constant.As mentioned, ε is the effective emissivity, which is the ratio of outgoing longwave radiation emitted at the top of the atmosphere to that emitted at the Earth's surface, and therefore represents the greenhouse effect.Thus, atmospheric effects are distilled into this single parameter, which we will use as our control between climate states.
The last terms in Eqs. ( 3) and ( 4) represent ice-ocean coupling by modeling heat transfer between the ocean box and its ice cover, a key component of our model.This includes the parameter D, having units of length, which parameterizes the under-ice exchange of energy with the ocean and which we refer to as the effective thermal boundary layer.We note that D is a key unconstrained parameter, which by default we set to D = 0.05 m, and explore the model's sensitivity to this parameter in Sect. 4.Here κ w is the seawater thermal conductivity, which we take to be constant, neglecting dependence on salinity and temperature.
The ocean heat transport in this model is then quantified as As mentioned, the modeling approach taken herein does not include an active atmospheric component for the sake of keeping the model as simple as possible while focusing on the ice-ocean coupling in a reduced model sense.In this spirit, the parameter ε plays the role of the net forcing of the atmosphere on the ocean-ice system, including any heat/moisture transport within the atmosphere, which we consider to be out of our system.

Sea ice component
We largely follow the "sea glacier" treatment of GP03 and later Li and Pierrehumbert (2011), with a few noted exceptions.The details of the rheology of sea glaciers are left to Appendix A. The equation for the ice thickness, by conservation of mass, is where h(θ, t) is ice thickness, v(θ, t) is meridional ice velocity, and M(θ, t) is the ice melting/production term.The equation for v is given by a Glen's flow law (GP03), where µ is a temperature dependent viscosity parameter accounting for the non-Newtonian rheology of the ice, the details of which are left to Appendix A. The ice melting/accumulation term M (θ, t) is a departure from the treatment of GP03.Ice melting or production can occur either from heat transferred through the ice from the surface, or from heat transferred through the ocean through the effective ice-ocean thermal boundary layer D, and is given by where L is the latent heat of fusion of ice.When M(θ, t) > 0, there is net accumulation of ice, and when M(θ, t) < 0, there is net melting.The first term on the right side accounts for heat transfer through the ice, assuming a linear temperature profile in the ice from the surface T s (θ ) to the base at freezing T f .The second term accounts for heat transfer with the ocean through the parameter D; an equivalent term appears in the energy budget for the ocean box temperatures in Eqs. ( 3) and (4).Since only average ocean box temperatures are computed by our model, to prevent an artificial and arbitrary jump in temperature across the box boundary from influencing the melting term, the step function surface temperature profile T ut (t), T up (t) is regularized to a smooth T 1,2 (θ, t) for use in Eq. ( 10).Note that our model only accounts for melting and freezing at the base of the ice, and there are no terms that model melting at the upper surface or accumulation due to evaporation/precipitation forcing.
The ice surface temperature T s (θ, t) is given by a primary radiative balance, as well as a term accounting for heat transfer through the ice.The (average annual) ice surface temperature T s (θ ) is given by where c i is the specific heat, α i the albedo, and κ i the thermal conductivity of ice.The last term of Eq. ( 11) accounts for heat transfer through the ice, as in Eq. ( 10).

Model setup
Settings for the parameters are listed in Table 2.The ocean component is run on a yearly time step, with the ice dynamics subcycled on a monthly time step.(This does not imply a seasonal cycle; rather we keep the insolation constant, as given in Eq. 7.) For each ocean time step, we solve the set of differential equations for box temperatures and salinities (Eqs.2-6), and ice surface temperature Eq. ( 11) using a simple forward Euler method.At each ice time step, we solve (Eq.8) using a second-order upwind scheme after solving for the velocity in Eq. ( 9).We discretize our latitudinal domain with 100 points, solving for ice thickness h and velocity v on staggered grids, and set v = 0 at the pole for the boundary condition in Eq. ( 9).We initialize the model in an ice-free state, so any ice is formed through the model by Eq. ( 10).Following the setup of Griffies and Tziperman (1995)

Results
As mentioned, we use the effective emissivity ε as a control between model climate states.We emphasize that, as we have no atmospheric heat transport or other atmospheric effects in our model other than this parameter to crudely account for greenhouse effects, we do not attribute physical meaning to this parameter value.With a choice of ε = 0.5, which is estimated to be a reasonable value for current climate (Voigt et al., 2011), the model remains in an ice-free planet with a thermally driven poleward circulation of ≈ 62.1770 Sv and associated heat transport ≈ 3.7332 PW.We note the circulation strength is in line with results referenced in Griffies and Tziperman (1995) that give an approximate meridional circulation strength of 20 Sv from the coupled ocean-atmosphere model of NOAA's Geophysical Fluid Dynamics Laboratory (GFDL) model, and the ocean heat transport is in line with estimates for the North Atlantic (≈ 1 PW; Trenberth and Caron, 2001).The salinities of the boxes quickly mix and converge to roughly the same value of 35 psu.The equatorial surface box temperature settles to T ut = 304.0K, whereas the other three boxes converge to a closer temperature range with T up = 289.3,T lp = 289.3,and T lt = 289.4K. Hence, we have a strong, thermally dominated poleward circulation in this simulation.Depending on the choice of the ε parameter (which controls the radiative balance and parameterizes any atmospheric effects), the model can reproduce present-day partial ice cover conditions, as well as Snowball Earth global ice cover conditions, as we will now describe.

Partial glaciation
Raising the effective emissivity from ε = 0.5 to ε = 0.7, we move to an equilibrium climate state with a small, stable ice cover.To determine the role of oceanic heat transport, we run the model with the circulation rate f set to zero for comparison against the full model run.Figure 2a and b show the ice thickness and ice velocity profiles with and without ocean circulation, and we see the ocean circulation is effective at reducing ice thickness as well as pulling the ice margin north.
Without the additional heat from the equatorial region moving poleward, the polar region remains cool, facilitating ice growth.
The ice cover is approximately 700 m thick at the pole in the full model run, much thicker than current sea ice cover; however, this is consistent with partial glaciation results from GP03.It is because the ice is this thick that viscous flow effects need to be considered.There is a strong response in the ice velocity, largely due to the thicker ice cover when ocean circulation is not included, resulting in stronger viscous flow.The ice velocity reaches its maximum just before the ice margin (approximately 2 km yr −1 in the full simulation).We also note that the ice surface temperature seen in Fig. 2c, calculated by Eq. ( 11), is in line with the air surface temperature forcing used in the experiments of GP03.In Fig. 2d we see the accumulation term becoming negative near the ice margin, indicating a region of net melting that stabilizes the ice margin.
The steady-state ocean circulation strength in this partial ice cover scenario is ≈ 44.44 Sv, which is closer to the numerical results of the GFDL model referenced in Griffies and Tziperman (1995) than the ice-free run.As with the ice-free runs, the box salinities quickly mix to the same value of approximately 35.97 psu, while the surface equatorial box temperature settles to T ut ≈ 281.6 K, and the other boxes mix to T up ≈ 271.63,T lp ≈ 271.65,T lt ≈ 271.71K.By increasing the effective emissivity ε, the model steady-state ice profile moves smoothly further equatorward, until a large ice cap instability threshold is reached.When ice appears south of this instability threshold, the entire planet is covered in ice and a Snowball Earth state is reached.

Global glaciation
To approximate the Neoproterozoic climate in our model, we lower insolation to 94 % of its current value, accounting for a weaker, younger sun (Voigt et al., 2011), and raise the geothermal heat flux from 0.05 to 0.08 W m −2 .Raising the effective emissivity from ε = 0.7 to ε = 0.83, we move from a climate state with a small, stable ice cover to global ice cover, a Snowball Earth, as shown in Fig. 3 progresses slowly towards equilibrium between ice thickness and geothermal heat, so simulations in this regime are run for 500 000 years to reach steady state in the ice component.In the comparative run with no ocean circulation, the system will not reach steady state, as there is no geothermal heat flux reaching the ice layer, and as a result there is a positive ice accumulation rate at all latitudes.In this scenario, the ocean circulation weakens to ≈ 16.5 Sv, with an associated heat transport of ≈ 6.27 × 10 −4 PW.The resulting ice thickness is relatively uniform with a polar maximum of 2 km, down to 1300 m near the Equator.The back-pressure term in the boundary condition for ice velocity in Eq. (A2) brings the ice velocity to zero at the Equator (Fig. 3b), the effect of which is to move the maximum velocity (≈ 1 km yr −1 ) to a location in the midlatitudes.
Examining the model results between these very different climate states of a small, stable ice cap at ε = 0.7 and Snowball Earth at ε = 0.83, we find a threshold where the model transitions.With an effective emissivity of ε = 0.82, the equilibrium climate state is near the large ice cap instability threshold, and we get a strong response from the ocean circulation.In Fig. 4, we show the results of a simulation with and without ocean circulation dynamics.We observe in Fig. 4a that, without oceanic heat transport, we get a Snowball Earth, but with oceanic heat transport, the ice line is held back from the Equator.As with the previous Snowball Earth state, the ocean circulation strength here, 17.5 Sv (poleward) is weaker than in the small ice cap simulation.However, even this weakened circulation, and thus weakened oceanic heat transport, is still enough to drive the climate into a snowball state if turned off, demonstrating strong sensitivity in this regime.

Model sensitivity to ocean and atmosphere parameters
We have already seen the model sensitivity to the effective emissivity parameter ε, which we used to drive the model through vastly different climate states in Sect.3.There are other key sensitivities the model exhibits, which we now discuss, namely in the parameterization of energy transfer between the ocean and ice components, as well as the scaling strength of ocean circulation.By the nature of the large ice cap instability, the threshold effective emissivity of ε = 0.82 seen in Sect.3.2 is sensitive to the parameter choices representing these key processes.We also explore the model's bistability in the small ice cap and large ice cap instabilities.

Ocean-ice energy transfer parameterization
The parameter D is responsible for parameterizing the transfer of energy between the ocean and ice systems, a so-called effective thermal boundary layer, and is the least constrained parameter in our model.It appears in equations for the surface ocean box temperatures (Eq.3), Eq. ( 4), and in ice melting/production (Eq.10).The role of D has competing effects in these two equations.For the ocean box temperatures, D appears in the term corresponding to energy loss due to the presence of the ice cover, and thus increasing D cools the ocean.However, this energy loss is balanced in the  There is a particularly strong response in the resulting melting/accumulation term which was used to constrain the value of D. Furthermore, the instability near the large ice cap threshold is reflected in nearby values of D to the default D = 0.05, resulting in global ice cover for ε = 0.82.Note that the discontinuity in ice accumulation rate for D = 0.1 is an artifact of the regularization of Eq. ( 10) as h approaches 0.
system by the ice melting/production term, where increasing D encourages melting.In the limiting case of D = 0, the ocean and basal ice would be forced to the same temperature, and increasing D further insulates the ocean from the ice, by slowing the heat transfer between the two component.sThere are also other feedbacks in the system, notably the indirect effect of D on the strength of the ocean circulation, and thereby oceanic heat transport, which we have already seen can strongly affect ice cover.
We explore the model's sensitivity to values of this parameter over two orders of magnitude in Fig. 5 in both the small ice cap regime (ε = 0.75) and near the large ice cap instability (ε = 0.82, insolation at 94 % current values, F geo = 0.08 W m −2 ), showing ice thickness profiles and ice melting/accumulation rate at the end of a 50 000-year run to equilibrium.In the small ice cap regime, increasing D steadily reduces ice thickness, though the ice margin remains in a small high-latitude range across changes in D, apart from the large end of the range, where the margin pushes further equatorward.The ocean circulation (poleward) also steadily reduces with increased D, from 32.9 Sv circulation with D = 0.01 down to 20.5 Sv circulation for D = 1.There is, however, a strong response in the ice melting/accumulation term, where large values of D yield small magnitude melting terms.This proved to be critical in melting excess ice, which we saw in our experiments with changing solar forcing discussed below in Sect.4.3.

Hydraulic constant
The default value of the hydraulic constant k = k 0 in Eq. ( 1) from Griffies and Tziperman (1995) results in a circulation strength of ≈ 32.7 Sv, in line with estimates of present-day mean meridional ocean circulation.This may not be appropriate for the Neoproterozoic, warranting an examination of the sensitivity of the model to changes to this parameter, explored in Fig. 6.We again show results for the same experimental setups as in Sect.4.1.Here we see a stronger effect on the ice margin in the small ice cap regime, where smaller hydraulic constants reduce heat transport from the tropics to the poles, thereby facilitating ice growth (note that circulation always remains poleward in this regime).In the large ice cap regime, we again see the sensitivity of the large ice cap instability, as deviating from the default circulation strength constant, in either direction moves the threshold emissivity away from ε = 0.82, and the resulting climate state is a Snowball Earth.

Model bistability
As discussed, a well-known feature of EBMs is hysteresis with respect to radiative forcing.To study this in our model, we set up experiments to investigate the hysteresis loop in the small ice cap instability and the large ice cap instability.In the small ice cap case, we began with conditions that resulted in a small ice cap, gradually increased radiative forcing by decreasing the emissivity ε until the ice cap melted away, and then increased ε to its starting value.With each Snowball Earth (ε = 0.6 −→ 0.9), then subsequent warming (ε = 0.9 −→ 0.6).A strong hysteresis loop is seen in the ice margin, as the model is unable to escape Snowball Earth even as radiative forcing is raised to levels associated with ice-free states, although strong enough increases in radiative forcing will eventually melt ice away from global ice cover.Circulation remains poleward, ever increasing through the Snowball Earth state manual radiative forcing increase.
incremental change in ε, the model is run for 50 000 years to ensure that the model is in equilibrium.A similar experiment was run for the large ice cap case, except with the forcing changes in cooling first, until Snowball Earth is reached, and then subsequent warming.The results of these experiments are shown in Fig. 7.The expected hysteresis loops in each scenario manifest themselves in the left column of Fig. 7.In the small ice cap case, The ice margin begins with ε = 0.7 near 60 • latitude, and following the red path through decreasing ε, we see the ice cap abruptly change from a margin at 80 • latitude, completely disappearing at ε = 0.55.As the forcing is then cooled through the blue path, the model remains ice-free until the appearance of a small ice cap near ε = 0.62, returning to its starting location at ε = 0.7.Looking at the response of the ocean circulation through the experiment, we see that there is a small response in the strength of circulation in the hysteresis loop.As alluded to in Sect.4.1, this behavior was actually used to constrain the value of D = 0.05, as larger values of D did not yield adequate melting rates to melt ice once it appeared, resulting in very large hysteresis loops.It is worth noting that the hysteresis loop in this model is present with only ocean heat transport, in the absence of an explicit atmospheric heat transport parameterization, which has previously been found to be necessary for such bistability (Held and Suarez, 1974).
In Fig. 7c, we see the familiar large ice cap instability hysteresis loop, where as the climate is cooled through increasing ε, the ice margin gradually increases until the large ice cap instability is reached near 20 • latitude, beyond which Snowball Earth is reached.As the climate is subsequently warmed, Snowball Earth is maintained, but notably the ocean circulation remains poleward, continuing to transport heat to the poles, and, in fact, with a strong increase in warming (up to ε = 0.4), Snowball Earth is escaped in the model.

Conclusions
We have presented a low-dimensional conceptual climate model consistent with elements of classical low-dimensional models that is able to reproduce both present-day, partial ice cover climate, and a Snowball Earth global ice cover climate.The radiative balance terms similar to those in EBMs produce states of ice cover consistent with classical EBM results of two stable solutions, a small ice cover and a global ice cover, and an unstable solution, a large but finite ice extent.A summary of these results is given in Table 1.
Our primary interest is investigating the role of dynamic ocean circulation in the initiation of Snowball Earth, particularly in the large ice cap instability.We find in parameter regimes where there is no ice or a small ice cap that the ocean circulation is expectedly thermally dominated and poleward.Discounting ocean circulation altogether allows for easier transition into global ice cover.We find, together with the effective emissivity which parameterizes the atmospheric component, that the energy transfer between the ice and ocean components plays a crucial role in determining the model's resulting climate state.
The results from our ice component are largely in line with those of GP03, in both ice thickness and position of the ice margin, despite some key differences in modeling approach.There is a notable difference in the ice velocity, however, in that our computed ice velocities are an order of magnitude less those reported in GP03.One possible reason for this is that the viscosity parameter µ, which GP03 are only required to calculate once due to a static surface temperature forcing, is recalculated in our model in response to the dynamic surface temperature in Eq. (11).Our steady-state surface temperature (shown in Fig. 4) is cooler than the forcing used in GP03 partial glaciation case, and thus our cooler surface temperature creates more viscous ice, slowing down ice flow.The ice thickness profiles in our Snowball Earth experiments vary smoothly by a few hundred meters from pole to Equator, in contrast to sea glacier models that obtain a more uniform thickness in snowball state, as in Li and Pierrehumbert (2011).Our results also qualitatively agree with that of Pollard and Kasting (2005)  The ocean in our Snowball Earth scenario has a considerably weaker circulation strength of 16 Sv than present-day estimates, and while this is not indicative of a stagnant ocean, it is not as strong as the circulation strengths of approximately 35 Sv that were achieved in the study of Ashkenazy et al. (2013).Their 2-D simulations allowed for resolution of both vertical mixing and horizontal eddies, and while they also did not have land in their simulations, they did have an underwater ocean ridge that had a highly localized and higher strength geothermal heat forcing corresponding to a spreading center which, together with strong vertical salinity profiles, drove the strong circulation.
The main conclusion we reach from this study is that ocean circulation and its associated heat transport play a vital role in determining the global climate state and ice cover.We have seen in the partial glaciation case that the ocean circulation severely inhibits ice growth, and we have seen in the nearglobal glaciation case that even in that state's severely weakened ocean circulation, a lack of oceanic transport leads to a drastically different Snowball Earth state.In particular, our study found that both heat flux between the ice and ocean and the value of the circulation constant that controls circulation strength (tuned to present-day conditions) played a crucial role in determining the global climate state.Moreover, we find bistability in the small and large ice cap instability regimes with regard to radiative forcing.
Our inclusion of a simple parameterization of under-ice ocean heat transport, mediating dynamic oceanic heat transfer not present in traditional EBM diffusive heat transport models, is here found to be crucial in determining the steadystate climate regime.This warrants further investigation in a GCM setting in order to examine the role played by thermal processes in the ice and ocean that account for heat transport in the sub-ice cover layer.To be more direct, we find that oceanic heat transport is crucial to understanding Snowball Earth initiation, that the ice cover affects it significantly, that the results are sensitive to the water-ice thermal coupling and the factors driving the circulation, and that it is therefore worthwhile using GCMs to investigate these factors in detail.
While atmospheric effects were largely neglected for simplicity and because our focus was on the role of oceanic heat transport in the Snowball Earth setting, including an atmospheric component would be the natural progression of this work, particularly a precipitation-evaporation component that would facilitate ice surface melting/accumulation, though it is worth noting that we were able to reproduce both present-day and Snowball Earth conditions without a precipitation-evaporation forcing.

Data availability
Model code has been uploaded to http://cims.nyu.edu/~comeau/model_code.tar.gz.It and the data are also freely available from the corresponding author.

Figure 1 .
Figure 1.Hemispheric four-box arrangement.Boxes 1 and 2 are the surface ocean boxes of depth d u , and Boxes 3 and 4 are the deep ocean boxes of depth d l .The water in Box i has (well-mixed) salinity S i and temperature T i .The boundary between polar and equatorial boxes is at latitude ζ .Ice cover sits atop the surface boxes with height h and the ice margin at η.The arrows between the boxes represent density-driven circulation f .

Figure 2 .
Figure 2. (a) Ice thickness h and (b) meridional velocity v profiles in partial ice cover scenario with effective emissivity ε = 0.7, with and without ocean circulation.Oceanic transport is seen to dramatically affect the ice predictions.(c) Ice surface temperature T s and (d) ice basal accumulation rate in partial ice cover scenario with effective emissivity ε = 0.7, with and without ocean circulation.The accumulation term stabilizes the ice margin, and the model produces reasonable values of ice surface temperatures.M is positive over the region of ice cover, indicating net accumulation of ice, and negative values of M beyond the ice margin indicate any ice present would be melted.

Figure 3 .
Figure 3. (a) Ice thickness h and (b) velocity v profiles with 94 % insolation, F geo = 0.08 W m −2 , and effective emissivity ε = 0.83, the Snowball Earth scenario.(c) Ice surface temperature T s , which is seen to not be impacted by ocean circulation.In (d), the simulation without ocean circulation cannot reach steady state due to geothermal heat not reaching the ice layer.

Figure 4 .
Figure 4. (a) Ice thickness h and (b) velocity v profiles with 94 % insolation and ε = 0.82, with and without ocean circulation.We note neglecting oceanic heat transport leads to drastically different global climate states.(c) Ice surface temperature T s and (d) ice basal accumulation M, which without oceanic heat transport is positive everywhere.

Figure 5 .
Figure5.Effect of D parameter on equilibrium climate state, small ice cap regime (ε = 0.75, top), and large ice cap regime (ε = 0.82, 94 % insolation, F geo = 0.08 W m −2 bottom).There is a particularly strong response in the resulting melting/accumulation term which was used to constrain the value of D. Furthermore, the instability near the large ice cap threshold is reflected in nearby values of D to the default D = 0.05, resulting in global ice cover for ε = 0.82.Note that the discontinuity in ice accumulation rate for D = 0.1 is an artifact of the regularization of Eq. (10) as h approaches 0.

Figure 6 .
Figure6.Effect of k 0 parameter on small ice cap regime (top) and large ice cap regime (bottom).Smaller circulation constants reduce heat transport from the tropics to the poles, thereby facilitating ice growth.The large ice cap instability is again reflected in the sensitivity to changes to the hydraulic constant, with the shown deviations resulting in Snowball Earth.

Figure 7 .
Figure7.Hysteresis experiments, radiative forcing.Top: small ice cap regime -forcing changes in ε, warming first (ε = 0.7 −→ 0.5), then cooling (ε = 0.5 −→ 0.7).Resulting ice margin (left) and circulation (right) are shown, demonstrating a hysteresis loop in the ice margin.Bottom: large ice cap regime -with insolation set to 94 % present values, forcing changes in ε, cooling first to achieve Snowball Earth (ε = 0.6 −→ 0.9), then subsequent warming (ε = 0.9 −→ 0.6).A strong hysteresis loop is seen in the ice margin, as the model is unable to escape Snowball Earth even as radiative forcing is raised to levels associated with ice-free states, although strong enough increases in radiative forcing will eventually melt ice away from global ice cover.Circulation remains poleward, ever increasing through the Snowball Earth state manual radiative forcing increase.

Table 1 .
Summary of results.Asterisk indicates simulations with insolation at 94 % of present-day values.Negative circulation values indicate surface poleward flow.

Table 2 .
Physical parameters used in simulations.See text for details on model initialization and settings of unconstrained parameters.
, with initial box temperatures T ut = 298 K, T up = T lp = T lt = 273 K and salinities S ut = 36.5 psu, S up = 34.5 psu, and S lp = S lt = 35 psu.Models that equilibrate to a climate state with no ice or a small ice cap typically reach an equilibrium climate steady state after 10 000 model years; however, when the model reaches a state of global ice cover, a longer period of about 100 000 years is needed before equilibrium is reached.www.earth-syst-dynam.net/7/937/2016/Earth Syst.Dynam., 7, 937-951, 2016 944 D. Comeau et al.: Ocean heat transport in Snowball Earth (e.g., ice thickness, vewww.earth-syst-dynam.net/7/937/2016/Earth Syst.Dynam., 7, 937-951, 2016 948 D. Comeau et al.: Ocean heat transport in Snowball Earth locity, ice accumulation rates, surface temperature) in the steady-state Snowball Earth scenario, with the notable exception of a sharper decline in tropical ice thickness, and while we are confident that adopting their more developed ice model would not significantly change the results of our model presented here in the Snowball Earth scenario (apart from a thinner tropical ice), it would certainly be of interest to study transient behavior and parameter ranges that lead to climate regimes other than Snowball Earth.