Water contamination in an energy balance model of climate with coupled iceline

: Budyko’s energy balance model with coupled iceline is a historic system of diﬀerential equations, which incorporates water albedo as a major determinant to climate trends. In this class of models, the position of the Earth’s iceline, and its coupling with temperature, play a very important role in determining long term climate patterns (which may range between extremes such as a water covered or an ice covered Earth, but also allow for stable states with intermediate ice-lines). Some extensions of this model have speciﬁcally focused on studying the contribution of greenhouse eﬀects to the temperature-iceline long term dynamics. Our model builds upon this existing work, and studies in a computational setup the potential eﬀects of global perturbations in ocean water purity. We found that small changes in physical parameters such as water salinity or chemical/biological contamination may lead to sizeable changes in long term climatic outcomes. Moreover, these eﬀects appear to be enhanced when occurring in combination, and even more so when in the presence of elevated greenhouse eﬀects.


Introduction
Earth's climate is constantly changing. It exhibits fluctuations and cycles, some of which occur naturally, as part of the intrinsic dynamics of a self-regulating system. There is, however, increasing evidence that some climate trends are triggered by external, human-induced factors. Over the past decade, the idea that the Earth's climate may be significantly and potentially irreversibly altered by human activities has become a scientifically documented and accepted theory. Climate change is the subject of ample data-driven studies and of massive modeling efforts. These aim to understand which factors and human activities have the most detrimental effect on the system's dynamics and regulation, and predict whether this regulation can be restored if some of the negative impact is eliminated.
One of the primary effects studies in conjunction with climate change is the green-house effect. As the density of CO 2 (and of other similar gases) has been gradually increasing in the Earth's atmosphere over the past century, the annual average Earth temperature has been steadily climbing. This correlation is not accidental: CO 2 accumulating in the atmosphere prevents the light reflected by the Earth's surface from escaping back into space, hence the trapped solar radiation leads to the consistent increase in Earth average temperature.
There are now many mathematical models based on global patterns in the ocean and atmosphere, that can be used to support and understand empirically patterns displaying the green-house effect. One of the first pure theoretical investigations was Budyko's model, an energy balance model (EBM) that provides a broad brushstroke view of the dynamics of climate (as reflected into Earth's average temperature evolution), based on the interaction of major components known to govern climate patterns (such as incoming solar radiation, planetary albedo).
Budyko's EBM captures in basic and beautifully simple form the the ice-albedo feedback mechanism: the more snow and ice there is, the more solar radiation is reflected back into space. When this occurs, Earth grows colder and this results in more snow.
The model assumes the annual average Earth surface temperature T (in • C) to be constant around latitude circles, and symmetric across the equator. Then T = T (t, y), as a function of time t and of the sine of the latitude (y = sin(θ)), is broadly governed by the following equation: where T is the mean temperature over all latitudes: The change in temperature is primarily driven by the amount of incoming solar energy (insolation). The proportion of solar energy reflected back into space is determined by the Earth's albedo α.
In reality the albedo is a local measure, that differs between land and water, and even with local physical properties of these. For simplification, the model assumes that the Earth is covered in water (extrapolating the fact that the planetary ocean covers a significantly large proportion of the Earth surface). Then the albedo only depends on latitude (via the position with respect to the fixed iceline η ∈ [0, 1]), and can only take two discrete values: that of water (below the iceline), and that of ice (above the iceline): More radiation is reflected by ice than water, hence α s > α w . Altogether, the average annual rate at which solar energy is absorbed by the Earth's surface can be then expressed as Qs(y)(1−α(y, η)), where Q is the annual mean insolation, and s(y) represents the distribution of this insolation over all latitudes. The second term represents the energy that leaves the system via outgoing longwave radiation (OLR), which decreases with the green-house effect (greenhouse gasses in the atmosphere block radiation from escaping). The term is represented as a linear approximation, in which lower values of the parameters A and B are tied to higher amounts of atmospheric CO 2 . The final term represents heat transfer between latitudes (also known as "convection").
Over time, different models have improved on the Budyko framework with more specific interactions. Our work builds upon a two-dimensional extension constructed and studied by McGehee and Widiasih [12,15] as well as by McGehee and Walsh [14]. The model introduces a moving iceline η, and reproduces the coupled behavior of the temperature-iceline system.
If we assumed, as in Budyko's model, that the iceline were stationary, the average temperature across that iceline would be T c ∼ −10 • C (the approximate glacier formation temperature). It reality, however, the stationary assumption is not realistic [11]: if the temperature T at the iceline latitude y = η is greater than T c , the ice line will recede towards the pole (η = 1) and if the temperature is below T c , the iceline will advance towards the equator. To capture this effect, Widiasih [15] introduced an additional equation that couples iceline dynamics with that of temperature. The rate of change of the ice line is considered proportional to the difference between the equilibrium temperature profile and the glacier-forming temeprature T c (with very large time scale 1/ε reflecting glacier dynamics).
Broadly speaking, this coupling captures the idea that, when ice is melting, more water is exposed, absorbing more of the Sun's energy. This causes the surface to warm up, and more ice to melt. Conversely, when the ice is advancing, more of the Sun's energy is reflected, causing the surface to cool, causing more ice to form.
Via a Legendre expansion and change of variables to find a two dimensional invariant subspace, McGehhe and Widiasih obtained a simplified system of two coupled differential equations, describing the temperature-iceline dynamic interactions [12]: where a quadratic approximation was used for the insolation distribution: s(y) = 1 + s 2 2 (3y 2 − 1), and Ω represents the amount of energy required to melt a square meter of ice. Note that the iceline η is measured within the range of 0 to 1. Once the iceline reaches the equator, η = 0, or the pole, η = 1, the equations stop accurately describing these dynamics, and the system is governed by new regulatory drives, corresponding to an ice covered or water covered Earth, respectively.
This system was studied by McGehee, Walsh and Widiasih [14,12], with a combination of analytical and numerical computations. It was established that, for a parameter range compatible with current satellite measured-values, there are two locally stable equilibria (one with high iceline, and one in the nonviable negative η range), and a saddle (corresponding to a lower η), as shown in Figure 1a. Subsequently, the Earth may stabilize to long term steady state with small ice caps, or may end up completely covered in either ice or water. One direction of the analysis targeted the dependence of the long-term behavior of the system on the level of atmospheric greenhouse gasses (as captured by the parameter A). Numerical simulations revealed a saddle-node bifurcation, suggesting that excessively low A values (high greenhouse effect) may lead to a water covered Earth, excessively high A values (very low greenhouse effect) may lead to an ice covered Earth, while intermediate As, compatible to current measured values, fall in a range where a steady high iceline is possible (see Figure 1b). Figure 1: Dynamic behavior of the EBM with coupled iceline as described in [14]. A. Phase plane showing the three coexisting equilibria, for A = 202: a locally stable equilibrium (blue dot) at (η * , w * ) ∼ (0.95, 5); a saddle equilibrium (red dot) at (η * , w * ) ∼ (0.25, −17.5), and a second locally stable equilibrium with η * < 0 (out of range). B. Bifurcation diagram with respect to A. The high iceline stable equilibrium branch collides with the pole η = 1 at A − = 198.7, and with the saddle equilibrium branch at the saddle node bifurcation point A * = 211.6. The water covered Earth prediction is represented by an orange line at η = 1, and the ice-covered Earth is shown as a purple line at η = 0. The other parameters were fixed to the nominal values in Table 1, with α = 0.32 and T c = −10 • C.
In this paper, we will build upon this system and focus in more detail on two additional aspects that may contribute to the temperature-iceline coupling -and which were not included in the original model, in favor of simplicity. One aspect relates to the effect of clouds. Recall that planetary albedo is the percent of the sunshine that is reflected back into space, and is captured in our current model by the parameters α w (above water surfaces) and α s (above ice). Globally, the albedo is ruled by clouds, the formation of which is temperature dependent. The higher the temperature w, the warmer the sea surface, hence the higher the evaporation, and the more water in the atmosphere (the higher the cloud cover[5]). Empirical data sets have documented (lagged) correlations between the two, in a variety of settings, and some studies have attempted to quantify this relationship, at different spacial and temporal scales [9,6,13]. For example Oreopoulos and Davies [13] suggest that the variation in albedo due to temperature changes may depend linearly on the change in cloud cover with surface temperature and on local cloud albedo changes (i.e., changes in the albedo of each cloud pixel at ISCCP resolution). However, determining the influence of surface temperature on water albedo in a planetary ocean setting is not trivial, since the coupling appears to also involve strong feedback (cloud albedo contributes to ocean temperature regulation, at a local scale). This makes it harder to separate and quantitatively infer the direct influence of one onto the other simply from data. A model setting can help us set testable predictions with respect to the effect of clouds on the overall systemic dynamics.
The second aspect relates to water quality, and its potential effect on the system's dynamics -via its impact not only the water albedo, but also on its freezing point. Due to the nature of the model, we will only discuss this in terms of global changes to water quality. We consider the potential of water pollution to produce changes to the chemistry of the water and to aquatic vegetation, which may in turn alter the absorption index, hence also the albedo. We also consider the fact that most water contaminants (such as industrial pollutants, but also natural salt) are not retained in ice when water freezes. Hence lower icelines (larger ice caps) lead to higher concentrations in water contaminants and potentially lower freezing points, slowing down further ice formation, and adding a regulatory mechanism against the progression to a completely frozen Earth.
The rest of the paper is organized as follows. In Section 2, we describe the new aspects that we incorporated into the existing Budyko model with coupled iceline, and the modeling steps taken to address them. In Section 3, we present and contextualize the results of our simulations. In Section 4, we compare our predictions with those of the original model, as well as with realistic predictions based on climate data. Finally, in Section 5, we describe the limitations of our model, and potential ideas to pursue in future work.

Modeling methods
We will be studying the Budyko model with coupled iceline, as described in Equations (1)-(2), with parameter values and ranges specifies in Table 1. Known consequences of water pollution include physical and chemical changes. Among the known and well-studied effects are changes in freezing point in conjunction with certain contaminants, and changes in optical properties (reflectiveness of the water surface). In the case of our model, these properties are respectively captured by two parameters: T c , the water temperature of glacier formation, and α w , the water albedo. To incorporate potential effects of water contamination on the evolution of the temperature-iceline system, we allowed these parameters to vary within the ranges specified in the table, while keeping the other parameters fixed to the nominal table values. We will examine the dependence on parameters via bifurcation diagrams, computed numerically with the Matlab [10] and Matcont [8] softwares (unless otherwise specified).
As in our previous analysis [7], we performed a change of scale (e.i., a change in units of our time   [14]. constants). Using a millennium (10 3 years ∼ 3.16 × 10 10 seconds) as our time unit, the original values of our time constants R = 4 · 10 9 , Ω = 1.5 · 10 11 , ε = 3.9 · 10 −13 become R = 0.1266, Ω = 0.474, ε = 0.01264. All simulations are performed with these transformed values, hence the convergence kinetics in the simulations occurs on a time scale measured in thousands of years.
Modeling the effects of varying the freezing point T c . In the original model, the temperature T c where glaciers form is considered as a fixed parameter. However, as pointed in the Introduction, the freezing temperature of water is known to depend on the water purity. Some contaminants (among which salt, one of the main components in planetary ocean water) have the ability to significantly lower the freezing point. Hence we deem it important to investigate the effects on the system of changing the value of T c around its nominal value in Table 1. For simplicity, we assume that the planetary water is well stirred, so the water quality is the same everywhere (assumption which we further discuss in the Limitations section).
Modeling the effects of varying the albedo α w . In the original model, the water albedo is considered a fixed parameter. However, there is ample evidence that the albedo varies widely with circumstances (such as presence of clouds), and that these variations can have significant contributions to the iceline system dynamics. We study the effects to the long-term prognosis of the system of perturbations in the water albedo values α w around the nominal value in Table 1. Intrinsically, the albedo is higher up north, where the surface is covered in ice. However, sun exposure is dramatically reduced as one gets close to the poles (winters are longer). Therefore, under normal circumstances, ice contributes less to the total planetary albedo (does not reflect much sunshine in an average year) Hence, in our analysis, we will keep α s fixed to its table nominal value.
Modeling how greenhouse effects modulate the effects of salinity and albedo. We will first perturb T c and α w independently, to observe the effects of each factor in isolation. We then aim to study their combined contributions, as well as how these are modulated by other aspects significant to the system, such as the green house effect (captured by the parameter A). The system's sensitivity and dynamic changes in response to perturbations in the green house effect A have been abundantly studied by McGehee and Walsh [14]. In our analysis, we will focus on understanding whether the impact of changing T c and α w on dynamics is increased if the water pollution is accompanied by an increased green house effect (lower A).
Modeling the coupling of water albedo with iceline dynamics. For our modeling purposes, we will work under the assumption that water albedo increases with cloud coverage, which is turn increases with surface temperature, and that this dependence is sigmoidal: The dependence relates the idea that, at low w values, temperature variations do not present with cloud production, hence do not lead to significant changes in albedo. As w rises, temperature variations produce increasingly large effects in albedo, crossing a window of high sensitivity (where the slope of the sigmoidal function is high). The value w 0 (the inflection point) is the temperature where the water albedo is the most sensitive to changes in temperature w. The coefficient k quantifies this maximum sensitivity, and the value of the albedo at the inflection point is A w /2. The effect then gradually saturates, to an asymptotic albedo value A w . In our analysis, we will fix A w = 0.64, so that the optimal value A w /2 = 0.32 (corresponding to our average value used in the previous analysis, in which the albedo was considered constant). This is consistent with the estimate that cloud cover likely adds in average ∼ 0.25 to the surface albedo []. We allowed the temperature value where cloud production is highest to vary between w 0 = 10 • C and w 0 = 20 • C (variability which may indirectly reflect the dependence of the optimal cloud production point on atmospheric conditions other than temperature). We allowed the sensitivity to vary from a small value of k = 0.001 to a ten-fold value of k = 0.01.
Modeling the coupling of water salinity with iceline dynamics. As mentioned earlier in the section, one of the physical aspects that trigger changes in freezing temperature T c is the concentration of water contaminants. On one hand, this concentration is controlled by humans, if one considers water pollutants sourced in human activities and various industries. On the other hand, this concentration broadly depends on the amount of water in the planetary ocean (since most foreign substances are not incorporated in glacier ice through the freezing process). To fix our ideas, we will concentrate in this paper on one compound: salt (sodium chloride, always present in ocean water). In this context: the less glaciers, the more planetary water, and the lower the concentration of salt. This suggests that water salinity, and as a consequence the freezing point T c depends tightly on the iceline dynamics: η increasing from 0 to 1, leads to less ocean water and increased salinity, thus lowering the glacier freezing point T c . We will express salinity as a mass concentration δ, expressing the mass of salt per liter (or kilogram) of ocean water. We incorporate this direct dependence into our model through a simple approximate computation, based on the measure of water salinity for the current icecap position.
For simplicity, we assume each Earth hemisphere (from η = 0 for η = 1) to be a flat disc (incorporating the curvature into a spherical model is not expected to qualitatively change the predictions). We call R the total radius of this "disk-Earth" (R = 10 4 km [3]), and r 0 the current ice cap radius, measured from the pole (r 0 ∼ R/3). The total volume of water (in liters) is then is the average depth of the planetary , and r is the icecap radius; this also expresses the mass of water in kilograms (since 1L ∼1kg of water).
We call m the total mass of salt in the planetary ocean (which is independent on the amount of water, hence on the iceline, since salt does not freeze into glaciers). With this setup, we can write the dependence of salinity on icecap size r: While m cannot be directly measured, we can estimate it from the current water salinity δ 0 = 35g of salt per L (or kg) of water [2], as follows: The effect of salinity δ on water freezing temperatures has been estimated empirically as (in other words, an increase of 5 part per thousand in salt concentration decreases the freezing temperature by ∼0.28 degrees, with an approximately linear dependence [4]. We can then express the dependence of the freezing temperature on icecap size as: Using equation (7) and differentiating equation (5) with respect to r, we get: Hence: From equation (6), it follows that: Notice that an increasingly higher iceline (dη/dt > 0, corresponding to receding polar ice caps, hence more water in the planetary ocean and lower salinity) leads to a higher T c (dT c /dt > 0). For this part of the analysis, we will consider T c as a coupled variable in the system, together with the iceline η and average temperature w. We will analyse the long-term dynamics of the 3-dimensional system formed of equations (3) and (11).
Modeling greenhouse effects on iceline dynamics with coupled salinity and albedo. We will finally investigate whether variations in the greenhouse effect, as captured in the model by the parameter A, may have a more significant contribution to the evolution of the extended system, when incorporating the coupled salinity aspect and the coupled albedo component [14].

Dependence of the system dynamics on water quality
For this analysis, our key parameters are T c (the glacier forming water temperature) and α w (the water albedo), since they directly depend on water quality. In this section, to fix our ideas, we allowed the system to operate under a moderate amount of atmospheric carbon, specified by A = 202. In Section 3.2, we study how perturbations to higher or lower green-house green house gasses may affect the behavior of the system. For all simulations, all other parameters were fixed to their nominal values in Table 1.  We studied changes in the behavior of the system in response to perturbations to the freezing temperature T c around its baseline value T c = −10 • C, reflecting the effects of water contaminants (including salt) on the chemical potential of the resulting fluid to freeze. All other parameters were fixed as specified at the start of Section 3.1. In addition, the water albedo was fixed to its baseline value α w = 0.32. Figure 2 shows the behavior of the system's equilibrium under perturbations to T c in the form of a bifurcation diagram. A locally stable and a saddle branch of the steady state meet at a saddle node point at T * c = −4.92. In addition, notice that the stable (blue) equilibrium branch hits the top η = 1 iceline bound for T + c = −11.75, and that the unstable (dotted red) branch hits the bottom η = 0 iceline bound at T − c = −18.4. Hence the system only has access to the stable equilibrium (η * , w * ) for values of T c in the interval [-11.75,-4.92], and otherwise approaches asymptotically either a water-covered Earth (for values of T c < T − c ), or an ice-covered Earth (for high values of T c > T * c ). More specifically, if the freezing temperature T c is higher than the saddle node bifurcation value T * c = −4.92, the long-term result is an ice-covered Earth with long-term temperature lower than −20 • C (represented in the figure by the purple line). As T c is lowered past the bifurcation point, the fate of the system is determined by its initial state. A combination of high enough iceline and temperature places the system in the attraction basin of the stable equilibrium, and the result is a partly ice-covered Earth (blue curve), with equilibrium iceline and temperature increasing with lowering T c . If the initial state is not in this attraction basin, the result is an ice-covered Earth (purple line).
If one continues to lower T c past the value T + c = −11.75 • C, the stable equilibrium is no longer reachable within the biological range (since η * > 1), and all initial conditions within its attraction basin will converge in reality to an ice-free Earth, with a high average temperature over 5 • C (orange line). The other initial conditions (marked by the downward arrow) converge to the ice-covered Earth (purple line). Lowering T c beyond the value T − c = −18.4 • C will seal the long-term fate of the system to the ice-free scenario, independently on the starting conditions.
Recall that the current measured value of T c is approximately −10 • C, placing the system within the interval [T + c , T * c ], and predicting an equilibrium with iceline η * = 0.95 (small ice caps) and temperature w * = 4.95 • C, corresponding roughly to our current living conditions. Perturbations to higher T c values produce a smooth decline in the iceline, until the bifurcation value T * c marks the crash to the ice-covered Earth. Perturbations to lower T c values produce a smooth retraction of the iceline, until the value T + c marks the transition to an ice-free Earth. Both these extreme scenarios are relative well-separated from our current situation (T c = −10), since the values of T * c and T + c are far enough to not be reachable via small perturbations in water quality. However, it is possible that changes in other system parameters (such as the water albedo α w , or the green house effects A) could change the position and shape of this bifurcation scheme to the point where these scenarios may become dangerously within reach. In the next sections, we investigate whether this may be the case.

Dependence on water albedo α w
We next studied changes in the behavior of the system in response to independently perturbing the water albedo parameter around its baseline value α w = 0.32, to capture the effects of water contaminants on the water surface reflective properties. All other parameters were fixed as specified at the start of Section 3.1. In addition, the glacier freezing temperature was fixed to its baseline value T c = −10 • C.
The evolution of the system's equilibrium as α w changes is illustrated as a bifurcation diagram in This results in the system only having access to the locally stable equilibrium for values of α w in the interval [α + w , α * w ]. When the albedo is not within this interval, the system approaches in the long-term either a water-covered or a snowball Earth. The water-covered Earth is guaranteed for α w < α − w , and the ice-covered Earth is guaranteed for α w > α * w ; either outcome is possible when , very small positive perturbations around this value may shrink the attraction basin of the stable equilibrium, and subsequently redirect some initial conditions originally converging to a high iceline, to instead converge to an ice-covered Earth. Larger positive perturbations, that transcend the critical bifurcation value α * , would seal the icecovered Earth as the only possible fate for the system. Negative perturbations past α − w would cause all initial states to converge to the ice-free Earth.

Simultaneous dependence on α w and T c
The behavior of the equilibrium with respect to changes in the freezing temperature T c may significantly depend on further context. Since our discussion revolves around water quality, we can focus on tracking how this behavior changes when the water albedo α w varies between 0 and 1. For different fixed values of α w , the equilibrium curve with respect to T c has qualitative the same  U-shape, with two branches (a top, stable one and a bottom, unstable one) meeting at a saddle node bifurcation. However, the position of the saddle node changes, both in terms of the parameter value T * c and of the coordinates (η * , w * ) corresponding to it. So do the value T c = T + c , where the high (stable) equilibrium branch intersects the top iceline bound η = 1, and the value T c = T − c , where the low (unstable) equilibrium branch intersects the bottom iceline bound η = 0. Figure 4a illustrates the evolution of the saddle node point (α * w , T * c ) (with the corresponding evolution of the saddle node coordinates in the phase plane (η, w) being shown in Figure 4b). This bifurcation curve was computed with the Matcont software. Notice that, as α * w increases and T * c decreases along the saddle node curve, the iceline value η * declines, allowing for the possibility of lower and lower icelines. Figure 4a also shows the evolution of the collision point (α + w , T + c ) of the equilibrium with the upper boundary η = 1 (orange curve) and the collision point (α − w , T − c ) of the point with the lower boundary η = 0 (purple curve). The saddle node curve was computed and illustrated directly using the Matcont continuation software. For a good approximate illustration of these two curves, the boundary points T + c and T − c were extracted for five equilibrium curves (for α w ∈ [0.22, 0.62], in 0.1 increments), and these sample points were interpolated by hand, as the orange and respectively purple curves.
Notice that the purple curve intersects the saddle node curve at α * w = 0.62, T * c = −31.3268. For values of α * w > 0.62, the the saddle node leaves the physical domain (which is very far from the current values of α w = 0.32 and T c = −10). The thin, elongated parameter locus enclosed by these three curves (R) represents the region of water purity where the system has access to a nontrivial longterm iceline 0 < η * < 1. This gives a good representation of the effects of perturbations in water purity (as captured by the values of its albedo and its freezing temperature) from the current state α w = 0.32, T C = −10, which is inside the region R. If the α w and T c are simultaneously increased, the system eventually exits R across the saddle node bifurcation curve, and lands in the realm of the frozen Earth prognosis. If either α w or T c are decreased sufficiently, the system hits the brown curve, and access to the nontrivial iceline is replaced by access to the water-covered Earth (with the possibility of frozen Earth remaining available to certain initial conditions). If T c is lowered while α w is increased, the system hits the purple curve, and the water covered Earth becomes the only possible outcome.

Compounding greenhouse effects to water quality
In this section, we revisit the results in Section 3.1, examining their dependence on the parameter A, which encodes the strength of the greenhouse effect. This will allow us to better understand whether changes in water quality enhance the effects of climate change to the Earth's long term dynamics.

Dependence on glacier forming temperature T c
We want to explore to what extent greenhouse effects modulate the dependence of the system's iceline dynamics on water purity (as reflected in the glacier forming temperature T c ). In our model, increased greenhouse effects correspond to smaller values of the key parameter A. Hence, in Figure 5, we illustrate the system's equilibrium behavior and bifurcations with respect to T c , for three different values of A, corresponding to lower (A = 210), medium (A = 202) and higher (A = 190) concentrations of greenhouse gas in the atmosphere. The bifurcation diagram for each different A value is shown as a different curve. The evolution of the iceline (η * ) and temperature (w * ) components of the equilibrium are shown separately in the left, and respectively in the right figure panel.
Notice that, for the A range represented here, all equilibrium curves exhibit a saddle node bifurcation (green square), marking the collision of a locally stable branch (blue curve) and a saddle branch (red curve) of the equilibrium. Also, as before, the stable equilibrium curve hits in Figure 5: Bifurcation of the system's equilibrium with respect to the glacier forming temperature T c , for different values of A. The panels show the evolution of the iceline component η * (left) and of the temperature component w * (right) of the equilibrium, as the parameter T c changes. Each different curve represents the evolution of the equilibrium for a different fixed value of the parameter A. In each case, a locally stable branch (blue solid curve) and a saddle branch of the equilibrium (red dotted curve) collide at a saddle node bifurcation point, marked as a green square. The point where the stable equilibrium curve hits the upper iceline boundary η = 1 is marked with a brown dot, and the iceline and temperature corresponding to this extreme, ice-free Earth scenario, are shown for each A value as orange lines. The point where the saddle equilibrium curve hits the lower iceline boundary η = 0 is marked with a purple dot, and the iceline and temperature corresponding to this ice-covered Earth scenario, are shown as violet lines. The other parameters were fixed to the values in Table 1. 1306. This means that freezing temperatures T c which would lead to a snow cover Earth under low carbon conditions (higher A) may lead to a high iceline, or even to a completely water-covered Earth for high carbon conditions (lower A). Also notice that the temperature w corresponding to a water-covered Earth (η = 1) increases with higher carbon conditions.
Recall that the current freezing point T c , leading to formation of glaciers, is T c = −10 • C. Notice that for a small interval of high A values around 210, T c = −10 is close to the saddle node T * c , hence a small variation in water purity negatively perturbing T c may lead to placing the system on the left instead of the right side of the saddle node bifurcation, leading to a dramatically different longterm outcome (high iceline versus ice covered Earth). When A is a little lower (close to 205), small negative perturbations of T c around the current value of T c = −10 • C may push a high iceline system into a regime where only the extremes (water covered or ice covered Earth) are possible, depending on initial conditions. For A even lower (close to 190), a small negative perturbation from T c = −10 • C may leave the water covered Earth as the only possible long term scenario. In each case, a locally stable branch (blue solid curve) and a saddle branch of the equilibrium (red dotted curve) collide at a saddle node bifurcation point, marked as a green square. The point where the stable equilibrium curve hits the upper iceline boundary η + = 1 is marked with a brown dot, and the iceline and temperature corresponding to this extreme, ice-free Earth scenario, are shown for each A value as orange lines. The point where the saddle equilibrium curve hits the lower iceline boundary η − = 0 is marked with a purple dot, and the iceline and temperature corresponding to this ice-covered Earth scenario, are shown as violet lines. The other parameters were fixed to the values in Table 1. All equilibrium curves exhibit a saddle node bifurcation (green square), where a stable (blue) and unstable (red) branches meet. The stable equilibrium branch hits the upper iceline boundary at α w = α + w (brown dots); the unstable equilibrium branch hits the lower iceline boundary at a parameter value α − w . As A decreases, the equilibrium curves shift to the right, as follows: for A = 190, we have α Notice that water albedos α w which would lead to a snow-covered Earth under low carbon conditions (higher A) may lead to a high iceline, or even to a completely water-covered Earth for high carbon conditions (lower A). Also notice that the temperature w corresponding to a water-covered Earth (η = 1) decreases with higher carbon conditions. Notice that, for low carbon (A close to 210), the current water albedo α w = 0.32 is between the values α + w and α * w , and close to both of these values. A small perturbation to higher albedo may lead to placing the system on the right instead of the left side of the saddle node bifurcation, leading to a dramatically different long term outcome (ice covered Earth versus high iceline). A small negative perturbation of α w may push a high iceline system into a regime where only the extremes (water covered or ice covered Earth) are possible, depending on initial conditions. Figure 7: Saddle-node bifurcation curves with respect to the albedo α w , for different values of the greenhouse effect A. The physical portion of each bifurcation curve (corresponding to iceline 0 ≤ η * ≤ 1) is shown in green; the portion corresponding to negative iceline is shown as a dotted pink continuation. The two scenarios are separated on each curve by a cyan dot. The pink diamond represents a cusp point, which occurs in the negative iceline range.

Simultaneous dependence on T c and α w
To understand the simultaneous dependence of the system's dynamics on freezing temperature and water albedo, a traditional approach is to separate the regions corresponding to different behaviors in the parameter plane (α w , T c ), and how this profile changes with different circumstances. In Section 3.1.3, our simulations identified a saddle-node curve for a fixed value of A = 202. This curve delimits a region of (α w , T c ) values that may lead to a high iceline, or even a water-covered Earth (to the left of the curve) from a region of (α w , T c ) in which only an ice-covered Earth is possible in the long-term (to the right of the curve). Our further computations show that is in fact representative for a larger range of A.
To illustrate the changes brought on by perturbations in A, Figure 7 shows the corresponding saddle-node curves in the plane (α w , T c ) for three different A values. While these curves are very similar in shape, they visibly shift to the right as A is decreased. This suggests, in line with our intuition, that an increased greenhouse effect (lower A) increases the system's long-term propensity for water versus ice.

Iceline dynamics with adjustable albedo
Our next step was to consider the effect of clouds on the dynamics of the coupled iceline-temperature system defined by equations (3). We did so by allowing the water albedo α w to vary as a function of temperature, as described in equation (4), thus simulating some basic aspects of the temperaturedependent mechanism contributing to cloud formation. The albedo dependence on w (via formation of clouds) is shaped by the two sigmoidal parameters: w 0 (the temperature where the albedo is most sensitive to temperature) and k (which modulates this maximum sensitivity). Since the values of w 0 and k are now known empirically (although physics and common sense can suggest realistic ranges) we studied how the behavior of the system depends on perturbations of these parameters within their ranges. Figure 8: Effect of introducing clouds on model's phase-plane dynamics. Top. Trajectories of the original system in the (η, w) phase plane may converge towards a negative iceline equilibrium (out of range), leading to a water-covered earth (blue curves), or may converge to a high iceline equilibrium (green curves, with the equilibrium shown as a green circle). The two behaviors are delimited in the phase plane by the stable curve of the saddle equilibrium (not shown), which separates the two attraction basins. Bottom. Trajectories of the system with clouds, for sigmoidal parameter values k = 0.005 and w 0 = 10. For comparison, the same initial conditions were labaled with the same numbers (1)-(4) in the two panels. Figure 8 illustrates some basic differences between the typical phase plane dynamics of the original system, and that of the system with clouds, for two fixed values of the sigmoidal parameters. The subtle dependence of α w on temperature does not seem to significantly affect (in this particular illustration) the number or position of the locally stable equilibria, yet there are important differences between the two phase planes. These stem from the fact that the dynamics around the saddle equilibirum has been altered, thus affecting the attraction basins of the two asymptotic outcomes, in a way that redirects initial conditions with lower icelines to evolve in time towards the high iceline, locally stable equilbrium. Figure 9 illustrates how the equilibrium curves with respect to T c and A are affected when the sigmoidal parameters are varied. In both cases, the curves were pushed to the left by increasing values of k (from k = 0.001 to k = 0.1) and with higher values of w 0 (from w 0 = 10 to w 0 = 20). This means that, for the same initial conditions, higher cloud formation sensitivity and peak temperatures may lead to pushing a system that would otherwise converge to a water-covered outcome, to a high iceline instead. Along the same lines, higher k and w 0 may push a system converging to a high iceline towards a lower iceline, or even to an ice-covered outcome instead.

Iceline dynamics with coupled salinity
We considered the system of three equations, with coupled variables defined by the iceline η, the temperature w, and the glacier freezing point T c , as described in equations (3) and (11).
We compared the long-term behavior of this system, which incorporates the effects of coupled salinity, with that of the original system (i.e., with T c viewed as a fixed parameter). The first notable difference when introducing the third equation is a first order degeneracy, leading to the system having a whole one-dimensional subspace of nonisolated equilibria (as shown in Figure 10). This degeneracy is associated with a zero eigenvalue along this curve (hence the equilibria are neutrally-attracting along the direction of the curve). The other two eigenvalues are either both negative (in which case the equilibrium locally attracts a two-dimensional subspace, hence we call it neutral-attracting), have opposite signs (neutral-saddle equilibria), or are order-two degeneracies that mark the transitions between these two behaviors. Figure 10 illustrates schematically how this situation affects the behavior of phase space trajectories (for two values of A, for comparison). In both cases, the equilibria corresponding to 0 < η * < 0.16 (thick blue portion of the curve) and η * > 0.58 (thick green portion) are netralattracting points, while those corresponding to 0.16 < η * < 0.58 (dotted red portion) are neutralsaddle. A few phase space trajectories are also shown, to suggest the shape of the attraction basins for the two neutral-attracting equilibrium portions: the colors of the trajectories match that of the outcome; the initial point is marked as a dot, and the equilibrium to which they converge asymptotically is marked as a circle with the corresponding color. Note that, while all trajectories were initiated at the same freezing temperature T c = −10, the value of T c changes in time as expected, along with the position of the iceline and the water salinity. For further comparisons, however, we will focus on the iceline-temperature coupling, hence we will consider two-dimensional projections onto the (η, w) phase plane. Figure 11 compares the behavior of the (η, w) coupled dynamics between the three-dimensional system with evolving salinity and the original system. For better comparison, the (η, w) phase plane was represented in both cases, with the original behavior shown in the top panels, and the new behavior shown in the bottom panels (which are 2D projections of the two panels in Figure 10). The left panels illustrate this comparison for A = 202, the right panels for A = 190.
Recall that the original (η, w) system with A = 202 exhibits three equilibria (as shown in Figures 10a): a locally attracting equilibrium with high iceline (green circle), broadly attracting initial conditions with higher icelines; a locally attracting equilibrium with negative iceline, causing the initial conditions with lower iceline to evolve into a water covered Earth; a saddle equilibrium (red circle), the unstable curve of which determines the separation between these two behaviors. The major difference in the system with coupled salinity (shown in Figures 10c) is that different initial conditions lead to different outcomes along the curve of nonisolated equilibria, with the possibility to converge to any of a continuum of different steady states with high iceline (green portion of the curve) or to different steady states with low icelines (blue portion of the curve). In particular, this new context allows for the possibility of the system to converge to low icelines η * > 0. This is a significant qualitative difference form the original system, which will be more  Table 1.
thoroughly interpreted in Section 4. Notice that the separation between these attraction basins is very similar (if not quite identical) to that in the original system, although the convergence to either of the two behaviors is now determined not by the position of one single saddle point, but by the neutral-saddle portion of the curve.
Recall that in the original system, where decreasing A pushes the locally attracting equilibrium to increasingly higher icelines (to the point where η > 1, leading to an ice-covered Earth), and pushes the saddle equilibrium to lower icelines, broadening the attraction basin of the locally attracting equilibrium (i.e., the range of initial states that in the long term will lead to high icelines, or to an ice-covered Earth versus a water covered Earth). When comparing between the behavior of the coupled salinity system for two different levels of greenhouse effect A, the effect is broadly similar, in the sense that lower A enhances convergence to higher icelines, by a dual effect on the position of the equilibria, and on the attraction basins. But this effect is enriched by the existence of an entire curve of equilibria, including those with low icelines. For example, an initial point that would lead to a blue trajectory for A = 202 (convergence to a low iceline or to a water-covered Earth) may lead to a green trajectory (converge to a high iceline) for A = 190. An initial condition which leads to a green trajectory under A = 202 may converge to a higher iceline, or even to an ice-covered Earth, for A = 190.

Iceline dynamics with coupled water salinity and albedo
In this section, we explore briefly a few interesting aspects of a combined model that considers simultaneously adjustable albedo and coupled salinity. Figure 12 illustrates the behavior of the system where adjustable albedo (as described and analyzed in Section 3.3) was added to the 3-dimensional model with coupled salinity (as analyzed in Section 3.4). The two panels illustrate the phase space dynamics for the same two levels of greenhouse effect considered in Figure 10, for the system with coupled salinity, but fixed albedo. Figure 11: Comparison between (η, w) phase plane dynamics between the original system and the system with coupled salinity, for two different levels of greenhouse effect: A = 202 (left) and A = 190 (right). Top. Phase plane for the original system. The high iceline locally attracting equilibrium is marked with a green circle (when in the 0 < η * < 1 range), and a few trajectories attracted to it are shown in green. A few trajectories attracted to the negative iceline equilibrium are shown in blue. The saddle point separating the two behaviors is shown as a red circle. When trajectories hit the boundary η = 0 or η = 1, the exit point is marked with a matching color triangle. Bottom. (η, w) phase plane for the system with coupled salinity (2D projections of the two panels in Figure 10). The system has a curve of nonisolated equilibria: for low iceline η * < 0.16 (left side, in blue) and for high iceline η * > 0.85 (right side, in green) the equilibria are neutral-attracting; for intermediate iceline 0.16 < η * < 0.85 (red dotted portion) they are neutralsaddle. A few trajectories that converge to high iceline equilibria are shown in green, and a few trajectories that converge to low iceline equilibria are shown in blue. When trajectories hit the boundary η = 0 or η = 1, the exit point is marked with a matching color triangle.
When comparing The two figures, notice that adjusting the albedo affects the separation between initial conditions evolving toward high versus low icelines (adjustable albedo enhances convergence to high icelines, hence some trajectories which are blue in Figure 10 are green in Figure 12). (This is not surprising, since in the original system, the effect of introducing adjustable albedo was to broaden the attraction basin of the high iceline equilibrium.) However, the size and distribution of Figure 12: (η, w, T c ) phase space dynamics for the model with adjustable albedo and coupled salinity, for two different levels of greenhouse effect: A = 202 (left) and A = 190 (right). In both cases, the system has a curve of nonisolated equilibria (not shown); low iceline neutral-attracting equilibria are separated from the high iceline neutral-attracting equilibria by the dynamics around the intermediate iceline, neutral-saddle equilibria. A few trajectories that converge to high iceline equilibria (shown as green circles) are shown in green, and a few trajectories that converge to low iceline equilibria (shown as blue circles) are shown in blue. Some trajectories hit the boundary η = 1, leading to a water-covered Earth. The simulations were performed for w 0 = 10 and k = 0.005; all other parameters were fixed at nominal values.
this effect in the phase space seems to depend on the values of A (notice that, while for the higher value of A = 202, the attraction basin of high iceline equilibria is increased, for the lower value A = 190, the trajectories are virtually unchanged between the two figure). Studying these patterns in more depth is one of the goals of our future work.

Discussion
Significant prior research on the energy balance model with dynamic iceline had explored the contribution of greenhouse effects to the evolution of the Earth's climate system. In our study, we analyzed extensions of the energy balance model with dynamic iceline, centered around the impact of water quality (in conjunction with greenhouse effects) on the system's long term dynamics. To do this, we first focused on understanding how the dynamics of the original EBM with coupled iceline depends on two key parameters, that may be considered to characterize different aspects of water purity: the glacier formation temperature T c , and the albedo α w . We analyzed the system's asymptotic behavior under perturbations of each of these two parameters independently.
When we examined the dependence on glacier forming temperature T c , we found that the system undergoes a saddle node bifurcation with respect to this parameter, so that access to a high iceline stable equilibrium is only possible for a relatively small interval for T c . For all values outside of this interval, the system either approaches a water-covered or an ice-covered Earth. This suggests that even a small content of water impurities (which alter the freezing point T c ) may compound to a dramatic long-term effect.
When we analyzed the system's dependence on water albedo α w , we also identified the presence of a saddle-node bifurcation, relatively close to the currently measured value of α w . This suggests that the contributions of water purity to the oceanic albedo augment the effect that small fluctuations in water purity may have to the Earth's climate in the long-term.
Once we established broadly the dependence of the system on α w and T c considered independently, we aimed to understand the effects of them being varied simultaneously (as a realistic scenario for changes in water quality), and under different levels of the greenhouse parameter A. First, for a level of A = 202 (close to the value estimated for our current greenhouse cover), we identified the region in the (α w , T c ) parameter plane where the system has access to a high iceline steady state. We noticed that our current values of α w and T c place us not only inside this regions, but also far enough from its boundaries, meaning that small perturbations are not sufficient to cause the system to leave this region through a dramatic phase transition. For example, when α w and T c increase simultaneously within this region, the equilibrium iceline decreases (leading to larger ice caps). However, we also noticed that the dependence on T c and α w (hence the position and shape of this stability region) are significantly modulated by the value of A, suggesting that this robustness of our current state to small perturbations may be lost at higher levels of the greenhouse effect. Broadly speaking, we found that higher greenhouse effects increase the system's long-term bias towards water versus ice. For example, we found that lower values for A result in the potential for the Earth to reach the water-covered state at lower temperatures.
To incorporate the fact that water albedo has been empirically tied to the presence of clouds (which in turn depend on temperature), we extended our model to give α w a phenomenological (sigmoidal) dependence on temperature. We found that this dependence affected primarily the shape of the attraction basins in the phase plane of the system, with cloud cover enhancing the tendency for high icelines, and more water. This is not surprising, since, in our modeling setup, higher temperatures lead to more clouds, which in turn lead to even higher temperatures. In reality, however, clouds are known to have a cooling effect on the Earth surface, which may tone down this feed-forward cycle (see the section on Limitations).
Since water contaminants (such as salt) impact the value of T c , and since the concentration of contaminants depends on the amount of planetary ocean water (i.e., on the position of the iceline), we expanded our model to incorporate T c as an iceline-dependent system variable (rather than as a fixed parameter). Since salt is a "contaminant" in all ocean water, we called this the "model with coupled salinity." A simple calculation based on empirical evidence lead to an equation relating the derivative dTc dt proportionally with the derivative of the iceline dη dt . This lead to a system having non-isolated equilibria ("neutral" in one direction), which opened up the potential for a few interesting phenomena. First, this scenario allows different initial conditions to converge to different (neutral-stable) outcomes, some of which have high iceline, and some which have low iceline. This already transcends the behavior of the original system (with fixed T c ), in which low icelines were not accessible asymptotically. Second, the curve portions with low versus high iceline neutral-stable equilibria are separated by a portion with neutral-unstable points, which governs to a large extent the attraction basins of the system. On one hand, this introduces a new "sensitivity" feature into the system: in order to know the long-term state of the system, one always needs to know precisely the initial state. While this is not equivalent to the traditional concept of sensitivity to initial conditions, it would still create predictability problems in a real-life system (related to the accuracy of measuring the initial state).
Our last analysis aimed to combine the adaptable albedo and the coupled salinity into a more comprehensive model. Not surprisingly, this model combined features of the two extensions considered separately. When exploring the potential outcomes to variations of the greenhouse effect in this system, we noted that some of the features of this extended model (such as the adjustable albedo) appear to be highlighted at some greenhouse levels, and not as relevant at other greenhouse levels.
Our future work will aim to better understand the contributions of the greenhouse effect, in combination with other environmental factors, to climatic dynamics. Overall, our analysis suggests the intriguing possibility that, in addition to the well-studied and known direct effects of greenhouse gasses accumulating in the Earth's atmosphere, there may be indirect aspects, coupled with the contribution of other factors pertinent to climate dynamics (such as water purity). Our results place increased emphasis on the idea that seemingly unrelated human-generated effects, each known to have a detrimental effect to the environment, may in fact compound effects, and have dramatic consequences to the future and survival of life on our planet.

Limitations and future work
There are a few serious limitations to our model, which should hence be considered only a first, proof of concept step towards addressing the potential climatic effects of water pollution in conjunction with greenhouse effects.
First, all our model variables and parameters are defined as "global" measures. In reality, the surface temperature w, the water albedo α w , the glacier forming temperature T c , all depend on local factors. A more refined model would have to take into consideration such heterogeneities.
Second, as mentioned earlier in the text, the assumption we made on clouds is that the cover increases with temperature. In fact, one of the difficulties of establishing precisely from empirical data the dependence of cloud formation on temperature raises precisely from the fact that the two are bidirectionally interconnected at the local level. For example, there is ample evidence that a thick cloud cover may cool off the Earth surface (depending on the type of clouds, on other atmospheric parameters like pressure or wind, on presence or absence of rain, etc). A more careful inspection of data relating cloud coverage and surface temperature may be able to suggest a finer and more realistic feedback coupling between these two components.
Other limitations relate to the assumption of a flat Earth for the salinity computation, and from the extrapolations made from empirical measures (e.g., T c may not depend linearly on salinity throughout the salinity range).