A non-autonomous framework for climate change and extreme weather events increase in a stochastic energy balance model

The purpose of this paper is twofold: to describe a scheme for a three-timescale non-autonomous dynamical system suitable for the description of climate change and to present a new example of an energy balance model (EBM) with spatial structure. This model fits into our abstract framework and may contribute to a conceptual explanation of the increased frequency of extreme events. The non-autonomous nature of this example is crucial, making its connection to the foundational scheme particularly interesting.

It is essential to emphasize two key aspects. First, the three-timescale scheme is designed to motivate the structure of the proposed EBM, including the adiabatic $CO 2$ variable, the presence of noise, and the system’s timescale. Second, the EBM in question is a simplified climate model, representing the lowest level of complexity among climate models. As such, it can only provide qualitative insights into climate dynamics. Therefore, we deliberately avoid making quantitative comparisons between our model’s results and real-world data.

From the viewpoint of climate change analysis, our goal is to provide a simplified model that demonstrates how temperature fluctuations increase with rising carbon dioxide (⁠ $CO 2$⁠) concentrations. Probably, the most dramatic example of climate change today is the increased frequency of extreme events. As pointed out by the IPCC Report, Ref. 1, this can result from both an increase in the mean temperature and an increase in the variance of temperature, as shown in Fig. 1, not only by the increase in the mean value but also by the increase in the variance of a climatic variable.

This phenomenon is not just speculative and can be observed by looking at real data, as can be done in Fig. 2(a). It shows the histogram of the daily mean temperature recorded in August in two different periods, from $1910$ to $1940$ in blue, and from $1993$ to $2023$ in red. In our opinion, being able to explain the joint increase in mean value and variance is a challenging and, at the moment, a completely open question. Indeed, a classical paradigm to explain the increase in the frequency of extreme events is given by a dynamical system approaching a bifurcation point, a system subject to random fluctuations, which are amplified near the bifurcation point.^2–8 But the main question then is whether the Earth, nowadays, is close to a bifurcation point or not, an issue which is speculated but not proven. Looking to the far past, it is rather clear that bifurcations took place connecting a moderate climate like our own today with glacial climates. But a clear indication that now we are close to a new bifurcation point in the direction of a warmer climate is missing. Except for certain data near the Tropics: localized in space at that latitude, there is experimental evidence of a potential bistability and, therefore, of the possibility of a fast transition to a warmer climate.⁹ But far from the Tropics, the data are different. Therefore, we have developed a space-dependent model that incorporates this difference. Seen globally, as an infinite dimensional dynamical system on the whole globe, the Tropics bistability does not add a new bifurcation to the model.

We focus on the class of EBMs, which give an elementary, yet useful, representation of Earth’s climate by capturing the fundamental mechanisms governing its behavior.^7,10–15 This type of models, which describes temperature evolution, assumes that it evolves according to the radiation balance of the budget, i.e., the difference in the radiation absorbed and emitted by the planet. Other key factors, such as insolation, atmospheric composition, and surface properties, are considered in EBMs in order to get a radiation budget as accurate as possible.

More in detail, we work with a Budyko–Sellers one-dimensional energy balance model (1D EBM), in which a diffusion term modeling meridional heat transport is added as a driver of temperature evolution.¹⁶ The novelty of our model consists in adding, for the first time, a particular phenomenon that may happen at the Tropics. Indeed, one of the key mechanisms to stabilize Earth’s temperature is the Planck feedback. It consists of an increase in outgoing longwave radiation (OLR) as surface temperature increases, thanks to the Stefan–Boltzmann law. But there exist cases, and our model highlights one of them, where this feedback can fail, as in the super-greenhouse effect (SGE).^9,17,18 This phenomenon is feedback between water vapor, surface temperature, and greenhouse effect. Once the sea surface temperature or greenhouse gas concentration reaches a certain level, the increase in absorbed thermal radiation from the surface due to augmented evaporation with rising sea surface temperature (SST) outweighs the concurrent elevation in OLR, causing OLR to decline as SST increases.

From the mathematical viewpoint, the EBM proposed here is a non-autonomous stochastic partial differential equation (SPDE) of reaction–diffusion type.^19,20 The non-autonomy stands in a slowly varying term corresponding to $CO 2$ concentration. We clarify the link between such a scheme and a more classical adiabatic approach, where the slowly varying term is kept constant and the long-term behavior of the corresponding autonomous dynamical system is investigated. It is worth pointing out that the timescale of our model is intermediate between the daily timescale of weather and the centurial timescale of climate. It is sometimes called macroweather, and it is of the order of a month to one year. For this reason, the abstract framework of a two-timescale separation of a non-autonomous dynamical system is introduced to describe the macroweather-climate dichotomy, interpreting the weather as a faster timescale than macroweather which is averaged out by looking at the temperature on time interval of months.

This paper is organized as follows. In Sec. II, we introduce the non-autonomous framework for weather, macroweather, and climate, with a particular focus on the latter two. We outline the scale separation between them and provide their tentative definitions. This first, abstract part culminates in a link between macroweather and climate, justifying an adiabatic approach to studying the former in the presence of a non-autonomous dependence, such as $CO 2$ concentration in the atmosphere. In Sec. III, we begin by recalling the key concepts of EBMs, starting with the zero-dimensional version in Sec. III A. We discuss why EBMs are suitable for describing macroweather and how they can explain the increase in global mean temperature (GMT) due to rising $CO 2$ concentrations. In Sec. III B, we introduce our new 1D EBM, detailing the parametrization of all the terms of the parabolic partial differential equation (PDE). Section III C describes the deterministic properties of the model, such as the existence of one or more steady-state solutions and their dependence on the $CO 2$ parameter. In Sec. III D, we discuss the stochastic extension of the model and its properties, including the existence of an explicit invariant measure. Section III E presents the main results of the second part. We demonstrate how our model predicts an increase in the variance of temperature fluctuations under $CO 2$ in the current climate configuration, without introducing new bifurcation points. We use both numerical tools and theoretical arguments to explain and understand this phenomenon. Finally, Sec. IV summarizes our findings, while  Appendixes A and  B detail the finite-difference scheme used for numerical simulations and some spectral properties needed to deduce the existence of the invariant measure for the stochastic problem.

In Sec. III, we introduce a stochastic EBM with suitable space dependence, and a slowly varying parameter corresponding to $CO 2$ concentration, in order to investigate the time change of fluctuations, as also discussed in Sec. I. There are three different timescales involved in this modeling. We could skip the first one by just mentioning Ref. 21 and restrict ourselves to two timescales only, but at the same time we want to insist on the non-autonomous structure of the modeling, hence it is convenient to enlarge the discussion a bit.

In the modeling we have in mind, there are three timescales, called

1. weather, where variations are visible at the timescale of hours/day [variables will be denoted by $V w ( t )$⁠, $T w ( t )$ etc.];

2. macroweather, where variations are visible at the timescale of months/year [variables will be denoted by $V m w ( t )$⁠, $T m w ( t )$ etc.];

3. climate, changing at the timescale of dozens of years (probability measures and their expected values characterize this level, still changing in time).

Subdivision and attribution of a precise timescale are not strict.

The parameter $τ= 1 ϵ$ corresponds to the typical reaction time of the slow variables, measured in the unitary time of $V$⁠. Appreciable changes of $T$ happen in a time of order $τ$⁠, at the weather scale.

With great simplification, we may think that $V$ collects the fluid dynamic variables (the fluid Velocity plus other related variables), which are very unstable and rapidly changing at the daily level, while $T$ represents Temperature. In this case, the unit of time at the weather level is of the order of hours, and $τ= 1 ϵ$ is of the order of a few months, hence, e.g., of order 100. On the contrary, the time change of $CO 2$ concentration, call it $τ CO 2$⁠, is of the order of a dozen of years, hence, e.g., of order 10 000 in the weather scale. With these figures, $q ( t )$ has a relaxation time of $100$ and $q ( ϵ t )$ of order $10000$⁠.

Let us heuristically describe the averaging approximation, which can be made rigorous under proper assumptions for suitable systems, see Ref. 22.

For our purposes below, adhering to Hasselmann’s proposal is essential, since our results are the consequence of random perturbations of a nonlinear system representing climate dynamics at the macroweather timescale, namely, a stochastic version of the EBM. Random perturbations are often accepted just based on the generic justification of an unknown coupling with other segments of the physical system (which at the end of the story is the reason also here, namely, the coupling with the fast variables) but Hasselmann’s proposal is a more precise explanation.

Concerning precisely the concept of climate, let us introduce some formalism. We again limit ourselves to stochastic differential equations (SDEs) on an Euclidean space $R d$ (e.g., the space-discretization of a stochastic partial differential equation, as in our main example below) but the concepts can be widely generalized, see, for instance, Ref. 25 for a non-autonomous abstract random dynamical system framework related to the weather-climate dichotomy (that we improve hereby introducing a third level, the macroweather). Call $Pr ( R d )$ the set of probability measures on Borel sets of $R d$⁠. Consider the SDE (1) on the full real line of time. Assume that $W ( t )$ is a $d$-dimensional two-sided Brownian motion defined on a Probability space $( Ω , F , P )$ with expectation $E$⁠. Assume that, for the Cauchy problem on the half line $[ t 0 , ∞ )$ with initial condition $T 0$ at time $t 0$⁠, with arbitrary $t 0$ and $T 0$⁠, at least weak global existence and uniqueness in law holds, and denote the solution by $T t 0 , T 0 ( t )$⁠, $t∈ [ t 0 , ∞ )$⁠; assume $T 0↦ T t 0 , T 0 ( t )$ is Borel measurable from $R d$ to $Pr ( R d )$ endowed with the weak convergence of measures.

Under suitable assumptions for the stochastic equations, there are results of existence (easy, in particular, relying on the existence of a compact global attractor) and also uniqueness (more difficult) for such invariant families. This part of the theory is in progress. When the invariant measure ${ μ t } t ∈ R$ is unique, we call it “the climate.”

If $q ( t )$ varies very slowly, we expect that also $μ t$ varies accordingly, and thus an adiabatic approach to the numerical computation of $μ t$ is reasonable, as already remarked above.

It is crucial to emphasize the following point. We believe that viewing climate as a time-dependent invariant measure, constructed in the pull-back sense, is not only a rigorous definition but also a physically meaningful one. Indeed, the current climate is the result of a long-term evolution that began in the distant past, where the dependence on the initial condition has been lost. This idea, together with the application to geophysical science of concepts from dynamical systems theory, developed in the 1990s,^26,27 began gaining traction approximately 15 years ago.^28,29

EBMs are elementary climate models where, in the simplest form, the temperature of the planet evolves according to the balance of the radiation absorbed and emitted by the Earth.^7,10–13 Their ability to capture the essential dynamics of Earth’s climate system while remaining computationally tractable makes EBMs valuable tools for understanding a wide range of climate phenomena, from the onset of ice ages to the impacts of greenhouse gas emissions on global temperatures.^30,31 There exists a spectrum of complexity for this kind of model, starting from the zero-dimensional (0D) case, moving to the one-dimensional (1D) case, and arriving at higher-dimensional models.¹⁶ Before delving into the details of our 1D EBM with space heterogeneous radiation balance, we illustrate (i) why an EBM has a macroweather timescale, and (ii) why an increase in $CO 2$ concentration in a stochastic 0D EBM leads to an increase in GMT, but not the variance of the solution.

In this section, we propose a 1D EBM with space-dependent radiation balance with local bistability in the outgoing radiation term. The new term does not lead to the addition of a new bifurcation in the model, even if results in a non-linear change in the global mean temperature with respect to $CO 2$ concentration, in comparison with a linear increase in the case of the non-space dependent emitted radiation case. However, by adding a noise component to the model, we detect an increase in fluctuations over time for those values of $CO 2$ concentration in which the non-linear behavior of GMT is triggered. We connect the increase in fluctuations over time, that we denote time variance, to a local (in space) notion of relaxation time. The latter indicator in a sense gives information about the local stability of a space point for, in our case, a global stable temperature configuration.

Third, we describe the modeling of the radiation emitted $R e$⁠, which is the main innovation of our model. Local bistability in the Outgoing Longwave Radiation (OLR) may arise in tropical regions due to a positive feedback loop involving surface temperature and moisture. This feedback, above a threshold of sea level pressure and GHG concentration, causes the atmosphere to become optically thick, thus reducing OLR.⁹ 

Lastly, the additive positive parameter $q$ models the effect of $CO 2$ concentration on the energy budget.^31,39 A higher $q$ leads to lower radiation emitted back to space from the surface. We emphasize the reason behind considering a constant value for the $CO 2$ parameter $q$⁠. Greenhouse gas concentration can be regarded as constant at a macroweather, i.e., monthly, timescale, by avoiding seasonality insertion since it evolves on a much larger timescale. For instance, the $CO 2$ concentration has increased by around 47% from 1850 to 2020, moving from 284 parts per million (ppm) to 412 ppm.⁴⁰ In a first approximation, we assume that the spatial distribution of greenhouse gases is uniform over the globe. This assumption, although no longer the state of the art, was widespread at the turn of the century and is based on the well-mixing property of GHG. This means that since most GHGs, such as $CO 2$⁠, have a large lifetime, many studies have been conducted using global average values for the spatial distribution.^1,39,41

In this section, we describe the deterministic properties of our model as the number of steady-state solutions and their stability. The new feature of our model is the rise of a strong non-linear increase in GMT with respect to the greenhouse parameter $q$⁠, in the proximity of the model configuration describing the actual climate of the Earth.

Note how the second part of the previous result gives qualitative information on the behavior of the GMT. Furthermore, if the diffusion coefficient is sufficiently large and the 0D EBM, obtained by removing the diffusion term and averaging in space the radiation balance, is bistable, it can be rigorously proven the existence of a second stable steady-state solution and a third unstable one.³² 

Given the non-linear nature of the problem, a rigorous demonstration of all properties of the model is not always possible. In this case, the problem can be overcome by using numerical simulations. In particular, for each fixed $q$⁠, we numerically simulated the solutions of Eq. (12). As $q$ varies, we obtained, according to the large literature on this kind of model, that there can be either 1 or 3 stationary solutions. In the former case, the solution is stable. In the latter case, two solutions are stable, one describing a “snowball” configuration $u S$ for Earth’s temperature with ice all over its surface. The other describes a warm climate $u W$⁠, similar to the one in which we are living. Additionally, an unstable solution $u M$⁠, whose average global mean temperature (GMT) lies between the GMT of $u S$ and $u S$⁠, also arises. See Fig. 4(a) for a graphical representation of the steady-state solutions in the bistable case.

In general, proving theoretically the results for such kinds of models is non-trivial due to the presence of nonlinear terms in the energy budget parametrization. The bifurcation diagram of the model in the $(q, u ¯ ∗)$ plane, where $u ∗$ denotes a solution of the elliptic PDE, and $u ¯ ∗= 1 2 ∫ − 1 1 u ∗(x)dx$ is its GMT, is depicted in Fig. 4(b). We highlight how the addition of a space-dependent OLR with tropics bistability does not introduce a new bifurcation. Thus, no new steady-state solutions are added or the stability of the already existing ones is altered. However, as highlighted in Fig. 4(c), a strong non-linear behavior of the GMT for the warm solution $u W$ appears for the value of $q$ in the neighborhood of $q≈11.3$⁠, a thing that does not happen if we remove the bistability in $R e$⁠.³² We aim to focus on that phenomenon, when a noise component, describing the weather effect, is added to the model. This is the main topic of Secs. III D and III E.

As discussed in the previous paragraphs, the EBM presented in this work corresponds to a macroweather timescale. However, as far as it has been introduced, it lacks in taking into account the effect of fast components of Earth’s system, such as atmospheric pressure, wind, and precipitation. In the literature, this has been achieved elementary by adding a stochastic term, such as white noise, to the radiation budget.^21,45–47 Note that this addition is necessary to describe the increase in the frequency of rare events. In fact, while a deterministic EBM is useful for obtaining information on global mean temperature, it is not able to investigate fluctuations around it due to all phenomena, such as weather, not included in the radiation balance of the model.

It is widely acknowledged that greenhouse gas emissions from human activities have led to more frequent and intense weather and climate extremes since the pre-industrial era, particularly temperature extremes.⁴⁰ Despite numerous definitions proposed to assess what constitutes an extreme event, there is currently no universally accepted definition.⁴⁹ One commonly used definition considers an extreme weather event as one that exceeds a predefined threshold for a climate variable.

As it is well understood that such occurrences become more likely as the variance of that climate variable increases, we consider the time variance as a proxy indicator of extreme weather events for our EBM setting. As explained in Sec. III D, an explicit formula for the invariant measure of the stochastic EBM (4) exists. However, due to the presence of a non-linear term inside the Gibbs factor in Eq. (18), obtaining theoretical information is challenging. Thus, we rely on numerical simulations to capture the behavior of the variance.

We observe two distinct behaviors of the variance. First, depending on $q$⁠, there is an increase in $σ t 2$⁠, peaking around $q≈11.3$⁠, followed by a decrease. This peak corresponds to the $q$ value that results in the largest increase in GMT, as shown in Fig. 4(c). Second, as expected, for a given value of $q$⁠, the spatial profile of $σ t 2$ is symmetric with respect to $x=0$⁠. Additionally, three local maximum points emerge one at $x=0$ and the others at $|x |≈0.8$⁠. Our interpretation is that $σ t 2$ identifies the highest fluctuating regions as those where the freezing water temperature is crossed (sub-arctic regions, $|x |=0.8$⁠), or where the SGE triggers bistability (tropical areas, $x=0$⁠).

To understand why, given $x∈[−1,1]$⁠, the map $q↦ σ t 2(x)= σ t 2 , ( q )(x)$ increases until around $q≈11.3$ and then decreases, we offer the following explanation. As $q$ increases, the temperature $u W(x)$ of the warm climate increases, as can be checked by numerical simulations and partially understandable from Proposition 1. This leads the temperature, especially in the tropical area, to approach the region where instability due to the SGE arises. Once the tropical temperature surpasses this region, the instability, and thus the variance, decreases. We support our statements as follows.

• Figure 6 shows how the steady-state solution $u W ( q )$ approaches unstable areas, defined by points $(x,u)$ where $∂ u R ( x , u ) > 0$⁠. For $q 1=11.21$⁠, $q 2=11.28$⁠, and $q 3=11.4$⁠, the map $(x,u)↦ ∂ uR(x, u W(x))$ is presented, with the curve $x↦(x, u W(x), ∂ uR(x, u W(x))$ depicted in red. When the time variance peaks (i.e., $q= q 2$⁠), the steady-state solution values around the equator are close to the local maximum of $∂ uR$⁠, approximately close to $( x M, u M)=(0,305)$⁠. For $q= q 1$ and $q= q 3$⁠, the equatorial steady-state value is smaller and larger than $u M$⁠, respectively.

• We consider the local mean stability indicator

  $γ ¯(q)= ∫ − 1 1γ(x)dx= ∫ − 1 1 ∂ uR(x, u W ( q )(x))dx.$

  Its plot is shown in Fig. 7(a). Although providing only average information, it exhibits behavior qualitatively similar to that of $σ t 2$ for each fixed value of $x$⁠, peaking around $q≈11.3$⁠.

Finally, we show the distribution of the solution $u(x,t)$ of the stochastic 1D EBM (4) for the fixed space point $x=0.6$ (similar results are observed for different space points). Figure 7(b) shows the distribution as a histogram, comparing two distinct values of $q$⁠. Blue depicts the results for $q=11.2$⁠, while red shows $q=11.3$⁠. From the plot, we deduce that our model captures the shift in the mean value and the increase in variance, correlating with the IPCC schematic climate change reproduction in Fig. 1 and weather observations in Modena in Fig. 2(a).

In the first part of this work, we have presented a non-autonomous framework to describe the scale separation between weather, macroweather, and climate. According to Hasselmann’s proposal, weather is conceptualized as arising from a set of deterministic equations, describing variables such as temperature on a fast timescale, typically from an hour to one day. Macroweather, on the other hand, encompasses variations that are neither as short-term as weather nor as long-term as climate. We assume that temperature on a macroweather timescale satisfies a stochastic equation, reflecting the most original aspect of Hasselmann’s work. Additionally, at the macro weather timescale, we include a non-autonomous term to model the atmospheric $CO 2$ concentration, which evolves on a timescale of years. Finally, we identify climate with the invariant measure arising from the stochastic equation at the macro weather level.

In the second part, we have introduced a new 1D EBM on a macroweather timescale, which can predict the increase in the number of extreme weather events associated with climate change. The novelty of our model relies on the parametrization of the non-linear space-dependent radiation budget. Indeed, we have inserted the presence of the SGE, a phenomenon typical of tropical areas that may lead to an instability in the OLR. Our model includes a non-autonomous term $q=q(t)$⁠, that in light of the first part, we have considered constant, that is the effect of the $CO 2$ concentration on the radiation balance. We have recalled the basic mathematical properties of our model, such as the existence of steady-state solutions, using also numerical simulations. Also, we have shown how the GMT of the steady-state solution increases with increasing $CO 2$ concentration. Third, we have the stochastic version of our 1D EBM, obtained by perturbing our model with an additive space-time white noise. Being in the class of stochastic reaction–diffusion SPDE, it is possible to write explicitly the invariant measure as a Gibbs one. However, the presence of the non-linearity arising from the radiation budget prevents us from deducing more information from the invariant measure.

The most important results of our work are presented in Sec. III E. There, we have exploited numerical simulations to study the changes in the invariant measure as the $CO 2$ concentration increases. To obtain this, we have performed numerical integration of the stochastic 1D EBM on a time interval of $500$ years and with the initial condition the warm steady-state solution $u W$ of the deterministic 1D EBM, changing only $q$ in different runs of the simulation. Given $q$⁠, for each space point $x∈[−1,1]$⁠, we associate the extreme weather event frequency with the time variance $σ t 2(x)$ of a trajectory of the stochastic 1D EBM at point $x$⁠. We explain the spatial behavior of $σ t 2(x)$⁠, which presents two local maximum points for $x≈±0.8$ and one for $x=0$⁠, combining heuristic reasoning with empirical indicators. The informal explanation is that the presence of the diffusion term forces the stable warm stationary solution, and the oscillations around it due to noise, to take on values that are not stable if we were to consider the ODE obtained by removing the diffusion term. These regions of temperature, which are locally unstable, correspond to the areas where the variance $σ t 2$ is highest. We also motivated the behavior of $σ t 2(x)= σ t 2 ( q )(x)$ with respect to $q$⁠, at fixed $x$⁠. The variance observed in this way increases to a maximum value and then decreases. This is explained by the increase in GMT due to $q$ and the presence of the diffusion term, leading the solution of the stochastic EBM to be increasingly in, and then out of, the region of instability due to the SGE.

Lastly, we are aware that our model, although based on physical laws, is largely phenomenological. In some parts of our model, we have included simplifications of convenience, such as in the perturbation of the diffusion function, which makes the problem non-singular and thus easier to deal with mathematically; in other parts, such as in the parametrization of the space-dependent and tropic-restored OLR, we have followed the principle of simplicity. Other choices would certainly have been possible, but the elementary nature of the model leads us to avoid choices that are too fine or complex.

Concerning open questions, an important one would be understanding the correct form of the noise term. We have assumed it to be additive and without correlation in space and time. This is an extreme simplification, which allows us to deduce certain properties of the invariant measure. But already Hasselmann in his seminal work proposed a different, i.e., multiplicative, type of noise. In addition to this, understanding how the properties of the model change if more physical processes are considered, such as advection, is a problem we would like to address in the future, as well as the extension of our work to a two-dimensional setting.

We acknowledge fruitful discussions with Paolo Bernuzzi and Giulia Carigi. The Modena temperature data were provided by the Geophysical Observatory of Modena, University of Modena and Reggio Emilia, Italy (www.ossgeo.unimore.it, Ref. 50) and are available upon request from the corresponding author. The research of G.D.S. is supported by the Italian National Interuniversity PhD course in sustainable development and climate change. The research of F.F. is funded by the European Union (ERC, NoisyFluid, No. 101053472). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

The authors have no conflicts to disclose.

G. Del Sarto: Conceptualization (equal); Methodology (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (lead). F. Flandoli: Conceptualization (equal); Methodology (equal); Writing – original draft (equal).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.11609953, Ref. 53.

In this section, we describe the numerical method adopted to approximate the solutions of the stochastic PDE (19). We used an implicit Euler-Maruyama method, which is a small modification of the semi-implicit Euler-Maruyama method presented in Ref. 51, Sec. 10.5.