The interaction between sea ice and salinity-dominated ocean circulation: implications for halocline stability and rapid changes of sea ice cover

Abstract

Changes in the sea ice cover of the Nordic Seas have been proposed to play a key role for the dramatic temperature excursions associated with the Dansgaard–Oeschger events during the last glacial. In this study, we develop a simple conceptual model to examine how interactions between sea ice and oceanic heat and freshwater transports affect the stability of an upper-ocean halocline in a semi-enclosed basin. The model represents a sea ice covered and salinity stratified Nordic Seas, and consists of a sea ice component and a two-layer ocean. The sea ice thickness depends on the atmospheric energy fluxes as well as the ocean heat flux. We introduce a thickness-dependent sea ice export. Whether sea ice stabilizes or destabilizes against a freshwater perturbation is shown to depend on the representation of the diapycnal flow. In a system where the diapycnal flow increases with density differences, the sea ice acts as a positive feedback on a freshwater perturbation. If the diapycnal flow decreases with density differences, the sea ice acts as a negative feedback. However, both representations lead to a circulation that breaks down when the freshwater input at the surface is small. As a consequence, we get rapid changes in sea ice. In addition to low freshwater forcing, increasing deep-ocean temperatures promote instability and the disappearance of sea ice. Generally, the unstable state is reached before the vertical density difference disappears, and the temperature of the deep ocean do not need to increase as much as previously thought to provoke abrupt changes in sea ice.

1 Introduction

The mechanism behind Dansgaard–Oeschger (D–O) events, the abrupt climate changes which occurred repeatedly during the last glacial period, is still debated. One of the hypothesized mechanisms involves a retreating North Atlantic sea ice cover (Broecker 2000; Gildor and Tziperman 2003; Masson-Delmotte et al. 2005; Li et al. 2005). Due to the ice-albedo feedback, the strong insulating effect, and impact on the atmospheric circulation, a changing sea ice cover in the Nordic Seas is a likely agent for the abrupt warming observed on Greenland during D–O events (e.g. Dansgaard et al. 1993; Severinghaus and Brook 1999). Thus, a retreating sea ice cover can explain the transition from cold stadial conditions to warmer, interstadial conditions on Greenland and in the northern North Atlantic (Li et al. 2005). In addition, a changing sea ice cover can explain several other features of the abrupt climate events, including precipitation changes (e.g. Peterson et al. 2000; Li et al. 2005, 2010), and changes in dust content (Mayewski et al. 1997). Here, we investigate the interactions between sea ice and a salinity-dominated ocean circulation in a conceptual model of the Nordic Seas to explore the possible role of sea ice in the abrupt warming associated with D–O events.

Internal instability in the oceanic thermohaline circulation, or the Atlantic Meridional Overturning Circulation (AMOC) has also been suggested as a mechanism for rapid climate change (Broecker et al. 1990; Tziperman 1997; Marotzke 2000; Ganopolski and Rahmstorf 2001). Broecker et al. (1990) first suggested internal oscillations in the volume transport of the AMOC as a mechanism for D–O events. The idea is that a strong circulation transports more heat from the South Atlantic to the North Atlantic, thereby warming the atmosphere at high latitudes. A weak, or an off state of the AMOC could lead to less heat import to the North Atlantic causing cold stadial conditions. Several modeling studies have simulated an AMOC with a strong and a weak state (e.g. Stocker and Wright 1991; Fanning and Weaver 1997; Manabe and Stouffer 1995; Stouffer et al. 2006). However, paleoclimatic reconstructions have shown active formation of North Atlantic Deep Water also during stadials (Yu et al. 1996). Latitudinal shifts in the North Atlantic Deep Water formation site is also hypothesized to account for the D–O events (Ganopolski and Rahmstorf 2001; Arzel et al. 2010; Colin de Verdiere and Raa 2010; Sevellec and Fedorov 2015). Freshwater forcing has typically been the agent to force the different states of the circulation, and the duration of the cold and the warm states generally depend on the duration of the freshwater forcing (Ganopolski and Rahmstorf 2002). On the other hand, Marotzke (1989), Winton (1993) and Winton and Sarachik (1993) found self-sustaining oscillations with strong and weak overturning due to increasing deep-ocean temperatures. Their idea is that freshwater input at the poles causes a strong, vertical salinity gradient (halocline) which suppresses deep-water formation. The stability of the high latitude halocline breaks down as diffusion from the south heats the deep ocean. Deep-water formation becomes active again, and warm, light water is transported to the high latitudes which then again stabilizes the water column. These deep-decoupling oscillations bear resemblances to the D–O stadial and interstadials (Timmermann et al. 2003). Seen in isolation, a changing AMOC gives a temperature response on Greenland which is too slow and too small compared to the observations. However, the addition of a responsive sea ice cover in the Nordic Seas would enhance the temperature response to small internal oscillations in the ocean.

Paleoclimatic reconstructions support a warming of the subsurface Nordic seas during cold stadial conditions (Dokken and Jansen 1999; Rasmussen and Thomsen 2004; Dokken et al. 2013). In particular, it has been proposed that the warming of the subsurface ocean during cold stadial conditions could explain the rapid retreat of sea ice in the Nordic Seas at the start of each interstadial (Dokken et al. 2013). From a marine sediment core from the Faeroe-Shetland channel, Dokken et al. (2013) showed systematic changes in the hydrography of the Nordic Seas during cold stadials and warm interstadials. During stadial conditions, the hydrography was similar to the Arctic Ocean today, with a cold, fresh surface layer above a warmer Atlantic layer. A sharp halocline marks the transition between the two different water masses, protecting the surface layer from the heat below. The authors proposed that the vertical density gradient and the halocline would gradually weaken due to a slow warming of the deep ocean. Eventually, as the vertical density difference disappeared, the stability of the water column would break down. When the halocline is gone, the Atlantic layer would mix to the surface, melting the sea ice cover. This scenario is dependent on a mechanism that gradually increases the temperature of the Atlantic layer.

A different scenario involves slow changes in the freshwater supply to the ocean. By using an analytical two-layer ocean model, Nilsson and Walin (2010) showed that the salinity-dominated (estuarine) circulation and halocline can break down even before the vertical density gradient vanishes. If this is also true for ice-covered conditions, which can entail additional feedbacks, the convective overturning proposed by Dokken et al. (2013) could occur at low freshwater inputs, even before the vertical density gradient disappears. This would be highly relevant for cold, dry glacial climates. An interesting question is at what point the halocline disappears given a constant temperature below the halocline.

However, the properties of thermohaline circulations are dependent on the relation between flow strength and density differences. In thermohaline models which include vertical mixing, the flow strength is sensitive to how the vertical velocity is represented (Lyle 1997; Huang 1999; Nilsson and Walin 2001). Assuming a fixed external energy supply to the small scale mixing, Nilsson and Walin (2010) represented the vertical velocity with a stratification-dependent vertical diffusivity. The result is a diapycnal velocity which decreases with increasing stratification and an ocean circulation which acts in the same way. If instead, diapycnal velocity is represented with a constant diffusivity, the ocean circulation increases with increasing stratification (e.g. Bryan 1986; Zhang et al. 1999; Nilsson and Walin 2001). As the classical box model by Stommel (1961) involves a flow which speeds up with increasing density difference, the latter representation bears resemblances to the typical Stommel model.

Thus, the stability of the salinity-dominated circulation depends on the assumed properties of the diapycnal mixing. In the conceptual model by Nilsson and Walin (2010) the authors found that there is a minimum freshwater input below which the halocline solution no longer exists. In models where the flow strength increases with density differences, there exists a solution for all freshwater forcing for the salinity-dominated circulation (Zhang et al. 1999; Longworth et al. 2005). If this is the case, other factors, such as increasing subsurface temperatures, are needed to initiate convective overturning. However, Longworth et al. (2005) and Guan and Huang (2008) showed that the presence of wind-driven, horizontal salt transport destabilizes the thermohaline circulation, and removes the equilibria at weak freshwater supplies. Therefore, additional factors may have the potential to introduce convective overturning due to low freshwater forcing in models with a constant diffusivity as well.

As the proposed mechanism by Dokken et al. (2013) involves a sea ice covered Nordic Seas, we add the sea ice model from Thorndike (1992) to the simple two-layer ocean model of Nilsson and Walin (2010). The aim of the study is to investigate how the presence of sea ice affects the dynamics and stability of the salinity-dominated circulation as the climatic conditions are varied. As the detailed features of the diapycnal mixing are not well known, we will in this study examine the two different vertical mixing representations discussed above. The experiments are divided into an “energy-constrained model”, where we retain the diapycnal velocity from Nilsson and Walin (2010), and a “constant-diffusivity model” where the diapycnal velocity is parametrized with a constant diffusivity. The vertical mixing is expected to change with a changing sea ice cover and climate. Therefore, even though the first case is more physically consistent, it is useful to study different representations of the vertical mixing.

Interactions between sea ice and thermally-dominated thermohaline circulation have been studied in conceptual models by for instance Yang and Neelin (1993, 1997) and Jayne and Marotzke (1999). Using a Stommel-type box model coupled to an energy-balance representation of the atmosphere, Jayne and Marotzke (1999) found that sea ice related feedbacks destabilized the thermally-dominated mode of the thermohaline circulation. In their model, feedbacks between sea ice, meridional temperature gradients and atmospheric moisture transport were of key importance, whereas brine rejection associated with sea ice formation played a more central role in the box model of Yang and Neelin (1993). However, less attention has been paid to how sea ice feedbacks affect a salinity-dominated thermohaline circulation. Stigebrandt (1981) developed a process-oriented two-layer model of the Arctic Ocean, which incorporates inflows and outflows, air-sea heat exchange, and sea ice dynamics. Interestingly, Stigebrandt found that the sea ice cover in the model vanished abruptly when the freshwater supply to the Arctic Ocean was reduced below some threshold value. In fact, Nilsson and Walin (2010) essentially used the same geostrophic flow representation as Stigebrandt (1981). However, Stigebrandt’s model contains additional complexities which makes it difficult to examine the sea ice related feedbacks in more detail.

In this study, we develop a conceptual model that is simple enough to admit analytical analyses of the key feedbacks. First, the two-layer ocean model by Nilsson and Walin (2010) and the simplified sea ice model from Thorndike (1992) will be presented in Sect. 2. We review the response of the model to freshwater forcing in the absence of sea ice in Sect. 3. Then, the response of the model to freshwater forcing in the presence of sea ice is presented in Sect. 4. We compare the stabilizing effect of sea ice in the energy-constrained model with the constant-diffusivity model in Sect. 4.2. In Sect. 5, we investigate the effect of increasing deep-ocean temperatures on the stability of the system. In Sect. 6, we discuss the impact of the results on the Dokken et al. (2013) hypothesis.

2 Model formulation

Here we will show how the conceptual two-layer ocean model from Nilsson and Walin (2010) can be coupled to the simplified sea ice model from Thorndike (1992). The two model components are presented below, and the coupled ocean–ice system can be seen in Fig. 1. Variables that are not defined in the text can be found in Table 1. The simplified ocean and ice models used here are well established and used in numerous studies. The novelty of this study is to examine the physics that results from their coupling. Useful background information on the physical motivation and derivation of two-layer ocean models can be found in e.g. Gnanadesikan (1999), Nilsson and Walin (2001) and Johnson et al. (2007). Bitz and Roe (2004) provide a concise derivation of the sea ice model from Thorndike (1992).

Table 1 Definitions of variables for the model

Fig. 1

A sketch of the coupled model. S and T are the surface layer salinity and temperature, respectively, while $S_A$ and $T_A$ are the salinity and temperature of the deeper, Atlantic layer. $M_D$ and $M_G$ are the diapycnal flow and geostrophic outflow, respectively. $F_{riv}$ denote the freshwater input while $F_{ice}$ denote the sea ice export. The surface layer depth is represented by H, the ocean heat flux by Q, and the sea ice thickness by h

The ocean model consists of a cold, light surface layer, and a warmer, homogeneous deeper layer with the same hydrography as the open ocean. The circulation is that of a salinity-dominated regime, in the sense that supply of freshwater at the surface establishes a stable stratification. The topography is that of a semi-enclosed basin (Fig. 1). The basin receives freshwater from river run-off, precipitation, and ice-import, together denoted as $F_{riv}$. The ocean flow consists of the diapycnal flow $M_D$, and the horizontal geostrophic outflow $M_G$. The diapycnal flow increases the depth of the light surface layer H, while $M_G$ thins the surface layer.

The subsequent derivations follow Nilsson and Walin (2010). By applying a time-dependent version of the Knudsen’s relations (Knudsen 1900), the conservation of volume and salinity can be written as

$$\begin{aligned} A\frac{dH}{dt}&= M_D - M_G + F_{net}, \end{aligned}$$

(1)

$$\begin{aligned} A\frac{dHS}{dt}&= S_AM_D - SM_G, \end{aligned}$$

(2)

where S is the salinity of the surface layer, $S_A$ is the salinity of the Atlantic layer, $F_{net} = F_{riv} - F_{ice}$ is the net freshwater input to the surface, $F_{ice}$ is the sea ice export, and A is the area of the Nordic Seas. Note that in steady-state, the above equations become the classical Knudsen’s relations.

By combining Eqs. 1 and 2, we get

$$\begin{aligned} AH\frac{d{\Delta } S}{dt}=-{\Delta } S M_D + (S_A-{\Delta } S)F_{net}, \end{aligned}$$

(3)

where ${\Delta } S = S_A-S$, and $S_A$ is a constant. This equation can be simplified further: In the open ocean, away from coastal regions, one typically finds that $S_A \sim 35$ and ${\Delta } S \sim 3$. Therefore, ${\Delta } S/S_A$ is significantly smaller than unity. This is widely used in oceanography to derive the Boussinesq approximation equations (Vallis 2006) where the surface freshwater flux is substituted by a virtual salt flux boundary condition and the freshwater flux is neglected in the volume conservation equation. We exploit this and derive simplified equations formally valid in the limit where ${\Delta } S/S_A\ll 1$. Note that the restriction ${\Delta } S/S_A\ll 1$ also implies that $F_{net}$ is small compared to the volume flows $M_D$ and $M_G$. This is consistent with the Boussinesq approximation. We can demonstrate this by examining the steady-state version of Eq. 3 in the limit ${\Delta } S/S_A\ll 1$, which yields ${\Delta } S M_D \sim S_A F_{net}$, or equivalently $F_{net} /M_D \sim {\Delta } S/S_A \ll 1$. From the steady-state volume conservation equation it follows that $M_D \sim M_A \gg F_{net}$. Therefore, when ${\Delta } S/S_A\ll 1$, we obtain the following approximate relations

$$\begin{aligned} A\frac{dH}{dt}&= -M_G + M_D, \end{aligned}$$

(4)

$$\begin{aligned} AH\frac{d{\Delta } S}{dt}&= -{\Delta } SM_D + S_AF_{net}. \end{aligned}$$

(5)

These equations yield a model that is simpler to analyse analytically, without changing the results qualitatively.

To derive representations of the flows $M_G$ and $M_D$, we rely on physical reasoning based on the classical thermocline scaling (e.g. Welander 1971, 1986). To begin with, the horizontal flow, related to $M_G$, is taken to be in geostrophic and hydrostatic balance, which leads to the thermal wind balance. Further, we assume that the outflow for our “Nordic Seas” domain occurs in a baroclinic boundary current attached to a meridional boundary; analogously to how the low-salinity Arctic surface water exits from the Nordic Seas in the East Greenland Current today (e.g. Nilsson et al. 2008). Thus, we assume that the basic north–south density variation is similar to the east–west density variation in the outflow region (see e.g. Gnanadesikan 1999; Park and Bryan 2000; Nilsson and Walin 2001; Johnson et al. 2007, for further discussion). By applying the thermal wind balance to compute the volume transport in a two-layer baroclinic frontal current, one obtains (e.g. Nilsson and Walin 2001)

$$\begin{aligned} M_G = \frac{g{\Delta }\rho H^2}{2f\rho _0}, \end{aligned}$$

(6)

where

$$\begin{aligned} {\Delta }\rho = \rho _0\beta {\Delta } S -{\Delta }\rho _T, \end{aligned}$$

(7)

and ${\Delta }\rho _T$ is the density difference due to the temperature difference between the two ocean layers.

Previous studies have shown that the resulting dynamics of the ocean model is sensitive to the assumed physical features of the diapycnal flow $M_D$ (e.g. Lyle 1997; Guan and Huang 2008; Nilsson and Walin 2001). Therefore, we examine two different representations of $M_D$ that yield different dynamical behaviours of the model. Following the classical thermocline scaling, we assume that the vertical stratification is governed by an advective-diffusive balance (Munk 1966). Conservation of buoyancy is now used to obtain the scaling of $M_D$. By taking the two-layer model variables H and ${\Delta }\rho$ to represent the vertical scale and density variation of the stratification, respectively, the advective–diffusive balance suggests the following relation (Welander 1971, 1986)

$$\begin{aligned} M_D = \frac{A\kappa }{H}, \end{aligned}$$

(8)

where $\kappa$ is a constant vertical turbulent diffusivity. However, a strong background stratification tends generally to suppress the turbulent diffusivity (e.g. Huang 1999). If the turbulent diffusivity is taken to be constrained by the energy input from tides and winds, one can argue for the following expression for $\kappa$ (Nilsson and Walin 2001)

$$\begin{aligned} \kappa = \frac{\epsilon }{g {\Delta }\rho }, \end{aligned}$$

(9)

where $\epsilon$ is the rate of work against the buoyancy force associated with the vertical mixing. In this formula, $\kappa$ decreases with increasing density stratification. Using the diffusivity formula (Eq. 9) in Eq. 8, we obtain the following alternative expression for the diapycnal flow

$$\begin{aligned} M_D = \frac{A\epsilon }{g{\Delta }\rho H}. \end{aligned}$$

(10)

We note that this is essentially the formula from Kato and Phillips (1969), describing the entrainment rate across a mixed layer with a density jump at its base. Kato and Phillips derived their formula by equating a constant mechanical energy input with the rate of increase in potential energy caused by the entrainment.

We are interested in the steady-state solutions where $M = M_D = M_G$ and ${\Delta } SM = S_AF_{net}.$ We get the steady-state flow and surface layer depth by combining Eq. 6 with Eqs. 8 or 10. For the energy-constrained model, we get;

$$\begin{aligned} M&= \left( \frac{A^2\epsilon ^2}{g2\rho _0f{\Delta }\rho }\right) ^{1/3}, \end{aligned}$$

(11)

$$\begin{aligned} H&= \left( \frac{A\epsilon 2\rho _0f}{g^2{\Delta }\rho ^2}\right) ^{1/3}. \end{aligned}$$

(12)

For the constant-diffusivity model;

$$\begin{aligned} M&= \left( \frac{A^2\kappa ^2g{\Delta }\rho }{2f\rho _0}\right) ^{1/3}, \end{aligned}$$

(13)

$$\begin{aligned} H&= \left( \frac{A\kappa 2\rho _0f}{g{\Delta }\rho }\right) ^{1/3}. \end{aligned}$$

(14)

2.1 Sea ice model

For the sea ice component, we follow Bitz and Roe (2004) and Thorndike (1992) which estimated equilibrium sea ice thickness h by assuming a quasi-steady balance between sea ice growth G and decay D. The sea ice thickness is estimated with only the main essential elements included; the net atmospheric flux incident at the top of the ice, and the heat flux Q from the ocean at the base of the ice. Making use of steady-state heat conduction through sea ice, and an atmosphere which is in thermal radiative equilibrium with the sea ice, we calculate G and D following Bitz and Roe (2004);

$$\begin{aligned} G(h)&= \frac{\tau }{L}\left[ \frac{A_i + BT_i(h)}{n_w} - \frac{F_{atm}}{2} - Q\right], \end{aligned}$$

(15)

$$\begin{aligned} D&= \frac{\tau }{L}\left[ -\frac{A_i}{n_s} + \frac{F_{atm}}{2} + Q + (1 - \alpha )F_{sw}\right], \end{aligned}$$

(16)

where the temperature of the ice $T_i$ is given by

$$\begin{aligned} T_i(h)&= \left( \frac{n_wh}{kn_w + Bh}\right) \left( -\frac{A_i}{n_w} + \frac{F_{atm}}{2}\right). \end{aligned}$$

(17)

$F_{atm}$ is the atmospheric meridional heat advection, $F_{sw}$ is the shortwave insolation, $A_i$ and B are the coefficients of the linearized Stefan Boltzmann’s law, Q is the heat flux from the ocean to the ice, $\alpha$ is the sea ice albedo, $\tau$ is one half year, L is the latent heat of fusion, k is the thermal conductivity, while $n_w$ and $n_s$ are the optical depths for winter and summer, respectively.

During summer, $T_i = 0\,^\circ$C is assumed, while $F_{sw} = 0$ during winter. $F_{sw}$ is defined as the annual input of shortwave radiation, distributed over the melting season. The growth and melting seasons are each estimated as half of the year.

Decay melt is independent of h, while $\partial G/\partial h < 0$. Thin ice grows faster than thick ice, and as a consequence thick ice is more sensitive to changes in the forcing (Bitz and Roe 2004). Bitz and Roe (2004) showed that the change in thickness increases approximately quadratically with sea ice thickness as thicker ice needs to thin more to establish a new equilibrium.

In this study, we keep the atmospheric and solar energy fluxes constant; see Bitz and Roe (2004) and Stranne and Björk (2012) for the effect of variations in $F_{sw}$ and $F_{atm}$, respectively. Therefore, variations of the steady-state sea ice thickness are controlled by the ocean heat flux, which is a function of the steady-state flow in the ocean;

$$\begin{aligned} Q = \frac{c_p\rho _0M{\Delta } T}{A}, \end{aligned}$$

(18)

where ${\Delta } T = T_A - T, T_A$ is the Atlantic layer temperature, and T is the temperature of the surface layer. A positive heat flux is directed upward into the ice.

Following Stranne and Björk (2012), we include a thickness dependent ice export E in the sea ice model,

$$\begin{aligned} E = \frac{2\tau A_{ex}h}{A}, \end{aligned}$$

(19)

where $A_{ex}$ is the areal export. As we study changes in the sea ice thickness on time scales longer than seasonal, we can write;

$$\begin{aligned} 2\tau \frac{dh}{dt}= G-D-E. \end{aligned}$$

(20)

The equilibrium sea ice thickness can be found by solving for $G(h)-D=E(h)$. Hence, Eqs. 15–17 and 19 can be combined to one equation which can be solved for the unknown h.

When the ocean-model is coupled to the sea ice-model, the ice-export transports freshwater out of the system. Therefore,

$$\begin{aligned} F_{net} = F_{riv} - F_{ice}, \end{aligned}$$

(21)

where¹

$$\begin{aligned} F_{ice}=\frac{(G-D)A}{2\tau }. \end{aligned}$$

(22)

3 Model response to freshwater in absence of sea ice

In this section, we study the features of the model in the absence of sea ice. When the diapycnal flow decreases with increasing stratification, the steady-state flow decreases with increasing density differences (Eq. 11). As freshwater inputs at the surface increase the density difference, the flow decreases with freshwater input in the energy-constrained case (Fig. 2c). The surface layer depth also decreases with freshwater forcing (Fig. 2b). Thus, in the limit of small freshwater input, the low salinity surface layer becomes very deep, implying that the upper limit of the warmer Atlantic water layer is displaced to greater depths; a feature which may have relevance for the glacial Arctic Ocean stratification (Jakobsson et al. 2010; Cronin et al. 2012).

Fig. 2

Salinity contrast ${\Delta } S$, surface layer depth H, steady-state flow M,  and ice thickness h as a function of freshwater input $F_{riv}$ in an energy-constrained model. An Atlantic water temperature of $1\,^\circ \hbox {C}$ is used. The results are shown both in the presence and absence of sea ice, specified in legend. Full and stippled lines represent the stable and unstable solutions, respectively. The stippled line ends where the vertical density difference vanishes

Nilsson and Walin (2010) showed that the salinity-dominated circulation does not have a solution below a threshold value of $F_{riv}$. The steady-state solution becomes unstable at small salinity contrasts (Fig. 2, stippled lines). To illustrate the underlying physics, we consider the feedbacks acting on a positive salinity perturbation, which decreases the vertical density difference. As we will show in Sect. 4.2, the mean flow always stabilizes, advecting away perturbations. However, the reduced vertical density difference is associated with a positive perturbation in the diapycnal flow (Eq. 10), which increases the salt transport to the upper layer. This further reduces the vertical salinity and density contrast, resulting in a positive feedback. For sufficiently small salinity contrasts, the positive feedback becomes larger than the negative feedback, and the system becomes unstable. In the more detailed stability analysis presented in the next section, it is shown that the flow becomes unstable when $\rho _o\beta {\Delta } S < 1.5{\Delta }\rho _T$. As this limit is approached, the halocline solution breaks down and the surface layer is suggested to extend toward great depths (Fig. 2b). Additional physics need to be incorporated into the model to study this regime.

In contrast to the energy-constrained model, the constant-diffusivity model does have a solution for all density differences (Fig. 3). By applying a constant vertical diffusivity, the flow increases with density differences (Eq. 13). Since $H \propto {\Delta }\rho ^{-1/3}$, the surface layer depth decreases with freshwater forcing also in this case (Fig. 3b). From Eq. 8 and 6, we see that $M_D$ and $M_G$ increases and decreases, respectively, with a reduced H. However, because $M_G$ increases with ${\Delta }\rho$ and this response dominates that of H, the flow increases with freshwater forcing (Fig. 3c).

Fig. 3

Salinity contrast ${\Delta } S$, surface layer depth H, steady-state flow M,  and ice thickness h as a function of freshwater input $F_{riv}$ in a constant-diffusivity model. An Atlantic water temperature of $1\,^\circ \hbox {C}$ is used. The results are shown both in the presence and absence of sea ice, specified in legend. Full and stippled lines represent the stable and unstable solutions, respectively. The stippled line ends where the vertical density difference vanishes

4 Model response to freshwater in presence of sea ice

4.1 The steady-state response of sea ice thickness to changes in freshwater supply

Here, we focus on the impact of freshwater supply on the steady-state sea ice thickness. The sea ice is linked to the freshwater supply via the heat flux from the ocean to the ice. The heat flux is proportional to M (Eq. 18). In essence, a stronger steady-state heat flux results in thinner sea ice. However, the thickness-dependent sea ice export affects the net freshwater forcing of the ocean. Therefore, the sea ice is dynamically coupled to the ocean circulation.

The energy-constrained model shows increasing sea ice thickness with increasing freshwater forcing (Fig. 2d). It is also evident that the steady-state sensitivity of the sea ice thickness depends on the strength of freshwater supply. At large freshwater inputs, the sea ice thickness is rather insensitive and large changes in the freshwater forcing are needed to change the sea ice thickness. On the other hand, for small freshwater inputs, only small changes in the freshwater forcing are needed to change the sea ice thickness drastically. It is worth noting that in the uncoupled sea ice model, thick sea ice is more sensitive to a given change in the forcing than thin ice (Bitz and Roe 2004). The same basic physics still applies here, which reflects the fact that thinner ice grows faster. However, the sensitivity of the ocean heat flux becomes so strong at low freshwater supplies that it dominates over the thickness-dependent effect of sea ice growth.

Figure 2a shows that for a given external freshwater input the salinity difference is lower in the presence of sea ice, i.e. the salinity of the upper layer is higher. This is a consequence of the sea ice export, which carries away part of the external freshwater input. For the energy-constrained model which only has steady-state solutions above a threshold value of the net freshwater input, the sea ice tends to extend the range of steady-states. However, the sea ice export actually limits the steady-states with respect to $F_{riv}$ (Fig. 2). Figure 4 illustrates how the net freshwater supply and freshwater exported by sea ice vary with the external freshwater input. For large inputs of $F_{riv}$, the role of the sea ice in the freshwater budget is small. However, for small inputs of $F_{riv}$, sea ice plays a leading order role for the freshwater budget.

Fig. 4

Net freshwater input $F_{net}$, and sea ice export $F_{ice}$ as a function of $F_{riv}$ in the energy-constrained case. An Atlantic water temperature of $1\,^\circ \hbox {C}$ is used. Full and stippled lines represent the stable and unstable solutions, respectively

In the constant-diffusivity case the sea ice response is essentially reversed. Here, the sea ice thins as the freshwater supply increases (Fig. 3d). For a sufficiently large freshwater input, the sea ice will vanish. This transition to an ice-free state is gradual and continuous, and presumably requires too large freshwater inputs to be relevant in the present context. More interestingly, Fig. 3 shows that there is no longer a steady-state solution for arbitrary small freshwater inputs when sea ice is present. The underlying physics are illustrated in Fig. 5: As the freshwater supply reduces, the sea ice thickness and export increase. Thus, a critical state is approached where the net freshwater supply vanishes. However, as can be inferred from Figs. 3 and 5, the steady-state solution branch becomes unstable well before the net freshwater supply reaches zero. A formal stability analysis is performed in the next section. However, the physics is straightforward and tied to a positive feedback between the salinity difference and the ice export. If the upper layer becomes more saline (i.e. ${\Delta } S$ decreases), the flow and the heat flux diminish. This causes the sea ice thickness and export to increase, which amplifies the negative perturbation in ${\Delta } S$. This positive feedback excludes solutions below a threshold value of $F_{riv}$ and introduces an unstable steady-state for which the ice-export nearly balances the freshwater input.

Fig. 5

Net freshwater input $F_{net}$, and sea ice export $F_{ice}$ as a function of $F_{riv}$ in the constant-diffusivity model. An Atlantic water temperature of $1\,^\circ \hbox {C}$ is used. Full and stippled lines represent the stable and unstable solutions, respectively

Thus, the sea ice thickness responds differently to freshwater forcing in the two cases. The sea ice thickness increases with freshwater forcing in the energy-constrained model and decreases with freshwater forcing in the constant-diffusivity model.

4.2 The effect of sea ice on ocean circulation

Because the presence of sea ice affects the ocean, it also changes the stability of the system. The following analysis shows the stabilizing effect of sea ice on a freshwater perturbation in both the energy-constrained and constant-diffusivity case.

4.2.1 Stability analysis

Here we examine the feedbacks between ocean circulation and sea ice. The main physical features can be studied with a simplified set of equations, where we neglect the time-rate of change in the equation for H (Eq. 4), but keep it in the salinity equation (Eq. 5). Thus, we assume that the halocline depth adjusts almost instantaneously to yield $0\approx -M_G+M_D$. As a result, we can write the strength of the circulation as $M=M({\Delta } S)$ where $M({\Delta } S)$ is the steady-state version of the circulation (Eqs. 11, 13). Whereas this approximation is not formally justifiable, it is reasonable as the approximation yields the correct stability criteria in the case without sea ice interaction (Nilsson and Walin 2010).

The simplified salinity equation (Eq. 5) now becomes

$$\begin{aligned} AH \frac{d {\Delta } S}{ d t }=- {\Delta } S M + S_A F_{riv} - S_A F_{ice}. \end{aligned}$$

(23)

We also assume that perturbations to the sea ice thickness (Eq. 20) are small, and linearize Eqs. 20 and 23 around a steady-state. We obtain

$$\begin{aligned}&A \overline{H} \frac{d {\Delta } S'}{ d t }=- {\Delta } S' \overline{M}- \overline{{\Delta } S} M' - S_A F_{ice}', \end{aligned}$$

(24)

$$\begin{aligned}&2\tau \frac{dh'}{dt}=G'-D'-E', \end{aligned}$$

(25)

where the overbar and primed variables denote steady-state and perturbation quantities, respectively.

Before performing the stability analysis, we provide a qualitative understanding of the three different terms on the right hand side of Eq. 24. The first term on the right-hand side represents a negative feedback due to the mean flow, acting to attenuate salinity perturbations.

The feedbacks from the remaining two terms are determined by the dependence of the flow and the ice export on the salinity perturbation. The flow perturbation is related to the salinity perturbation as

$$\begin{aligned} M'= \left( \frac{\partial \overline{M}}{\partial \overline{{\Delta } S}} \right) {\Delta } S'. \end{aligned}$$

(26)

This makes the second term on the right-hand side of Eq. 24 a positive feedback for the energy-constrained case where $\left( \frac{\partial \overline{{M}}}{\partial \overline{{\Delta } S}}\right) <0$. On the other hand, for the constant-diffusivity case, the feedback becomes negative as $\left( \frac{\partial \overline{{M}}}{\partial \overline{{\Delta } S}}\right) >0.$

We now consider the feedback from the last term in Eq. 24; the sea ice export dependence on the salinity perturbation. Equation 22 relates $F_{ice}$ and $G-D$. The growth and decay depends indirectly on the salinity via the ocean heat flux which is proportional to M (Eq. 18). Therefore,

$$\begin{aligned} F_{ice}'= -\frac{c_p\rho _0{\Delta } T}{L}M'=-\frac{c_p\rho _0{\Delta } T}{L}\left( \frac{\partial \overline{M}}{\partial \overline{{\Delta } S}} \right) {\Delta } S'. \end{aligned}$$

(27)

Equation 27 shows that a positive flow perturbation, associated with increased ocean heat flux, yields negative perturbations in the ice export. We thus find that the sea ice export results in a positive feedback in the constant-diffusivity model where a positive perturbation in ${\Delta } S$ is associated with a positive perturbation in the flow. This reduces the sea ice export, which increases the net amount of freshwater retained in the surface layer. The reverse is true for the energy-constrained model where the sea ice export results in a stabilizing feedback.

As we now have a qualitatively feel for the different responses of Eq. 24 to salinity perturbations, we move on to the proper linear stability analysis to find out whether the equilibrium solutions are stable or not. The equations governing small perturbations in sea ice thickness and salinity differences, Eqs. 24 and 25, can be written

$$\begin{aligned} \frac{d}{dt}\left( \genfrac{}{}{0.0pt}{}{H'}{{\Delta } S'}\right) = \left( \genfrac{}{}{0.0pt}{}{a \ \ b}{c \ \ d}\right) \left( \genfrac{}{}{0.0pt}{}{H'}{{\Delta } S'}\right) \end{aligned}$$

(28)

where a, b, c,  and d are the coefficients of the stability matrix $A_s$, given by the partial derivatives of the right hand side of Eqs. 24 and 25. To perform a stability analysis which is valid for both mixing representations considered in this study, we write

$$\begin{aligned} \overline{M}=c\overline{{\Delta } \rho }^{\gamma }, \end{aligned}$$

(29)

where c is a constant, $\gamma = -1/3$ for the energy-constrained case, and $\gamma = 1/3$ for the constant-diffusivity case. This gives,

$$\begin{aligned} a= & {} \left( \frac{\partial (G-E)}{\partial h}\right) \frac{1}{2\tau }, \ \ \ \ b=-\frac{c_p{\Delta } T\rho _0}{L}\frac{\overline{M}\gamma \beta \rho _0}{A\overline{{\Delta }\rho }},\nonumber \\ c= & {} -\left( \frac{\partial G}{\partial h}\right) \frac{S_A}{2\tau \overline{H}}, \ \ \ d=-\frac{\overline{M}}{A\overline{H}}\left[ 1 + \frac{\gamma \beta \rho _0}{\overline{{\Delta }\rho }}\left( \overline{{\Delta } S} - S_A\frac{c_p{\Delta } T\rho _0}{L} \right) \right], \end{aligned}$$

(30)

where a results from $\partial D/\partial h =0$ (Eq. 16), and b results from $\partial E/\partial Q =0$ (Eq. 19). c is found given that only the last term in Eq. 24 is a function of h, and d is found making use of Eqs. 26 and 27.

An equilibrium solution is stable provided that the real parts of the eigenvalues of the stability matrix are negative. We are looking for solutions on the form exp(xt). Therefore, the eigenvalues can be found by solving det($A_s-Ix)=0$. This gives us

$$\begin{aligned} x_{1,2}=\frac{a+d}{2} \pm \left( \left( \frac{a+d}{2}\right) ^2 + bc - ad\right) ^{1/2}. \end{aligned}$$

(31)

The eigenvalues are negative and an equilibrium solution is stable provided that $bc-ad<0$ and $a+d<0$. The first constraint, $bc-ad<0$, gives

$$\begin{aligned} \overline{{\Delta } S} > \frac{{\Delta }\rho _T}{\beta \rho _0(1+\gamma )} + \frac{\gamma S_A\lambda }{1+\gamma }, \end{aligned}$$

(32)

where

$$\begin{aligned} \lambda \equiv \left\{ \begin{array}{l l}-\frac{c_p {\Delta } T\rho _0}{L} \left( \frac{\partial E}{ \partial h } \right) \left( \frac{\partial (G-E)}{ \partial h } \right) ^{-1} &{} \text {for} \quad h>0 \\ 0 &{} \text {for} \quad h=0. \end{array} \right. \end{aligned}$$

(33)

The second constraint, $a+d<0$, gives

$$\begin{aligned} \overline{{\Delta } S} > \frac{{\Delta }\rho _T(1+\phi )}{\beta \rho _0(1+\gamma +\phi) } + \frac{\gamma S_A}{1+\gamma + \phi }\frac{c_p{\Delta } T\rho _0}{L}, \end{aligned}$$

(34)

where

$$\begin{aligned} \phi = -\frac{A\overline{H}}{\overline{M}2\tau }\left( \frac{\partial (G-E)}{\partial h}\right). \end{aligned}$$

(35)

Note that $\frac{\partial (G-E)}{\partial h}<0,$ and that the stability criteria are not analytical solutions as $\lambda =\lambda ({\Delta } S)$ and $\phi =\phi ({\Delta } S)$.

The parameter $\lambda$ controls the sea ice feedback in this model. It depends mainly on the ice thickness, which is indirectly linked to the ocean state which is determined mainly by the salinity contrast (Fig. 6). The parameter $\lambda$ attains a maximum value of $c_p {\Delta } T\rho _{0}/ L$ when $\left( \frac{\partial \overline{E}}{ \partial \overline{h} } \right) \gg \left( \frac{\partial \overline{G}}{ \partial \overline{ h} } \right)$. For $T_A=1^\circ$C, that value is 0.0372. Due to the thermodynamics of sea ice, the rate of sea ice growth decreases for thick ice, and $\lambda$ approaches its maximum for very thick sea ice. If sea ice export is insensitive to sea ice thickness, $\lambda$ is zero.

Fig. 6

The parameter $\lambda S_A$ as a function of the salinity contrast ${\Delta } S$ for the a energy-constrained and b constant-diffusivity model

The parameter $\phi$ is essentially the ratio of the time scales for the response of the salinity and the sea ice. In the limit of fast sea ice response, $\phi$ approaches infinity, and the second stability criterion (Eq. 34) yields $\overline{{\Delta } S}>{\Delta }\rho _T/\beta \rho _0$ which is always satisfied. In the limit of slow sea ice response, $\phi =0$. In this case, the second stability criterion is stricter when $\gamma >0$. The first criterion applies for $\gamma <0$. For both representations of the diapycnal flow considered in this study, the first criterion determines where the solution becomes unstable.

For the energy-constrained case where $\gamma <0$, the sea ice feedback from the last term in Eq. 32 is negative, confirming sea ice as a negative feedback. As the presence of sea ice stabilizes the circulation, while flow perturbations destabilize, it is interesting to see whether the presence of sea ice compensates for the destabilizing effects in the energy-constrained model. We get the stability criterion for the energy-constrained model $\overline{{\Delta } S_E}$ by inserting $\gamma = -1/3$ in Eq. 32;

$$\begin{aligned} \overline{{\Delta } S_E} > \frac{3{\Delta }\rho _T}{2\rho _0\beta } - \frac{S_A\lambda }{2}. \end{aligned}$$

(36)

The critical salinity contrast where the vertical density difference vanishes $\overline{{\Delta } S_T}$ is given by $\overline{{\Delta } S_T} = {\Delta }\rho _T/(\beta \rho _o).$ For $\overline{{\Delta } S} > S_A\lambda$ the critical salinity $\overline{{\Delta } S_E}$ is reached before the vertical density difference disappears. As a consequence, the equilibrium becomes unstable in the energy-constrained model before the vertical density difference disappears for high Atlantic water temperatures and thin sea ice. This is always true for the parameters used in this study (Fig. 7a).

Fig. 7

Stability criteria for the a energy-constrained and b constant-diffusivity model for given Atlantic water temperatures $T_A$ and salinity contrasts ${\Delta } S$. Shaded areas represent the stable parts of the system. Dotted, black lines mark the vertical density difference ${\Delta }\rho$

For the constant-diffusivity model where $\gamma >0$, Eq. 32 shows that the sea ice feedback is positive; the presence of sea ice increases the salinity difference where the solution becomes unstable. The stability criterion for the constant-diffusivity model $\overline{{\Delta } S_S}$ becomes;

$$\begin{aligned} \overline{{\Delta } S_S} = \frac{3{\Delta }\rho _T}{4\rho _0\beta } + \frac{S_A\lambda }{4}. \end{aligned}$$

(37)

For $\overline{{\Delta } S} < S_A\lambda$ the salinity contrast $\overline{{\Delta } S_S}$ is reached before $\overline{{\Delta } S_T}$. Hence, the equilibrium solution becomes unstable before the vertical density difference disappears. Concluding, the constant-diffusivity model becomes unstable before the vertical density difference disappears for low Atlantic water temperatures, thick sea ice, or large ice exports. As h decreases with $F_{riv}$, this occurs at small freshwater supplies (Fig. 7b). Given the parameters in Table 1, and a freshwater input of 0.5 Sv, the system becomes unstable before the vertical density difference disappears for Atlantic water temperatures below $2\,^\circ \hbox {C}$.

5 Quasi-steady response to increasing Atlantic water temperatures

We have now established how the system responds to freshwater perturbations in the presence of sea ice. As the hypothesis from Dokken et al. (2013) involves a destabilization of the vertical stratification through increasing subsurface temperatures, we here investigate how increasing the Atlantic water temperature affects the steady-state solutions. Increasing Atlantic water temperatures decreases the vertical density difference and increases the heat flux from the ocean to the ice. As these two variables are indirectly coupled, the response of the steady-state solution is not straightforward, and will be examined in the following sections.

5.1 Energy-constrained case

First, we study the dynamic response of the steady-state flow to increasing $T_A$. We recall that M intensifies with decreasing density differences in the energy-constrained case (Eq. 11). The stronger ocean flow advects more saline water into the surface layer, thereby reducing the salinity contrast (Fig. 8, inner loop). In addition, the value of the critical salinity contrast increases with Atlantic water temperatures (Eq. 36). The dynamic response to increasing Atlantic water temperatures is therefore to decrease the stability of the system by driving the system toward the critical salinity contrast where the steady-state solution becomes unstable.

Fig. 8

The isolated effect of Atlantic water temperatures $T_A$ on salinity differences ${\Delta } S$. The thermodynamic response (outer loop) is valid for both models, while the dynamic response is only valid for the energy-constrained model (inner loop). ${\Delta }\rho$ is the density difference, M the steady-state flow, Q the heat flux from the ocean to the ice, h the sea ice thickness and $F_{ice}$ the sea ice export.

Figure 9a shows instead that the salinity difference predominantly increases with $T_A$ for a given freshwater forcing. This is linked to the thermodynamic sea ice response and the following decrease in sea ice thickness. As the heat flux from the ocean to the ice is proportional to the temperature contrast between the two ocean layers, ${\Delta } T = T_A - T$, and T is taken to be fixed, ${\Delta } T$ increases linearly with $T_A$. Therefore, the sea ice thins when the Atlantic water temperature increases (Fig. 9d). The following melting of sea ice (or decrease in sea ice export) freshens the upper layer, thereby enhancing ${\Delta } S$ (Fig. 8, outer loop).

Fig. 9

Salinity contrast ${\Delta } S$, surface layer depth H, steady-state flow M,  and ice thickness h as a function of Atlantic water temperature $T_A$ in an energy-constrained model. $F_{riv}$ is specified in legend (Sv). The salinity contrast where the density difference disappears is shown by the grey line. Full and stippled lines represent the stable and unstable solutions, respectively

When sea ice is present, the dynamic and thermodynamic effects of increasing $T_A$ have opposing effects on the density contrast. The sea ice melt buffers the system’s response to increasing $T_A$. The result is that ${\Delta }\rho$ is nearly constant, and as a consequence, the steady-state flow is close to invariant to changes in $T_A$ (Fig. 9c). Therefore, the heat flux from the ocean to the ice is mainly determined by ${\Delta } T$ in the energy-constrained case. However, when the sea ice thins sufficiently, the effect of melting is too small to balance the increase in $T_A$ with respect to the density contrast. At that point, ${\Delta }\rho$ decreases with $T_A$. Therefore, as h approaches zero, M increases sharply with $T_A$. As a consequence, the salinity contrast decreases with Atlantic water temperatures when h is small (Fig. 9a, d). At this point, high Atlantic water temperatures destabilize the system. The sea ice feedback, which is controlled by $\lambda$ (Eq. 33), is not strong enough to balance the destabilizing effect of the flow.

Due to the destabilizing effect of the flow, the steady-state solution becomes unstable at high Atlantic water temperatures (Fig. 9). For low freshwater inputs, this occurs at salinity contrasts where sea ice still is present. As a consequence, the sea ice abruptly disappears (Fig. 9d). For higher freshwater inputs, the steady-state solution turns unstable as the sea ice disappears due to the removal of the stabilizing effect of the sea ice (not shown). For all values, the system turns unstable before the vertical density difference disappears (Fig. 7).

5.2 Constant-diffusivity case

The thermodynamic response of sea ice thickness to increasing Atlantic water temperatures is independent of the mixing representation. The increase in the temperature contrast increases the ocean heat flux, and hence decrease the sea ice thickness in the constant-diffusivity model as well (Fig. 8, outer loop). The consequent melting from the sea ice (or decrease in ice-export) leads to an increase in the salinity contrast.

On the other hand, as the flow slows down with increasing $T_A$ in the constant-diffusivity case (Eq. 13), the dynamic response differs from the response of the energy-constrained case. A weaker flow enhances ${\Delta } S$, thereby stabilizing the system. However, at low temperatures and freshwater inputs, the steady-state flow actually speeds up with increasing Atlantic water temperatures (Fig. 10c). The reason why is related to the thick sea ice which is more sensitive to changes in the forcing than thin ice. The sea ice decreases rapidly with increasing heat flux. Therefore, the strong sea ice melt leads to an increase in the salinity contrast which dominates the changes in the temperature contrast. As a consequence, the density difference increases and hence the steady-state flow speeds up. However, as the melting effect exceeds the anomalous advective salt transport, the positive flow perturbation freshens the upper layer. Therefore, this intermediate increase in flow with $T_A$ does not change the stabilizing effect of the Atlantic temperatures; the salinity contrast still increases with Atlantic water temperatures (Fig. 10a).

Fig. 10

Salinity contrast ${\Delta } S$, surface layer depth H, steady-state flow M,  and ice thickness h as a function of Atlantic water temperature $T_A$ in a constant-diffusivity model. $F_{riv}$ is specified in legend (Sv). The salinity contrast where the density difference disappears is shown by the grey line. Full and stippled lines represent the stable and unstable solutions, respectively

In the constant-diffusivity case, the sea ice thickness response depends on the competing effects of the dynamic and thermodynamic responses to increasing $T_A$. The heat flux is proportional to ${\Delta } TM$, and the two terms have opposing effects on the heat flux as the steady-state flow mainly decreases with $T_A$ (Fig. 10c). Interestingly, we get a local minimum in sea ice thickness with increasing $T_A$ (Fig. 10d). This occurs when the decrease in M balances the increase in ${\Delta } T$, i.e., when the dynamic and thermodynamic effects are equal. When the decrease in the steady-state flow is large enough to dominate changes in the temperature contrast, the heat flux decreases with increasing $T_A$. The result is an increase in h with increasing Atlantic water temperatures at large $T_A$ (Fig. 10d).

The change in the response of the sea ice to increasing Atlantic water temperatures changes the stability properties of the system. At low freshwater inputs the steady-state solution becomes unstable at low Atlantic water temperatures (Fig. 10, blue line, and Fig. 7). This is due to the thick sea ice and large ice-export which removes most of the freshwater introduced to the surface layer. By increasing $T_A$, steady-state solutions are introduced as the sea ice melt stabilizes the system. However, at higher Atlantic water temperatures, where sea ice thickens with $T_A$, the sea ice growth acts destabilizing. At high Atlantic water temperatures, the ice export dominates the freshwater input and flow perturbations, and the system becomes unstable. Note that the salinity contrast still increases at high temperatures due to a weaker flow (Fig. 10a). However, due to thicker sea ice and a larger $\lambda$, so does the salinity contrast where the steady-state solution becomes unstable (Eq. 37). Interestingly, the system becomes unstable at both high and low Atlantic water temperatures for low freshwater inputs (Fig. 10, blue line). For high freshwater inputs, the system stays stable until the vertical density difference disappears (Fig. 7).

We have shown that the system can turn unstable at high Atlantic water temperatures for both mixing representations, depending on the freshwater input. What happens at this point goes beyond the physics contained in the present version of our model. Still, it is fair to assume that the sea ice suddenly disappears as a new regime is established.

6 Discussion

The dynamical effect of sea ice on a freshwater perturbation is seen to depend on the representation of vertical mixing. In a system where the diapycnal flow increases with density differences (constant-diffusivity case), sea ice destabilizes against a freshwater perturbation (Fig. 11, red part) and introduces unstable solutions. If the diapycnal flow decreases with density differences (energy-constrained case), the presence of sea ice stabilizes the system (Fig. 11, blue part). In this case, the presence of sea ice extends the range of stable steady-state solutions. However, the stabilizing effect is not strong enough to allow for stable solutions for arbitrarily small freshwater inputs.

Fig. 11

The stabilizing effect of sea ice on a freshwater perturbation. Left/red part of the figure represents the constant-diffusivity case, while the right/blue part represents the energy-constrained case. $F_{riv}$ is the freshwater input at the surface, ${\Delta } S$ the salinity contrast, M the steady-state flow, Q the heat flux from the ocean to the ice, h the sea ice thickness and E the sea ice export

As the salinity-dominated circulation state abruptly becomes unstable for threshold values of freshwater input and deep-ocean temperature, the model suggests that the halocline and sea ice cover in the Nordic Seas could break down abruptly in response to a slow gradual change of the climatic conditions. The systematic stratification changes between cold stadials and warm interstadials as observed in the Nordic Seas by e.g. Dokken and Jansen (1999), Rasmussen and Thomsen (2004) and Dokken et al. (2013) can be explained by a reduction in the freshwater supply to the Nordic Seas. This is highly relevant for a cold, glacial climate with a weaker hydrological cycle. Generally, we expect reduced input of freshwater to the Nordic Seas during cold, stadial conditions. In this case, only a small change in the freshwater supply is needed to initiate large and abrupt changes in both sea ice and ocean stratification (Fig. 2). This can be contrasted to the traditional freshwater hosing experiments where a large additional freshwater forcing of about 1 Sv to the North Atlantic is applied to maintain cold stadial conditions with extensive sea ice cover (e.g. Manabe and Stouffer 1995; Stouffer et al. 2006). In this case, an abrupt termination of the freshwater forcing is required to trigger a transition to a warm interstadial with weak stratification and reduced sea ice (Liu et al. 2009).

The abrupt changes in sea ice and ocean stratification at small freshwater supplies can occur independently of a change in subsurface temperature. However, the simple model also transitions to an unstable solution with a removal of sea ice when the temperature of the subsurface Atlantic water is increased (Fig. 9d). This is similar to the hypothesized transition from cold stadial states to warm interstadials triggered by warm Atlantic water as in Dokken et al. (2013). Interestingly, the unstable state is reached before the vertical density difference disappears. This is mainly true for the energy-constrained case. For the constant-diffusivity case, the stable solutions disappear for both high and low Atlantic water temperatures when the freshwater supply is low. At higher freshwater supplies, the system stays stable until the vertical density difference disappears. However, we recall that the constant-diffusivity model may be a less realistic representation of the system and is not thought of as a good analogue for the Nordic Seas/Arctic.

For the energy-constrained case, at low freshwater inputs, only a small increase in temperature is needed to destabilize the equilibrium state, and as a consequence, sea ice abruptly disappears. As both increasing Atlantic water temperatures and decreasing freshwater inputs promote an unstable system, the removal of the halocline as proposed by Dokken et al. (2013), can occur with smaller changes in Atlantic water temperatures than previously thought.

Our results are highly dependent on the use of a thickness dependent sea ice export. If we instead assume that the ice-export is constant and independent of the ice thickness (not shown), then the sea ice related feedbacks on the flow and the salinity stratification vanish. The equilibrium states become unstable at the same salinity contrast as when sea ice is absent. However, as the ice-export reduces the salinity contrast for a given freshwater input, the system with stability-dependent mixing still becomes unstable at a larger freshwater input in the presence of sea ice.

As the present conceptual model is highly simplified and essentially one-dimensional, it neglects several processes that could be of importance for the dynamics of the halocline and sea ice. In particular, the model does not represent effects of wind forcing on either the ocean circulation or the sea ice export. Possible impacts of ice mechanics are also ignored in the parametrization of the sea ice export, which simply is taken as proportional to the sea ice thickness. Note further that even though the sea ice model crudely represents an annual growth and melt cycle, effects of seasonal changes in the freshwater storage due to growth and melt of sea ice are not accounted for in the model. Adding representations of seasonality and natural variability can obviously change the stability range of the steady-state solutions examined here. Reductions in air-sea momentum exchange for high sea ice concentrations could also affect the stability properties of the system. In the presence of sea ice, the energy input from the atmosphere to oceanic mixing via e.g. internal wave generation should diminish (Rainville and Woodgate 2009). Thus, the energy supply to the vertical mixing $\epsilon$ would hence decrease. In the energy-constrained model this would introduce a feedback between the strength of the flow and the sea ice concentration, with ensuing impacts on the stability of the system.

The results from our simple model suggest a mechanism for how a reduction in freshwater input to the Nordic Seas could initiate large changes in the sea ice cover. Due to the instability at small freshwater supplies, small changes in the freshwater input or increases in subsurface temperatures could terminate the salinity-dominated mode with sea ice, and allow for a less stratified regime without sea ice. This could aid in explaining why the sea ice and hydrography of the Nordic Seas were highly variable during the last glacial cycle. Our simple model only accounts for the retreat of the sea ice cover, and does not provide any mechanism for the re-appearance of sea ice in the Nordic Seas.

7 Conclusion

The main results from the examination of the conceptual model of sea ice-ocean-circulation feedbacks are:

• The presence of sea ice stabilizes against a freshwater perturbation when the vertical velocity is represented with a constant energy-supply to the diffusivity.

• The presence of sea ice destabilizes against a freshwater perturbation when the vertical velocity is represented with a constant diffusivity.

• For sufficiently weak freshwater supply and irrespective of the representation of the vertical mixing, the salinity-dominated circulation becomes unstable.

• The sea ice is highly sensitive to changes in subsurface Atlantic water temperatures.

• The results from the simple conceptual model suggest that during cold glacial conditions, with reduced input of freshwater to the Nordic Seas, relatively small changes in freshwater or Atlantic water temperature could have triggered abrupt transitions in sea ice cover.