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

Gianmarco Del Sarto  Scuola Normale Superiore, Pisa, Italy  University School
for Advanced Studies IUSS Pavia, Pavia, Italy  Franco Flandoli  Scuola Normale
Superiore, Pisa, Italy

###### Abstract

We develop a three-timescale framework for modelling climate change and
introduce a space-heterogeneous one-dimensional energy balance model. This
model, addressing temperature fluctuations from rising carbon dioxide levels
and the super-greenhouse effect in tropical regions, fits within the setting
of stochastic reaction-diffusion equations. Our results show how both mean and
variance of temperature increase, without the system going through a
bifurcation point. This study aims to advance the conceptual understanding of
the extreme weather events frequency increase due to climate change.

##  1 Introduction

The purpose of this paper is twofold: to describe a scheme for 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 increase in the frequency of extreme events. The
non-autonomous nature of this example is crucial, making its connection to the
foundational scheme particularly interesting.

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 ($\text{CO}_{2}$) concentrations. Probably the most
dramatic example nowadays of climate change is the increased frequency of
extreme events. As pointed out by the IPCC Report [HDG+01], this can result
from both an increase in the mean temperature and an increase in the variance
of temperature, as shown in Figure 1, not only by the increase of the mean
value but also by the increase in the variance of a climatic variable.

Figure 1: Schematic showing the effect on extreme temperatures when (a) the
mean temperature increases, (b) the variance increases, and (c) when both the
mean and variance increase for a normal distribution of temperature.
Reproduced from the IPCC Report on Climate Change 2001: The Scientific Basis
[HDG+01, Figure 2.32]

This phenomenon is not just speculative and can be observed by looking at real
data, as can be done in Figure 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 of 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 [DSvN+08, SBB+09,
AWVC12, LL+12, GL20]. 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: localised 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 [DG18]. But far from the Tropics, the data are
different. Therefore we have developed a space-dependent model which
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 behaviour [Bud69, Sel69, Nor75, Ghi76, Día97, GL20, CLM+23].
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 solar 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 modelling meridional heat transport
is added as a driver of temperature evolution [NK17]. 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 of 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)
[EWV14, BC16, DG18]. This phenomenon is a feedback between water vapour,
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.

(a)

(b)

Figure 2: (a) Daily temperature recorded in August in Modena, Italy. The blue
histogram shows the values for the time period $1910-1940$ (blue), and the red
histogram refers to $1993-2023$ (red). Data provided by the Osservatorio
Geofisico di Modena www.ossgeo.unimore.it. (b) Observational Outgoing Longwave
Radiation (OLR) dependence on Sea Surface Temperature (SST), and
$\text{SMART}^{1}$ OLR output for various humidity values. The red dashed line
is the mean. Reproduced from [DG18, Figure 2a], with permission from Wiley.

111SMART is a dataset of OLR measurements used in [DG18].  From the
mathematical viewpoint, the EBM proposed here is a non-autonomous stochastic
partial differential Equation (SPDE) of reaction-diffusion type [DPZ96,
DPZ14]. The non-autonomy stands in a slowly varying term corresponding to
$\text{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 behaviour 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 Section 2, 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
$\text{CO}_{2}$ concentration in the atmosphere. In Section 3, we begin by
recalling the key concepts of EBMs, starting with the zero-dimensional version
in Section 3.1. We discuss why EBMs are suitable for describing macroweather
and how they can explain the increase in global mean temperature (GMT) due to
rising $\text{CO}_{2}$ concentrations. In Section 3.2, we introduce our new 1D
EBM, detailing the parametrization of all the terms of the parabolic partial
differential equation (PDE). Section 3.3 describes the deterministic
properties of the model, such as the existence of one or more steady-state
solutions and their dependence on the $\text{CO}_{2}$ parameter. In Section
3.4, we discuss the stochastic extension of the model and its properties,
including the existence of an explicit invariant measure. Section 3.5 presents
the main results of the second part. We demonstrate how our model predicts an
increase in the variance of temperature fluctuations under $\text{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, Section 4 summarizes our findings, while
Appendix A details 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.

##  2 The non-autonomous scheme with three timescales

In Section 3 below we introduce a stochastic EBM with suitable space
dependence, and a slowly varying parameter corresponding to $\text{CO}_{2}$
concentration, in order to investigate the time-change of fluctuations, as
discussed also in the Introduction. There are three different timescales
involved in this modelling. We could skip the first one by just mentioning
[Has76], and restrict ourselves to two timescales only, but at the same time
we want to insist on the non-autonomous structure of the modelling, hence it
is convenient to enlarge the discussion a bit.

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

  * •

weather, or meteo, where variations are visible at the timescale of hours/day
(variables will be denoted by $V_{w}\left(t\right)$, $T_{w}\left(t\right)$
etc.);

  * •

macroweather, where variations are visible at the timescale of months/year
(variables will be denoted by $V_{mw}\left(t\right)$, $T_{mw}\left(t\right)$
etc.);

  * •

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.

###### Remark 1

The reader will realize that, for the purpose of Section 3, we could start
from Subsection 2.5 Let us then explain why we think that Subsections 2.1-2.4
are also very important. As already said, a main purpose of this work is to
identify possible explanations for an increase in variance, when a parameter
changes. In a sentence, the two main ingredients are a noise in the system
equations and a suitable nonlinearity which amplifies the variations of the
noise in a different way for different values of the parameter. The noise,
then, is crucial. Subsections 2.1-2.4 are devoted to explaining its origin.

###  2.1 The weather timescale

Following [Has76], it is natural to model the weather scale by deterministic
equations, ordinary equations for simplicity of notation (but the ideas are
the same for PDEs); randomness can be introduced but it is not strictly
necessary, except maybe for a description of the uncertainty about initial
conditions and parameters, not included in the present discussion. Following
[Has76], we distinguish the main physical variables in fast and slow ones,
according to a system of the form

$\displaystyle\partial_{t}V_{w}$ $\displaystyle=f\left(V_{w},T_{w}\right)$
$\displaystyle\partial_{t}T_{w}$ $\displaystyle=\epsilon
g\left(V_{w},T_{w},q\left(\epsilon t\right)\right)$ $\displaystyle
q\left(t\right)\text{ slowly varying}$

for a small $\epsilon>0$. Here $V_{w}$ changes in unitary time (corresponding
to hours/day), and $T_{w}$ varies very slowly (monthly, say). In addition, the
slow variable is influenced by a slow time-change of structure, described by
the time-varying parameters $q\left(\epsilon t\right)$. The function
$q\left(t\right)$ is assumed to be slowly varying, hence $q\left(\epsilon
t\right)$ is super-slowly varying (from here the three timescales arise).

The parameter $\tau=\frac{1}{\epsilon}$ 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 $\tau$, at the meteo-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 meteo level is of the order
of hours, and $\tau=\frac{1}{\epsilon}$ is of the order of a few months, hence
e.g. of order 100. On the contrary, the time-change of $\text{CO}_{2}$
concentration, call it $\tau_{\text{CO}_{2}}$, is of the order of a dozen of
years, hence e.g. of order 10000 in the meteo scale. With these figures,
$q\left(t\right)$ has a relaxation time of $100$ and $q\left(\epsilon
t\right)$ of order $10000$.

###  2.2 Macroweather timescale for $T$

Then we change the scale and set

$\widetilde{T}_{mw}\left(t\right):=T_{w}\left(\frac{t}{\epsilon}\right)$

so that we observe variations of $\widetilde{T}_{mw}\left(t\right)$ in unitary
time. We call this the macroweather timescale. It holds

$\partial_{t}\widetilde{T}_{mw}\left(t\right)=g\left(V_{w}\left(\frac{t}{\epsilon}\right),\widetilde{T}_{mw}\left(t\right),q\left(t\right)\right).$

Let us recall that, at this timescale, $q$ varies very slowly. Let us look for
a simplification of this equation, where $V_{w}$ does not appear anymore.

###  2.3 The averaging approximation

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

At the integral level, we have

$\widetilde{T}_{mw}\left(t\right)-\widetilde{T}_{mw}\left(t_{0}\right)=\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),\widetilde{T}_{mw}\left(s\right),q\left(s\right)\right)ds.$

If $t-t_{0}$ is small, let us use the reasonable approximation

$\displaystyle\sim\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),\widetilde{T}_{mw}\left(t_{0}\right),q\left(t_{0}\right)\right)ds$
$\displaystyle=\left(t-t_{0}\right)\frac{1}{t-t_{0}}\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),\widetilde{T}_{mw}\left(t_{0}\right),q\left(t_{0}\right)\right)ds.$

Then, if $\epsilon<<t-t_{0}$, we heuristically invoke an ergodic theorem and
approximate

$\sim\left(t-t_{0}\right)\int
g\left(v,\widetilde{T}_{mw}\left(t_{0}\right),q\left(t_{0}\right)\right)\nu_{\widetilde{T}_{mw}\left(t_{0}\right)}\left(dv\right)$

where $\nu_{\tau}\left(dv\right)$ is invariant for

$\partial_{t}V=f\left(V,\tau\right).$

Setting

$\overline{g}\left(\tau,q\right)=\int
g\left(v,\tau,q\right)\nu_{\tau}\left(dv\right)$

we may write

$=\left(t-t_{0}\right)\overline{g}\left(\widetilde{T}_{mw}\left(t_{0}\right),q\left(t_{0}\right)\right)$

and then again approximate it to

$\sim\int_{t_{0}}^{t}\overline{g}\left(\widetilde{T}_{mw}\left(s\right),q\left(s\right)\right)ds.$

Hence we get the simplified model

$\partial_{t}\overline{T}_{mw}\left(t\right)=\overline{g}\left(\overline{T}_{mw}\left(t\right),q\left(t\right)\right).$

###  2.4 Hasselmann’s proposal

However, in our case, this simplification is not realistic. If $\epsilon$ is
of order $\frac{1}{100}$, then $t-t_{0}$ is of order one, because we need the
validity of the approximation

$\frac{1}{t-t_{0}}\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),a,b\right)ds\sim\int
g\left(v,a,b\right)\nu_{a}\left(dv\right).$

But on a time of order one, we observe variations of
$\widetilde{T}_{mw}\left(t\right)$, we said above, hence the approximation

$\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),\widetilde{T}_{mw}\left(s\right),q\left(s\right)\right)ds\sim\int_{t_{0}}^{t}g\left(V_{w}\left(\frac{s}{\epsilon}\right),\widetilde{T}_{mw}\left(t_{0}\right),q\left(t_{0}\right)\right)ds$

is not so strict (on the contrary, it is excellent for the
$q\left(s\right)\sim q\left(t_{0}\right)$ approximation).

We need to keep fluctuations into account at the macroweather scale. A
phenomenological way (Hasselmann proposal) is to replace the model above by

$dT_{mw}\left(t\right)=\overline{g}\left(T_{mw}\left(t\right),q\left(t\right)\right)dt+\sqrt{\epsilon}\sigma\left(T_{mw}\left(t\right),q\left(t\right)\right)\circ
dW\left(t\right)$

for a suitable ”volatility” $\sigma$ (Stratonovich integral $\circ$ looks more
appropriate). In [Has76], heuristic justifications are given, inspired for
instance the random displacements of a Brownian particle in a fluid of
molecules (which on their own are subject to a deterministic fast dynamics,
coupled with the slow deterministic dynamic of the bigger particle). The
Bremen school on Random Dynamical Systems and other research groups explored
for some time rigorous justifications for this proposal, but a final answer is
not known, see for instance [Arn01, Kif01]. However, the observation of
temperature time series at the timescale of month-year clearly shows some form
of stochasticity and thus Hasselmann’s proposal looks very appealing.

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 more
precise explanation.

Let us however advise the reader that we shall start, in our example below,
from a stochastic PDE for the temperature macroweather-scale, given a priori,
not derived precisely from the weather scale as described above. We want to
concentrate on the consequences of particular nonlinearities. Our model will
have the simplified form

$dT_{mw}\left(t\right)=\overline{g}\left(T_{mw}\left(t\right),q\left(t\right)\right)dt+\sqrt{\epsilon}\sigma
dW\left(t\right)$

with constant $\sigma$.

###  2.5 Macroweather and climate

As announced at the end of the last subsection, our investigation starts from
an equation of the form (stochastic differential equation or SPDE)

$\displaystyle dT\left(t\right)$
$\displaystyle=\overline{g}\left(T\left(t\right),q\left(t\right)\right)dt+\sigma
dW\left(t\right)$ (1) $\displaystyle q\left(t\right)\text{ slowly varying}$

where we skip the subscripts but keep in mind that it is a macroweather model,
we have skipped the factor $\sqrt{\epsilon}$ but we shall choose a small
diffusion coefficient $\sigma$, and the nonlinear function $\overline{g}$ will
be chosen by means of typical arguments related to EBM’s. The slowly varying
function $q\left(t\right)$ will describe the effect, in the model, of slowly
varying $\text{CO}_{2}$-concentration, appreciated on a timescale of dozens of
years.

The climate is a collection of statistical information from the time series of
this model. If it were an autonomous system ($q\left(t\right)$ equal to a
constant), we would invoke invariant measures. Due to the time-change in
$\overline{g}$, we have to use the formalism of time-varying invariant
measures. However, at the simulation level, we approximate this time-varying
system by an adiabatic system parametrized by a parameter $q$:

$dT\left(t\right)=\overline{g}\left(T\left(t\right),q\right)dt+\sigma
dW\left(t\right)$

and investigate its invariant measures, parametrized by $q$. The slow change
of statistics for the true non-autonomous system is mimicked by the change of
statistics when the parameter $q$ is changed.

Concerning precisely the concept of climate, let us introduce some formalism.
We again limit ourselves to stochastic differential equations (SDEs) on an
Euclidean space $\mathbb{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 [FPT22] for a non-
autonomous abstract random dynamical system framework related to the weather-
climate dichotomy (that we improve here by introducing a third level, the
macroweather). Call $\Pr\left(\mathbb{R}^{d}\right)$ the set of probability
measures on Borel sets of $\mathbb{R}^{d}$. Consider the SDE (1) on the full
real line of time. Assume that $W\left(t\right)$ is a $d$-dimensional two-
sided Brownian motion defined on a Probability space
$\left(\Omega,\mathcal{F},\mathbb{P}\right)$ with expectation $\mathbb{E}$.
Assume that, for the Cauchy problem on the half line
$\left[t_{0},\infty\right)$ 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}}\left(t\right)$, $t\in\left[t_{0},\infty\right)$; assume
$T_{0}\mapsto T^{t_{0},T_{0}}\left(t\right)$ is Borel measurable from
$\mathbb{R}^{d}$ to $\Pr\left(\mathbb{R}^{d}\right)$ endowed with the weak
convergence of measures.

For all $t_{0}<t$ introduce the Markov semigroup $\mathcal{P}_{t_{0},t}$
mapping $\Pr\left(\mathbb{R}^{d}\right)$ into $\Pr\left(\mathbb{R}^{d}\right)$
defined by the identity

$\int_{\mathbb{R}^{d}}\phi\left(y\right)\left(\mathcal{P}_{t_{0},t}\nu\right)\left(dy\right)=\mathbb{E}\int_{\mathbb{R}^{d}}\phi\left(T^{t_{0},T_{0}}\left(t\right)\right)\nu\left(dT_{0}\right)$

for every $\nu\in\Pr\left(\mathbb{R}^{d}\right)$ and every bounded continuous
test function $\phi$ on $\mathbb{R}^{d}$. One can prove that

$\displaystyle\mathcal{P}_{r,t}\mathcal{P}_{s,r}$
$\displaystyle=\mathcal{P}_{s,t}$ $\displaystyle\mathcal{P}_{s,s}$
$\displaystyle=Id.$

Moreover, one can link the Markov operator for equation (1) to the Fokker-
Planck equation

$\partial_{t}f+\operatorname{div}\left(\overline{g}\left(\cdot,q\right)f\right)=\frac{\sigma^{2}}{2}\Delta
f$

but we do not stress the rigorous results here.

With these notations, a time-varying invariant measure is a family

$\left\\{\mu_{t}\right\\}_{t\in\mathbb{R}}\subset\Pr\left(\mathbb{R}^{d}\right)$

such that

$\mathcal{P}_{s,t}\mu_{s}=\mu_{t}\qquad\text{for every }s\leq t\text{.}$

Under suitable assumptions for the stochastic equations, there are results of
existence (easy) and also uniqueness (more difficult) for such invariant
families. This part of the theory is in progress. When the invariant measure
$\left\\{\mu_{t}\right\\}_{t\in\mathbb{R}}$ is unique, we call it ”the
climate”.

When uniqueness does not hold or it is not known, the idea could be to look
for families $\mu_{t}$ not only invariant but also with additional properties
of interest for physical sciences or other reasons. A typical one could be a
pull-back version of ”convergence to equilibrium”:

$\mu_{t}=\lim_{s\rightarrow-\infty}\mathcal{P}_{s,t}\lambda$ (2)

where $\lambda$ is a ”natural” measure, like a rotation invariant centered
Gaussian measure on $\mathbb{R}^{d}$.

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

##  3 Energy balance models

###  3.1 The macroweather timescale and the global mean temperature increase
due to $\text{CO}_{2}$ concentration

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 [Bud69, Sel69, Nor75, Ghi76, GL20]. 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 [DF11,
BDvdH22]. 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 [NK17]. 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 of $\text{CO}_{2}$ concentration in a stochastic 0D EBM leads to an
increase of GMT, but not the variance of the solution.

First, we consider a 0D EBM for the global mean temperature $T=T(t)$ which
given a positive initial condition $T_{0}$, is an ODE of the form

$\begin{split}C\frac{dT}{dt}&=Q_{0}\beta+q-A-B\cdot(T-273),\\\
T(0)&=T_{0}.\end{split}$ (3)

In this model, the radiation emitted $R_{e}$ by the planet is assumed,
according to Budyko empirical radiation formula [Bud69], of the form

$R_{e}(T)=A+B\cdot(T-273)-q,$

where $A,B$ are positive constants that can be derived by a best-fit estimate
with real data observations; $C>0$ denotes the heat capacity per square meter,
while $Q_{0}$ and $\beta$ are respectively the global mean radiation and
coalbedo. Lastly, the additive parameter $q>0$ the effect of $\text{CO}_{2}i$
concentration on the radiation balance [NK17, BDvdH22, DSBFK24]. Denoting by
$T_{*}$ the unique stable fixed point of the model, i.e.

$T_{*}=\frac{Q_{0}\beta-A+q}{B}+273,$

the solution of (3) is given by

$T(t)=T_{*}+\left(T_{0}-T_{*}\right)e^{-t/\tau_{0}},$

where $\tau_{0}=C/B$ is the relaxation time, i.e. the timescale at which a
deviation from the equilibrium temperature is reabsorbed. Considering an all-
land planet, it is reasonable to consider the heat capacity as half the heat
capacity at constant pressure of the column of dry air over a square meter, as
pointed out in [NK17], leading to

$C=5\cdot 10^{7}\text{ J}\text{ K}^{-1}\text{ m}^{-2}.$

On the other hand, satellite data suggest [NK17, GLN93]

$B=1.90\text{ W }\text{m}^{2}\text{ K}^{-1}.$

This leads to a relaxation time of the other of one month, as

$\tau_{0}=\frac{C}{B}=2.63\cdot 10^{7}\text{ s}\approx 30\text{ days}.$

It is worth pointing out that the hypothesis of an all-land planet is a huge
simplification of reality. Indeed, the Earth’s system has various components
capable of storing heat efficiently, each with its own unique capacity [GL20].
The value for $C$ we have considered corresponds to the capacity of the
atmospheric column. But, even considering a planet with a mixed-layer only
ocean, the heat capacity would be $60$ times larger ([NK17]), i.e. in the
order of a few years, remaining thus in the macro weather timescale.

Second, we force the model with a stochastic noise modelling the effect of
fast terms with respect to the slow radiation balance terms of equation (3).
Denote by $(W_{t})_{t\geq 0}$ a Brownian motion, and consider the SDE given
by:

$\begin{split}CdT_{t}&=\left(Q_{0}\beta+q-A-BT_{t}\right)dt+\sigma dW_{t},\\\
T(0)&=T_{0},\end{split}$ (4)

where $\sigma>0$ is the noise intensity. Denoting by
$\tilde{A}:=Q_{0}\beta+q-A,$ the solution can be explicitly written, using a
variation of parameters technique ([Bal17]) as

$T_{t}=T_{0}e^{-t/\tau_{0}}+\frac{\tilde{A}}{\tau_{0}}\left(1-e^{-t/\tau_{0}}\right)+\sigma\int_{0}^{t}e^{-(t-s)/\tau_{0}}dW_{s}.$

The solution is a Gaussian process with a mean value

$\mathbb{E}\left[T_{t}\right]=T(0)e^{-t/\tau_{0}}+\tilde{A}\left(1-e^{-t/\tau_{0}}\right),$
(5)

and variance

$\begin{split}Var(T_{t})&=Var\left(\frac{\sigma}{C}\int_{0}^{t}e^{-(t-s)/\tau_{0}}ds\right)=\frac{\sigma^{2}}{C^{2}}\int_{0}^{t}e^{-2(t-s)/\tau_{0}}ds\\\
&=\frac{\sigma^{2}}{2BC}\left(1-e^{-2t/\tau_{0}}\right).\end{split}$ (6)

Further, the stochastic EBM in (4) has a unique gaussian invariant measure
$\nu\sim\mathcal{N}(\mu_{\nu},\sigma^{2}_{\nu})$, whose mean $\mu_{\nu}$ and
variance $\sigma^{2}_{\nu}$, can be obtained taking the limit for the time
that tends to infinity in (5) and (6), leading to

$\mu_{\nu}=\tilde{A}=\frac{Q_{0}\beta+q-A}{B},\qquad\sigma^{2}_{\nu}=\frac{\sigma^{2}}{2BC}.$

Thus, a change in the $\text{CO}_{2}$ concentration leads to a change in the
mean value of the climate, i.e. the invariant measure, but not in its
variability. Usually, the variability increase results from a critical
transition, such as a saddle-node bifurcation, which arises in the model. This
is the concept of bifurcation-induced tipping point, which leads to the
critical slowing down behaviour of the system, resulting in an increase of
variance and autocorrelation close to the bifurcation point [DSvN+08, SBB+09,
AWVC12, LL+12]. However, the presence of a bifurcation point for the global
scale climate is questionable, and the presence of bistable regimes for
climate components, sometimes called tipping elements such for Atlantic
meridional overturning circulation or polar ice sheets, is localized in space.
For all these reasons, in the next section, we will describe a new one-
dimensional model with a local in-space change in the non-linear term that is
able to explain the variance increase.

###  3.2 Model formulation

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 w.r.t. $\text{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 of fluctuations over time for those values
of $\text{CO}_{2}$ concentration in which the non-linear behaviour 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.

The main characteristics of a 1D EBM are that the temperature $u=u(x,t)$,
depending on time $t$ and space $x=sin(\phi)$, where $\phi\in[-\pi/2,\pi/2]$
denotes latitude, are that it evolves according to the diffusion of heat, and
the planet energy balance [Bud69, Sel69, Ghi76, CLM+23]. Our model, which is
an extension of the one considered in [DSBFK24], assumes that the temperature
$u$ satisfies the non-linear parabolic partial differential equation (PDE)

$\begin{split}C_{T}\partial_{t}u&=\partial_{x}\left(\kappa(x)u_{x}\right)+R_{a}(x,u)-R_{e}(x,u;q),\quad
x\in[-1,1],t\geq 0\\\ u(x,0)&=u_{0}(x),\quad x\in[-1,1],\\\
u_{x}(-1,t)&=u_{x}(1,t)=0,\quad t\geq 0.\end{split}$ (7)

The heat capacity $C_{T}=5\cdot 10^{-7}\text{ J}\text{ K}^{-1}\text{ m}^{-2}$
is considered uniform over the whole planet, which we assume is an all-land
planet, as the one presented at the beginning of Section 3. The PDE is a
reaction-diffusion type [Tem97, Smo12]. In the following, we are going to
describe the terms governing its time evolution, while the values of the
constants appearing in the parametrizations can be found in Table 1.

In scientific modelling, heat diffusion is often depicted assuming the planet
as a thin shell. This leads to a non-constant diffusion coefficient of the
form $\kappa(x)=D\cdot(1-x^{2})$, with $D>0$, where the term $1-x^{2}$ arises
due to the spherical setting [NK17, Chapter 5]. However, this choice
introduces difficulties for the mathematical treatment, resulting in
degenerate problems [Flo14, CMM20]. For instance, proving the existence of
steady-state solutions requires the use of weighted Sobolev spaces. To address
this issue, we introduce a simplifying perturbation that removes the
singularity at the border [DSBFK24]. We consider the diffusion function given
by:

$\kappa(x)=D\cdot(1-x^{2})+\delta(x),$

with $D=0.3$ and $\delta\in C^{\infty}(-1,1)$, $\delta(x)=0$ if
$\absolutevalue{x}\leq\eta$, $\delta$ even and non-decreasing in $(0,1)$. The
radiation absorbed by the planet, denoted as $R_{a}$, is the product of a
spatially dependent solar radiation function
$Q_{0}(x)=\hat{Q}_{0}\cdot(1-x^{2})$, where $\hat{Q}_{0}>0$, and a
temperature-dependent co-albedo $\beta=\beta(u)$. The co-albedo
$\beta(u)=1-\alpha(u)$, where $\alpha$ is parameterized by a smooth, non-
increasing, bounded function [BDvdH22, DSBFK24]

$\alpha(u)=\alpha_{1}+\frac{\alpha_{2}-\alpha_{1}}{2}\cdot\left[1+\tanh{\left(K\cdot\left(u-u_{ref}\right)\right)}\right],$
(8)

with $0<\alpha_{1}<\alpha_{2}<1$, $\alpha_{1}$ denotes ice albedo,
$\alpha_{2}$ denotes water albedo, $u_{ref}=275\text{ K}$, and $K$ is a
parametrization constant.

Third, we describe the modelization 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 [DG18].

To model this effect, we choose a space-dependent, latitude-symmetric Outgoing
Longwave Radiation (OLR) of the form:

$R_{e}(x,u;q)=q+\left|x\right|R_{e}^{pl}(u)+\left(1-\left|x\right|\right)R_{e}^{eq}(u),$
(9)

where $R_{e}^{pl}$ and $R_{e}^{eq}$ denote respectively the OLR at the pole
and the equator. The convex combination between $R_{e}^{pl}$ and $R_{e}^{eq}$
Figure 3 shows the space-dependent $OLR$ $R_{e}$, using different colours to
represent distinct latitude points.

Figure 3: Representation of the Outgoing Longwave Radiation (OLR)
$R_{e}=R_{e}(x,u;q)$ in (9) for $q=0$. Different colours are used to represent
the dependence on the space. At the tropics, the OLR presents the bistability
due to the super-greenhouse effect shown in Figure 2. The bistable regimes
progressively disappear as the space point moves to the poles, where the
Stefan-Boltzmann law in (10) is applied to parametrize OLR.

We define them as:

$R_{e}^{pl}(u)=\varepsilon_{0}^{pl}\sigma_{0}u\absolutevalue{u}^{3},$ (10)

where $\varepsilon_{0}$ is the emissivity constant, $\sigma_{0}$ is the
Stefan-Boltzmann constant. On the other hand, at the poles, we reproduce the
super-greenhouse effect by considering:

$R_{e}^{eq}(u)=R_{e}^{pl}(u)g(u)+(1-g(u))\varepsilon_{0}^{eq}\sigma_{0}u\absolutevalue{u}^{3},$
(11)

where $g$ is a smooth transition function of the form

$g(u)=\frac{1}{1+e^{\left(u-u_{ref}^{SGE}\right)K^{SGE}}},$

where $u_{ref}^{SGE}=303.2$ K and $K^{SGE}=0.36$ are respectively a reference
temperature and a constant in the transition velocity involved in the SGE
parametrization. Further, $\varepsilon_{0}^{eq}<\varepsilon_{0}^{pl}$ is an
emissivity constant at the equator taking into account the OLR reduction due
to the SGE. Note that the OLR at the poles and at the equator in (10)-(11) has
been defined, just for mathematical convenience, in the physical meaningless
range of negative Kelvin temperature.

Lastly, the additive positive parameter $q$ models the effect of
$\text{CO}_{2}$ concentration on the energy budget [MMS01, BDvdH22]. 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
$\text{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
$\text{CO}_{2}$ concentration has increased by around $47\%$ from $1850$ to
$2020$, moving from $284$ parts per million (ppm) to $412$ ppm [oCCI23]. 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
$\text{CO}_{2}$, have a large lifetime, many studies have been conducted using
global average values for the spatial distribution [MHSS98, MMS01, HDG+01].

Table 1: Parameters and constants appearing in the 1D EBM (7). Symbol | Meaning | Value
---|---|---
$D$ | Diffusivity constant | $0.45$
$\hat{Q}_{0}$ | Mean solar radiation | $341.3\;\text{W}\,\text{m}^{-2}$
$\varepsilon_{0}^{pl}$ | Emissivity at the poles | $0.61$
$\varepsilon_{0}^{eq}$ | Emissivity at the equator | $0.478$
$\sigma_{0}$ | Boltzmann’s constant | $5.67\cdot 10^{-8}\text{W}\,\text{m}^{-2}\,\text{K}^{-1}$
$\alpha_{1}$ | Ice albedo | $0.7$
$\alpha_{2}$ | Water albedo | $0.289$
$K$ | Constant rate - albedo parametrization | $0.1$
$K^{SGE}$ | Constant rate - SGE parametrization | $0.1$
$u_{ref}$ | Reference temperature - albedo parametrization | $275\;\text{K}$
$u_{ref}^{SGE}$ | Reference temperature - SGE parametrization | $303.2\;\text{K}$
$C_{T}$ | Heat capacity | $5\cdot 10^{7}\;\text{J}\,\text{m}^{-2}\text{K}^{-1}$

###  3.3 Deterministic properties of the model

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 of GMT with respect to the
greenhouse parameter $q$, in the proximity of the model configuration
describing the actual climate of the Earth.

In general, in dynamical system theory, huge information are given by the
study of the steady-state solutions of the model, which are, in a sense, the
long-time behaviour solutions. More in detail, four our model consists of the
non-negative solutions $u=u(x)$ of the following elliptic PDE:

$\begin{split}0&=\left(\kappa(x)u^{\prime}\right)^{\prime}+R_{a}(x,u)-R_{e}(x,u;q),\\\
0&=u^{\prime}(-1)=u^{\prime}(1).\end{split}$ (12)

A common way to prove the existence of at least one stable steady-state
solution involves a variational approach by studying the minimization problem

$\inf\left\\{F_{q}(u)\mid u\in H^{1,2}(-1,1),\;u\geq 0\right\\},$ (13)

where:

$F_{q}(u):=\int_{-1}^{1}\mathcal{R}(x,u(x))\,dx+\frac{1}{2}\int_{-1}^{1}\left[\kappa(x)u^{\prime}(x)\right]^{2}\,dx,$
(14)

with $\partial_{u}\mathcal{R}=R_{e}-R_{a}$ and $H^{1,2}(-1,1)$ denoting the
Sobolev space of order $1$ and exponent $2$ [Nor75, NHPW79, Nor90, Bre11].
Applying the results [DSBFK24, Theorem 1 and Theorem 1], it is possible to
gain information about the existence and uniqueness of a steady state
solution.

###### Proposition 2

  1. (i)

There exists an unique minimizer $u_{*}\in C^{\infty}([-1,1])$ of the
variational problem (13). Further, $u_{*}$ solves the elliptic problem (12)
and it is stable for the dynamics of the 1D EBM (7).

  2. (ii)

The map $q\mapsto\frac{1}{2}\int_{-1}^{1}u_{*}(x)ds$ is non-decreasing.

Note how the second part of the previous result gives qualitative information
on the behaviour 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 [DSBFK24, Theorem 2].

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 (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 Figure 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,\bar{u}_{*})$
plane, where $u_{*}$ denotes a solution of the elliptic PDE, and
$\bar{u}_{*}=\frac{1}{2}\int_{-1}^{1}u_{*}(x)\,dx$ is its GMT, is depicted in
Figure 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 Figure 4(c), a strong non-linear behaviour
of the GMT for the warm solution $u_{W}$ appears for the value of $q$ in the
neighbourhood of $q\approx 11.3$, a thing that does not happen if we remove
the bistability in $R_{e}$ [DSBFK24]. Our aim is to focus on that phenomenon,
when a noise component, describing the weather effect is added to the model.
This is the main topic of the next section.

(a)

(b)

(c)

Figure 4: (a) Steady-state solutions for the 1D EBM (7) for $q=11.3$. Solid
lines denote stable solutions, dotted lines denote unstable solutions. The
snowball solution $u_{S}$ is plotted in blue, the middle solution $u_{M}$ in
red, and the warm solution $u_{W}$ in yellow. (b) Bifurcation diagram in the
$(q,\bar{u}_{*})$ plan for the 1D EBM (7), where $u_{*}$ denotes a steady-
state solutions and $\bar{u}_{*}=\frac{1}{2}\int_{-1}^{1}u_{*}(x)dx$ is its
global mean temperature (GMT). The $S$-shaped bifurcation diagram is
characterized by the two classical saddle-node bifurcations around $q\approx
7$ and $q\approx 38$, and a non-linear (with respect to $q$) increase in GMT
around $q\approx 22$. (c) Zoom of the bifurcation diagram around $q\approx
11.3$ and $u_{W}$.

###  3.4 Stochastic properties of the model

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
literature, this has been achieved in an elementary way by adding a stochastic
term, such as white noise, to the radiation budget [Has76, Imk01, DLV09,
DD22]. 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.

We consider $H=L^{2}(-1,1)$ and the stochastic EBM given by:

$\begin{split}du_{t}&=\left[Au_{t}+R(x,u_{t})\right]dt+\sigma dW_{t},\\\
u_{\mid_{t=0}}&=u_{0}\end{split}$ (15)

where $(W_{t})_{t}$ denotes a cylindrical Wiener process on $H$, $u_{0}\in H$
is a non-negative initial condition and $A:D(A)\to L^{2}(-1,1)$ is the
operator

$\begin{split}D(A)&=\left\\{u\in
H^{2}(-1,1)\;\mid\;u^{\prime}(-1)=u^{\prime}(1)=0\right\\}\\\
Au&=\left(\kappa(x)u^{\prime}\right)^{\prime}.\end{split}$ (16)

We would like to apply the invariant measure theory for gradient systems to
deduce the existence, and uniqueness of an invariant measure and its explicit
formula. This can not be applied to the operator $A$ since it is invertible on
$L^{2}(-1,1).$ For this reason, we consider the operator

$\tilde{A}:=\lambda Id-A,$ (17)

where $\lambda>0$ is a positive constant. By exploiting the Sturm-Liouville
theory, it can be proved that $-\tilde{A}$ is a self-adjoint, negative
definite operator, with eigenvalue $0<\lambda_{1}<\lambda_{2}<\cdots$, such
that the trace of $Tr\left[\left(-\tilde{A}\right)^{\beta-1}\right]<+\infty$,
for some $\beta\in(0,1)$. See Appendix A.2 for a sketch of the proof of these
facts. In this way, it is possible to deduce the following result [DP06,
DPZ14, DSBFK24].

###### Proposition 3

Let $E=C([-1,1]).$ Then, there exists an unique $\mathbb{P}$-a.s. $E$-valued
mild solution of the SPDE (15). Further, there exists an unique Gibbs
invariant measure $\tilde{\nu}$, and $\tilde{\nu}\ll\tilde{\mu}$ with explicit
formula

$\tilde{\nu}(du)=\frac{1}{Z}\exp\left(-\frac{2}{\sigma^{2}}\tilde{I}(u)\right)\tilde{\mu}(du),$
(18)

where
$\tilde{\mu}\sim\mathcal{N}\left(0,-\frac{\sigma^{2}}{2}\tilde{A}^{-1}\right)$
and

$\tilde{I}_{q}(u)=\int_{-1}^{1}\mathcal{R}(x,u(x))dx-\frac{\lambda}{2}\norm{u}_{2}^{2}.$

We remark that the invariant measure $\tilde{\nu}$ is the object, in our
context of stochastic EBM, that in Section 2.5 has been taken as the
definition of climate. Then, it is worth pointing out the link between the
invariant measure $\tilde{\nu}$ and the functional $F_{q}$ building up the
variational problem (13). Indeed, at least formally, we can write the Gaussian
measure $\tilde{\mu}$ as

$\tilde{\mu}(du)=\frac{1}{Z_{1}}\exp\left(-\frac{1}{2}\langle\mathcal{Q}^{-1}u,u\rangle\right)"du",$

where $Z_{1}$ is a normalization constant,
$\mathcal{Q}=-\frac{\sigma^{2}}{2}\tilde{A}^{-1}$ is the covariance operator
of $\tilde{\mu}$, $\langle\cdot,\cdot\rangle$ denotes the scalar product on
$H=L^{2}(-1,1)$ and $"du"$ is the formal notation for the Lebesgue measure on
$H$. By an integration of parts, we deduce

$-\frac{1}{2}\langle\mathcal{Q}^{-1}u,u\rangle=-\frac{2}{\sigma^{2}}\left(\frac{\lambda}{2}\norm{u}_{2}^{2}-\frac{1}{2}\langle\kappa(x)u^{\prime},u^{\prime}\rangle\right).$

Hence, substituting back the formal expression for $\tilde{\mu}$ in (18), we
conclude

$\tilde{\nu}(du)\propto\exp\left(-\frac{2}{\sigma^{2}}F_{q}(u)\right)"du".$

From this expression, we can deduce two important facts. First, the invariant
measure of the stochastic EBM is concentrated on the global minimum points of
the functional $F_{q}$. Second, the study of $\tilde{\nu}$ and its spread
depending on $q$ is as difficult as understanding how the functional $F_{q}$
changes. For this reason, in the next section, we will use numerical
simulation to investigate how $\tilde{\nu}$ changes, depending on the
adiabatic parameter $q$, as the $\text{CO}_{2}$ lies in the critical interval
in Figure 4(c).

###  3.5 Variance and extreme weather events increase

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 [oCCI23]. Despite
numerous definitions proposed to assess what constitutes an extreme event,
there is currently no universally accepted definition [SDM08]. 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 Section 3.4, 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 (18), obtaining theoretical information is
challenging. Thus, we rely on numerical simulations to capture the behaviour
of the variance.

Specifically, given a fixed value of $q>0$ and the warm steady-state solution
$u_{W}=u_{W}^{(q)}$, we numerically integrate the stochastic PDE:

$\begin{split}C_{T}\partial_{t}u&=\partial_{x}\left(\kappa(x)u_{x}\right)+R_{a}(x,u)-R_{e}(x,u;q)+\sigma
dW_{t},\quad(x,t)\in[-1,1]\times[0,T],\\\ u(x,0)&=u_{W}(x),\quad
x\in[-1,1],\\\ u_{x}(-1,t)&=u_{x}(1,t)=0,\quad t\in[0,T].\end{split}$ (19)

The simulation runs for $T=500$ years to capture the properties of the
invariant measure around the warm climate $u_{W}$, and we chose a noise
intensity $\sigma=0.2.$ The finite difference method applied to simulate the
equation is detailed in Appendix A.1. Here, we describe what we mean by
variance, its properties detected by numerical experiments, and our
observations.

Given a space point $x\in[-1,1]$ and a realization $\omega\mapsto u(\omega)$
of the solution of the equation (19), we consider the variance of the process
$t\mapsto u_{t}(x,\omega)=u(x,t)$. We denote the numerical approximation of
the solution by $U=(u_{ij})_{ij}$, where $u_{ij}=u(x_{i},t_{j})$, and
$(x_{i})_{i=1,...,n}$ and $(t_{j})_{j=1,...,m}$ represent the spatial and
temporal meshes, respectively, over the domain $[-1,1]$ and $[0,T]$. The time-
variance is calculated as

$\sigma^{2}_{t}(x_{i})=\frac{1}{m}\sum_{j=1}^{m}\left(u_{ij}-\bar{u}_{i}\right)^{2},$

where $\bar{u}_{i}=\frac{1}{m}\sum_{j=1}^{m}u_{ij}$. Figure 5(a) illustrates
the time-variance indicator, with different colours representing different
values of the parameter $q$.

We observe two distinct behaviours of the variance. First, depending on $q$,
there is an increase in $\sigma^{2}_{t}$, peaking around $q\approx 11.3$,
followed by a decrease. This peak corresponds to the $q$ value that results in
the largest increase in GMT, as shown in Figure 4(c). Second, as expected, for
a given value of $q$, the spatial profile of $\sigma^{2}_{t}$ is symmetric
with respect to $x=0$. Additionally, three local maximum points emerge one at
$x=0$ and the others at $|x|\approx 0.8$. Our interpretation is that
$\sigma^{2}_{t}$ 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 support this interpretation, it is useful to introduce the local stability
indicator

$\gamma(x):=\partial_{u}R(x,u_{W}(x)),$

where $R(x,u):=R_{a}(x,u)-R_{e}(x,u;q)$. If the model were a simple Ordinary
Differential Equation (ODE) without the coupling diffusion term, $\gamma(x)$
would be the eigenvalue obtained by linearizing the reaction term around the
stable steady-state temperature $u_{W}(x)$. Positive values of $\gamma(x)$
indicate local instability, while negative values indicate stability if the
diffusion term is removed. Figure 5(b) presents the local stability indicator
$\gamma$, with colours denoting different values of the parameter $q$.
Notably, there is a clear relationship between time-variance and the local
stability indicator: regions with high time-variance $\sigma^{2}_{t}$ coincide
with areas of positive $\gamma$. This explains the spatial profile of
$\sigma^{2}_{t}$ for a fixed value of $q$.

(a)

(b)

Figure 5: Indicators for the stochastic 1D EBM (19) over the time interval
$[0,T]$, with $T=500\text{ yrs}$. (a) Time-variance indicator
$\sigma^{2}_{t}$. Different colours represent different values of $q$ for
which the stochastic 1D EBM is simulated. (b) Local stability indicator
$\gamma$, using the same colour scheme as in subfigure (a).

The reason for $\gamma(x)>0$ at some points is due to the presence of the
diffusion term. If $\kappa\equiv 0$, then $u_{W}$ is not a minimizer of
$F_{q}$, but instead solves the problem

$u_{W}(x)=\operatorname*{arg\,min}_{u\geq 0}\mathcal{R}(x,u),$ (20)

where $\partial_{u}\mathcal{R}:=R_{e}-R_{a}$. Here,
$\partial_{u}\mathcal{R}_{\mid(x,u_{W}(x))}=0$, and since $u_{W}(x)$ is a
minimum point, it follows
$\partial^{2}_{u}\mathcal{R}_{\mid(x,u_{W}(x))}=\partial_{u}R_{\mid(x,u_{W}(x))}\leq
0$. But with $\kappa(x)>0$, $u_{W}(x)$ no longer satisfies (20), but instead
is a minimizer of $F_{q}$, where the diffusion term appears in the second term
on the right-hand side of (14). This term, which intuitively minimizes the
temperature gradient, allows for some points to fluctuate around a mean value
that would be unstable if $\kappa\equiv 0$.

To understand why, given $x\in[-1,1]$, the map
$q\mapsto\sigma^{2}_{t}(x)=\sigma^{2,(q)}_{t}(x)$ increases until around
$q\approx 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
2. 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.

  1. (i)

Figure 6 shows how the steady-state solution $u_{W}^{(q)}$ approaches unstable
areas, defined by points $(x,u)$ where $\partial_{u}R(x,u)>0$. For
$q_{1}=11.21$, $q_{2}=11.28$, and $q_{3}=11.4$, the map
$(x,u)\mapsto\partial_{u}R(x,u_{W}(x))$ is presented, with the curve
$x\mapsto(x,u_{W}(x),\partial_{u}R(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 $\partial_{u}R$, 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.

  2. (ii)

We consider the local mean stability indicator

$\bar{\gamma}(q):=\int_{-1}^{1}\gamma(x)dx=\int_{-1}^{1}\partial_{u}R(x,u_{W}^{(q)}(x))dx.$

Its plot is shown in Figure 7(a). Although providing only average information,
it exhibits behaviour qualitatively similar to that of $\sigma^{2}_{t}$ for
each fixed value of $x$, peaking around $q\approx 11.3$.

(a) $q=11.21$

(b) $q=11.29$

(c) $q=11.4$

Figure 6: Plot of the map $(x,u)\mapsto\partial_{u}R(x,u)$ and the curve
$x\mapsto(x,u_{W}(x),\partial_{u}R(x,u_{W}^{(q)}(x))$ in red, for different
values of $q$.

(a)

(b)

Figure 7: (a) Average local stability indicator
$\bar{\gamma}(q)=\int_{-1}^{1}\partial_{u}R(x,u_{W}^{(q)}(x))dx$. It
increases, peaks around $q\approx 11.3$, and then decreases. (b) Histogram of
the distribution of the solution of $t\mapsto u(\bar{x},t)$ for the stochastic
EBM (19) for $q=11.2$ (blue) and $q=11.3$ (red), at $\bar{x}=0.6$.

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 Figure 1 and weather
observations in Modena in Figure 2(a).

##  4 Conclusions

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 $\text{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 in 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 $\text{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 $\text{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 Section 3.5. There, we
have exploited numerical simulations to study the changes in the invariant
measure as the $\text{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\in[-1,1]$ we
associate the extreme weather event frequency with the time variance
$\sigma^{2}_{t}(x)$ of a trajectory of the stochastic 1D EBM at point $x$. We
explain the spatial behaviour of $\sigma^{2}_{t}(x)$, which presents two local
maximum points for $x\approx\pm 0.8$ and on 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 $\sigma^{2}_{t}$ is highest. We also motivated the
behaviour of $\sigma^{2}_{t}(x)=\sigma^{2}_{t}{(q)}(x)$ concerning 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 parameterisation 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, 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.

##  Appendix A Appendix

###  A.1 Numerical scheme

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 [LPS14, Section 10.5].

First of all, it is worth recalling that the numerical experiments presented
in Section 3.5 are all performed for a fixed value of $q$, and using
$u_{W}=u_{W}^{(q)}$ as initial condition of the parabolic stochastic problem
(19). The steady-state solution $u_{W}$ solves the elliptic problem (12). To
numerically approximate it, we have applied the same finite difference scheme
described in [DSBFK24, Appendix A].

Second, we move to describe the numerical method for the parabolic problem. We
denote by $R(x,u;q)=R_{a}(x,u)-R_{e}(x,u;q)$ the non-linear radiation budget,
and the underline notation to denote a vector (e.g.
$\underline{y}\in\mathbb{R}^{4}).$ We consider two uniform meshes for the
spatial domain $[-1,1]$ and the time domain $[0,T]$, i.e.

$x_{0}=-1<x_{1}<\cdots<x_{n}=1,\quad x_{i}=-1+i\Delta x,\quad
i=0,...,n,\quad\Delta x=\frac{2}{n},$

and

$t_{0}=0<t_{1}<\cdots<t_{m}=1,\quad t_{j}=j\Delta t,\quad
j=0,...,m,\quad\Delta t=\frac{T}{m}=\frac{500}{m}.$

The number of points in the space and time mesh, respectively $n+1$ and $m+1$,
are chosen in a way that $\Delta x=\Delta x=0.01$. Then, the solution to the
problem can be approximated by considering the system

$\begin{split}C_{T}\frac{u_{i,j+1}-u_{i,j}}{\Delta
t}&=\frac{u_{i-1,j+1}\kappa_{i-\frac{1}{2}}-u_{i,j+1}(\kappa_{i-\frac{1}{2}}+\kappa_{i+\frac{1}{2}})+u_{i+1,j+1}\kappa_{i+\frac{1}{2}}}{\Delta
x^{2}}\\\ &+R(x_{i},u_{i,j+1})+\sqrt{\frac{\Delta t}{\Delta x}}\sigma
z_{i,j}\quad 0\leq i\leq n,\,0\leq j\leq m-1,\\\
\frac{u_{1,j}-u_{-1,j}}{2\Delta x}&=\frac{u_{n+1,j}-u_{n-1,j}}{2\Delta
x}=0\quad j=1,...,m,\\\ u_{i,0}&=u_{W}(x_{i}),\quad i=0,...,n,\end{split}$

where $u_{-1,j},u_{n+1,j}$ are ghost points, $\left(z_{i,j}\right)_{i,j}$ is a
collection of independent identically distributed (i.i.d.) normal random
variables, and
$u_{i,j}=u(x_{i},t_{j}),\;\kappa_{i\pm\frac{1}{2}}=\kappa(x_{i\pm\frac{1}{2}}),\;x_{i\pm\frac{1}{2}}:=x_{i}\pm\Delta
x/2.$ Since $u_{1,j}=u_{-1,j}$ and $u_{n+1,j}=u_{n-1,j}$, the previous system
of equation can be rewritten as

$\left(I_{n+1}-rA\right)\underline{u}^{j+1}=I_{n+1}\underline{u}^{j}+\frac{\Delta
T}{C_{T}}\underline{R}(\underline{u}^{j+1})+\sqrt{\frac{\Delta t}{\Delta
x}}\sigma\underline{Z}^{j},$

where $I_{n+1}$ denotes the $(n+1)\times(n+1)$ identity matrix ,
$\underline{u}^{j}=\left(u_{0,j},\cdots,u_{n,j}\right)^{T}$, $r=\frac{\Delta
t}{C_{T}\Delta x^{2}}$,
$\underline{R}(\underline{y})=\left(R(x_{0},y_{1}),\cdots,R(x_{n},y_{n+1})\right)^{T}$,
$\left(\underline{Z}^{j}\right)_{j}$ denotes a set of i.i.d.
$\mathcal{N}(\underline{0},I_{n+1})$ normal random vectors, and
$A\in\mathbb{R}^{(n+1)\times(n+1)}$ is the tridiagonal matrix

$A=\begin{bmatrix}-\frac{\kappa_{-1/2}+\kappa_{1/2}}{2}&\frac{\kappa_{-1/2}+\kappa_{1/2}}{2}&&&&\\\
\kappa_{-1/2}&-(\kappa_{-1/2}+\kappa_{1/2})&\kappa_{1/2}&&&\\\
&\ddots&\ddots&\ddots&&&\\\
&&&\kappa_{n-1/2}&-(\kappa_{n-1/2}+\kappa_{n+1/2})&\kappa_{n+1/2}\\\
&&&&\frac{\kappa_{n-1/2}+\kappa_{n+1/2}}{2}&-\frac{\kappa_{n-1/2}+\kappa_{n+1/2}}{2}\end{bmatrix}.$

Thus, given the numerical approximation $\underline{u}^{j}$ at time $t_{j}$,
to advance the scheme at time $t_{j+1}$ it is needed to solve the previous
non-linear system of algebraic equations. To do this, we apply the Newton-
Raphson Method (NRM), that we recall in the following. We consider the map

$F^{j}(\underline{y})=\left(I_{n+1}-rA\right)\underline{y}-\underline{u}^{j}-\frac{\Delta
T}{C_{T}}\underline{R}(\underline{y})-\sqrt{\frac{\Delta t}{\Delta
x}}\sigma\underline{Z}^{j}.$

To approximate the vector $\underline{y}\in\mathbb{R}^{n+1}$ such that
$F^{j}(\underline{y})=0$, the NRM consider the following sequence

$\begin{cases}\underline{y}^{k+1}&=\underline{y}^{k}-J_{F^{j}}(\underline{y}^{k})^{-1}F^{j}(\underline{y}^{k}),\\\
\underline{y}^{0}&=\underline{u}^{j},\end{cases}$

where $J_{F^{j}}\in\mathbb{R}^{(n+1)\times(n+1)}$ is the Jacobian matrix of
$F^{j}$ with respect to $\underline{y}.$ Since the noise intensity
$\sigma=0.2$ is small, we expect that the choice
$\underline{y}^{0}=\underline{u}_{j}$ is a good initial guess of the solution.
The iteration of the NRM are stopped when
$\norm{\underline{F}^{j}(\underline{y}^{k}}\leq 10^{-10}$.

###  A.2 Spectral properties of the operator $\tilde{A}$

To apply the invariant measure theory, the operator $\tilde{A}=\lambda Id-A$
defined in Section 3.4 should satisfies the following assumptions [DP06,
Section 11]:

  1. (i)

$\tilde{A}$ is self-adjoint and there exists $\omega>0$ such that

$\langle\tilde{A}x,x\rangle\leq-\omega\absolutevalue{x}^{2},\quad x\in D(A).$

  2. (ii)

There exists $\beta\in(0,1)$ such that
$Tr\left[(-\tilde{A})^{\beta-1}\right]<+\infty$.

We denote by $(\lambda_{n})_{n}$ the eigenvalues of $-A$, where

$\begin{split}D(A)&=\left\\{u\in
H^{2}(-1,1)\;\mid\;u^{\prime}(-1)=u^{\prime}(1)=0\right\\}\\\
Au&=\left(\kappa(x)u^{\prime}\right)^{\prime}.\end{split}$ (21)

Since $A$ is a Sturm-Liouville regular operator, classical results assure the
existence of the real eigenvalues $\lambda_{n}$, and furthermore, it can be
proved

$\lambda_{1}<\lambda_{2}<\cdots\cdots<\lambda_{n}<\cdots\to+\infty.$

First, to check hypothesis (i) it is thus sufficient to prove that

$\lambda_{n}\geq 0,\quad\forall n.$

Indeed, by definition of eigenvalue, there exists an eigenfunction $v_{n}\in
D(A)$ such that

$Av_{n}=-\lambda_{n}v_{n}.$

Multiplying the previous identity by $v_{n}$ and integrating over the domain
$[-1,1]$, we get

$\int_{-1}^{1}\left[\kappa(x)v^{\prime}_{n}(x)\right]^{\prime}v_{n}(x)dx=-\lambda_{n}\int_{-1}^{1}v_{n}(x)^{2}dx.$

Performing an integration by parts on the left-hand side, we deduce

$-\int_{-1}^{1}\kappa(x)\left(v_{n}^{\prime}(x)\right)^{2}dx=-\lambda_{n}\int_{-1}^{1}v_{n}(x)^{2}dx.$

In conclusion,

$\lambda_{n}=\frac{\int_{-1}^{1}\kappa(x)\left(v_{n}^{\prime}(x)\right)^{2}dx}{\int_{-1}^{1}v_{n}(x)^{2}dx},$

and our claim follows from the fact that $\kappa(x)>0$ on $[-1,1]$.

Second, the hypothesis (ii) is a consequence of the following asymptotic
estimate for regular Sturm-Liouville problems [FP94, Section 4]: there exist
$B>0$ and $n_{0}>0$ such that

$\left[\frac{\left(n-\frac{1}{2}\right)\pi}{B}\right]^{2}<\lambda_{n}<\left[\frac{\left(n+\frac{1}{2}\right)\pi}{B}\right]^{2},\quad\forall
n\geq n_{0}.$

## Code and data availability

This work does not include any externally supplied code. All material in the
text and figures was produced by the authors using standard mathematical and
numerical analysis tools. The code for the numerically simulation of the 1D
EBM is available at Zenodo (https://doi.org/10.5281/zenodo.11609953). The
Modena temperature data were provided by the Geophysical Observatory of
Modena, University of Modena and Reggio Emilia, Italy (www.ossgeo.unimore.it)
and are available upon request from the corresponding author.

## Acknowledgments

We acknowledge fruitful discussions with Paolo Bernuzzi and Giulia Carigi. 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.

## References

  * [Arn01] Ludwig Arnold.  Hasselmann’s program revisited: the analysis of stochasticity in deterministic climate models.  In Peter Imkeller and Jin-Song von Storch, editors, Stochastic Climate Models, pages 141–157, Basel, 2001. Birkhäuser Basel.
  * [AWVC12] Peter Ashwin, Sebastian Wieczorek, Renato Vitolo, and Peter Cox.  Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1166–1184, Mar 2012.
  * [Bal17] Paolo Baldi.  Stochastic Calculus.  Springer International Publishing, 2017.
  * [BC16] Tom Beucler and Timothy W. Cronin.  Moisture-radiative cooling instability.  Journal of Advances in Modeling Earth Systems, 8(4):1620–1640, 2016\.
  * [BDvdH22] Robbin Bastiaansen, Henk A Dijkstra, and Anna S von der Heydt.  Fragmented tipping in a spatially heterogeneous world.  Environmental Research Letters, 17(4):045006, mar 2022.
  * [Bre11] Haim Brezis.  Functional analysis, Sobolev spaces and partial differential equations, volume 2.  Springer, 2011.
  * [Bud69] Mikhail I Budyko.  The effect of solar radiation variations on the climate of the earth.  tellus, 21(5):611–619, 1969.
  * [CLM+23] P. Cannarsa, V. Lucarini, P. Martinez, C. Urbani, and J. Vancostenoble.  Analysis of a two-layer energy balance model: Long time behavior and greenhouse effect.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(11):113111, 11 2023.
  * [CMM20] Piermarco Cannarsa, Martina Malfitana, and Patrick Martinez.  Parameter determination for energy balance models with memory.  Mathematical Approach to Climate Change and its Impacts: MAC2I, pages 83–130, 2020.
  * [DD22] Gregorio Díaz and Jesús Ildefonso Díaz.  Stochastic energy balance climate models with legendre weighted diffusion and an additive cylindrical wiener process forcing.  Discrete and Continuous Dynamical Systems - S, 15(10):2837–2870, 2022.
  * [DF11] Dietmar Dommenget and Janine Flöter.  Conceptual understanding of climate change with a globally resolved energy balance model.  Climate Dynamics, 37(11–12):2143–2165, March 2011.
  * [DG18] M. Dewey and C. Goldblatt.  Evidence for radiative-convective bistability in tropical atmospheres.  Geophysical Research Letters, 45(19):10,673–10,681, 2018.
  * [Día97] Jesús Ildefonso Díaz.  On the mathematical treatment of energy balance climate models.  In Jesús Ildefonso Díaz, editor, The Mathematics of Models for Climatology and Environment, pages 217–251, Berlin, Heidelberg, 1997\. Springer Berlin Heidelberg.
  * [DLV09] J.I. Díaz, J.A. Langa, and J. Valero.  On the asymptotic behaviour of solutions of a stochastic energy balance climate model.  Physica D: Nonlinear Phenomena, 238(9):880–887, 2009.
  * [DP06] Giuseppe Da Prato.  An Introduction to Infinite-Dimensional Analysis.  Springer Berlin Heidelberg, 2006.
  * [DPZ96] G. Da Prato and J. Zabczyk.  Ergodicity for Infinite Dimensional Systems.  London Mathematical Society Lecture Note Series. Cambridge University Press, 1996.
  * [DPZ14] Giuseppe Da Prato and Jerzy Zabczyk.  Stochastic Equations in Infinite Dimensions.  Cambridge University Press, April 2014.
  * [DSBFK24] G. Del Sarto, J. Bröcker, F. Flandoli, and T. Kuna.  Variational techniques for a one-dimensional energy balance model.  Nonlinear Processes in Geophysics, 31(1):137–150, 2024.
  * [DSvN+08] Vasilis Dakos, Marten Scheffer, Egbert H. van Nes, Victor Brovkin, Vladimir Petoukhov, and Hermann Held.  Slowing down as an early warning signal for abrupt climate change.  Proceedings of the National Academy of Sciences, 105(38):14308–14312, 2008.
  * [EWV14] Kerry Emanuel, Allison A. Wing, and Emmanuel M. Vincent.  Radiative-convective instability.  Journal of Advances in Modeling Earth Systems, 6(1):75–90, 2014\.
  * [Flo14] Giuseppe Floridia.  Approximate controllability for nonlinear degenerate parabolic problems with bilinear control.  Journal of Differential Equations, 257(9):3382–3422, 2014.
  * [FP94] C.T. Fulton and S.A. Pruess.  Eigenvalue and eigenfunction asymptotics for regular sturm-liouville problems.  Journal of Mathematical Analysis and Applications, 188(1):297–340, 1994.
  * [FPT22] Franco Flandoli, Umberto Pappalettera, and Elisa Tonello.  Nonautonomous attractors and young measures.  Stochastics and Dynamics, 22(02):2240003, 2022.
  * [FW12] Mark I. Freidlin and Alexander D. Wentzell.  Random Perturbations of Dynamical Systems.  Springer Berlin Heidelberg, 2012.
  * [Ghi76] Michael Ghil.  Climate stability for a sellers-type model.  Journal of Atmospheric Sciences, 33(1):3 – 20, 1976.
  * [GL20] Michael Ghil and Valerio Lucarini.  The physics of climate variability and climate change.  Rev. Mod. Phys., 92:035002, Jul 2020.
  * [GLN93] Charles E Graves, Wan-Ho Lee, and Gerald R North.  New parameterizations and sensitivities for simple climate models.  Journal of Geophysical Research: Atmospheres, 98(D3):5025–5036, 1993.
  * [Has76] K. Hasselmann.  Stochastic climate models part i. theory.  Tellus, 28(6):473–485, 1976.
  * [HDG+01] J. T. Houghton, Y. Ding, D. J. Griggs, M. Noguer, P. J. van der Linden, X. Dai, K. Maskell, and C. A. Johnson, editors.  Climate Change 2001: The Scientific Basis.  Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change, 2001.
  * [Imk01] Peter Imkeller.  Energy balance models — viewed from stochastic dynamics.  In Peter Imkeller and Jin-Song von Storch, editors, Stochastic Climate Models, pages 213–240, Basel, 2001. Birkhäuser Basel.
  * [Kif01] Yuri Kifer.  Averaging and climate models.  In Stochastic climate models, pages 171–188. Springer, 2001.
  * [LL+12] T. M. Lenton, V. N. Livina, V. , E. H. van Nes, and M. Scheffer.  Early warning of climate tipping points from critical slowing down: comparing methods to improve robustness.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1185–1204, 2012.
  * [LPS14] Gabriel J. Lord, Catherine E. Powell, and Tony Shardlow.  An Introduction to Computational Stochastic PDEs.  Cambridge Texts in Applied Mathematics. Cambridge University Press, 2014\.
  * [MHSS98] Gunnar Myhre, Eleanor J. Highwood, Keith P. Shine, and Frode Stordal.  New estimates of radiative forcing due to well mixed greenhouse gases.  Geophysical Research Letters, 25(14):2715–2718, July 1998.
  * [MMS01] Gunnar Myhre, Arne Myhre, and Frode Stordal.  Historical evolution of radiative forcing of climate.  Atmospheric Environment, 35(13):2361–2373, 2001.
  * [NHPW79] Gerald R. North, Louis Howard, David Pollard, and Bruce Wielicki.  Variational formulation of budyko-sellers climate models.  Journal of Atmospheric Sciences, 36(2):255 – 259, 1979.
  * [NK17] Gerald R. North and Kwang-Yul Kim.  Energy Balance Climate Models.  Wiley, August 2017.
  * [Nor75] Gerald R. North.  Theory of energy-balance climate models.  Journal of Atmospheric Sciences, 32(11):2033 – 2043, 1975.
  * [Nor90] Gerald R. North.  Multiple solutions in energy balance climate models.  Global and Planetary Change, 2(3):225–235, 1990.
  * [oCCI23] Intergovernmental Panel on Climate Change (IPCC).  Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change.  Cambridge University Press, 2023.
  * [SBB+09] Marten Scheffer, Jordi Bascompte, William Brock, Victor Brovkin, Stephen Carpenter, Vasilis Dakos, Hermann Held, Egbert Nes, Max Rietkerk, and George Sugihara.  Early-warning signals for critical transitions.  Nature, 461:53–9, 10 2009.
  * [SDM08] David B Stephenson, HF Diaz, and RJ Murnane.  Definition, diagnosis, and origin of extreme weather and climate events.  Climate extremes and society, 340:11–23, 2008.
  * [Sel69] William D. Sellers.  A global climatic model based on the energy balance of the earth-atmosphere system.  Journal of Applied Meteorology and Climatology, 8(3):392 – 400, 1969.
  * [Smo12] Joel Smoller.  Shock waves and reaction—diffusion equations, volume 258.  Springer Science & Business Media, 2012.
  * [Tem97] Roger Temam.  Infinite-Dimensional Dynamical Systems in Mechanics and Physics.  Springer New York, 1997.


# Steno AI at SemEval-2023 Task 6: Rhetorical Role Labeling of Legal Documents
using Transformers and Graph Neural Networks

Anshika Gupta, Shaz Furniturewala, Vijay Kumari, Yashvardhan Sharma
BITS Pilani,Pilani,Rajasthan
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

A legal document is usually long and dense requiring human effort to parse it.
It also contains significant amounts of jargon which make deriving insights
from it using existing models a poor approach. This paper presents the
approaches undertaken to perform the task of rhetorical role labelling on
Indian Court Judgements as part of SemEval Task 6: understanding legal texts,
shared subtask A Modi et al. (2023). We experiment with graph based approaches
like Graph Convolutional Networks and Label Propagation Algorithm, and
transformer-based approaches including variants of BERT to improve accuracy
scores on text classification of complex legal documents.

##  1 Introduction

Rhetorical Role Labelling for Legal Documents refers to the task of
classifying sentences from court judgements into various categories depending
on their semantic function in the document. This task is important as it not
only has direct applications in the legal industry but also has the ability to
aid several other tasks on legal documents such as summarization and legal
search. This task is still in it’s early stages, with huge scope for
improvement over the current state-of-the-art.

To facilitate automatic interpretation of legal documents by dividing them
into topic coherent components, a rhetorical role corpus was created for Task
6, sub-task A of The International Workshop on Semantic Evaluation Modi et al.
(2023). Several applications of legal AI, including judgment summarizing,
judgment outcome prediction, precedent search, etc., depend on this
classification.

##  2 Related Works with Comparison

The predominant technique used in Rhetorical Role Labeling over large datasets
is based on the use of transformer-based models like LEGAL-BERT Chalkidis et
al. (2020) and ERNIE 2.0 Sun et al. (2020), augmented by various heuristics or
neural network models. The accuracy of these approaches has remained low over
the years. The results are summarized in Table 1.

The dataset Parikh et al. (2022) used to implement the above approaches is
relatively small, consisting only of a few hundred annotated documents and 7
sentence classes.

Model | F1 score
---|---
`LEGAL-BERT` | 0.557
`LEGAL-BERT + Neural Net` | 0.517
`ERNIE 2.0` | 0.505
Table 1: Summary of related works on the task of rhetorical role labelling on
legal text. Parikh et al. (2022)

##  3 Dataset

The dataset Kalamkar et al. (2022) is made up of publicly available Indian
Supreme Court Judgements. It consists of 244 train documents, 30 validation
documents and 50 test documents making a total of 36023 sentences.
For every document, each sentence has been categorized into one of 13 semantic
categories as follows:

  1. 1.

PREAMBLE: The initial sentences of a judgement mentioning the relevant parties

  2. 2.

FAC: Sentences that describe the events that led to the filing of the case

  3. 3.

RLC: Judgments given by the lower courts based on which the present appeal was
made to the present court

  4. 4.

ISSUE: Key points mentioned by the court upon which the verdict needs to be
delivered

  5. 5.

ARG_PETITIONER: Arguments made by the petitioner

  6. 6.

ARG_RESPONDENT: Arguments made by the respondent

  7. 7.

ANALYSIS: Court discussion of the facts, and evidence of the case

  8. 8.

STA: Relevant statute cited

  9. 9.

PRE_RELIED: Sentences where the precedent discussed is relied upon

  10. 10.

PRE_NOT_RELIED: Sentences where the precedent discussed is not relied upon

  11. 11.

Ratio: Sentences that denote the rationale/reasoning given by the Court for
the final judgement

  12. 12.

RPC: Sentences that denote the final decision given by the Court for the case

  13. 13.

None: A sentence not belonging to any of the 12 categories

##  4 Proposed Techniques and Algorithms

We try several different approaches for the task at hand. All our models use
LEGAL-BERT as their base, and use various methods for further processing and
refining of results.

The LEGAL-BERT family of models is a modified pretrained model based on the
architecture of BERT Devlin et al. (2019). The variant used in this paper is
LEGAL-BERT-BASE, a model with 12 layers, 768 hidden units, and 12 attention
heads. It has a total of 110M parameters and is pretrained for 40 epochs on a
corpus of 12 GB worth of legal texts.

This model was fine-tuned on the task dataset for 2 epochs with a learning
rate of 1e-5 using the Adam optimizer and Cross entropy loss

###  4.1 Direct Classification of CLS tokens

First, we used the default classifier of LEGAL-BERT to find the first set of
predictions, to establish a baseline for our further experiments. Our next
step used the CLS tokens extracted from the final hidden layer of this trained
model.

Similar to the methodology of Gaoa et al.(2020) and Furniturewala et al.(2021)
we utilised the CLS tokens from LEGAL-BERT for further classification models.
This CLS token is a 768-dimensional semantic feature that represents BERT’s
understanding of the text input. It is a fixed embedding present as the first
token in BERT’s output to the classifier and contains all the useful extracted
information present in the input text.

We tried directly applying various multi-layer neural networks to the
extracted CLS tokens. These two models served as a baseline to assess the
efficacy of our methods.

Figure 1: Extracting CLS Tokens (Furniturewala, 2021)

###  4.2 Graph-Based Approaches

We implemented classificaton systems based on graph architectures. We modeled
the data into a graph using cosine similarity on the CLS tokens generated by
LEGAL-BERT. An edge was created between two sentences if and only if their CLS
tokens had cosine similarity greater than 0.5, with the cosine similarity
acting as edge weight. The threshold was included to minimize the presence of
noise-heavy edges in the graph.

$\cos({\bf x},{\bf y})=\frac{\sum_{i=1}^{n}{{\bf x}_{i}{\bf
y}_{i}}}{\sqrt{\sum_{i=1}^{n}{({\bf x}_{i})^{2}}}\sqrt{\sum_{i=1}^{n}{({\bf
y}_{i})^{2}}}}$ (1)

The cosine similarity between two nodes, X and Y, is defined in equation (1),
where x and y are the CLS tokens for nodes X and Y respectively, and n is the
length of the CLS token, i.e. 768 in this case. The function for the final
adjacency matrix is defined equation (2).

$A_{XY}=\begin{cases}\cos({\bf x},{\bf y})&{if\hskip 5.69054pt\cos({\bf
x},{\bf y})>0.5}\\\ 0&{otherwise}\end{cases}$ (2)

On this graph, we performed the label diffusion algorithm Zhou et al. (2003),
to establish a graph-based baseline for our system. Random walk label
diffusion assigns labels to an unlabeled node using the average of it’s
neighbours, weighted by their distance from the node.

$\centering{F^{t+1}=\alpha\cdot P\cdot F^{t}+(1-\alpha)*Y}\@add@centering$ (3)
$\centering{P=D^{-1/2}\cdot A\cdot D^{-1/2}}\@add@centering$ (4)
$\centering{F^{*}=(1-\alpha)*(I-\alpha P)^{-1}\cdot Y}\@add@centering$ (5)

Figure 2: GCN Architecture Kipf and Welling (2016)

To implement it, we combined the train and validation label array, one-hot
encoded it and masked the validation labels. We then used equation (5) to
generate predictions for each sentence. Here P is the normalised adjacency
matrix, Y is the array of one-hot encoded labels, $\alpha$ is a hyper-
parameter, D is the degree matrix, and Z is the array of predicted labels.

The matrix P is obtained via equation (4), normalizing the adjacency matrix A
using the square root inverse of the degree matrix D. For our experimentation,
we used $\alpha=0.5$.

Furthermore, we used a two-layer Graph Convolution Network (GCN) Kipf and
Welling (2016) to perform classifications on the data. Inspired by the
methodology of BERTGCN Lin et al. (2021), we used the LEGAL-BERT embeddings of
each sentence as the node representation for our graph, and then performed
graph convolutions on it.

The GCN architecture uses trainable weights to identify the optimal weightage
that each neighbour of each node should have on its label. The use of two
layers allows us to incorporate the context of one-hop neighbours into the
label of a particular node.

$\displaystyle Z$ $\displaystyle=f(X,A)$ (6)
$\displaystyle=softmax(\hat{A}\cdot ReLU(\hat{A}XW^{(0)})W^{(1)})$ (7)

We used equation (7) to predict the labels of the validation set. Here, Â
represents the symmetrically normalized adjacency matrix, X is the feature
vector which in this case is the LEGAL-BERT embeddings of the nodes, $W^{i}$
is the matrix of trainable weights in layer $i$.

The calculations required for this approach were extremely computationally
expensive, so we were not able to train the model on the entire training set
on a V100 server. We used half of the training documents for graph building
and the prediction of labels. However, the LEGAL-BERT embeddings were
generated by fine-tuning the model on all training documents.

###  4.3 Context-Based LEGAL-BERT

Our final approach was a Context-Based LEGAL-BERT. We cleaned each sentence by
removing all stopwords (such as ’a’, ’an’, ’the’) present using the NLTK
library. Then we created a 5 sentence input corresponding to any given input
by concatenating its two preceeding sentences and its two succeeding sentences
in order. These 5 sentences were separated using LEGAL-BERT’s separater token
</s>. Sentences at the beginning or end of a document were padded using a
string of <pad> tokens.

These 5 sentence inputs were then tokenized using LEGAL-BERT’s tokenizer and
fed into the model using the baseline parameters. We used the default
classifier to perform classification on these context-based inputs.

##  5 Results

We trained the models and tested them on the validation set. The accuracy
scores have been reported in Table 2.

We see that the performance of these models is significantly better than the
previous attempts at this problem. The improvement of the results of
previously studied models can be attributed to the increase in dataset size,
along with other changes in the structure of the task.

However, our Context-based LEGAL-BERT approach outperforms the other
frameworks by a significant margin. This exhibits that the context of each
sentence is critically important in determining its label, and that we are
successful in incorporating the context of each sentence into its
representation.

We saw that graph-based approaches did not significantly improve performance
compared to the current state-of-the-art models. However, it is important to
note that we were unable to run the Graph Convolution Network using the entire
train dataset due to compute constraints.

Despite such constraints, there might be other reasons for the mediocre
performance of graph-based models. One possible reason is that the
representation of the sentences used for building the model was not able to
capture information necessary to make better predictions. This also explains
how the Context-based LEGAL-BERT performed so much better - it improved the
quality of sentence representation, successfully capturing a wider range of
features pertaining to the task at hand.

Model | Accuracy
---|---
`LEGAL-BERT` | 65.56%
`LEGAL-BERT + Classifier` | 67.15%
`Graph Label Diffusion` | 66.34%
`GCN` | 67.42%
`Context-based LEGAL-BERT` | 71.02%
Table 2: Summary of results obtained by the models on validation dataset

##  6 Conclusion and Future Work

In this paper, we tried several different techniques to perform a sentence
classification task on legal documents. Through our experiments, we show that
incorporating context into the CLS tokens of sentences offers a significant
improvement of 5.5 percentage points over LEGAL-BERT.

Moreover, through our experiments on graph-based models, we show that
improving the CLS tokens results in a better classification, compared to the
regular CLS tokens used in a variety of different ways. The Context-based
LEGAL-BERT model was not only more accurate but also less resource intensive.

For future improvements on these models, we could try the Graph Convolutional
Network approach on the complete dataset. We could also try the various
methods of classification, such as a custom neural network or label diffusion,
on the context-based CLS tokens.

Moreover, we could further try to incorporate more sentences for context of
each target sentence. This would require the use of a long-former model, since
the total number of tokens passed into the model will increases.

## References

  * Chalkidis et al. (2020) Ilias Chalkidis, Manos Fergadiotis, Prodromos Malakasiotis, Nikolaos Aletras, and Ion Androutsopoulos. 2020.  LEGAL-BERT: The muppets straight out of law school.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 2898–2904, Online. Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Furniturewala et al. (2021) Shaz Furniturewala, Racchit Jain, Vijay Kumari, and Yashvardhan Sharma. 2021.  Legal text classification and summarization using transformers and joint text features.
  * Gaoa et al. (2020) Jiaming Gaoa, Hui Ninga, Zhongyuan Han, LeiLei Kongb, and Haoliang Qib. 2020.  Legal text classification model based on text statistical features and deep semantic features.
  * Kalamkar et al. (2022) Prathamesh Kalamkar, Aman Tiwari, Astha Agarwal, Saurabh Karn, Smita Gupta, Vivek Raghavan, and Ashutosh Modi. 2022.  Corpus for automatic structuring of legal documents.  In _Proceedings of the Thirteenth Language Resources and Evaluation Conference_ , pages 4420–4429, Marseille, France. European Language Resources Association.
  * Kipf and Welling (2016) Thomas Kipf and Max Welling. 2016.  Semi-supervised classification with graph convolutional networks.
  * Lin et al. (2021) Yuxiao Lin, Yuxian Meng, Xiaofei Sun, Qinghong Han, Kun Kuang, Jiwei Li, and Fei Wu. 2021.  BertGCN: Transductive text classification by combining GNN and BERT.  In _Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021_ , pages 1456–1462, Online. Association for Computational Linguistics.
  * Modi et al. (2023) Ashutosh Modi, Prathamesh Kalamkar, Saurabh Karn, Aman Tiwari, Abhinav Joshi, Sai Kiran Tanikella, Shouvik Guha, Sachin Malhan, and Vivek Raghavan. 2023.  SemEval-2023 Task 6: LegalEval: Understanding Legal Texts.  In _Proceedings of the 17th International Workshop on Semantic Evaluation (SemEval-2023)_ , Toronto, Canada. Association for Computational Linguistics (ACL).
  * Parikh et al. (2022) Vedant Parikh, Upal Bhattacharya, Parth Mehta, Ayan Bandyopadhyay, Paheli Bhattacharya, Kripa Ghosh, Saptarshi Ghosh, Arindam Pal, Arnab Bhattacharya, and Prasenjit Majumder. 2022.  Aila 2021: Shared task on artificial intelligence for legal assistance.  FIRE ’21, page 12–15, New York, NY, USA. Association for Computing Machinery.
  * Sun et al. (2020) Yu Sun, Shuohuan Wang, Yukun Li, Shikun Feng, Hao Tian, Hua Wu, and Haifeng Wang. 2020.  Ernie 2.0: A continual pre-training framework for language understanding.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(05):8968–8975.
  * Zhou et al. (2003) Dengyong Zhou, Olivier Bousquet, Thomas Lal, Jason Weston, and Bernhard Schölkopf. 2003.  Learning with local and global consistency.  In _Advances in Neural Information Processing Systems_ , volume 16. MIT Press.