The tipping times in an Arctic sea ice system under influence of extreme events

Arctic sea ice variations are important indicators of climate changesMin et al. (2008); Notz and Marotzke (2012). Satellite observations have revealed a substantial decline in September Arctic sea ice extent since the late 1970s Scheweiger, Wood, and Zhang (2019). Several studies Budyko (1969); Sellers (1969); North, Gahalan, and Jr (1981); Crowley (2000); Zhang, Rothrock, and Steele (2000) have shown that the generation and melting of sea ice are affected by energy flux involving seasonal variations in solar radiation, thermodynamics, and heat transport in the atmosphere and ocean. Due to the complexity of the overall system, it is necessary to use mathematical models to determine how these physical processes interact in order to provide explanations for both observed and possible future behaviors of the time evolution of Arctic sea ice.

In the context of energy flux balance models, Budyko Budyko (1969) and Sellers Sellers (1969) have recognized the advantages of simple deterministic theories of climate that provide a clear assessment of stability and feedbacks Moon and Wettlaufer (2017). Afterwards, a deterministic single-column energy flux balance model has been proposed for the Arctic Ocean by Eisenman and Wettlaufer Eisenman and Wettlaufer (2009). Subsequently, this model has been improved by including additional physical mechanisms, such as the influence of clouds Abbot, Silber, and Pierrehumbert (2011); Hill, Abbot, and Silber (2016), in the time-dependent terms of the equation.

There are two reasons to include random variables and/or processes in a deterministic model: (i) to provide a statistical model for differences between the output of the deterministic model and observational data, such as in time series models (for example, ARMA models) of observational data Gao et al. (2016), and (ii) to circumvent the problem of accounting for environmental influences that are too complex to include in the deterministic equations Kampen (1976).

External disturbances that fluctuate rapidly relative to slowly responding state variables are often treated as stochastic processes with very short correlation times. For example, Hasselman Hasselman (1976) showed that the evolution of slowly varying “climate variables” subject to short-time, fluctuating “weather variables” can be approximated by diffusion processes. In these models, climate variability, driven by short time-scale fluctuations, may be large or small depending on the degree of stabilizing feedback. A common feature of observations is that ice extent exhibits Gaussian noise structure on annual to biannual time scales Agarwal and Wettlaufer (2018). Consequently, a stochastic Arctic sea ice model with Brownian motion has been considered. This stochastic model can be used to quantify how white noise impacts the potential transition between the ice-covered state and the ice-free oneMoon and Wettlaufer (2017); Chen et al. (2019). Recently, increasing attention has been given to extreme events that trigger the variations and evolutions of Arctic sea ice. Extreme weather events, such as heatwaves, droughts, floods, hurricanes, blizzards and other events, that occur rarely and unpredictably, can have a significant impact on climate changes and human survival. Meanwhile, it is reported that extreme events of the type that occur in weather and climate fluctuations are more appropriately modeled as realizations of power law, or heavy tailed, distributions Farazmand and Sapsis (2017); Selmi et al. (2016). These characteristics of extreme events, described above, have strong non-Gaussianity, which cannot be described by general Gaussian noise. A Lévy process is thought to be an appropriate model for such non-Gaussian fluctuations, with properties such as intermittent jumps and heavy tail. Furthermore, researchers have indicted that the presence of $\alpha$-stable Lévy noise could imply that the underlying mechanisms for abrupt climatic changes are single extreme events Ditlevsen (1999); Zheng et al. (2020). We therefore consider an Arctic sea ice model under influence of $\alpha$-stable non-Gaussian Lévy noise in the following study.

In view of the rapid decrease in summer sea ice extent in the Arctic Ocean during the past decade Livina and Lenton (2013), a potential tipping point in summer sea ice has been considered in recent studies of dynamical systems. The term tipping point commonly refers to a critical threshold at which a tiny perturbation can qualitatively alter the state or development of a system Lenton et al. (2008). Tipping points associated with bifurcations or induced by noise are studied in a simple global energy balance modelAshwin et al. (2012); Sutera (1981); Lucarini and Bódai (2017, 2019).

In this paper, we use the maximal likely trajectory to determine tipping times for the most probable transitions from a perennially ice-covered state to a seasonally ice-free one, and from a seasonally ice-free state to a perennially ice-free one, respectively. We expect that the most probable tipping time serves as a valid indicator to predict when Arctic sea ice begins to melt, and to estimate the transition time to the most devastating state corresponding to when the Arctic sea ice melts completely.

The structure of this paper is as follows. We present our methods, and introduce an Arctic sea ice model driven by $\alpha$-stable Lévy process in Section II. Then we conduct numerical experiments to investigate the impact of non-Gaussianity and the greenhouse effect on the tipping times for transitions in this stochastic Arctic sea ice model in Section III. Concluding remarks appear in Section IV. In Appendix A we provide a brief description of the symmetric $\alpha$-stable Lévy process. Appendix B contains definitions and nominal values of parameters used in our calculations.

In this section we define the maximal likely trajectory for a stochastic dynamical system driven by $\alpha$-stable Lévy noise. The definition is based on the solution of its associated nonlocal Fokker-Planck equation which is presented below. We then introduce a model for stochastic Arctic sea ice subject to extreme weather events which, in turn, are modeled by $\alpha$-stable Lévy noise.

A. Nonlocal Fokker-Planck equation for the probability density

Consider the following scalar nonautonomous stochastic differential equation (SDE):

where $X(t)$ is an $\mathbb{R}$-valued stochastic process. Here $L^{\alpha}_{t}$ is a symmetric $\alpha$-stable Lévy process, with $\alpha\in(0,2]$, defined on the probability space $(\Omega,\mathcal{F},P)$ (See Appendix A). The positive noise intensity is $\epsilon$. The nonautonomous drift term $f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ satisfies a Lipschitz condition to ensure the existence and uniqueness of the solution of the SDE $({\ref{e1}})$. A symmetric $2$-stable process is a Brownian motion, which is a Gaussian process. When Lévy index $\alpha\in(0,2)$, the $\alpha$-stable Lévy process is a jump process. A brief introduction to the $\alpha$-stable Lévy process is given in Appendix A.

For $x_{0}\in\mathbb{R}$, we suppose that the nonautonomous SDE $(\ref{e1})$ has a unique strong solution, and the probability density for this solution exists and is strictly positive. The probability density function $p(x,t)\triangleq p(x,t|x_{0},0)$ of the solution process $X(t)$ driven by a non-Gaussian $\alpha$-stable Lévy process satisfies the following nonlocal Fokker-Planck equation Duan (2015); Sun and Duan (2012):

with initial condition $p(x,0)=\delta{(x-x_{0})}$.

When $\alpha=2$, the SDE (1) is driven by Brownian motion. In this case the density function $p(x,t)$ is the solution of the Fokker-Planck equation Duan (2015):

We use the same initial condition for equation (3) that we used for equation (2), namely $p(x,0)=\delta{(x-x_{0})}$.

In the present paper, we use the “punched-hole" trapezoidal numerical algorithm of Gao et al.Gao, Duan, and Li (2016) to find the solution of the nonlocal Fokker-Planck equation (2).

B. The maximal likely trajectory

When it comes to trajectory for a nonautonomous SDE, one likely option is to plot sample solution orbits from an initial state, mimicking deterministic phase portraits. However, each sample solution trajectory is an “outcome” of a trajectory, which could hardly offer useful information for understanding dynamics. So we consider the maximal likely trajectory, which is determined by the maximizes of the probability density function $p(x,t)$, at every time $t$.

Given $X(0)=x_{0}$, the maximal likely state Zeitouni and Dembo (1987, 1988) $x_{m}(t)$ of the stochastic system (1) at each time $t\in[0,T_{f}]$ is defined by

Here the maximizer $x_{m}(t)$ for $p(x,t)$ indicates the most probable location of these orbits at time $t$.

Thus, for the point $t_{i}$, based on equations $(\ref{e2})$ and $(\ref{e3})$, we get the maximal likely state $x_{m}(t_{i})$ by computing the maximum of $p(x,t_{i})$. We connect this series of $\{(x_{m}(t_{i}),~{}i=1,2,\ldots\}$ to get the maximal likely trajectory. The distances between the mesh points $\{t_{i},~{}i=1,2,\ldots\}$ need to be sufficiently small to get a good approximation of the maximal likely trajectory. Note that the maximal likely trajectory $x_{m}(t)$ is not a solution of nonautonomous SDE (1).

C. Arctic sea ice model

We consider a model for Arctic sea ice established by Eisenman and Wettlaufer Eisenman and Wettlaufer (2009). This model is based on energy balance model where the energy per unit surface area, $E$ (with units $Wm^{-2}yr$). When $E<0$, the layer is interpreted to consist entirely of ice. Conversely, when $E\geq 0$, the layer is interpreted to consist entirely of water, that is, the layer is in its ice-free state. In this model $E$ satisfies the following nonautonomous differential equation Eisenman and Wettlaufer (2009):

where

and

Here $\alpha(E)$ is the surface albedo, which is a central control element involved in transition dynamics. The fraction $(1-\alpha(E))F_{S}(t)$ is the amount of absorbed short-wave radiation by the ice albedo $\alpha_{i}$ and the ocean albedo $\alpha_{ml}$. The deterministic term $\alpha(E)$ describes the energy flux balance at the atmosphere ice (ocean) interface where we calculate the surface temperature $T(t,E)$. Specifically, the state variable $E$ has the physical interpretation that the energy is stored in sea ice as latent heat when the ocean is in the ice-covered state (i.e. $E<0$) or in the ocean mixed layer as sensible heat when the ocean is in the ice-free state (i.e. $E\geq 0$). The term $\Delta F_{0}$ represents the reduction in outgoing long-wave radiation due to increased greenhouse gas forcing levels. Incident surface short-wave radiation $F_{S}(t)$ and basal heat flux $F_{B}$ are specified at central Arctic values Maykut and Untersteiner (1971). The final term $\nu_{0}\mathit{R}(-E)$ in equation $(\ref{e4})$ is the fraction of sea ice pushed by wind out of the Arctic each year. $\mathit{R}(-E)$ ensures this term is zero when there is no sea ice, where

The seasonally varying parameters $F_{0}(t)$ and $F_{T}(t)$, which are used to determine the surface energy flux, have values computed by using a atmospheric model Eisenman and Wettlaufer (2009). Nominal values and additional details about the model parameters are described in Appendix B (Table (1)).

For the deterministic Arctic sea ice model $(\ref{e4})$, Eisenman and WettlauferEisenman and Wettlaufer (2009) have constructed a Poincaré map to obtain the periodic solutions, as shown in Fig. 1. Throughout this paper, the “lower” and “upper” stable periodic solutions of the deterministic system are denoted by $E^{*}_{l},E^{*}_{u}$, respectively, and the “middle” metastable one is denoted by $E^{*}_{m}$. In order to facilitate the numerical calculation, we write the deterministic model using the scale transformation $E/100$.

In the case of Arctic sea ice, a fixed point corresponds to a steady state solution of the yearly cycle of ice growth and retreat Abbot, Silber, and Pierrehumbert (2011). In this paper we use the method used by Hill, Abbot, and Silber Hill, Abbot, and Silber (2016) to construct Poincaré maps, which indicate the energy changes from 1 January in one year ($E_{n}(0)$) to 1 January in the next year ($E_{n+1}(0)=f(E_{n}(0))$) for a range of initial energies. Fixed points $E^{*}_{n}(0)$ of the Poincaré map occur when $E_{n+1}(0)=f(E^{*}_{n}(0))=E^{*}_{n}(0)$, and they correspond to the periodic solutions of the system with the same periodicity as the forcing (i.e., annual). For different greenhouse gas forcing, $\Delta F_{0}$, we use a Poincaré map to get the fixed points. As shown in Fig. 2, for $\Delta F_{0}\in[10,16)$, there is only one stable fixed point, corresponding to the perennially ice-covered state. For $\Delta F_{0}\in[16,20]$, there are two stable fixed points (the perennially ice-covered state and the perennially ice-free one) and one unstable intermediate fixed point. For $\Delta F_{0}\in(20,23]$, there are two stable fixed points: the seasonally ice-free state and the perennially ice-free state. For $\Delta F_{0}\in(23,25]$, there is only one stable fixed point, which is the perennially ice-free state.

D. Stochastic Arctic sea ice model

The dynamical system $(\ref{e4})$ is a deterministic model. Even though it successfully captures the seasonal cycle of Arctic sea ice thickness and predicts the nature of the transitions as the greenhouse gas forcing $\Delta F_{0}$ increases. However, a realistic feature of the dynamic behavior of the ice cover is its variability due to internal fluctuations and external forcing Fetterer and Untersteiner (1998). Following Ditlevsen Ditlevsen (1999), we consider the following nonautonomous differential equation driven by a scalar symmetric $\alpha$-stable Lévy process with $1\leqslant\alpha\leqslant 2$:

where the term is $f(E,t)=\left(1-\alpha(E)\right)F_{S}(t)-F_{0}(t)-F_{T}(t)T(t,E)+\Delta F_{0}+F_{B}+\nu_{0}\mathit{R}(-E)$ corresponds to the right side of equation (5). The positive quantity $\varepsilon$ denotes the noise intensity. Thus, equation (6) represents the Arctic sea ice model subject to extreme events as modeled by a non-Gaussian $\alpha$-stable Lévy process. This is a pure jump process when the Lévy index $\alpha$ satisfies $1\leq\alpha<2$. It is known that a pure $\alpha$-stable Lévy process has smaller jumps with higher jump probabilities as $\alpha$ approaches $2$. In this stochastic system, the noise intensity and Lévy index can be regarded as the intensity and the frequency of the extreme events. The special case $\alpha=2$ corresponds to the usual Gaussian process, which is used to model “normal" atmospheric fluctuations.

In this section, we analyse how the noise intensity $\epsilon$, Lévy index $\alpha$ and greenhouse gas forcing $\Delta F_{0}$ affect the maximal likely trajectory of the nonautonomous SDE $(\ref{e8})$. We determine the tipping times for transitions from the perennially ice-covered state to the seasonally ice-free one, and from the seasonally ice-free state to the perennially ice-free one. In the following, we choose one century (i.e. $T=100$) as the computational terminal time. Computations are carried out over a time interval of one century, that is, $T=100$ years. Following Eisenman and Wetlafer Eisenman and Wettlaufer (2009) we use the terminology “summer” sea ice to refer to the annual minimum of Arctic sea-ice area and thickness and “winter” sea ice to refer to the annual maximum.

A. Effect of $\alpha$-stable Lévy process for the weakened greenhouse effect level $\Delta F_{0}=19$

For the relatively weak greenhouse gas level $\Delta F_{0}=19$, Fig. 3 illustrates the maximal likely trajectory from the perennially ice-covered state to the seasonally ice-free state for Lévy index $\alpha=1.50$ and noise intensity $\epsilon=0.0450$. The blue curve indicates that the stochastic system maintains a perennially ice-covered state, which is enlarged in Fig. 3. The red curve indicates that the system has entered an unstable regime of operation. During the computational time interval of one century, the system begins to experience an ice-free state for a small fraction of time each year. Each year thereafter, the length of time in the fractionally ice-free state increases (see Fig. 3) until the system enters a stable regime of operation in which the length of time in the fractionally ice-free state remains constant, as shown in Fig. 3.

Based on the behavior discussed in the previous paragraph, we refer to the time when the system makes a first transition from a perennially ice-covered state to a state in which a fraction of each year is spent in an ice-free state as the early-warning tipping time, and denote it by $T_{1}$. After $T_{1}$, the ice-covered period is decreasing until the Arctic appears in ice-free state in summer. $T_{1}$ could be regarded as an early-warning signal of anomalous Arctic sea ice, and it helps us to predict the approximate time when Arctic sea ice will disappear in summer. If $T_{1}=100$, the system never appears ice-free and it remains in a perennially ice-covered state during our simulation cycle (one century).

We wish to understand the qualitative behavior of the dependence of the early-warning tipping time $T_{1}$ on the Lévy index $\alpha\in[1,2]$ and the noise intensity $\epsilon$. This dependence is shown by the graphs in Fig. 4. We find that $T_{1}$ decreases as $\epsilon$ increases. This agrees with the corresponding result for Gaussian noise. For $\epsilon\in[0.051,0.09]$, $T_{1}$ increases as $\alpha$ increases. However, when $\epsilon$ is greater than $0.09$, $T_{1}$ is insensitive to changes in $\alpha$. To determine the smallest value of $\epsilon$ that leads to a transition away from a perennially ice-covered regime within one century, draw a horizontal line in Fig. 4 at a level slightly below $T_{1}=100$. We see that the minimum noise intensity required for such a transition increases as $\alpha$ increases. For $\alpha=1.25,1.5,1.75$, the minimum values of $\epsilon$ are approximately $0.03528,0.043096,0.051$, respectively. For Gaussian noise, the minimum required noise intensity is clearly much larger.

We know that the $\alpha$-stable Lévy process has smaller jumps with higher jump probabilities for larger values of $\alpha$ ($1\leqslant\alpha<2$). This means that the frequency of extreme events increases when Lévy index $\alpha$ approaches $2$, and the intensity of extreme weather events increases as Lévy noise intensity increases. These results on early-warning tipping times show that increased extreme events intensity lead to the earlier melting of the sea ice. For the small intensity of extreme events, the increased frequency of extreme events has a positive effect on predicting the time for sea ice to melt in summer.

Once the sea ice begin to melt in summer, the system becomes unstable, and with evolution of a fraction of each year, the system will reach the seasonally ice-free state. In order to study the effect of non-Gaussian $\alpha$-stable noise on the time when ice-free state appeared in summer, we will introduce the ice-free time $T_{2}$ to represent the period of the year when the Arctic Ocean is ice-free as shown in Fig. 3. It means that the larger $T_{2}$ is, the longer the a period of time for ice-free state in one year. If the climate becomes warmer further, it will be easier to appear a perennially ice-free state for the Arctic Ocean. It is worth pointing out that the Arctic sea ice system remains in the ice-covered state all year round if the ice-free time $T_{2}=0$.

Fig. 4 demonstrates the dependence of the ice-free time $T_{2}$ on the Lévy index $\alpha\in[1,2)$ and the noise intensity $\epsilon$. We find that $T_{2}$ becomes longer as $\epsilon$ increases under influence of both Gaussian and non-Gaussian Lévy noise. Another observation is that there is an intersection point near $\epsilon=0.065$. When $\epsilon$ is larger thane the intersection point, the ice-free time is longer as the value of $\alpha$ increases. However, the ice-free time $T_{2}$ has opposite behavior when $\epsilon$ is smaller than $\epsilon=0.065$.

These results imply that one has to consider both the value of $\alpha$ and $\epsilon$ when we examine a fraction time of each year in an ice-free state. The ice-free period increases as the intensity of extreme events increases. Meanwhile, we see that the decreasing intensity and frequency of extreme events will be effective in reducing ice free times in summer. These results show the importance for studying extreme weather events.

B. Effect of $\alpha$-stable Lévy process for the enhanced greenhouse effect level $\Delta F_{0}=21$

For this larger value of $\Delta F_{0}=21$, the maximal likely trajectory of the nonautonomous SDE $(\ref{e8})$ from the seasonally ice-free state to the perennially ice-free state for Lévy index $\alpha=1.75$ and noise intensity $\epsilon=0.0250$ is shown in Fig. 5. The blue curve represents the stochastic system concentrated on the seasonally ice-free state, as enlarged in Fig. 5. The red curve means that the system has experienced an unstable ice-free state, the length of time in the fractionally ice-free state increases (Fig. 5) until the system remains in a perennially ice-free metastable states (as $E_{m}>0$ throughout the year), as shown in Fig. 5.

To analyze how the system of Arctic sea ice shifts from the seasonally ice-free state to the perennially ice-free state under noise, we will continue to use the early-warning tipping time $T_{1}$ to denote the first time when the seasonally ice-free state changes, as shown in Fig. 5. After $T_{1}$, a fraction of ice-free time increases until the winter sea ice vanishes to the perennially ice-free state. Once this time $T_{1}$ is reached, we should develop and deploy adaptive strategies, and take a more pre-emptive, precautionary policy approach to prevent the situation from getting worse.

Fig. 6 shows that the early-warning tipping time $T_{1}$ presents a ladder descending trend with the increased noise intensity $\epsilon$ for non-Gaussian Lévy noise with $\alpha=1.25,~{}1.50,~{}1.75$. Additionally, for different $\alpha$, we find that the range of Lévy noise intensity that enables the system to shift to the perennially ice-free state is different. For example, for $\alpha=1.25,1.75$, the corresponding transition range of $\epsilon$ is $[0.0211,0.058]$ and $[0.0245,0.0897]$, respectively. On the other hand, for the Gaussian noise, the behaviour of $T_{1}$ agrees in general with the results in previous pure jump case. However, compared with non-Gaussian Lévy noise, Gaussian noise requires the stronger noise intensity and have the larger range of noise intensity for transition to the seasonally ice-free state.

Furthermore, we find that $T_{1}$ obviously appears earlier for the enhanced greenhouse level $\Delta F_{0}=21$ by comparing with the weakened greenhouse level $\Delta F_{0}=19$, as shown in Fig. 6 and 4. For example, keeping $\alpha=1.25$ and $\epsilon=0.04$, the early-warning tipping time close to $35$ and $7$ for $\Delta F_{0}=19$ and $21$, respectively. This means that the Arctic sea ice will melt much earlier under influence of the enhanced greenhouse effect.

For $\Delta F_{0}=21$, the ice-free time $T_{2}\equiv 1$, because the system keep in the perennially ice-free state. In this case, the Arctic Ocean is ice-free all over the year. Therefore, the ice-free time $T_{2}$ can not capture the dynamical behaviour in this case. Next, we will introduce the other quantity, which can help us to predict the approximate time when the thickest ice sheet of the Arctic region in a year will disappear. The disaster-happening tipping time T3 will be introduced to express the last time when the Arctic Ocean has the winter sea ice, and after that time, the Arctic Ocean remains ice-free throughout the year, as shown in Fig. 5. The time $T_{3}$ is a signal of disaster that may happen in Arctic sea ice with unimaginable consequences–from the loss of polar bear habitat to possible increases of weather extremes at mid-latitudes Kerr (2012). We hope that the $T_{3}$ never shows up.

Fig. 6 shows the numerical results of the disaster-happening tipping time $T_{3}$ effect of the noise intensity $\epsilon$ and Lévy index $\alpha=1.25,~{}1.50,~{}1.75$. We find that $T_{3}$ shows a U-shaped changing trend with the increasing of value of $\epsilon$ for the case of non-Gaussian Lévy noise. In the beginning, $T_{3}$ decreases rapidly as $\epsilon$ increases. Then, the effect of the increasing $\epsilon$ on $T_{3}$ is relatively weak. Finally, $T_{3}$ will turn a corner when the curve reaches a nadir, after this inflection point, $T_{3}$ rapidly grows with increasing of $\epsilon$. Furthermore, the value of $\alpha$ is larger, the growth rate of $T_{3}$ is slower. This means that it will delay the time for the appearance of the perennially ice-free state when the intensity of extreme events increases to a certain degree. Meanwhile, the stronger the frequency of extreme events is, the slower the delay time will be. A possible explanation for these inflection points could be the interaction between the ice albedo feedback and the greenhouse effect. We note that for the case of Gaussian noise, relations between the disaster-happening tipping time $T_{3}$ and the value of noise intensity are similar to the non-Gaussian Lévy noise, while the appearance of inflection point requires the stronger noise intensity.

C. Combination of Lévy noise intensity $\epsilon$ and Lévy index $\alpha$ trigger transition

For $\Delta F_{0}=19$ and a fixed value of $\alpha\in[1,2)$, Fig. 7 shows the minimum $\epsilon$ that can trigger transition from the perennially ice-covered state to the seasonally ice-free state, which is marked with red solid points. For example, for $\alpha=1.0,1.75$, the minimum noise intensity that can trigger the transition is $\epsilon\approx 0.03,0.051$, respectively. The larger $\alpha$ is, the larger noise intensity is needed. We collect all the combinations of $\alpha$ and $\epsilon$ that can induce the transitions for Arctic sea ice system in the region I .

Similarly, the combinations of $\alpha$ and $\epsilon$ leading to the transition from the seasonally ice-free state to the perennially ice-free state can be obtained for $\Delta F_{0}=21$, as shown in Fig. 7. For example, for $\alpha=1.5$, the minimum value of noise intensity is $\epsilon\approx 0.02114$ and the maximum value of noise intensity is $\epsilon\approx 0.055$. As shown in Fig. 7, the blue curve and red curve denote an upper bound and a lower bound of $\epsilon$, which correspond to the maximum and minimum noise intensities inducing the transition, respectively. These two curves constitute the region II where the system can shift to a perennially ice-free state from a seasonally ice-free state .

Arctic sea ice in recent summers shows the record lows in ice extent Stroeve et al. (2005); Walker (2016). It has also been thinning at a remarkable rate over the past few decades Bitz and Roe (2004); Muller-Stoffels and Wackerbauer (2011). To insight into the Arctic sea ice variations under influence of extreme weather events, we propose a nonautonomous Arctic sea ice model under the non-Gaussian $\alpha$-stable Lévy noise. We use the maximal likely trajectory, based on the numerical solution of the nonlocal Fokker-Planck equation (2), to obtain the tipping times for the stochastic Arctic sea ice model (6). The early-warning tipping time $T_{1}$ and the disaster-happening tipping time $T_{3}$ are used to predict the time when the Arctic Ocean may appear in ice-free state in summer and in winter, respectively.

By numerical experiments, we find that the tipping times $T_{1}$ and $T_{3}$ depend strongly on the Lévy index $\alpha$, noise intensity $\epsilon$ and greenhouse effect $\Delta F_{0}$. For example, for the enhanced greenhouse level $\Delta F_{0}=21$, the early-warning tipping time $T_{1}$ decreases with the increased intensity of $\alpha$-stable Lévy noise, in agreement with the corresponding results for the weakened greenhouse level $\Delta F_{0}=19$, but $T_{1}$ comes earlier for $\delta F_{0}=21$. This implies that Arctic sea ice will melt much earlier in summer under influence of the enhanced greenhouse effect and increased noise intensity. On the other hand, the disaster-happening tipping time $T_{3}$ shows a U-shaped changing trend under non-Gaussian Lévy noise with the increase of value of $\epsilon$. Another observation from the results is that for the ice-free time $T_{2}$, we see that the decreased intensity and frequency of extreme events will be effective in reducing the length of time in the fractionally ice-free state.

Comparing with the non-Gaussian Lévy noise, we discover that a larger noise intensity is needed to induce a transition between metastable states in Gaussian case (see Fig.4 and Fig.6). Finally, we identify the combinations of $\alpha$ and $\epsilon$ that trigger the transitions from one state to the other one in the Arctic sea ice system under non-Gaussian Lévy noise for both the weakened and the enhanced greenhouse gas levels.

Our work provides a theoretical framework for studying and predicting the Arctic sea ice variations under the influence of extreme events.

It is well known that Brownian motion has the stationary, the independent increments and almost surely continuous sample paths. A Lévy process is a non-Gaussian process with heavy-tailed statistical distribution and intermittent bursts. The stable distribution $S_{\alpha}(\delta,\beta,\lambda)$ is determined by the following four indexes: Lévy index $\alpha\in(0,2)$, scale parameter $\delta\in[0,\infty)$, skewness parameter $\beta\in[-1,1]$ and shift parameter $\lambda\in(-\infty,\infty)$ Duan (2015).

A scalar symmetric $\alpha$-stable Lévy process $L^{\alpha}_{t}$ is defined via the following properties:

• (i)

  $L^{\alpha}_{0}=0$, almost surely (a.s.);

• (ii)

  $L^{\alpha}_{t}$ has independent increments: For every natural number $n$ and positive time instants with $t_{0}<t_{1}<\cdots<t_{n}$, the random variables $L^{\alpha}_{t_{1}}-L^{\alpha}_{t_{0}}$, $L^{\alpha}_{t_{2}}-L^{\alpha}_{t_{1}}$, $\cdots$, $L^{\alpha}_{t_{n}}-L^{\alpha}_{t_{n-1}}$ are independent;

• (iii)

  $L^{\alpha}_{t}$ has stationary increments: $L^{\alpha}_{t}-L^{\alpha}_{s}$ and $L^{\alpha}_{t-s}$ have the same stable distribution $S_{\alpha}((t-s)^{\frac{1}{\alpha}},0,0)$, for $t>s$;

• (iv)

  $L^{\alpha}_{t}$ has stochastically continuous sample paths: $L^{\alpha}_{t}\rightarrow L^{\alpha}_{s}$ in probability, as $t\rightarrow s$.

An $\alpha$-stable Lévy process $L^{\alpha}$ has the following “heavy tail” estimate Zheng et al. (2020):

i.e., the tail estimate decays in power law. The constant $C_{\alpha}$ is a gamma function depending on $\alpha$.

For a symmetric $\alpha$-stable Lévy process $L^{\alpha}_{t}$, the skewness index $\beta$ is equal to $0$. The paths of $L^{\alpha}_{t}$ have countable jumps, and the jump measure $\nu_{\alpha}(dy)$ depends on $\alpha$,

For $0<\alpha<1$, $L^{\alpha}_{t}$ has larger jumps with lower jump probabilities, while for $1<\alpha<2$, it has smaller jumps with higher jump frequencies. The special case for $\alpha=2$ corresponds to a Brownian motion, and it is a Cauchy process for $\alpha=1$. More information on $\alpha$-stable Lévy process, see references Saltzman (2002); Duan (2015).