Amplified cooling of Snowball Earth from a salt-albedo feedback

It is believed that the atmospheric circulation on Snowball Earth produced a net ablation zone exposing bare sea ice. Under sufficiently low temperatures, salt begins to precipitate out of sea ice, forming a lag deposit of crystals with high albedo as the ice sublimates. This could have resulted in a salt-albedo feedback that has not previously been included in modeling studies of Snowball Earth. We implement a salt-albedo feedback in a simple climate model and show that, once initiated, this mechanism could have intensified global cooling in the initial phase of Snowball Earth. Our results suggest that salt precipitation may have played a role in shaping the early climate of Snowball Earth.

1 Introduction

Modeling studies indicate a net ablation zone at lower latitudes on Snowball Earth, where annual mean sublimation exceeds precipitation (Pierrehumbert, 2005; Pierrehumbert et al., 2011; Abbot et al., 2011, 2013). This would have exposed sublimating bare sea ice under extremely cold temperatures for thousands of years during the global glaciation events of the late Neoproterozoic era (715–635 Ma). At sufficiently low temperatures, the salt in seawater becomes saturated and precipitates in crystal form (Light et al., 2009; Carns et al., 2015). In the net ablation zone on Snowball Earth, precipitated crystals at the surface could have accumulated a lag deposit as the ice sublimated (Light et al., 2009, 2016; Carns et al., 2015, 2016). Field studies (Carns et al., 2015) and laboratory experiments (Light et al., 2016; Carns et al., 2016) have shown that salt crystals have a high albedo, which can result in a salt-albedo feedback mechanism (Light et al., 2016; Shields and Carns, 2018). The effects of this salt-albedo feedback have not previously been implemented in a climate model for Snowball Earth.

2 Methods

We implemented a simplified salt-albedo feedback parameterization in a diffusive one-dimensional energy balance model (EBM) and solved it using analytical methods. The EBM describes the zonal mean surface temperature of the Earth, and its key features are meridional heat distribution through diffusion, the Planck response and a nonlinear albedo feedback. This is a well-known model that has been extensively studied in the literature (e.g. North, 1975; North et al., 1981). Using spherical coordinates, the model may be formulated as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }C{\displaystyle \frac{\partial T}{\partial t}}-D{\displaystyle \frac{\mathrm{1}}{\sin x}}{\displaystyle \frac{\partial }{\partial x}}\left(\sin x{\displaystyle \frac{\partial T}{\partial x}}\right)+BT\\ {\displaystyle }& {\displaystyle }=Qs\left(x\right)(\mathrm{1}-a(x,T\left)\right)-A\end{array}\end{array}$$

where x is the co-latitude, ranging from x=0 at the North Pole to $x=\mathit{\pi }/\mathrm{2}$ at the Equator. Here C is the heat capacity of the atmosphere–ocean column, Q is solar irradiance and the diffusion coefficient D scales the meridional heat transport. Solar radiation is distributed across latitudes according to the function s(x), and energy is radiated into outer space through the parametrization A+BT (Budyko, 1969). Additionally, the solution is subject to no-flux boundary conditions (North, 1975).

We adapted the albedo function a(x,T(x)) to include both an ice-albedo and a salt-albedo feedback mechanism. Following the standard approach for EBMs, we let latitudes with a zonal mean temperature below a specific isotherm T_s have an ice/snow cover albedo value denoted α₂. This will be separate from the albedo value of ice-free latitudes denoted α₁. The extent of the net ablation zone on Snowball Earth is denoted by the latitude r₁. We let the snow/ice albedo take on a lower value α₃ equatorward of r₁ to reflect that of bare sea ice. This dynamic, known as the Jormungand mechanism, has been shown to produce Waterbelt states, i.e., a Snowball Earth with a tropical water belt, in several climate models (Abbot et al., 2011; Braun et al., 2022; Walsh, 2018; Samuelsberg and Jakobsen, 2025b). However, it is uncertain whether the Snowball Earth events could have been a Waterbelt state (Hoffman et al., 2017; Braun et al., 2022). We let the surface take a salt albedo value α₄ for temperatures below some eutectic temperature T_e in the bare sea ice region. The eutectic temperature refers to the temperature at which the salt in seawater becomes solid and starts to precipitate, and is about −36 °C for the entire seawater system (Marion et al., 1999; Light et al., 2016; Vancoppenolle et al., 2019). A lag deposit of precipitated salt crystals is expected to form at the surface when the sea ice sublimates (Light et al., 2016). This dynamic only operates in the net ablation zone. We initially set the eutectic temperature to $T_{\mathrm{e}}=-\mathrm{36}$ °C. The albedo of hydrohalite salt crystals is extremely high at α₄=0.93 (Light et al., 2016), and we used this value throughout this study. Other model parameters are given in Table 1. Mathematically, the proposed albedo function takes the form:

$$\begin{array}{}\text{(2)}& a(x,T)=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{1}},& T>T_{\mathrm{s}}\\ & \text{(ice-free latitudes)}\\ {\mathit{\alpha }}_{\mathrm{2}},& T\le T_{\mathrm{s}},x<r_{\mathrm{1}}\\ & \text{(ice/snow cover; poleward of}{r}_{\mathrm{1}})\\ {\mathit{\alpha }}_{\mathrm{3}},& T_{\mathrm{e}}<T\le T_{\mathrm{s}},x\ge r_{\mathrm{1}}\\ & \text{(bare sea ice; equatorward of}{r}_{\mathrm{1}})\\ {\mathit{\alpha }}_{\mathrm{4}},& T\le T_{\mathrm{e}},x\ge r_{\mathrm{1}}\\ & \text{(salt deposit)}\end{array}\right.\end{array}$$

Here r₁ refers to the latitude in spherical coordinates corresponding to the geographic latitude (in degrees north) listed in Table 1.

Table 1 Model parameters in the EBM. A₀, B, Q and s(x) are from Abbot et al. (2011). The critical isotherm for ice formation T_s and the heat capacity C are from North et al. (1981). D is take from Pierrehumbert et al. (2011). The eutectic temperature T_e and the albedo of salt crystals α₄ are from Light et al. (2016). The albedo over ice-free ocean α₁ is from Braun et al. (2022), the snow albedo value α₂ is from Abbot et al. (2011) and α₃ is a plausible bare sea ice albedo based on the measurements in Dadic et al. (2013).

Equilibrium solutions to Eq. (1), with the albedo in Eq. (2), were found using boundary integrals. Here we briefly outline the main procedures of the method and how it is adapted to handle the salt-albedo feedback parametrization.

We non-dimensionalized Eq. (1) and introduced the dimensionless, equilibrium temperature field $u\left(x\right)=\mathrm{1}/T_{\mathrm{s}}T\left(x\right)$, where T_s is the ice line temperature. In Samuelsberg and Jakobsen (2025a), it was shown that u(x) satisfy the following relation on some interval (x₁,x₂) in the domain:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }u\left(x\right)={\int }_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}\mathrm{d}s\sin sg(s,x)f(s,u)\\ {\displaystyle }& {\displaystyle }+\left\{g(s,x)\sin s{\displaystyle \frac{\partial u}{\partial x}}\left(s\right)-u\left(s\right)\sin s{\displaystyle \frac{\partial g}{\partial x}}(s,x)\right\}{|}_{s=x_{\mathrm{1}}}^{s=x_{\mathrm{2}}},\end{array}\end{array}$$

where g(x,s) is a Green's function associated with the problem and

$$\begin{array}{}\text{(4)}& f(x,u)={\displaystyle \frac{Q}{T_{\mathrm{s}}D}}s\left(x\right)(\mathrm{1}-a(x,u\left(x\right)\left)\right)-{\displaystyle \frac{A}{T_{\mathrm{s}}D}}.\end{array}$$

Here the albedo function a(x,u(x)) takes dimensionless temperatures and the input is scaled to give the same albedo response as the function in Eq. (2). We used the Green's function from Samuelsberg and Jakobsen (2025a). The primary objective for any method for solving diffusive, one-dimensional EBMs with a step-function albedo response at the ice line, is to locate the ice line. In the model, the defining property of the ice line latitude, denoted x_s, is that T(x_s)=T_s, and it represents the transition between snow/ice and open ocean. We let the salt-albedo feedback take on a similar temperature response associated with the eutectic temperature T_e. The salt-albedo feedback can therefore introduce a salt line, which is the transition latitude between salt cover and bare sea ice within the bare sea ice region equatorward of r₁. We let x_e denote the salt line latitude. Mathematically, the defining property of the salt line is that T(x_e)=T_e. We located the ice line and the salt line by developing a set of boundary integral equations. Details are left out here, but the approach was the same as in Samuelsberg and Jakobsen (2025a) and Samuelsberg and Jakobsen (2025b).

To understand the effects of the salt-albedo feedback mechanism on the system, we drew bifurcation diagrams for the equilibrium solutions to Eq. (1). Figure 1 shows a bifurcation diagram obtained using the model parameters from Table 1. Following Abbot et al. (2011), we let the bifurcation parameter be A, which controls the outgoing longwave radiation. The x axis in Fig. 1 show a deviation in the bifurcation parameter from the reference value A₀ given in Table 1, such that radiative forcing $\mathrm{\Delta }A=A_{\mathrm{0}}-A$. This represents forcing in the model as a reduction in A results in lower outgoing energy, e.g. due to increased atmospheric greenhouse gas concentrations. We chose to display the equatorial temperature on the y axis of the bifurcation diagram. In the literature on EBMs, bifurcation diagrams typically show ice line latitude on the y axis. However, in this study, we focused on alternative Snowball Earth states, and all states below the point b₂ in the bifurcation diagram in Fig. 1 are fully ice-covered. To differentiate these states, we represented them using their equatorial temperature. Note that this projection partially obscures the small ice cap instability, and the branch labeled “High-latitude ice line” in Fig. 1 is not directly connected to the “Hothouse”-branch. Stability was determined using the slope stability theorem (Cahalan and North, 1979), which applies to the classical two-ice-line solutions of Eq. (1), and through simulations performed with the finite-difference code from Samuelsberg and Jakobsen (2025b). The transitions shown in Fig. 1 were verified using finite-difference simulations.

3 Results

The salt-albedo feedback gave rise to new climatic states not previously seen in modeling studies (Fig. 1). It is well established that the EBM for an idealized aquaplanet has three stable states: the ice-free hothouse state, the high-latitude ice line state and the Snowball Earth state (North et al., 1981). We prescribe a lower bare sea ice albedo in the net ablation zone, which gives the EBM Waterbelt solutions as well (Abbot et al., 2011; Samuelsberg and Jakobsen, 2025b). The salt-albedo feedback most notably gives rise to two new, stable, coexisting Snowball Earth states: a Snowball Earth with a salt deposit (<b₆ in Fig. 1) and a Snowball Earth without a salt deposit (b₂−b₃). The latter is characterized by global ice cover and equatorial bare sea ice. The salt-deposit Snowball Earth state is likewise fully ice covered, but with a salt deposit on the bare sea ice. In addition, several unstable states appeared, with two smaller salt deposits on the bare sea ice (b₃−b₄ and b₅−b₆ in Fig. 1) or with more than two deposits (b₄−b₅). We will not focus on these states here.

Figure 1 Bifurcation diagram for the EBM with a salt-albedo feedback (Eq. 1 with the albedo function in Eq. 2). Radiative forcing on the x axis is defined as $\mathrm{\Delta }A=A_{\mathrm{0}}-A$, where A₀ is the reference value given in Table 1. Black solid lines indicate stable equilibrium states and dashed lines indicate unstable states. Stable climate states are labeled. The branch terminating at b₆ extends beyond the plotted frame. Likewise, the branch originating in b₁ extends beyond the frame. Gray solid arrows show transitions in the Snowball Earth hysteresis. The dashed arrow indicates the unfeasible transition from a warm climate to a Snowball Earth without salt deposits. Model parameters are given in Table 1. The gray line between point b₄ and b₅ represents a truncation where Snowball Earth states with multiple salt deposits are omitted, the situation here is analogous to the complicated ice distributions reported in Samuelsberg and Jakobsen (2025b). These solutions are unstable.

The presence of two new Snowball Earth states raises the question of which state is more likely to represent the climate of the Snowball Earth events in the Late Neoproterozoic. To assess this, we adopted an approach previously applied to assess Waterbelt states as an explanation to Neoproterozoic geology (Abbot et al., 2011; Braun et al., 2022). Namely, that the state had to be accessible from warm climates (Fig. 2) and had to exist over a large range of atmospheric CO₂ concentrations. This is to determine which of the two new states is consistent with the well-established hysteresis loop associated with the Snowball Earth events (Hoffman et al., 1998, 2017). That is, from a warm climate to Snowball Earth and back again following a substantial increase in atmospheric CO₂ concentrations. The peculiar aspect of the Neoproterozoic Snowball Earth events is the lack of evidence for ice sheets in 1.5 billion years leading up to these global glaciations, known as the Proterozoic glacial gap (Hoffman et al., 2017). Thus, the climate system must have rested in some warm, more or less ice-free equilibrium state prior to these events. In the EBM, such states exist in the upper branch of the bifurcation diagram in Fig. 1 and the bifurcation point b₁ is associated with the onset of Snowball Earth. For the model parameters used in Fig. 1, the Snowball Earth without a salt deposit is not accessible because a parameter drift below b₁ drives the system into the salt-deposit state. Therefore, the transition from a warm climate to a Snowball Earth without a salt deposit is not possible here.

Figure 2 Diagram showing the accessibility of the Snowball Earth without a salt deposit for different values of the critical temperatures T_s and T_e. The lines represent the boundary between accessible and inaccessible regimes of the parameter space for different values of the bare sea ice extent r₁=20° N, r₁=15° N and r₁=10° N. Above the lines the Snowball Earth without a salt deposit is inaccessible from warm climates. Other model parameters are given in Table 1.

The Snowball Earth without a salt deposit is inaccessible for plausible model parameters. The accessibility of the Snowball Earth state without a salt deposit depends most prominently on the eutectic temperature T_e and the ice formation temperature T_s. Figure 3 shows a selection of bifurcation curves with different model parameters. A transition into the salt-deposit Snowball Earth from a warm, ice-free state is favored by higher diffusion rates (Fig. 3a) and larger albedo differences α₂−α₁ (Fig. 3c), whereas lower diffusions and albedo differences prevent a direct transition (pink curves in Fig. 3a and c). Assuming a direct transition to Snowball Earth, the parameters controlling whether such a transition sends the system into a Snowball Earth with or without a salt deposit are T_e, T_s and r₁ (see Fig. 3b and d). The curves in the (T_e,T_s)-plane in Fig. 2 mark the boundary between the accessible and inaccessible regimes of the parameter space for r₁=20° N, r₁=15° N and r₁=10° N. The extent of the net ablation zone, here denoted by the latitude r₁, varies among general circulation models simulating Snowball Earth atmospheres (Abbot et al., 2013), therefore we investigate the sensitivity to this parameter. The Snowball Earth without a salt deposit is only accessible for large bare sea ice regions and extremely low eutectic temperatures, suggesting it is unlikely that Snowball Earth existed without salt deposits.

Figure 3 Bifurcation diagram for the EBM in Eq. (1) with the albedo function in Eq. (2) and: (a) Different diffusion rates D. The numerical value D=0.312 gives the temperate solutions to the EBM a modern-day temperature field, and is taken from North et al. (1981). The value D=0.25 is taken from Pierrehumbert et al. (2011) and represents a cold climate diffusion. A lower diffusion rate with numerical value D=0.13 was also explored. (b) Different combinations of critical isotherms T_s and T_e. (c) Different albedo values α₁ and α₂. (d) Different bare sea ice latitudes r₁. (e) Different salt crystal albedo values α₄. Other model parameters are given in Table 1. Solid lines indicate stable equilibrium solutions to Eq. (1) and dashed lines indicate unstable solutions. Solutions with multiple salt lines (see bifurcation points b₄ and b₅ in Fig. 1) are omitted here.

Based on this study, it is not possible to determine how long a salt deposit lasted on Snowball Earth. We observe that the range of atmospheric CO₂ concentrations supporting Snowball Earth without a salt deposit in the model is limited. The physical interpretation of this is that there is a limited window where the temperature on Snowball Earth is high enough to avoid activating the salt-albedo feedback. However, the model includes no ice dynamics, and it is widely recognized that sea glaciers must have formed and flowed equatorwards on Snowball Earth (Goodman and Pierrehumbert, 2003; Goodman, 2006; Pierrehumbert et al., 2011). Thus, the ice within the net ablation zone would have become progressively fresher over time as meteoric ice (ice formed by precipitation) moved in from higher latitudes. This makes it difficult to assess the duration and extent of the salt deposit for the Snowball Earth events. Goodman (2006) found that only a narrow zonal band of marine ice (ice formed by basal freezing) remained near the Equator over million-year timescales. Consequently, the climatic effects of a salt deposit probably became less important towards the later stages of these glaciations. Furthermore, the model has no explicit CO₂ input. The EBM is therefore not the optimal modeling framework to determine conditions for deglaciation, and we omit discussion of the range of atmospheric CO₂ concentrations associated with the alternative Snowball Earth states.

Figure 1 illustrates why assessing the duration of the surface salt deposit is important for Snowball Earth. In the model, terminating the salt-deposit Snowball Earth requires a much higher atmospheric CO₂ concentration than the classical Snowball Earth without a salt deposit. The deglaciation threshold for the salt-deposit Snowball Earth is the bifurcation point b₆ in Fig. 1, and the radiative forcing required to trigger the melting is substantially higher than for the state without a salt deposit, b₂. Perhaps the greatest challenge to the Snowball Earth hypothesis is explaining how deglaciation could occur within a realistic range of atmospheric CO₂ concentrations consistent with the duration of these events. Several efforts (Abbot and Pierrehumbert, 2010; Li and Pierrehumbert, 2011; Wu et al., 2021; Abbot et al., 2012; Abbot, 2014) have been made to lower the CO₂ deglaciation threshold in climate models to align with geochemical proxy data (Bao et al., 2008, 2009). If a salt deposit lasted for a significant time, the salt-albedo feedback may have provided an additional cooling element to Snowball Earth.

Our results indicate that the salt-albedo feedback could have played an important role in the initial phase of Snowball Earth by rapidly lowering temperatures due to the high albedo of salt crystals. Particularly since the full seawater system contains numerous species that precipitate, and initial precipitation of, e.g., mirabilite and hydrohalite occurs around −8 and −23 °C, respectively (Light et al., 2009), suggesting that some precipitation can be expected above $T_{\mathrm{e}}=-\mathrm{36}\mathrm{°}\mathrm{C}$. We note that our analysis only includes stationary states, and widespread glaciation must occur before a lag deposit of salt crystals forms. However, it is possible that the salt-albedo feedback promoted favorable conditions for the initiation of Snowball Earth by keeping sea ice albedo high as the Earth cooled. All states with partial salt deposits are unstable in the model. Once the initial salt deposits form, the system collapses into the Snowball Earth state with fully salt-covered deep tropics, given a warm pre-glacial climate. Moreover, it does not matter at what temperature this process begins, the system ultimately ends up in the fully salt-covered state (Fig. 2, for $T_{\mathrm{e}}>-\mathrm{36}\mathrm{°}\mathrm{C}$). Although the causes of the Snowball Earth initiation remain uncertain, our results show that it is difficult to avoid activating the salt-albedo feedback shortly after global glaciation sets in.

There are some important caveats to our results. First, it is uncertain whether a lag deposit of salt will form and persist on Snowball Earth. The simplified parametrization used here neglects ice dynamics and may therefore provide an overly crude representation of this phenomenon. Accumulation of surface salt crystals ultimately depends on the rates of sublimation and the downward migration of liquid brine in the ice (Light et al., 2009), processes that are not explicitly included in this model. The spatial and temporal extent of saline marine ice, where salt precipitation could occur during Snowball Earth, is also uncertain (Goodman, 2006). Furthermore, our salt parametrization relies on the assumption that water evaporated in the tropics is largely transported poleward before falling as precipitation, as consistently found in GCM simulations of Snowball Earth atmospheres (Pierrehumbert, 2005; Pierrehumbert et al., 2011; Abbot et al., 2013). This leads to a tropical net ablation zone, exposing bare sea ice and subsequently the precipitated salt crystals. Second, the strength of the salt-albedo feedback is uncertain. This feedback may be weaker in a more comprehensive model where processes such as cloud and wind dynamics are included. Braun et al. (2022) found that clouds can mask the albedo effect of bare sea ice, and clouds may likewise mask the salt crystals. Surface dust (Abbot and Pierrehumbert, 2010) may also reduce the albedo by partially covering salt crystals over time. Moreover, wind-driven redistribution of salt crystals (Light et al., 2016) into regions with net precipitation may ensure larger areas of salt-free bare sea ice in the net ablation zone. We used a high salt crystal albedo value in this study, however a lower value merely shifts the location of the bifurcation points b₄−b₆ in Fig. 1 (see Fig. 3e). More research is needed to assess the strength of the salt-albedo feedback on Snowball Earth, and this study should be seen as an initial investigation into the possible impacts of a salt-albedo feedback.

4 Conclusions

The initial phase of Snowball Earth was likely marked by rapid cooling associated with both a salt-albedo feedback mechanism and the well-known ice-albedo feedback. We have shown that a salt-albedo feedback introduces two coexisting Snowball Earth states in a simple climate model, one with a lag deposit of salt crystals and one without, where the former is significantly colder. We find that the warmer state is unlikely to participate in the initial stage of the Snowball Earth hysteresis, suggesting that the colder state is consistent with Neoproterozoic geology. Although the persistence of the colder Snowball Earth state remains uncertain, our results highlight salt precipitation as an important physical process that warrants further research in future modeling studies of Snowball Earth.

Appendix A: Accessibility of Snowball Earth without a salt deposit

Using boundary integrals, we derived an expression for the accessibility of Snowball Earth without a salt deposit. Specifically, we found an expression for when b₁ and b₃ from Fig. 1 take on the same A value, which represents the boundary between accessible and inaccessible regimes of the model parameter space.

Initially, we found the A value associated with the bifurcation point b₁ in Fig. 1. An expression for the hothouse solution to the stationary, dimensionless form of Eq. (1) can be obtained by inserting x₁=0 and $x_{\mathrm{2}}=\mathit{\pi }/\mathrm{2}$ in Eq. (3),

$$\begin{array}{}\text{(A1)}& u\left(x\right)=\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\mathrm{d}s\sin sg(s,x)f_{\text{hothouse}}\left(s\right),\end{array}$$

where

$$\begin{array}{}\text{(A2)}& f_{\text{hothouse}}\left(x\right)={\displaystyle \frac{Q}{T_{\mathrm{s}}D}}s\left(x\right)(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{1}})-{\displaystyle \frac{A}{T_{\mathrm{s}}D}}.\end{array}$$

Here the boundary values vanish from Eq. (A1) because of the no-flux boundary conditions (North, 1975), i.e.,

$$\begin{array}{}\text{(A3)}& \underset{x\to \mathrm{0}}{lim}\sin x{\displaystyle \frac{\partial u}{\partial x}}\left(x\right)=\underset{x\to \mathit{\pi }/\mathrm{2}}{lim}\sin x{\displaystyle \frac{\partial u}{\partial x}}\left(x\right)=\mathrm{0}.\end{array}$$

The bifurcation point b₁ in Fig. 1 is the A value where the temperature at the pole take the value for ice formation, T(0)=T_s. Let u_s denote the dimensionless ice formation temperature and u_e denote the dimensionless eutectic temperature. Inserting u(0)=u_s into Eq. (A1) we get

$$\begin{array}{}\text{(A4)}& u_{\mathrm{s}}=\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}\mathrm{d}s\sin sg(s,\mathrm{0})f_{\text{hothouse}}\left(s\right).\end{array}$$

We took a similar approach to determine the A value associated with the bifurcation point b₃ in Fig. 1. The Snowball Earth without a salt deposit terminates when the temperature at the bare sea ice latitude r₁ reaches the eutectic temperature, T(r₁)=T_e. An expression for the Snowball Earth solution without salt deposits was found by inserting x₁=0 and x₂=r₁ in Eq. (3),

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}u\left(x\right)& =\underset{\mathrm{0}}{\overset{r_{\mathrm{1}}}{\int }}\mathrm{d}s\sin sg(s,x)f_{\text{Snowball Earth wo salt deposit}}\left(s\right)\\ & +g(r_{\mathrm{1}},x)\sin r_{\mathrm{1}}{\displaystyle \frac{\partial u}{\partial x}}\left(r_{\mathrm{1}}\right)-u\left(r_{\mathrm{1}}\right)\sin r_{\mathrm{1}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{1}},x),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& f_{\text{Snowball Earth wo salt deposit}}\left(x\right)\\ & =\left\{\begin{array}{ll}\frac{Q}{T_{\mathrm{s}}D}s\left(x\right)(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{2}})-\frac{A}{T_{\mathrm{s}}D}& x<r_{\mathrm{1}}\\ \frac{Q}{T_{\mathrm{s}}D}s\left(x\right)(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{3}})-\frac{A}{T_{\mathrm{s}}D}& x\ge r_{\mathrm{1}}\end{array}\right..\end{array}\end{array}$$

Given the value $\partial u/\partial x\left(r_{\mathrm{1}}\right)$, we may find the A value associated with the bifurcation point b₃ by inserting u(r₁)=u_e in Eq. (A5). In order to determine $\partial u/\partial x\left(r_{\mathrm{1}}\right)$, we assumed hemispheric symmetry and extend the analysis beyond the equator, and applied Eq. (3) to the bare sea ice region (r₁,r₂), where $r_{\mathrm{2}}=\mathit{\pi }-r_{\mathrm{1}}$,

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }u\left(x\right)={\int }_{r_{\mathrm{1}}}^{r_{\mathrm{2}}}\mathrm{d}s\sin sg(s,x)f_{\text{Snowball Earth wo salt deposit}}\left(s\right)\\ {\displaystyle }& {\displaystyle }+g(r_{\mathrm{2}},x)\sin r_{\mathrm{2}}{\displaystyle \frac{\partial u}{\partial x}}\left(r_{\mathrm{2}}\right)-u\left(r_{\mathrm{2}}\right)\sin r_{\mathrm{2}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{2}},x)\\ {\displaystyle }& {\displaystyle }-g(r_{\mathrm{1}},x)\sin r_{\mathrm{1}}{\displaystyle \frac{\partial u}{\partial x}}\left(r_{\mathrm{1}}\right)+u\left(r_{\mathrm{1}}\right)\sin r_{\mathrm{1}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{1}},x).\end{array}\end{array}$$

As the solution is symmetric about the equator we must have $u\left(r_{\mathrm{1}}\right)=u\left(r_{\mathrm{2}}\right)=u_{\mathrm{e}}$ and $\partial u/\partial x\left(r_{\mathrm{1}}\right)=-\partial u/\partial x\left(r_{\mathrm{2}}\right)$ at the bifurcation point b₃. Inserting this and x=r₁ into Eq. (A7), we solve for $\partial u/\partial x\left(r_{\mathrm{1}}\right)$ and get

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\displaystyle \frac{\partial u}{\partial x}}\left(r_{\mathrm{1}}\right)={\displaystyle \frac{\mathrm{1}}{g(r_{\mathrm{2}},r_{\mathrm{1}})\sin r_{\mathrm{2}}+g(r_{\mathrm{1}},r_{\mathrm{1}})\sin r_{\mathrm{1}}}}\\ & [-u_{\mathrm{e}}+\underset{r_{\mathrm{1}}}{\overset{r_{\mathrm{2}}}{\int }}\mathrm{d}s\sin sg(s,r_{\mathrm{1}}\left)f_{\text{Snowball Earth wo salt deposit}}\right(s)\\ & -u_{\mathrm{e}}\sin r_{\mathrm{2}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{2}},r_{\mathrm{1}})+u_{\mathrm{e}}\sin r_{\mathrm{1}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{1}},r_{\mathrm{1}})]\end{array}\end{array}$$

Inserting u(r₁)=u_e and Eq. (A8) into Eq. (A5) we got an expression for the A value associated with the bifurcation point b₃,

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& u_{\mathrm{e}}=\underset{\mathrm{0}}{\overset{r_{\mathrm{1}}}{\int }}\mathrm{d}s\sin sg(s,r_{\mathrm{1}})f_{\text{Snowball Earth wo salt deposit}}\left(s\right)\\ & -u_{\mathrm{e}}\sin r_{\mathrm{1}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{1}},r_{\mathrm{1}})\\ & +{\displaystyle \frac{g(r_{\mathrm{1}},r_{\mathrm{1}})\sin r_{\mathrm{1}}}{g(r_{\mathrm{2}},r_{\mathrm{1}})\sin r_{\mathrm{2}}+g(r_{\mathrm{1}},r_{\mathrm{1}})\sin r_{\mathrm{1}}}}\\ & \times [-u_{\mathrm{e}}+\underset{r_{\mathrm{1}}}{\overset{r_{\mathrm{2}}}{\int }}\mathrm{d}s\sin sg(s,r_{\mathrm{1}}\left)f_{\text{Snowball Earth wo salt deposit}}\right(s)\\ & -u_{\mathrm{e}}\sin r_{\mathrm{2}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{2}},r_{\mathrm{1}})+u_{\mathrm{e}}\sin r_{\mathrm{1}}{\displaystyle \frac{\partial g}{\partial x}}(r_{\mathrm{1}},r_{\mathrm{1}})].\end{array}\end{array}$$

We solved Eq. (A4) for A and inserted this into Eq. (A9). The resulting expression is omitted for brevity, but it represents an expression for the boundary between the accessible and inaccessible regimes of the parameter space for the Snowball Earth without salt deposits, as shown in Fig. 2.

Appendix B: Supplementary figure

Figure B1 shows the solar radiation, albedo, and temperature for the two Snowball Earth states: with and without a salt deposit.

Figure B1 (a) Temperature profiles for Snowball Earth with a salt deposit (black) and Snowball Earth without a salt deposit (blue) at radiative forcing level ΔA=45 in Fig. 1. Other model parameters are given in Table 1. (b) The ice/snow albedo functions ($a(x,T)\forall T<T_{\mathrm{s}}$) for Snowball Earth with a salt deposit (black) and Snowball Earth without a salt deposit (blue), together with the latitudinal distribution of Neoproterozoic solar radiation Qs(x) (red).