Optimal control of polar sea-ice near its tipping points (2024)

EBM: bifurcations and tipping points

The EBM was developed in the classic papers of Budyko and Sellers^37,38. Here, we use a standard one-dimensional diffusive EBM³²; it is derived by considering the heat balance on an infinitesimal strip centered at latitude θ. In our EBM, the zonally-averaged, annual mean surface temperature T evolves as:

where \(x=\sin \theta\) is the latitudinal coordinate, and subject to boundary conditions imposing no heat flux at the pole and equator. For simplicity, the system is assumed to be symmetric about the equator and therefore the spatial domain considered is 0° ≤ θ ≤ 90°. The ocean mixed layer heat capacity is given by C and is assumed to be a constant. The latitude-dependent incoming annual mean shortwave (SW) radiation is Q_SW(x), the outgoing longwave radiation (OLR) term is aT + b₀³⁹ and the external forcing applied to the EBM to mimic the effect of increasing/decreasing GHG concentration is given by F(t). The ice-albedo α(T) is temperature dependent and is assumed to be a smoothed step function with α = 0.32 for T ≫ −10 °C and α = 0.62 for T ≪ −10 °C; these asymptotic albedo values are based on the Budyko formulation in ref. ³⁷. A more detailed discussion of the EBM used in this study can be found in the “Energy balance model” section, and the EBM parameter values are given in Table 1.

The sharp transition in ice-albedo across the ice edge, taken to be at T = −10 °C, results in a strong positive nonlinear feedback with important consequences for the dynamics of this model. The EBM with albedo feedback is known to exhibit multistability⁴⁰, which is defined as the existence of multiple possible time-invariant states for a given set of parameters. Continuous changes in system parameters produce qualitative changes in the stability characteristics of the EBM. These changes are known as bifurcations, and the existence of these bifurcations in the EBM is what gives rise to multistability (for certain values of those system parameters) More specifically, all the bifurcations are of saddle-node type, characterized by the appearance or disappearance of a pair of steady-state solutions^41,42 as a system parameter is varied. Saddle-node bifurcations cause catastrophic shifts and are associated with hysteresis behavior; many other tipping points of the Earth system such as AMOC and West Antarctic Ice Sheet are also hypothesized to be associated with saddle-node bifurcations^5,43.

In the EBM, the system parameter creating these bifurcations is \(F\left(t\right)\) (i.e., the level of GHG forcing). The runaway feedback from the strongly nonlinear albedo during a global cooling event can lead to a sudden transition to the so-called “snowball” Earth⁴⁴ once the ice cover crosses a critical latitudinal threshold. This tipping point is commonly referred to as a large icecap instability.

A second tipping point is a small icecap instability (SICI)^40,45 which leads to irreversible loss of icecaps smaller than a characteristic length scale associated with diffusive energy transport. The bifurcation curve corresponding to the SICI (dotted line) in Fig. 1b represents the steady solutions to the EBM for different values of the control parameter F (refer to the “Numerical simulations” section for computational details). Three different branches of steady solutions can be noted: the top branch, with the ice-line at the pole, corresponds to the ice-free stable solutions; the bottom branch consists of stable solutions with a finite icecap; and the middle branch represents the intermediate unstable solutions. The yellow diamonds represent the two saddle-node bifurcation points associated with SICI. The right bifurcation point represents the limit past which the EBM has no finite icecap steady solutions. Similarly, the left bifurcation point represents the end of the ice-free branch. The small icecap instability has been a subject of intense scrutiny over the past decades owing to its implications for a rapid, irreversible transition to a perennially ice-free Arctic under greenhouse warming^7,46.

Response to greenhouse warming

In order to understand the response of the EBM to greenhouse warming in the vicinity of the SICI, we subject it to a spatially uniform, linear ramp-up in the GHG forcing parameter F such that the increase in forcing is 5 W/m² at the end of 150 years (≈0.033 W/m²/year) and is then brought back linearly to F = −3 W/m² over the next 150 years. The starting state is the steady solution corresponding to F = −3 W/m².

In the presence of the forcing defined above, the evolution of the temperature profile and ice-line, indicated by the black dashed line, is shown in Fig. 1a. The irreversibility associated with SICI manifests as asymmetry in the responses of mean temperature and ice-line during the ramp-up and ramp-down phases. Initially, the icecap recedes linearly; however, around year 140, there is an almost abrupt disappearance of sea-ice cover, which does not recover until year 270.

The hysteresis curve (solid line) corresponds to the time-varying solutions in response to the applied forcing and is shown in Fig. 1b along with the bifurcation diagram for this model. It does not follow the bifurcation curve (dotted line) exactly because the rate at which the applied forcing varies is too rapid for the system to adjust to steady state. Consequently, the bifurcation diagram only provides an estimate of when the system is likely to tip. The tipping time in real-world systems is determined by the interplay of the short-term variability in the climate system (modeled as noise to the long-term variability of climate) and system response time. As the system approaches the tipping threshold, the restorative forces that damp perturbations weaken (critical slowdown^25,47), and the susceptibility to noise increases—a strong perturbation may therefore cause the system to tip well before the forcing threshold⁴⁸. On the other hand, if the external noise levels are low, the tipping time will be determined by the system response time. Owing to the system inertia, the tipping threshold may be exceeded temporarily without triggering rapid transition, referred to as “overshooting”^35,36,49.

In Fig. 1b, we define the hysteresis width ΔF as the horizontal distance between the bifurcation points (yellow diamonds) and can be used as a measure of the degree of irreversibility and instability⁵⁰. The tipping size Δθ is defined to be the vertical distance between the bifurcation points and can be interpreted as the critical size for receding icecaps, past which the system transitions rapidly to an ice-free state in a warming climate. During the ramp-up (red line), the system initially remains close to the stable finite-ice branch but as it approaches the tipping point, there is accelerated loss of ice and the system quickly transitions to the ice-free branch, i.e., ice-line retreats to the pole. During much of the ramp-down stage (blue line), the system remains on the ice-free branch and it is only once it crosses the second saddle-node bifurcation that the system recovers sea-ice. The hysteretic behavior implies that once the tipping point is crossed, modest efforts to reduce forcing alone will not suffice to reverse sea-ice loss.

Control around the tipping point

An effective control strategy in this context is one that can arrest or even reverse the crossing of the small icecap tipping threshold. Such a control strategy must be amenable to the strong nonlinearity near the critical threshold and adaptable to the different regimes of the system resulting from multistability.

The control mechanism chosen for our study is an artificially introduced planetary albedo-based control u(x, t), that can vary based on latitude. Note that estimates of u(x, t) can be used to design any general solar radiation management methods that locally modify the planetary albedo such as stratospheric sulfate aerosol injection⁵¹, marine cloud brightening⁵² and surface albedo modification⁵³. Once we solve for the optimal albedo control, the optimal external forcing required for a general control strategy can be obtained as uQ_SW. The controlled EBM is then given by:

Note that there exists a control u_e(x, t) that exactly cancels any forcing F(x, t) at all times t:

For a spatially uniform forcing, the requisite exact control u_e increases monotonically from the equator to the pole. Solar insolation drops as we move away from the equator, so at higher latitudes, a larger fraction of the incident radiation needs to be reflected to offset the same forcing (see SI for details).

Exact control, being passive, does not alter the stability characteristics of the EBM and therefore is ineffective once the sea-ice extent drops below the critical value (tipping size). Additionally, exact control does not take into account the cost associated with the deployment of control, so it is not cost-efficient.

Optimal control

We first formulate a cost function that accounts for the deviation from some target temperature profile T_f(x) and the scale of intervention u(x, t) at every time t and spatial location x. Here, the cost function \({\mathcal{C}}\) is chosen as:

where T_f(x) is the desirable state to be attained and the deviation from the target profile and the control are suitably normalized by \(\bar{T}\) and \(\bar{u}\), respectively (see section “Variational approach for optimizing control”). The relative prioritization of control cost and temperature deviation, which will depend on the climate policy goals, will be reflected in the relative values of the weights w_u and w_T. The objective can then be written as:

where \({\mathcal{L}}[u,T]\) is the total cost functional, and t_f is the time period over which the control is to be applied. We also need a constraint on u to ensure that the net albedo remains less than one for all x and t—i.e.:

This constraint is imposed implicitly by restricting the weight ratio w_T/w_u. As the relative penalty on the control decreases, the fraction of the cost associated with control in Eq. (3), increases. Therefore, by keeping w_T/w_u under some upper bound λ_max, the constraint can be satisfied. This upper bound depends only on the EBM parameters and can be determined numerically.

The principles of variational calculus can be employed to derive the Euler–Lagrange equation³⁴ corresponding to the minimization problem in Eq. (4), and is given by:

Therefore, at every time t, the optimal control u can be computed using the instantaneous T from (2) and the local sensitivity to control, \({T}_{u}:= \frac{\partial T}{\partial u}\). Consequently, to compute u(x, t), T_u(x, t) must be known and this can be done by deriving an additional equation for the time evolution of T_u from the EBM. This system of PDEs for the controlled EBM is then solved by integrating forward in time. A detailed derivation of the optimal control can be found in the “Variational approach for optimizing control” section and the computational methods are described in the “Numerical simulations” section.

The weight ratio \(\frac{{w}_{T}}{{w}_{u}}\) is a free tuning parameter to balance the cost of control and the rate of approach to the target temperature profile. This will be explored in detail in a subsequent section.

In what follows, we consider scenarios relevant to two plausible tipping pathways—one in which the system tips prematurely and the other in which tipping is delayed due to system inertia or overshooting.

Case 1

The effects of internal variability and other sources of noise on the control are neglected here for simplicity. However, while examining the possible tipping scenarios, one must take into account the role of external noise in determining the tipping time—a system that is “close” to its bifurcation point may tip prematurely if the noise intensity is strong enough to cross the potential barrier height associated with the attractor basin for the steady solution^47,49, i.e., the system may transition to the ice-free regime before the GHG forcing exceeds the tipping threshold due to the destabilizing effects of noise.

In order to probe this tipping pathway, three different starting scenarios are considered for optimal control. Figure 2a shows the starting states on the three branches of the bifurcation diagram associated with SICI, all of which correspond to F = 0 W/m². Two of these states are stable (ice-free and finite-ice); the ice-free state represents a plausible final state following a perturbation-induced transition. We also choose a point on the unstable equilibrium line as it provides a convenient reference point. It could, for example, be a point on a dynamic non-equilibrium trajectory that undergoes a premature transition. The unstable starting state allows us to investigate the efficacy of the optimal control in stabilizing and reverting the system to a stable finite-icecap state—even in the presence of GHG forcing. The chosen target state (filled-in triangle) corresponds to F = −0.5 W/m² on the finite-icecap solution branch. The EBM is subject to a spatially uniform linear ramp-up and ramp-down forcing (up to 3.5 W/m² and back to zero) over a period of 140 years, mimicking the ramp-up of CO₂ forcing at a rate of 1%/year followed by a ramp-down at the same rate to the original CO₂ forcing level (Fig. 2b) that represents the gradual lowering of CO₂ to pre-industrial levels through the lowering of emissions, carbon capture, etc.⁵⁴. The model itself does not require us to specify the mechanism by which the forcing level changes, so it is sufficient for our purposes here that there exist plausible examples of such mechanisms. The weight ratio w_T/w_u is set to be 0.1 for which the applied albedo control is effective and remains physical.

For each of the three scenarios, the root-mean-square (RMS) of control at time t, a measure of the cost of that control (see Eq. (3)):

is compared in Fig. 2b.

The response of the EBM with control can be roughly divided into two segments. The first stage (t ≤ 35 years) is dominated by the initial transition via control to the target state; the differences in the requisite control amongst the three starting conditions are manifested mainly during this time interval. The time scale for this adjustment is set by the ocean heat capacity C and the weight ratio w_T/w_u which constrains the growth of control. Reducing w_T/w_u lowers the magnitude of control u during the adjustment to the target state.

Significant differences can be noted in the cost associated with the initial adjustment for the three scenarios. The area under the RMS control curve in Fig. 2b, during the adjustment phase for the unstable starting state is approximately twice that for the finite-ice starting case and the corresponding control cost, u_cost for the first 20 years is nearly four times that of the unstable case. From Fig. 2c, it is clear that most of this increased cost can be attributed to the transition from the unstable branch to the stable finite-ice branch. For the initial state in the ice-free regime, the maximum of RMS control for the adjustment is about four times that of the control required for the finite-ice initial state. However, u_rms does not tell us the whole story: the maximum local albedo control required for the stable finite-ice adjustment is 0.01 while those for the unstable and stable ice-free states are approximately 0.03 and 0.1, respectively, as seen in the three sub-panels in Fig. 2c. Therefore, the restoration of sea-ice starting from the ice-free branch requires a large external albedo imposed over a sizeable latitudinal zone and this will likely have significant ramifications for the global climate. Another characteristic of note is that for the case starting from a stable ice-free state, there is a sharp increase in applied control once the ice reappears at around year 10 and lasts until the target ice-line is achieved. The optimal control given by (5) depends on both the deviation from the ideal temperature profile and the sensitivity of temperature to control T_u. The reappearance of ice bolsters the sensitivity to control T_u since the local positive ice-albedo feedback now works in synergy with the external albedo u to reduce ΔT_f and results in the sharp increase of u near the ice-line to take advantage of the increased sensitivity there. Note that the time scale associated with ice reappearance is set by the choice of weight ratio w_T/w_u, which regulates the rate at which the applied control increases, and the heat capacity C, which modulates the temperature response rate.

The second stage corresponds to the stabilized stage when the optimal control is identical for all three scenarios and only works to counteract the ramp-up and ramp-down GHG forcing to maintain the target state. A key feature of the optimal control in both stages is the spatial localization around the ice-line: the largest values of u are along the ice-line. This is in contrast to the exact control, which requires maximum albedo change at the poles. This observation suggests that spatially localized control interventions may be adequate for reversing the rapid loss of sea-ice.

Case 2

In this section, we explore intervention scenarios that account for the possible overshoot of the tipping threshold due to system inertia, which is regulated by the heat capacity C in the EBM. To understand the dependence of the control requirement on the system response time, three out-of-equilibrium initial conditions are considered on a time-varying trajectory of the EBM. The system is allowed to evolve from the steady state at F = −3 W/m² subject to ramp-up forcing (0.05 W/m²/year) and three different intervention times are chosen as shown in Fig. 3a. The first corresponding to F = 0.25 W/m² that is just before the tipping threshold (F = 0.45 W/m²). The second initial condition is chosen at F = 1 W/m², which is in the overshoot phase. The final one is chosen at F = 2.25 W/m², by which time the transition to the ice-free regime is complete. The target state is chosen to be the same as in case 1.

The variation in u_rms for the three intervention times are shown in Fig 3b. Similar to case 1, the major differences between the three intervention times occur during the adjustment to the target state. All three states evolve identically after the target state is attained. Therefore, the appropriate measure to compare the cases is the additional control requirement associated with the delayed intervention (∫Δu_rmsdt) and is defined here as the area enclosed by the closed loop consisting of the u_rms curve for delayed intervention (green or red) and the pre-tipping (blue) control curve. As the time of intervention is delayed, this additional control requirement increases. In comparison with the overshoot case, ∫Δu_rmsdt for the ice-free starting state is significantly larger.

To further explore the role of the overshoot in determining the control requirement, the above analysis is repeated for a range of intervention times, 0 ≤ t_i ≤ 40 (in years), and for a range of heat capacities, 1.4 ≤ C ≤ 8.2 (in W year/m²/K). The additional control requirement (∫Δu_rmsdt) as a function of heat capacity and intervention time is shown in Fig. 4. The time (t = 9 years) at which the ramp-up forcing reaches the tipping threshold is indicated by a dash-dotted line. For a given heat capacity, as the intervention time exceeds a certain limit, an almost step change in the control requirement can be noted. The dashed line in Fig. 4 represents a simple quadratic fit for this limiting intervention time past which there is a steep rise in the control requirement. Physically, this corresponds to the onset of delayed tipping in the system. If the intervention is delayed beyond this time, the system rapidly loses sea-ice and becomes ice-free. The extra time afforded by the system response time scale or the “overshoot shoulder” (t_o) is defined as the distance between the two lines. As the system response becomes sluggish, there is a larger leeway in the intervention time, but this is accompanied by a steeper change in the control requirement past the overshoot shoulder. Therefore, even though increased system inertia delays tipping, that same system inertia significantly increases the control costs once the system is past the overshoot shoulder; the cost increase, in and of itself, is not surprising, but the suddenness of that increase beyond a certain time point is unexpected.

Tuning weight ratio for optimal control

In this formulation, and given the parameter values in Table 1, the weight ratio w_T/w_u has a bounded operational range. For values of w_T/w_u less than λ_min ≈0.04, the control, u does not scale up fast enough to effectively control the climate. For values of w_T/w_u greater than λ_max ≈6.7, the control u can produce a net albedo greater than one (which is unphysical).

For the three initial conditions described above, the cost of control:

varies with weight ratio within the operational range, and is shown in Fig. 5. From Fig. 5, it is clear the control costs for all three states increase when the deviation from the target temperature profile is given more weight than the applied control, as one would expect. The differences in costs between the three initial states also increase with larger w_T/w_u values—especially near the upper limit. Also noteworthy is the almost constant cost for the pre-tipping scenario (blue curve). The weak sensitivity of the total control cost for the case of the pre-tipping control implies that one may prioritize, and thus hasten, the attainment of the climate target without significantly increasing the cost of control. The results in Fig. 2 compare the control costs for the three scenarios corresponding to w_T/w_u = 0.1, which is close to the lower limit of the operational range (which is approximately 0.04, as mentioned previously) and therefore provides an approximate lower bound for the control required around the SICI.

Tipping point characteristics

Hysteresis width ΔF is defined as the width of the hysteresis loop and has units of W/m² as shown in Fig. 1. Hysteresis width can be used to quantify the irreversibility associated with a tipping point and may be used to characterize the behavior of a tipping point. The hysteresis width is mainly determined by the nonlinearity in the co-albedo and the diffusivity parameter D in the EBM. The latter can be tuned to vary the hysteresis width associated with SICI⁵⁰.

To understand how control costs scale with hysteresis width, we first compute the control cost of restoring the finite-ice stable state from the ice-free initial state at a distance of 0.1ΔF from the SICI tipping point (similar to the red and blue dots in Fig. 2a, which are at a distance of approximately 0.2ΔF) and then compare control costs as D and the associated ΔF systematically vary. This study is carried out at constant GHG forcing in order to isolate and analyze the cost of adjustment from the ice-free regime to the finite-ice regime. Intriguingly, the control costs (with w_T/w_u = 0.5) do not scale monotonically with hysteresis width (Fig. 6). Rather, the costs initially decrease as the hysteresis width decreases with increasing D. Then, upon reaching a minimum at around D = 0.4, the costs start to increase. This behavior can be understood once the tipping size Δθ (see Fig. 1) associated with the SICI—the distance between the latitude of the tipping point and the pole—is taken into account. The tipping size provides a measure of how much ice growth is required before the system will self-adjust to the stable finite-ice branch. The tipping size is known to increase in an EBM with increasing diffusivity (see ref. ⁴⁰) and this is responsible for the increase in costs at higher diffusivities. Significant uncertainties exist in our estimates of critical thresholds and irreversibility for tipping elements^3,4, so it is imperative to consider the effect of these attributes when evaluating system control characteristics and reversibility around a tipping point. We also carried out an analogous study that considers the scenario in Case 2 subject to a linear ramp-up in GHG forcing and these qualitatively similar results are shown in the SI. Additionally, the control sensitivity to the nonlinearity in the co-albedo is also explored in the SI.