# Extreme sensitivity and climate tipping points

Peter Ashwin  Anna S. von der Heydt  Peter Ashwin
Centre for Systems, Dynamics and Control,
Department of Mathematics, University of Exeter,
Exeter, EX4 4QF UK<EMAIL_ADDRESS>A. S. von der Heydt
Institute for Marine and Atmospheric Research,
Department of Physics, also Centre for Complex Systems Studies,
Utrecht University, Princetonplein,
Utrecht, The Netherlands<EMAIL_ADDRESS>
###### Abstract

A climate state close to a tipping point will have a degenerate linear
response to perturbations, which can be associated with extreme values of the
equilibrium climate sensitivity (ECS). In this paper we contrast linearized
(‘instantaneous’) with fully nonlinear geometric (‘two-point’) notions of ECS,
in both presence and absence of tipping points. For a stochastic energy
balance model of the global mean surface temperature with two stable regimes,
we confirm that tipping events cause the appearance of extremes in both
notions of ECS. Moreover, multiple regimes with different mean sensitivities
are visible in the two-point ECS. We confirm some of our findings in a
physics-based multi-box model of the climate system.

Keywords Climate sensitivity; Tipping point; Energy balance model; Stochastic
climate model

##  1 Introduction

The equilibrium climate sensitivity (ECS) is widely used as a measure for
expected future global warming. Following Charney’s definition [7], the ECS is
the increase in global mean surface temperature (GMST) per radiative forcing
change after the fast-acting feedback processes in the Earth System reach
equilibrium. Fast-acting means here that those processes are faster than the
time-horizon for global mean temperature evolution that interests us,
typically taken to be 100 years [32].

The value of ECS remains not very well constrained, as the expected warming
per doubling of atmospheric CO2 still contains a considerable uncertainty of
1.5 – 4.5∘C [22]. In fact, this range has not changed much since first ECS
estimates based on energy balance arguments, despite enormous developments in
climate modelling, and improved observational methods [23]. In particular,
large temperature changes and dangerous climate change as a consequence of
increased atmospheric CO2 cannot be excluded. For example, climate
observations from the instrumental period have not narrowed down the range of
expected climate change mainly because of uncertainties in quantification of
the forcing [33]. Recently there have been concerted attempts to constrain ECS
using emergent constraints in climate models [9].

Some sources of uncertainty for ECS lie in the classical measurement or model
uncertainty, although in particular for observations the quantification of the
applied forcing generally contains the largest uncertainties. Furthermore, by
going back in time further than the instrumental period (e.g. last millennium,
glacial cycles or even millions of years of palaeoclimate data) the
uncertainty in both the forcing and the global mean temperature response
becomes significant [23]. There has been a debate as to the cause of the
uncertainty and extremes of the ECS distribution, in particular the long tail
towards high sensitivity values. On the one hand, [31] suggest it is an
inevitable consequence of nonlinear transformation of normally distributed
feedbacks that appear in the denominator when calculating ECS. On the other
hand, [41] suggest they are a sign of ‘tipping points’ owing to nonlinearities
in the system - this has generated a lively debate. In this paper we highlight
(a) the notion of ECS can usefully be generalised to a truly nonlinear
geometric notion: the two-point sensitivity and (b) the distribution of ECS
values is a valuable tool for characterising both state-dependent (feedback)
dynamics and tipping of the climate system.

The climate system exhibits both internal and forced variability on many
timescales. The consequence is that any ‘equilibrium’ is only relative to
fixing part of the feedback processes that are internal to the climate system,
in particular the ‘slower’ part. This requires an assumption that a time scale
separation (into fast and slow processes) exists and the time scale of
interest sits between fast and slow. For climate model simulations and
observations of the last century there might be a time horizon where this is a
reasonable assumption [33], but as we include palaeoclimate data and model
simulations into the estimate of ECS this assumption needs to be carefully
evaluated. In particular, methods to estimate ECS from palaeoclimate data or
models differ from those of (short) climate model simulations; the latter
generally derive ECS from the decay of the energy imbalance at the top of the
atmosphere induced by a instantaneous doubling or quadrupling of CO2 [18];
palaeoclimate reconstructions instead make the assumption that the
reconstructed climate is in (energetic, short time scale) equilibrium and
compare different of these ‘equilibria’ to each other for estimating ECS.
Without compensating for slow feedback processes, palaeoclimate records give
the so-called Earth System Sensitivity (ESS) that includes the effect of slow
processes and boundary conditions (e.g. geography, vegetation and land ice)
[29]. If estimates of these slow processes are available then ECS can be
estimated from the ESS under an assumption of time scale separation [32, 38].

Note that the ECS is usually thought of as a linearized response of the GMST
to perturbations in the radiative balance of the earth. Next to incoming
(short-wave) and outgoing (long-wave) radiation, feedback processes in the
climate system play an important role in determining the ECS. In their sum
these feedback processes tend to enhance ECS (net positive feedback) and
associated time scales vary from fast to very slow. Examples of fast feedback
processes include cloud feedbacks, water vapour feedback and sea-ice
processes. The strength of each of these feedback processes depends on the
background (long-term mean) climate state [40] and it is therefore not
surprising that the sum of the fast feedback processes varies over time in
particular when considering climate states far back in time and under very
different boundary conditions. For example, from palaeoclimate records
together with ice sheet modelling, it has been found that ECS varies
considerably between glacial and interglacial states [24, 15, 25]. Both
present-day and Palaeogene climate model simulations suggest state-dependence
of ECS due to feedback processes [6, 30]. For example, Transient Climate
Response (TCR) considers the deep ocean warming as a slow process.

Abrupt climate shifts have occurred in the past climate system and therefore
seem likely to occur in the future for a variety of reasons [27, 35]. In the
present climate system, potential tipping elements have been identified some
of which may have a considerable impact on future values of GMST [27, 14].
Even those tipping elements that have little affect on the GMST may cause
significant regional damage and/or contribute to global mean climate change by
triggering cascades of transitions involving other tipping elements [11].
Across such an abrupt transition there is a breakdown of the assumption of a
linear response to perturbations, suggesting that the ECS does not adequately
represent the temperature response to radiative perturbations [39, 4]. In
practice, when deriving ECS from palaeoclimate time series, which include
abrupt transitions, these shifts may lead to extreme values of ECS.

In this paper we show that more general notions of ECS can be useful in
understanding the response of a climate state to changes in radiative forcing
- in addition to an ‘instantaneous’ linearised notion of ECS we explore
‘incremental’ and ‘two-point’ climate sensitivities that are distributions
related to dynamic properties of the climate system: they characterise the
geometry of the dynamics and are not simply estimators of a ‘mean ECS’. In
fact the distributions of sensitivities reflect the intrinsic uncertainty due
to climate system dynamics. The paper is organized as follows: Section 2
introduces these notions of ECS and relates them to the underlying climate
dynamics. We illustrate these concepts using a global energy balance model in
Section 3. In particular, we relate properties of extremes of the ECS to the
presence of tipping points and multistability. Section 4 finishes with a
discussion of conclusions and some challenges for the future. Appendix A
extends results of [38] and examines extremes of this sensitivity associated
with tipping points in a more realistic physics-based multi-box model of the
glacial cycles by Gildor and Tziperman [17].

##  2 Sensitivities and the climate attractor

In order to understand variability, abrupt transitions and response to
perturbations we consider the climate system as a high-dimensional multiscale
complex dynamical system whose evolving trajectories form a climate attractor.
The ECS can be defined on this attractor and regimes or states may be
identified where a linear approximation of the response may be reasonable.
Tipping points visible in the GMST will show up as large but occasional shifts
between different ‘climate regimes’ of the attractor, or indeed different
attractors. We visualise the attractor by projection onto climate observables
relevant for determining ECS, i.e. the GMST $T$ and the radiative forcing $R$
per unit area [38]. Consider the energy balance model

$c_{T}\frac{dT}{dt}=R_{\mathrm{forcing}}+R_{\mathrm{slow}}+R_{\mathrm{fast}}-R_{\mathrm{OLW}},$
(1)

where the left hand side represents the rate of change of the global mean
surface temperature $T$ (with specific heat capacity $c_{T}$) and on the right
hand side $R_{\mathrm{forcing}}$ is the (external) radiative forcing
(including changes in CO2), $R_{\mathrm{slow}}$ ($R_{\mathrm{fast}}$) is the
radiative perturbation due to all slow (fast) feedback processes within the
climate system and $R_{\mathrm{OLW}}$ is the outgoing longwave radiation,
respectively. Following the formalism of [32], the specific climate
sensitivity is

$S_{\mathrm{forcing},\mathrm{slow}}=\frac{\Delta T}{\Delta
R_{\mathrm{forcing}}+\Delta
R_{\mathrm{slow}}}\approx\frac{dT}{d(R_{\mathrm{forcing}}+R_{\mathrm{slow}})},$
(2)

which equals the Charney sensitivity $S$ if $\Delta R_{\mathrm{slow}}$ is the
sum of all slow feedback processes contributing to the ECS (and under the
assumption of time scale separation). In practise, only some of the slow
processes are accessible from palaeoclimate records (e.g. only land ice), in
which case the specific climate sensitivity is only an approximation of the
Charney sensitivity [32] (e.g. $S_{[CO_{2},LI]}$ is the specific climate
sensitivity considering only land ice changes as slow feedback). This ECS
gives a linear prediction for change in temperature:

$T^{\prime}=T+S_{\mathrm{forcing},\mathrm{slow}}\,\left(\Delta
R_{\mathrm{forcing}}+\Delta R_{\mathrm{slow}}\right).$ (3)

For a specific energy balance model including regime shifts we can explicitly
calculate ECS for the different regimes, see section 3. We note that several
other authors have highlighted the need to improved notions of ECS: this
includes [8] who propose to use a measure-based approach to understand climate
sensitivity and [12] who consider conditional climate sensitivities
constrained by temperature, coupled with resilience measures for switching to
other regimes.

###  2.1 Observation of the climate attractor

We consider the climate system as a high dimensional dynamical system that
evolves along trajectories $x(t)$ according to a smooth flow

$x(t)=\varphi_{t}(x_{o})$ (4)

where $x\in X$ represents the instantaneous state of the climate system in
some high dimensional state space and $\varphi_{t}(x_{0})$ evolves the initial
state $x_{0}$ along by a time $t$. The global mean temperature
$T:X\rightarrow\mathbb{R}$ and radiative forcing $R:X\rightarrow\mathbb{R}$
are considered to be observables of the underlying dynamical system on $X$. We
assume that the dynamics of $x$ are stationary, i.e. that there is a natural
probability measure $M$ on $X$ such that typical trajectories $x(t)$ of (4)
satisfy

$\lim_{t\rightarrow\infty}\frac{1}{t}\int_{s=0}^{t}p(x(s))\,ds=\int_{x}p(x)\,dM(x).$
(5)

for typical $x_{0}$ and any integrable observable $p:X\rightarrow\mathbb{R}$,
i.e. the long-time average of $p$ can be computed using an ergodic hypothesis,
by averaging over the measure $M$ in phase space. This implies that, for any
open set $\mathcal{A}\subset X$, the long-term average proportion of time a
typical trajectory spends in $\mathcal{A}$ is $M(\mathcal{A})$. For small
enough perturbations, linear response theory suggests a linear change in mean
observables: see for example [28].

In [38] it is supposed there is a stationary measure $\mu$ of points in the
$(\Delta R_{[CO_{2},LI]},T)$-plane according to how often they are visited
over asymptotically long times, i.e. for any measurable subset
$\mathcal{A}\subset\mathbb{R}^{2}$ we define

$\mu(\mathcal{A}):=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{s=0}^{t}\chi_{\mathcal{A}}(\Delta
R_{[CO_{2},LI]}(s),T(s))\,ds.$ (6)

for typical initial condition, where $\chi_{\mathcal{A}}$ is the indicator
function, $\chi_{\mathcal{A}}(\Delta R,T)=1$ if $(\Delta R,T)\in\mathcal{A}$
and $=0$ otherwise. Note that applying (5) with $p(x)=\chi_{\mathcal{A}}(x)$
(where $\chi_{\mathcal{A}}(x)=1$ if $x\in{\mathcal{A}}$ and $0$ otherwise)
gives

$\mu({\mathcal{A}})=M(\\{x~{}:~{}(R_{[CO_{2},LI]}(x),T(x))\in{\mathcal{A}}\\}).$
(7)

In other words, the measure $\mu$ is simply a projection of a natural measure
$M$ on the ‘climate attractor’ onto the two observables $(R_{[CO_{2},LI]},T)$.
In general we will consider throughout

$\Delta R_{[CO_{2},LI]}=R_{[CO_{2},LI]}-\tilde{R}_{[CO_{2},LI]}$

i.e. the change in radiative forcing relative to some fixed reference level
$\tilde{R}_{[CO_{2},LI]}$ usually chosen as the level during the pre-
industrial climate.

###  2.2 Incremental and two-point sensitivities

In order to predict the temperature at some fixed future time-horizon $\Delta
t$ in response to a change in radiative forcing, we consider the quotient (2)
for fixed changes in time. In this case we can view the distribution of what
we call incremental sensitivities as the spread of trajectories from the
current estimated values of instantaneous $(\Delta R,T)$, where from now on we
write the radiative forcing corrected for slow feedback as $R$. On the other
hand we can consider a time-independent choice of pairs of points on the
climate attractor to obtain two-point sensitivities.

Let us assume the current state at time $t=0$ of the climate system is given
by a measure $\sigma_{0}$ on some high dimensional phase space $X$ (that
projects onto a point A in Fig. 1). Note that this will always be a measure
rather than a point because of lack of knowledge of sub-grid parametrized
processes (e.g. [36]) but it will project onto the current values $(\Delta
R_{0},T_{0})=(\Delta R(x),T(x))$ for all $x$ in the support of $\sigma_{0}$.
As time progresses, this state will spread to give a measure at time $t$ that
is

$\sigma_{t}({\mathcal{A}})=\sigma(\varphi_{-t}({\mathcal{A}}))$

for any ${\mathcal{A}}\subset X$ (Fig. 1 shows a trajectory in black and
others from the ensemble starting at $A$ in grey). The incremental sensitivity
for a time interval $\Delta t$ is then

$S^{\Delta t}_{0}(x)=\frac{T(\varphi_{\Delta t}(x))-T_{0}}{\Delta
R(\varphi_{\Delta t}(x))-\Delta R_{0}}$ (8)

with distribution

$\mathbb{P}(S_{0}^{\Delta t}\in{\mathcal{A}})=\sigma(\\{x~{}:~{}S^{\Delta
t}_{0}(x)\in{\mathcal{A}}\\}).$

Over long time, if there is decay of correlations and mixing of trajectories
on the climate attractor [37, 36, 28] then $\sigma_{\Delta t}\rightarrow M$ in
the weak sense as $\Delta t\rightarrow\infty$, and so we expect the
distribution of long-term incremental sensitivities for $\Delta
t\rightarrow\infty$ become time-independent for typical trajectories within
the attractor:

$\mathbb{P}(S_{0}^{\infty}\in{\mathcal{A}})=M(\\{x~{}:~{}S^{\infty}_{0}(x)\in{\mathcal{A}}\\}).$
(9)

where

$S_{0}^{\infty}(x)=\frac{T(x)-T_{0}}{\Delta R(x)-\Delta R_{0}}.$

Note that (7) means that the distribution of long-term sensitivities starting
at $(\Delta R_{0},T_{0})$ can be written in terms of the geometry of the
projected measure $\mu$

$\mathbb{P}(S_{0}^{\infty}\in{\mathcal{A}})=\mu(\\{(\Delta
R_{1},T_{1})~{}:~{}S^{\infty}_{0,1}\in{\mathcal{A}}\\})$ (10)

where we define the two-point sensitivity as

$S^{\infty}_{0,1}=\frac{T_{1}-T_{0}}{\Delta R_{1}-\Delta R_{0}}.$ (11)

The distribution of long-term incremental sensitivities (10) for a generic
choice of the initial climate state suggests [38] a time-independent notion of
climate sensitivity that can be found by picking pairs of points
$(R_{0,1},T_{0,1})$ independently distributed according to $\mu$ and
evaluating (2).

This means that for any ${\mathcal{A}}\subset\mathbb{R}$ we can use $\mu$ to
assign a probability to the sensitivity being in ${\mathcal{A}}$:

$\mathbb{P}(S^{\infty}_{0,1}\in{\mathcal{A}}):=\mu\times\mu\left(\left\\{(\Delta
R_{0},T_{0}),(\Delta
R_{1},T_{1})~{}:~{}S^{\infty}_{0,1}\in{\mathcal{A}}\right\\}\right).$ (12)

with $S^{\infty}_{0,1}$ defined as in (11). This gives, in some sense, a
maximal set of possibilities for the sensitivities in that it compares the
observables $T$ and $\Delta R$ over all possible time points and possible
trajectories of the system. This is comparable to the conditional climate
sensitivity of [12] except rather than dividing into regimes, they restrict to
deviations of temperature at most $\delta_{T}$ from $T_{0}$. In the case that
the sensitivity is fixed at $S_{0}$, note that $S^{\infty}_{0,1}$ is a Dirac
$\delta$-distribution centred at $S_{0}$.

Figure 1: Schematic diagram to demonstrate the instantaneous, incremental and
two-point sensitivities for an ensemble of trajectories starting at point $A$
that evolves on the climate attractor. The equilibrium of a mean energy
balance model is shown as a red curve. For a specific trajectory (shown in
black) in the ensemble the slopes of $AC$ and $AD$ correspond to incremental
sensitivities (for fixed $\Delta t$) or two-point sensitivities (for varying
$\Delta t$). The instantaneous sensitivity is the slope of the tangent to the
closest solution in the equilibrium model at $B$. Note that the cold regime
$T\leq T_{\mathrm{thr}}$ and the warm regime $T>T_{\mathrm{thr}}$ have
different asymptotic (instantaneous) sensitivities corresponding to slopes
$EB$ and $DF$ respectively.

###  2.3 Sensitivities and climate regimes

By partitioning the climate attractor into a number of regimes, we can
condition the sensitivities on staying within a regime, or undergoing a
transition between regimes. By making an optimal partition of the attractor
projected into $(\Delta R,T)$ space we can hope to find localised
distributions of sensitivities for pairs in the same regime. As in [38] we
consider these sensitivities conditional on climate regime by partitioning
$\mu$ into two distributions

$\mu=\mu_{C}+\mu_{W}$

corresponding to being in a cold (C) or warm (W) state. In our case we set

$\displaystyle\mu_{C}({\mathcal{A}})$
$\displaystyle=\mu({\mathcal{A}}\cap\\{(\Delta R,T)~{}:~{}T\leq T_{thr}\\}),$
$\displaystyle\mu_{W}({\mathcal{A}})$
$\displaystyle=\mu({\mathcal{A}}\cap\\{(\Delta R,T)~{}:~{}T>T_{thr}\\})$

for any measurable ${\mathcal{A}}\subset\mathbb{R}^{2}$ and some threshold
temperature $T_{thr}$. As in [38] we define distributions of conditional
sensitivities by

$\mathbb{P}(S^{WW}\in{\mathcal{A}}):=\mu_{W}\times\mu_{W}\left(\left\\{(\Delta
R_{0},T_{0}),(\Delta
R_{1},T_{1})~{}:~{}S^{\infty}_{0,1}\in{\mathcal{A}}\right\\}\right)$ (13)

Conditional sensitivities for changes of regime are for example

$\mathbb{P}(S^{CW}\in{\mathcal{A}}):=\mu_{C}\times\mu_{W}\left(\left\\{(\Delta
R_{0},T_{0}),(\Delta
R_{1},T_{1})~{}:~{}S^{\infty}_{0,1}\in{\mathcal{A}}\right\\}\right)$ (14)

The distribution of sensitivities (12) is then a sum of the four conditional
sensitivities

$\mathbb{P}(S^{\infty}_{0,1}\in{\mathcal{A}})=\mathbb{P}(S^{CC}\in{\mathcal{A}})+\mathbb{P}(S^{CW}\in{\mathcal{A}})+\mathbb{P}(S^{WC}\in{\mathcal{A}})+\mathbb{P}(S^{WW}\in{\mathcal{A}}).$

Moreover, (14) means that we have a symmetry

$\mathbb{P}(S^{CW}\in{\mathcal{A}})=\mathbb{P}(S^{WC}\in{\mathcal{A}}).$

The distributions of $S^{CW}$ and $S^{WC}$ correspond to choices of pairs
across the two regimes: these distributions are associated with ‘tipping
between regimes’. Even though the two-point sensitivities may measure states
very far apart in time, we find that extreme values of the sensitivity are
usually associated with choice of points from two different regimes.

##  3 Sensitivity and tipping in climate models

To illustrate the notions of instantaneous and two-point sensitivities, we
consider a conceptual energy balance model: a more complex model is briefly
discussed in Appendix A. We consider a variant of the Budyko-Ghil-Sellers
energy balance model [5, 16, 34] for GMST. This model builds on [12, 41] and
has multiple regimes with state-dependent sensitivity in each. It is a special
case of (1) for global mean surface temperature $T(t)$ with atmospheric CO2
concentration $C(t)$ as a parameter:

$c_{T}\,\frac{dT}{dt}=F(T,C):=\left[Q_{0}(1-\alpha(T))+A\ln\left(\frac{C}{C_{0}}\right)-\epsilon(T)\sigma
T^{4}\right]$ (15)

For this equation, $Q_{0}$ represents the solar input modulated by the
temperature-dependent albedo $\alpha(T)$. The change in radiative forcing due
to atmospheric greenhouse gases is

$\Delta R_{[CO_{2}]}=A\ln(C/C_{0}).$

where $A=5.35$ Wm-2 is the direct forcing effect of CO2 and $C_{0}$ represents
pre-industrial CO2 levels. Finally, the outgoing long wave radiation $\sigma
T^{4}$ is modified by a temperature-dependent emissivity $0<\epsilon(T)<1$.

We consider a temperature-dependent emissivity decreasing from one plateau to
a lower one because of changes in water vapour and cloud feedbacks. There are
other choices [41], but for simplicity we assume here

$\epsilon(T)=\epsilon_{1}+\frac{\epsilon_{2}-\epsilon_{1}}{2}\left[1+\tanh\left(\frac{T-T_{0}}{T_{\epsilon}}\right)\right].$

where $\epsilon_{1}$ and $\epsilon_{2}$ are the limit emissivities for low and
high temperatures respectively, $T_{0}$ is the threshold and $T_{\epsilon}>0$
corresponds to the range of temperatures over which there is variation (see
Figure 2b).

Note that water vapour feedback is sometimes included in the CO2 term,
resulting in an additional constant and modified $A$ [20, 12]. Here we
separate radiative forcing due to CO2 and temperature-dependent water
feedbacks in emissivity. As in [12], we assume that the albedo varies with
temperature due to changes in land-ice feedback processes: we assume there are
threshold temperatures $T_{1}<T_{2}$ associated with changes of albedo
$\alpha(T)$ and define a function

$\Sigma(T)=\frac{(T-T_{1})}{T_{2}-T_{1}}H(T-T_{1})H(T_{2}-T)+H(T-T_{2})$ (16)

that switches from $0$ for $T<T_{1}$ to $1$ for $T>T_{2}$: $H(T)$ is
approximately a Heaviside unit step function and we use a smooth approximation

$H(T)=(1+\tanh(T/T_{\alpha}))/2$

as in [12]. As in that paper, we write the albedo

$\alpha(T)=\alpha_{1}(1-\Sigma(T))+\alpha_{2}\Sigma(T)$

so that it changes smoothly from a higher albedo $\alpha_{1}$ in the presence
of more ice surface ($T<T_{1}$) to a lower $\alpha_{2}$ in the presence of
more ocean surface ($T>T_{2}$), see Figure 2a. Note that [12] consider a
global transition from ice-covered to ocean-covered earth - here we model a
large but regional change in ice cover with a smaller contrast in global
albedo between the two states; our choice of parameters might be more
realistic for albedo variations between glacial and interglacial states.

Figure 2: Behaviour of the equilibrium energy balance model (15) with
parameters as in Table 1. (a) temperature-dependence of albedo $\alpha(T)$ and
(b) emissivity $\epsilon(T)$; (c) CO2 and (d) radiative forcing levels
necessary to give temperature equilibria, corresponding to (18) and (19),
respectively. Note the region of multistability, and temperature-dependence of
the sensitivity corresponding to slopes in the bottom right figure.

We add a stochastic term to (15) that represents unresolved subgrid processes
with a fixed amplitude $\eta_{T}$:

$c_{T}\,dT=F(T,C)dt+\eta_{T}dW_{T}.$ (17)

The motivation for noise on the $T$ variable is to model internal variability
of climate heat transport processes. The parameters listed in Table 1 are
used, except where specified. Note that the deterministic equilibria of (15)
are at $F(T,C)=0$, which gives

$C=\Gamma(T):=C_{0}\exp\left[\frac{\epsilon(T)\sigma
T^{4}-Q_{0}(1-\alpha(T))}{A}\right].$ (18)

From (18), this means we have equilibria at

$\Delta R_{[CO_{2}]}=A\ln(\Gamma(T)/C_{0})=\epsilon(T)\sigma
T^{4}-Q_{0}(1-\alpha(T))$ (19)

Figure 2 illustrates temperature dependence of albedo and emissivity as well
as the resulting equilibrium forcing $\Delta R=A\ln(\Gamma(T)/C_{0})$ needed
to give this temperature. Note there is a unique equilibrium for each $T$, but
not necessarily for each $C$: as discussed in [12, 41] there are three
branches of equilibria for a range of $C$: for the parameters used there is
bistability in the region

$-1.744~{}Wm^{-2}<\Delta R<3.004~{}Wm^{-2},\
202~{}\mathit{ppm}{}<C<490~{}\mathit{ppm}{},$ (20)

denoted using the red lines in Figs. 3 and 5. We can define the instantaneous
sensitivity111This is referred to as local slope sensitivity in [24]. as
$S=1/\lambda$, where

$\lambda=\frac{d}{dT}\Delta
R_{[CO_{2}]}=[\epsilon^{\prime}(T)T+4\epsilon]\sigma
T^{3}+Q_{0}\alpha^{\prime}(T)$ (21)

is the total feedback factor in this model: $S$ corresponds to the slope of
the tangent of the equilibrium (non-stochastic) model (see Fig. 1, point B).

The sensitivities on the stable branches differ due to both nonlinearity of
black body radiation and change in emissivity. Due to the choice of albedo and
emissivity changes in this model, $\alpha^{\prime}$ is nonzero only in the
bistable regime and $\epsilon^{\prime}$ is nonzero only in the temperature
range where it varies ($T_{0}-T_{\epsilon}/2\leq T\leq T_{0}+T_{\epsilon}/2$):
see Figure 2. State-dependence between glacial and interglacial states has
been detected in estimates of specific climate sensitivities from different
palaeo-data, suggesting lower sensitivity during cold periods than during warm
periods (e.g. [40, 15] who estimate a close approximation of the Charney
sensitivity and find warm (interglacial) climate states to be about 60% more
sensitive than cold (glacial) states).

We can compute the curvature of (19) as

$\frac{d^{2}}{dT^{2}}\Delta
R_{[CO_{2}]}=(\epsilon^{\prime\prime}T^{2}+8\epsilon
T+12\epsilon^{\prime})T^{2}-Q_{0}\alpha^{\prime\prime}$

Note that this is small except near the folds at $T\approx T_{1}$ and
$T\approx T_{2}$ with maximum absolute value $\frac{d^{2}}{dT^{2}}\Delta R$ of
order $T_{\alpha}^{-1}$. This confirms that the saddle-node becomes non-smooth
in the limit $T_{\alpha}\rightarrow 0$. The second derivative gives the size
of the quadratic correction $a$ in Zaliapin et al. [41]; very large values of
the slope near the saddle nodes correspond to the run-away climate observed in
[4].

For the model (17), the atmospheric CO2 concentration is a parameter for the
energy balance dynamics. We explore this by considering a ‘wandering’ CO2
profile such that $\gamma(t):=\ln(C(t))$ undertakes a Brownian motion with
growth in variance $\eta_{\gamma}$ per unit time between reflecting limit
values. More precisely, we consider CO2 dynamics governed by soft reflecting
boundary conditions at $\ln C_{min}$ and $\ln C_{max}$:

$d\gamma=K\theta(\gamma)\,dt+\eta_{C}\,dW_{\gamma}$ (22)

where the noise in the $C$ variable (with fixed amplitude $\eta_{C}$)
represents variability in CO2 forcing. We assume

$\theta(\gamma):=H(\ln C_{min}-\gamma)(\ln C_{min}-\gamma)+H(\gamma-\ln
C_{max})(\ln C_{max}-\gamma)$ (23)

and we use parameters

$K=10^{-7}~{}s^{-1},~{}C_{min}=10^{2}~{}\mathit{ppm},~{}~{}C_{max}=~{}10^{3}~{}\mathit{ppm},~{}\eta_{C}=2\times
10^{-6}~{}s^{-1/2}.$ (24)

Clearly there are common causes of variability of temperature and CO2 and so
in general there will be strong correlations between the noise terms $W_{T}$
and $W_{\gamma}$; for convenience we assume here that they are uncorrelated.
In most studies of climate sensitivity, carbon cycle processes are not treated
as feedbacks but rather as forcing (external to the system); this assumption
corresponds in our model to CO2 driving temperature changes with no direct
impact of temperature on CO2. The parameter $\eta_{C}$ determines the
timescale of wandering of the CO2: we consider cases where this is slower
than, or comparable to the the timescale of evolution of $T$.

Figure 3(b) shows a time series for a typical simulation of (17) with
wandering CO2 (22) and parameters as in Table 1, while the corresponding time
series of $T$ is shown in Figure 3(a); we see as $C$ crosses thresholds (for
the bistable regime as calculated in (20)) the state of the system tips
between warm and cold states. Global mean temperature $T$ vs $C$ and the
relative radiative forcing $\Delta R_{[CO_{2}]}$ are shown in Figure 2 (c,d)
for the non-stochastic model and in Figure 4 for the stochastic model,
respectively. Observe the region of bistability around $\Delta R_{[CO_{2}]}=0$
(Figure 4b) that switches rapidly between cool high albedo and warm low albedo
states via saddle-node bifurcations. There is approximate linearity away from
these tipping points, but with different mean slopes.

$\begin{array}[]{|c|rl||c|rl|}\hline\cr\alpha_{1}&0.52&-&\alpha_{2}&0.47&-\\\
T_{1}&278&K&T_{2}&288&K\\\ \epsilon_{1}&0.53&-&\epsilon_{2}&0.39&-\\\
A&5.35&Wm^{-2}&T_{0}&288&K\\\ T_{\alpha}&5&K&T_{\epsilon}&20&K\\\
C_{0}&280&\mathit{ppm}{}&Q_{0}&342&Wm^{-2}\\\ C_{T}&5\times
10^{8}&Jm^{-2}K^{-1}&\sigma&5.67\times 10^{-8}&Wm^{-2}K^{-4}\\\
\eta_{T}&5\times 10^{-6}&Ks^{-1/2}&\eta_{C}&2\times 10^{-6}&s^{-1/2}\\\
\hline\cr\end{array}$ Table 1: Parameters for the energy balance model (15,17)
adapted from [12] to include state-dependent emissivity. Note that [12]
consider a global transition and so use different values: $\alpha_{1}=0.7$
$\alpha_{2}=0.2$, $T_{1}=263~{}K$, $T_{2}=293~{}K$, $T_{\alpha}=0.273K$,
$A=20.5~{}Wm^{-2}$ and have an additional constant $150~{}Wm^{-2}$.

Figure 3:  Time series for the energy balance model (1) with wandering CO2
(22). (a) Global mean temperature $T(t)$ and (b,c) atmospheric CO2 $C$ in
original (b) and logarithmic coordinates (c). The red lines indicate the
limits imposed on the randomly wandering CO2. See text for details. Figure 4:
(a) Temperature versus (a) atmospheric CO2 concentration and (b) $\Delta R$ of
radiative forcing by CO2, for the time series in Figure 3. The red line
corresponds to the steady solution of (15) with dependence on CO2.

###  3.1 Extreme sensitivities, early warning signals and tipping between
regimes

The energy balance model (17) with wandering CO2 (22) gives a framework in
which one can test correlation between extreme values of sensitivity and
tipping between climate regimes, as well as testing other possible early
warning signals for a tipping event. Indeed, we find a rise in instantaneous
sensitivity seems to act as a good precursor in cases of slow variation of
CO2.

Figure 5(left) shows the variation of instantaneous sensitivity and two early
warning signals for tipping between regimes for the wandering variation of CO2
concentration $\ln C(t)$ relative to the timescale of evolution of $T$, with
$\eta_{C}=2\times 10^{-6}~{}s^{-1/2}$; this means that CO2 variations are
comparable in timescale to $T$ fluctuations. By considering the nearest
equilibrium point on $C=\Gamma(T)$ for fixed $C$ (18), we evaluate the
instantaneous sensitivity $S$ using (21) and plot this in middle panel.
Observe there is a qualitative change in $S$ before and after tipping events.
There is a clear precursor and then singularity as the tipping point is
crossed. Note that in this case the instantaneous sensitivity depends only on
the current CO2 level and the nearest branch - fluctuations in $T$ around the
branch do not affect $S$, while fluctuations in $C$ do. There are also
apparent ‘false alarms’: for example, the fluctuations of $S$ around 3 kyr and
18 kyr. Figure 5(left, panels d,e) show detrended estimates of $sd(t)$ and
$AR1(t)$ [42] using moving averages with length $\tau=500$yrs. Neither show
any clear precursors before tipping events. By contrast, Figure 5(right) shows
early warning signals for tipping between regimes for slower variation of CO2
concentration with $\eta_{C}=5\times 10^{-7}~{}s^{-1/2}$, where $T$ evolves
faster than $C$. Unlike the case on the left, this slower switching gives a
visible precursor of increasing $AR1$. In both cases there are increasing
fluctuations of instantaneous sensitivity $S$. Note however, that the
instantaneous $S$ we consider here uses the model equations, and hence will be
more difficult to access from complex model realisations or observations.

Figure 5: Examples of tipping events for $T$ for the energy balance model
(15,17) with wandering CO2 (22) with (left) parameters as in Figure 3, i.e.
with $(\eta_{T},\eta_{C})=(5\times 10^{-6},2\times 10^{-6})$ and (right)
slower variation of $C$ $(\eta_{T},\eta_{C})=(5\times 10^{-6},5\times
10^{-7})$. (a,b) time series of $\Delta R$ and $T$. The red lines on (a) show
the locations of the saddle-node bifurcations that bound the region of
bistability (20). (c) shows estimated instantaneous sensitivity from the
nearest equilibrium of (17). Note the gradual rise and fluctuations in $S$ on
approach to the tipping point, and the two levels of $S$ corresponding to the
differing sensitivities of the two stable branches. (d,e) show standard
deviation and AR1 coefficient: note that the AR1 coefficient seems to have
predictive power only for the right column.

###  3.2 Two-point sensitivities and tipping

An approximation of the stationary density of the global attractor for the
system (17,22) is shown in Figure 6. We classify the system regime as one of:

$\left\\{\begin{array}[]{lr}\mbox{Warm (W)}&\mbox{ if }T>T_{\mathrm{thr}}\\\
\mbox{Cold (C)}&\mbox{ if }T\leq T_{\mathrm{thr}}\end{array}\right.$

where we choose a threshold $T_{\mathrm{thr}}=10~{}C$ between the two stable
branches: see Figure 1. We simulate a single very long trajectory ($5\times
10^{5}$ years) of the energy balance model (15,17) with wandering CO2 and use
this to create a density plot of the climate attractor projected onto $T$ vs
$\Delta R$, as shown in Figure 6. This is used to consider the two-point
sensitivities and probabilities of tipping as in Figure 7. Panels (a-d) are
computed by sampling incremental sensitivities (8) from points for increments
up to $20~{}kyr$. Panels (e-h) are computed by sampling $10^{7}$ pairs of
points from the distribution in Figure 6 using the two-point sensitivity (11).
We observe:

  * •

There is good qualitative agreement between the incremental sensitivities
averaged over long delays and the two-point sensitivities sampled
independently from the attractor. Indeed, the autocorrelation of the
timeseries for $T$ (not shown) has substantially decayed and has its first
zeros around $20kyr$.

  * •

High probability of tipping (see Figure 7(b,d)) corresponds mostly to extremes
of $S$ that may be positive or negative $S$.

  * •

Within the W and C regimes, the sensitivities are closely clustered but have
different means for the W and C state. We can estimate these using average
temperatures and (21) as $S\approx 0.79~{}K[Wm^{2}]^{-1})$ for the W and
$S\approx 0.55$ for the C state, respectively.

  * •

Note that there are relatively low-probability ‘shoulders’ of the
distributions within-regime. These are due to the classification of regimes
also including states that are in transitions: although the system is in
transition, both points are still classified as the same regime.

  * •

A hysteresis-like effect may appear even in the absence of bistability. If the
dynamic change in a system parameter is fast enough then an apparent
hysteresis may arise due to different lags for a rising or falling parameter:
see for example [19]. However, Figure 4 suggests this is not the case here.

Figure 6: Density of $T$ versus $\Delta R=A\ln(C/C_{0})$ of radiative forcing
by CO2 for a longer time series similar to that in Figure 4(b). Observe the
two regions with approximately linear relationship (cf Figure 1): the colour
scale indicates density [arbitrary units]. Figure 7: (a,c,e,g) Conditional
two-point sensitivities and (b,d,f,h) probabilities of tipping from warm
(a,b,e,f) and cold (c,d,g,h) states, for the energy balance model (17) with
wandering CO2 (22). (a-d) are computed using a range of delays up to
$20~{}kyr$ while (e-h) are computed using $10^{7}$ independently sampled pairs
of points from the distribution in Figure 6. Note that sensitivities outside
the horizontal range are rounded to the last bin.

###  3.3 Sensitivities in the absence of tipping

Considering the same model (15,17) with wandering CO2 (22), we use different
parameters to contrast the results in the previous sections with cases where
there is no bistability. In particular we consider parameters as in Table 1
except for:

  * •

Default albedo contrast: $\alpha_{1}=0.52$, $\alpha_{2}=0.47$ (i.e. also as in
Table 1).

  * •

Low albedo contrast: $\alpha_{1}=0.50$, $\alpha_{2}=0.48$.

  * •

No albedo contrast: $\alpha_{1}=\alpha_{2}=0.495$.

Figure 8 shows the low and no albedo contrast cases, comparing to the default
case Figure 5(left). For the low albedo contrast case there is no longer a
region of bistability, but there is nontrivial variation of the instantaneous
sensitivity along the attractor. For the no albedo contrast case, the
instantaneous sensitivity is close to constant. The projection of the climate
attractor into the $(\Delta R,T)$ plane is shown in Figure 9(a-c), while the
corresponding distribution of two-point sensitivities in Figure 9(d-f).
Observe the presence of non-unimodal distributions for (d,e) associated with
regions with different two-point sensitivity, and clear skewness and tails
again associated with the geometry of the measure in (a,b). Note the higher
average in (e) corresponds to there being only a single regime in (b) that
runs over a wide range of temperatures.

In physical terms, the skewness (and long tails) in (d,e) originate from the
state-dependence and nonlinearity of feedbacks (i.e. non-constant feedback
factors). The bistability of the two regimes with different feedbacks gives
the two peaks in the distribution of Figure 9(d). However, Figure 9(e) still
has two peaks: these originate from state dependence on the same attractor
(Figure 9(b). For this low-albedo-contrast case, there is no ‘tipping’ but we
still find very non-Gaussian distribution of $S$ that comes from
nonlinearities in the system that, in this case, do not produce tipping. Note
that only in the no albedo contrast case (c) is there a plausible fit to
Gaussian.

Figure 8: Response of global mean temperature to the same realisation of
wandering CO2 for the energy balance model (17) with (left) low albedo
contrast and (right) no albedo contrast; (a-e) as in Figure 5. In both cases
there is no region of bistability. Observe there are large but bounded
fluctuations of the instantaneous sensitivity in case (left), indicating
state-dependency but no tipping. In case (right) there is comparatively little
fluctuation of sensitivity; compare with Figure 5(left) for the default
parameters with a region of bistability. Figure 9: (a-c) Plots of $T$ vs
$\Delta R$ for a long timeseries of the energy balance model (17) with
wandering CO2, parameters as in Figure 6. (d-f) plots of distributions of two-
point sensitivities corresponding to (a-c) respectively. (a,d) correspond to
the default scenario with bistability, (b,e) to the low albedo contrast and
(c,f) to the no albedo contrast scenario; cf Figure 8. Observe the presence of
two peaks (corresponding to state dependent sensitivities on the two branches)
and a large mass of extremes for (d). (e) shows a range of state-dependent
variation but relatively few extremes; while (f) resembles a normal
distribution. Observe that (d) can be decomposed into the conditional
sensitivities shown in Figure 7.

##  4 Discussion

We demonstrate that state-dependence and the presence of tipping points
produces signatures in the distribution of instantaneous and two-point notions
of ECS. We explore this using a global energy balance model where state-
dependence and multistability originate from the dependence of both albedo and
emissivity on temperature.

For the deterministic version of our model (15) with fixed CO2 the changes in
albedo mean there can be bistability between regimes, while the changes in
emissivity contribute to different sensitivities within these regimes. The
distribution of ECS comes from several sources - nonlinearities that result in
tipping points and/or state-dependence of the feedbacks and sub-grid
variability that we model here as stochastic perturbations. Such regime-
dependent sensitivity and extremes associated with tipping points are also
visible in the more complex Gildor and Tziperman model [17, 38], as outlined
in Appendix A.

For the stochastic model (17) with wandering CO2, regime-dependent sensitivity
is visible as differences in slope of the stable regimes for the $T$ vs
$\Delta R$ plots (see Figure 4(b)). The densities of the stable regimes for
the $T$ vs $\Delta R$ plots (see Figure 6) show varying slopes and so
conditional two-point sensitivities for the two regimes (see e.g. Figure 7(a)
and (c)) can have peaks at different sensitivities. We compare several notions
of sensitivity. These are the instantaneous sensitivity associated with the
slope of the equilibrium branch, incremental sensitivity associated with a
fixed delay, and two-point sensitivity that compares arbitrary points on the
climate attractor. The presence of tipping points gives extremes of
sensitivities in that (i) there are large fluctuations of instantaneous
sensitivities for the nearest equilibrium just before a tipping point (see
Figure 5); and (ii) in the distributions of conditional two-point
sensitivities that cross regimes (see Figures 7(d) and (f)). It is remarkable
that the two-point sensitivities are so informative, given that they compare
points that may be very far from each other in time.

There remains much work to be done to understand the relation between and
limitations of these notions of ECS, and indeed ECS calculated in other ways,
for example from instantaneous doubling of CO2: this will involve transient
non-equilibrium processes due to ocean thermal inertia. Note that determining
ECS from palaeoclimate time-series [39], the two-point notion clearly has an
advantage that we are not limited by the time-resolution of the time series.

The energy balance model can be criticised as being very simple and hard to
parametrize in terms of the various physical processes that contribute to
albedo and emissivity. Moreover, we consider CO2 in (17) purely as a forcing
term which ignores known land surface and ocean processes where temperature is
known to affect CO2 balance. However, the model is complex enough to confirm
that extremes in ensembles of computed climate sensitivities can indicate
nearby tipping points. Computations presented in Appendix A confirm this
picture in a more complex box-model for the glacial cycles, where the CO2 is
modelled dynamically.

###  4.1 Future perspectives

When (and how) extremes of sensitivity can be effective precursors of a
tipping event will depend on a number of factors. In particular, the timescale
of dynamics of the climate response needs to be faster than the timescale of
changes in forcing. Figure 5(left) shows that as CO2 variability is rapid,
this results in tipping points with little precursor visible in changes to
$AR1$, though it is visible in the instantaneous sensitivity. There may be
‘rate-induced’ tipping points [1, 2] that appear when the timescale of the
forcing interacts with that of the system. The size of the region of effective
nonlinearity can also vary. For (17) this is affected by the temperature
scales $T_{\alpha}$ and $T_{\epsilon}$ over which the albedo and emissivity
changes occur. Note also that although tipping points do give rise to extremes
in the distribution of ECS, extremes do not necessarily indicate a tipping
point.

Translating these results to a time-dependent setting and to more complex
models will be difficult: the possible states in the ‘climate attractor’ and
the associated invariant measure $\mu$ is harder to define in the presence of
non-stationary temporal variation of forcing, or for large recurrence times
and a variety of nonlinear multiscale processes. In such cases, interpreting
transitions as tipping points is a challenge; nonetheless, the palaeoclimate
record does show a variety of large and sudden transitions [26]. For example,
ice core/ocean core records indicate repeated sudden changes in (regional)
surface temperature associated with glacial cycles [20] or Dansgaard-Oeschger
[13] events as well as global transitions, for example the greenhouse-icehouse
transition at the Eocene-Oligocene transition [10]. Although glacial cycles
can be found in models such as [17] as relaxation oscillation with clear
regimes, for climate reconstruction data these regimes are not so clear (e.g.
[24, 15]).

Recent work [35] suggests we are at a crossroads in terms of the future earth
system state. On the one hand, looking at the palaeoclimate record for the
last 1 million years suggests that we are overdue descent into an ice age. On
the other hand, comparison of anthropogenic CO2 emissions with the palaeo
record suggest the next tipping point may be to much warmer ‘hothouse’ earth.
A better understanding of improved indicators such as two-point ECS and what
they say about the climate response to changes in greenhouse gases, together
with a better understanding of hothouse earth climate states that may have
existed in the past (e.g. the Palaeocene climate [21, 3]) should help our
understanding and guide future generations in their need to avoid dangerous
climate change.

### Acknowledgements

We thank Peter Cox, Henk Dijsktra, Valerio Lucarini, Michael Ghil and Georg
Gottwald for interesting comments and suggestions concerning this work. We
thank the EU ITN network CRITICS and CliMathNet for opportunities to discuss
and develop this research. AvdH thanks the program of the Netherlands Earth
System Science Centre (NESSC), financially supported by the Dutch Ministry of
Education, Culture and Science (OCW). This project is TiPES contribution #7:
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 820970.

### Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] Peter Ashwin, Clare Perryman, and Sebastian Wieczorek.  Parameter shifts for nonautonomous systems in low dimension: bifurcation-and rate-induced tipping.  Nonlinearity, 30(6):2185, 2017.
  * [2] 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, 2012.
  * [3] Michiel L J Baatsen, Anna S von der Heydt, M Huber, Michael A Kliphuis, Peter K Bijl, Appy Sluijs, and Henk A. Dijkstra.  Equilibrium state and sensitivity of the simulated middle-to-late Eocene climate.  Clim. Past Disuss., in review, 2018.
  * [4] Jonah Bloch-Johnson, Raymond T. Pierrehumbert, and Dorian S. Abbot.  Feedback temperature dependence determines the risk of high warming.  Geophysical Research Letters, 42:4973–4980, 2015.
  * [5] Mikhail I. Budyko.  The effect of solar radiation variations on the climate of the earth.  tellus, 21(5):611–619, 1969.
  * [6] Rodrigo Caballero and M Huber.  State-dependent climate sensitivity in past warm climates and its implications for future climate projections.  Proceedings of the National Academy of Sciences, 110(35):14162–14167, 2013.
  * [7] J. Charney.  Carbon Dioxide and Climate: A Scientific Assessment.  The National Academies Press, Washington, DC, 1979.
  * [8] Mickaël D. Chekroun, Eric Simonnet, and Michael Ghil.  Stochastic climate dynamics: Random attractors and time-dependent invariant measures.  Physica D: Nonlinear Phenomena, 240(21):1685 – 1700, 2011.
  * [9] Peter M. Cox, Chris Huntingford, and Mark S. Williamson.  Emergent constraint on equilibrium climate sensitivity from global temperature variability.  Nature, 553(7688):319, 2018.
  * [10] H.K Coxall, P.A. Wilson, H. Pälike, C.H. Lear, and J. Backman.  Rapid stepwise onset of antarctic glaciation and deeper calcite compensation in the pacific ocean.  Nature, 433:53–57, 2005.
  * [11] Mark M. Dekker, Anna S. von der Heydt, and Henk A. Dijkstra.  Cascading transitions in the climate system.  Earth System Dynamics, 9(4):1243–1260, 2018.
  * [12] H.A. Dijkstra and J.P. Viebahn.  Sensitivity and resilience of the climate system: A conditional nonlinear optimization approach.  Communications in Nonlinear Science and Numerical Simulation, 22(1-3):13–22, 2015.
  * [13] Peter D. Ditlevsen, Mikkel S. Kristensen, and Katrine K. Andersen.  The recurrence time of dansgaard–oeschger events and limits on the possible periodic component.  Journal of Climate, 18(14):2594–2603, 2005.
  * [14] Sybren Drijfhout, Sebastian Bathiany, Claudie Beaulieu, Victor Brovkin, Martin Claussen, Chris Huntingford, Marten Scheffer, Giovanni Sgubin, and Didier Swingedouw.  Catalogue of abrupt shifts in Intergovernmental Panel on Climate Change climate models.  Proceedings of the National Academy of Sciences of the United States of America, 112(43):E5777–86, 2015.
  * [15] Tobias Friedrich, Axel Timmermann, Michelle Tigchelaar, Oliver Elison Timm, and Andrey Ganopolski.  Nonlinear climate sensitivity and its implications for future greenhouse warming.  Science Advances, 2(11):e1501923–e1501923, 2016.
  * [16] Michael Ghil.  Climate stability for a sellers-type model.  Journal of the Atmospheric Sciences, 33(1):3–20, 1976.
  * [17] H. Gildor and E. Tziperman.  A sea ice climate switch mechanism for the 100-kyr glacial cycles.  Journal of Geophysical Research: Oceans, 106(C5):9117–9133, May 2001.
  * [18] Jonathan M. Gregory.  A new method for diagnosing radiative forcing and climate sensitivity.  Geophysical Research Letters, 31(3):147–4, 2004.
  * [19] Mátyás Herein, János Márfy, Gábor Drótos, and Tamás Tél.  Probabilistic concepts in intermediate-complexity climate models: A snapshot attractor picture.  Journal of Climate, 29(1):259–272, 2016.
  * [20] Andrew McC. Hogg.  Glacial cycles and carbon dioxide: A conceptual model.  Geophysical research letters, 35(1), 2008.
  * [21] M Huber and Rodrigo Caballero.  The early Eocene equable climate problem revisited.  Clim. Past, 7(2):603–633, 2011.
  * [22] IPCC.  Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change.  Cambridge University Press, 2013.
  * [23] Reto Knutti and Gabriele C. Hegerl.  The equilibrium sensitivity of the earth’s temperature to radiation changes.  Nature Geoscience, 1:735–743, 2008.
  * [24] Peter Koehler, Lennert B. Stap, Anna S. von der Heydt, Bas de Boer, Roderik S.W. Van De Wal, and J. Bloch-Johnson.  A State-Dependent Quantification of Climate Sensitivity Based On Paleodata of the Last 2.1 Million Years.  Paleoceanography, 32:1102–1114, 2017.
  * [25] P. Köhler, B. de Boer, A.S. von der Heydt, L.B. Stap, and R.S.W. van de Wal.  On the state dependency of the equilibrium climate sensitivity during the last 5 million years.  Climate of the Past, 11(12):1801–1823, 2015.
  * [26] Timothy M. Lenton.  Early warning of climate tipping points.  Nature climate change, 1(4):201, 2011.
  * [27] Timothy M. Lenton, Hermann Held, Elmar Kriegler, Jim W. Hall, Wolfgang Lucht, Stefan Rahmstorf, and Hans Joachim Schellnhuber.  Tipping elements in the Earth’s climate system.  Proceedings of the National Academy of Sciences of the United States of America, 105(6):1786–1793, 2008.
  * [28] Valerio Lucarini.  Revising and extending the linear response theory for statistical mechanical systems: Evaluating observables as predictors and predictands.  Journal of Statistical Physics, 173:1698–1721, 2018.
  * [29] Daniel J. Lunt, Alan M. Haywood, Gavin A. Schmidt, Ulrich Salzmann, Paul J. Valdes, and Harry J. Dowsett.  Earth system sensitivity inferred from Pliocene modelling and data.  Nature Geoscience, 3(1):60–64, 2010.
  * [30] Patrik L. Pfister and Thomas F. Stocker.  State-Dependence of the Climate Sensitivity in Earth System Models of Intermediate Complexity.  Geophysical Research Letters, 44(20):10,643–10,653, 2017.
  * [31] Gerard H. Roe and Marcia B. Baker.  Why is climate sensitivity so unpredictable?  Science, 318(5850):629–632, 2007.
  * [32] Eelco J. Rohling, Appy Sluijs, Henk A. Dijkstra, Peter Köhler, Roderick SW van de Wal, Anna S von der Heydt, David J Beerling, André Berger, Peter K Bijl, Michel Crucifix, et al.  Making sense of palaeoclimate sensitivity.  Nature, 491(7426):683, 2012.
  * [33] Stephen E Schwartz.  Determination of earth’s transient and equilibrium climate sensitivities from observations over the twentieth century: strong dependence on assumed forcing.  Surveys in geophysics, 33(3-4):745–777, 2012.
  * [34] William D. Sellers.  A global climatic model based on the energy balance of the earth-atmosphere system.  Journal of Applied Meteorology, 8(3):392–400, 1969.
  * [35] Will Steffen, Johan Rockström, Katherine Richardson, Timothy M. Lenton, Carl Folke, Diana Liverman, Colin P. Summerhayes, Anthony D. Barnosky, Sarah E. Cornell, Michel Crucifix, Jonathan F. Donges, Ingo Fetzer, Steven J. Lade, Marten Scheffer, Ricarda Winkelmann, and Hans Joachim Schellnhuber.  Trajectories of the earth system in the anthropocene.  Proceedings of the National Academy of Sciences, 115(33):8252–8259, 2018.
  * [36] Alexis Tantet, Valerio Lucarini, Frank Lunkeit, and Henk A. Dijkstra.  Crisis of the chaotic attractor of a climate model: a transfer operator approach.  Nonlinearity, 31(5):2221–2251, apr 2018.
  * [37] Alexis Tantet, Fiona R. van der Burgt, and Henk A. Dijkstra.  An early warning indicator for atmospheric blocking events using transfer operators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(3):036406, 2015.
  * [38] Anna S. von der Heydt and Peter Ashwin.  State dependence of climate sensitivity: attractor constraints and palaeoclimate regimes.  Dynamics and Statistics of the Climate System, 1(1):1–21, February 2017.
  * [39] Anna S. von der Heydt, Henk A. Dijkstra, Roderik S.W. Van De Wal, Rodrigo Caballero, Michel Crucifix, Gavin L. Foster, M. Huber, Peter Koehler, E. J. Rohling, Paul J. Valdes, Peter Ashwin, Sebastian Bathiany, T. Berends, L. van Bree, Peter D. Ditlevsen, Michael Ghil, Alan M. Haywood, Joel Katzav, Gerrit Lohmann, J. Lohmann, Valerio Lucarini, A. Marzocchi, H. Pälike, I. Ruvalcaba Baroni, D Simon, Appy Sluijs, L.B. Stap, A. Tantet, J.P. Viehbahn, and Martin Ziegler.  Lessons on Climate Sensitivity From Past Climate Changes.  Current Climate Change Reports, 2(4):148–158, 2016.
  * [40] Anna S. von der Heydt, Peter Koehler, Roderik S.W. Van De Wal, and Henk A. Dijkstra.  On the state dependency of fast feedback processes in (paleo) climate sensitivity.  Geophysical Research Letters, 41(18):6484–6492, 2014.
  * [41] I. Zaliapin and M. Ghil.  Another look at climate sensitivity.  Nonlinear Processes in Geophysics, 17(2):113–122, 2010.
  * [42] Xiaozhu Zhang, Christian Kuehn, and Sarah Hallerberg.  Predictability of critical transitions.  Phys. Rev. E, 92:052905, Nov 2015.

##  Appendix A Appendix: Tipping in a multi-box climate model with ocean
biogeochemistry

To investigate the notion of two-point sensitivity in a more complex and
physics based earth system model, we explore here the Gildor-Tziperman (2001)
model [17] with Milaknovich forcing and biogeochemistry in the ocean. In this
model the ocean and atmosphere consists of 4 latitudinal (zonally averaged)
boxes each, with two layers in the ocean and one in the atmosphere. Within the
boxes a variety of thermodynamic quantities (eg temperature $T$) and species
(eg ocean salinity $S$ and atmospheric CO2) are modelled. In addition there is
dynamic land ice (slowly evolving) and sea ice (rapidly evolving) as well as
fluxes that join these boxes. The system is sufficiently complex to allow
modelling of the Pleistocene ice-age oscillations of land-ice in response to
Milankovich forcing. The glacial-interglacial cycles appear in this model as
internally generated self-sustained oscillations, which are then modified by
the Milankovich forcing. More details are given in [17, 38].

Figure 10 shows projections of a long trajectory (500 kyr) on the climate
attractor of this more complex climate model onto the plane of global mean
temperature $T$ against (a) $\Delta R_{[CO_{2}]}$ (b) $\Delta R_{[LI]}$ and
(c) $\Delta R_{[CO_{2},LI]}$ [38]. Observe that all three projections clearly
show two climate regimes, a lower ‘cold’ state (corresponding to large amounts
of sea ice and $T<12.28C$) and an upper ‘warm’ state (corresponding to
$T>12.28C$). The projection on the combined radiative forcing of CO2 and the
slow feedback in land ice changes $\Delta R_{[LI]}$ shows a clear slope and
hence ’mean’ sensitivities in both regimes (Figure 10c). Following the
formalism of [32], the slopes in (a) should reflect the Earth System
sensitivity, while (c) should give a good approximation of the Charney or
equilibrium climate sensitivity. Note that in this model there is only one
slow feedback to correct for, namely the land-ice albedo feedback. When
projecting onto the $\Delta R_{[CO_{2}]}$ plane (Figure 10a) the cold regime
appears very diffuse and with very high (or sometimes negative) Earth System
sensitivities, suggesting that in the cold branch local climate dynamics are
not entirely determined by CO2. Similarly, when projecting only on land ice
changes $\Delta R_{{LI}}$ (Figure 10b) the warm branch appears more diffuse,
i.e. CO2 dynamics seem to be more important than land ice dynamics on this
branch.

Figure 10: Plots of $T$ vs $\Delta R$ for the Gildor-Tziperman model [17]; see
[38]. (a) projection onto the $(\Delta R_{[CO_{2}]},T)$-plane, which should
reflect the Earth System Sensitivity only considering CO2 as forcing; (b)
projection onto the $(\Delta R_{[LI]},T)$-plane, considering only the slow
land-ice albedo feedback as forcing; (c) projection onto the $(\Delta
R_{[CO_{2},LI]},T)$-plane, considering both CO2 and the slow feedbacks (in
this model there is only the land-ice albedo feedback slow) as forcing.
Following the formalism of [32] the slopes in this graph should reflect a good
approximation of the Charney ECS.

The left column of Figure 11 (a-d) shows the distribution of two-point climate
sensitivities for $R_{[CO_{2},LI]}$ conditional on regime: this is comparable
to Figure 7(a-d) in that one can observe (i) clearly localised distributions
of ECS in (a,c) conditional on remaining within the W or C regime, (ii) a
broader distribution in (c): this seems to be associated with the curvature of
the C regime branch in Figure 10c), (iii) a clear association of tipping from
W to C (b) or from C to W (d) being associated with extreme sensitivities.
Note that for this model there is no energy balance model available and so it
is not possible to compute the instantaneous sensitivity. For comparison we
show in the right column of Figure 11 (e-h) the same distributions for the
Earth System Sensitivities (ESS), which are not compensated for the slow
feedback (in this model the land-ice albedo feedback). Observe that for both
regimes there is a much broader distribution for the ESS (e and g) than for
the ECS in (a) and (c). In particular in the cold regime, earth system
sensitivities (g) are much higher than equilibrium sensitivities (c) because
the land-ice albedo feedback is very strong in the cold (land-ice covered)
states.

Figure 11: Conditional sensitivities and probabilities of tipping for glacial
cycles simulations with the Gildor-Tziperman model [17] (as in Figure 7 for
the energy balance model). (a-d) Conditional two-point sensitivities
compensating for slow feedbacks (i.e. reflecting ECS) and probabilities of
tipping (i.e. from the distribution shown in Figure 10c). (e-h) Conditional
earth system sensitivities ESS (not compensated for slow feedbacks) and
probabilities of tipping (i.e. from the distribution shown in Figure 10a).


11institutetext: Viny Saajan Victor 22institutetext: Fraunhofer ITWM,
Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany, 22email:
<EMAIL_ADDRESS>33institutetext: Manuel Ettmüller
44institutetext: Fraunhofer ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern,
Germany 44email<EMAIL_ADDRESS>55institutetext: Andre
Schmeißer 66institutetext: Fraunhofer ITWM, Fraunhofer-Platz 1, 67663
Kaiserslautern, Germany 66email<EMAIL_ADDRESS>77institutetext: Heike Leitte 88institutetext: Technical University of
Kaiserslautern, Gottlieb-Daimler-Straße 47, 67663 Kaiserslautern, Germany
88email<EMAIL_ADDRESS>99institutetext: Simone Gramsch
1010institutetext: Fraunhofer ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern,
Germany 1010email<EMAIL_ADDRESS>
# Machine Learning Based Optimization Workflow for Tuning Numerical Settings
of Differential Equation Solvers for Boundary Value Problems

Viny Saajan Victor  Manuel Ettmüller  Andre Schmeißer  Heike Leitte and Simone
Gramsch

###### Abstract

Several numerical differential equation solvers have been employed effectively
over the years as an alternative to analytical solvers to quickly and
conveniently solve differential equations. One category of these is boundary
value solvers, which are used to solve real-world problems formulated as
differential equations with boundary conditions. These solvers require certain
numerical settings to solve the differential equations that affect their
solvability and performance. A systematic fine-tuning of these settings is
required to obtain the desired solution and performance. Currently, these
settings are either selected by trial and error or require domain expertise.
In this paper, we propose a machine learning-based optimization workflow for
fine-tuning the numerical settings to reduce the time and domain expertise
required in the process. In the evaluation section, we discuss the
scalability, stability, and reliability of the proposed workflow. We
demonstrate our workflow on a numerical boundary value problem solver.

##  0.1 Introduction

Differential equations are one of the most effective mathematical tools for
understanding and predicting the behavior of dynamical systems in nature,
engineering, and society. The study of differential equations entails learning
how to solve them and interpret the solutions. In the field of differential
equations, a boundary value problem (BVP) is a type of differential equation
with a set of additional constraints known as boundary conditions. The
solution to a boundary value problem must satisfy the boundary conditions.
Despite the fact that analytical methods for solving boundary value problems
yield exact answers, they become difficult to apply to complex problems.
Hence, several numerical methods have become popular, leveraging the
development of computing capabilities. These methods provide approximate
solutions that have sufficient accuracy for engineering purposes. The boundary
value problem solvers developed based on the numerical methods usually consist
of settings that describe method-specific properties such as error tolerance,
mesh size, etc. These settings have an impact on the solvability and
performance of the solvers, and the impact varies depending on the problem. In
order to obtain qualitatively good solutions within acceptable solution times,
fine-tuning these settings is required. However, these settings are now tuned
manually by trial and error, which is ineffective due to the large parameter
space of the settings and the complex interactions between them.

Machine learning (ML) techniques have achieved remarkable success in recent
years due to their ability to identify patterns and structures in the given
data. Since ML models can predict in real-time, they have been deployed in a
variety of industries. In earlier work, we showed how ML-based quality
optimization is effective in the technical textile industry victor2022visual .
In this paper, we propose a two-stage optimization workflow based on ML to
make the process of tuning the numerical settings of boundary value problem
solvers more effective. In the first stage of the workflow, we design an ML
pipeline that maps the numerical settings to solver performance. The solver
performance metric includes solvability status and computational cost
statistics of the solver, such as the number of evaluations and mesh points
required by the solver and the corresponding residual error obtained from the
solver. In the second stage of the workflow, we use the trained ML model to
predict the influence of the numerical settings on the success and performance
of the solver. We use this knowledge to propose a multi-criteria optimization
strategy for fine-tuning the settings, thereby reducing the domain expertise
required in the process.

##  0.2 Foundations: The test bench for boundary value problem solvers

In collaborations, we often work together with domain experts in the technical
textile industry such as mathematicians or process engineers. They want to
model and simulate a variety of fiber formation processes such as polymer melt
spinning ettmuller2021flow  or the production of glass wool arne2011fluid .
These processes are typically modeled as a system of differential equations
that require a numerical solver. In the special case of ordinary differential
equations, the bvp4c solver from MATLAB is frequently used when dealing with
boundary value problems kierzenka2001bvp . The bvp4c solver is a finite
difference method based on a collocation method that provides a continuous
solution that is fourth-order accurate in the interval of integration. The
integration interval is thereby subdivided into smaller intervals using a mesh
of points. The subdivision scheme is adaptive, so further points can be both
added and removed. The boundary conditions and the collocation conditions
imposed on all of the sub-intervals create a global system of algebraic
equations. These equations are solved using the simplified Newton method to
arrive at a numerical solution. For the formulation of the ML problem
statement and demonstration of the proposed approach, we used a self-
implemented C++ version of the solver above, which uses a Newton-Armijo method
instead of a simplified Newton method for solving the system of non-linear
equations.

Figure 1: Input and output parameters of BVP solver.

The solvability status and the computational resources required by the BVP
solver are mainly affected by the type of differential equations, boundary
conditions, the initial guess of the solution, and the solver-specific options
as shown in Figure 1. In our workflow, we keep the initial guess for the
solver constant, as we are interested in analyzing the effect of solver-
specific options. We also want to analyze the correlation between the optimal
choice of numerical settings and the complexity of the boundary value problem.
For this, we selected ten different reference problems for our analysis. Out
of all the test problems defined in soetaert2010package , we have chosen four
linear (1, 3, 4, and 7) and six non-linear (19, 20, 22, 23, 24, and 33)
problems, since we expect to gain more insights from the more complex
differential equations.

##  0.3 Machine learning based optimization workflow for tuning the BVP solver
settings

In this section, we elaborate on the proposed workflow for tuning the
numerical settings of the BVP solvers. The work consists of the machine
learning process and the optimization process. The machine-learning process is
used to map the numerical settings of the BVP solvers to their performance.
This pipeline comprises a binary classifier and a multiple-output regression
model. The binary classifier is used to predict the solvability status of the
boundary value problem based on the numerical settings. The settings predicted
to be solvable by the classifier are then fed into the multiple-output
regression model which forecasts solver performance. Both machine learning
models rely on the training data collected for the reference boundary value
problems. The best models are chosen based on the accuracy obtained from the
testing data. The chosen models are used to solve the multi-objective
optimization problem to obtain the desired combination of numerical settings
that optimize the quality criteria. We discuss the details of feature
selection, data collection, model selection, training, and testing of ML
models in the following subsections.

Figure 2: Workflow of machine learning-based optimization for tuning the
numerical settings of BVP solvers. Different stages are colored corresponding
to the section numbers in this chapter.

###  0.3.1 Feature Selection

Our feature selection process involved discussion with the domain experts to
incorporate their experience in choosing the features based on the type of
boundary value problems under inspection. The selected features are also
confirmed using exploratory data analysis and the machine learning feature
selection process. In the former, the correlations of the features within
themselves and to the output features are considered, and in the latter,
methods such as feature elimination and dimensionality reduction are utilized.
We found a good agreement between the expert-based and ML-based feature
selection processes.

#### Input features

We identified nine input features that affect the solvability and performance
of the boundary value problem solvers. As discussed before, the BVP solver
under consideration internally uses the Newton-Armijo method to solve the
system of algebraic equations. Therefore, the input features consist of the
generic BVP solver settings and the specific Newton-Armijo solver settings.
The selected input features are as follows.

  1. 1.

_Test Case Type_ is the type of the boundary value problem chosen from the
test bench explained before. This is a categorical parameter where each
category represents a problem from soetaert2010package .

  2. 2.

_Newton Tolerance_ handles the termination of the Newton solver. If the
residual of the equation system is below the _NewtonTolerance_ the algorithm
is considered converged.

  3. 3.

_Newton Critical Tolerance_ is an additional tolerance that is considered when
the solver does not reach the more restrictive _NewtonTolerance_. If the
residual is below the _NewtonCriticalTolerance_ after the maximum number of
iterations it is still considered as converged, otherwise as diverged.

  4. 4.

_Newton Maximum Iterations_ is the maximum number of iterations for the Newton
solver.

  5. 5.

_Newton Armijo Probes_ is the maximum Armijo probes used in the Newton-Armijo
algorithm.

  6. 6.

_Maximum Grid Points_ is the maximum number of mesh points allowed when
solving the boundary value problem.

  7. 7.

_Add Factor_ is the feature that controls the addition of new mesh points
based on the difference between computed residual in the interval and the
relative tolerance of the solver.

  8. 8.

_Remove Factor_ is the feature that controls the removal of mesh points if it
is too fine based on the difference between the computed residual in the
interval and the relative tolerance of the solver.

  9. 9.

_Use Collocation Scaling_ is the boolean variable that specifies whether the
inner (non-boundary) conditions are scaled regarding the number of mesh
points.

#### Output features

The output features represent the solvability status of the boundary value
problem and the statistics that affect the time taken to solve the problem and
the resources required. The four chosen output features are:

  1. 1.

_Success Status_ is the boolean feature that specifies the solvability status
of the BVP solver.

  2. 2.

_Number of ODE Evaluations_ is the number of times the solver calls the right-
hand side function of the differential equation in the algorithm.

  3. 3.

_Number of Grid Points_ is the number of grid points of the numerical
solution.

  4. 4.

_Maximum Residuum_ is the maximum residual of the numerical solution.

###  0.3.2 Data Collection

As discussed before, a BVP solver based on MATLAB’s bvp4c was used to collect
the data for training the machine learning models. We set realistic ranges for
the input features based on discussions with the domain experts. Latin
Hypercube Sampling was used as the sampling technique for creating the input
data as this technique distributes samples evenly over sample space. Based on
the boundary value problems under investigation and the domain expertise we
chose either linear scale or logarithmic scale for the sampling. Table 1 shows
the ranges and scales used for each input feature for sampling. The generated
input data are fed to the solver to generate the output labels. We generated
100,000 data points for each of the ten reference problems.

Table 1: Ranges and scales used for sampling input features. The table also shows the default values of the features. Input Feature |  Feature Range |  Default Value |  Sampling Scale
---|---|---|---
Maximum Grid Points |  $[100,10000]$ |  $1000$ |  Linear
Newton Critical Tolerance |  $[${10}^{-12}$,100]$ |  $1$ |  Logarithmic
Newton Armijo Probes |  $[1,10]$ |  $2$ |  Linear
Newton Maximum Iterations |  $[1,100]$ |  $4$ |  Linear
Newton Tolerance |  $[${10}^{-12}$,${10}^{-2}$]$ |  ${10}^{-12}$ |  Logarithmic
Add Factor |  $[1,1000]$ |  $100$ |  Logarithmic
Remove Factor |  $[0,2]$ |  $0.5$ |  Linear
Use Collocation Scaling |  _{ True, False}_ |  _True_ |  Discrete
Test Case Type |  Ten discrete problem types |  — |  Discrete

###  0.3.3 Model Selection

Based on our problem statement, we evaluated several binary classification and
multi-output regression algorithms. The selected algorithms are classified
into five major categories which are support vector machine, neighbors-based,
boosting, bagging, and artificial neural net-based algorithms. One best
algorithm from each category is selected based on the accuracy of the test
dataset. The description of the chosen categories is listed as follows.

  1. 1.

Support Vector Machine (SVM) noble2006support  is a supervised learning
algorithm that finds the best hyperplane to separate the features into
different domains. They solve the non-linear problems with the help of kernels
which implicitly map the input features into a high dimensional space.

  2. 2.

Boosting algorithms are based on the gradient boosting technique, in which the
classification model is constructed as an ensemble of weak models. It is a
sequential process, where each subsequent model attempts to correct the errors
of the previous model. We selected the light gradient boosting machine (LGBM)
algorithm ke2017lightgbm  that includes sparse optimization, parallel
training, multiple loss functions, regularization, and early stopping
techniques to provide an accurate and generalized model.

  3. 3.

Bagging algorithms are based on ensemble techniques, where the dataset is
randomly sampled into subsets and each of the subsets is trained using
different models. The predictions are calculated by combining individual
subset predictions. Random Forest (RF) breiman2001random  is selected in this
category which is a bagging or bootstrap aggregation technique based on
decision trees.

  4. 4.

Neighbors-based algorithms work on the fundamental idea of finding the $n$
training samples closest in distance to a new point and assigning the label
from these. Euclidean distance is the most common metric used to measure
distance. These are non-parametric algorithms that do not make strong
assumptions about the form of the mapping function. Hence, they frequently
work in scenarios where the decision border is not regular. We selected the
K-nearest neighbor algorithm fix1989discriminatory , cover1967nearest , that
works based on the $K$ nearest neighbors of each data point.

  5. 5.

Artificial neural network (ANN) is a fully connected feed-forward network that
makes use of an aggregation of functions to understand and convert input data
that is presented in one form into the desired output. A traditional neural
network’s fundamental building blocks are nodes, which are arranged into
layers. To get the final forecast, the input characteristics are transferred
through these layers using a series of non-linear operations. For binary
classification, we used a multilayer perceptron (MLP) of two hidden layers
with 100 neurons each. For regression, we used an MLP comprising six hidden
layers with 250 neurons each.

###  0.3.4 Training and Testing

We divided the dataset into a $80\text{\,}\mathrm{\char 37\relax}$ training
data and a $20\text{\,}\mathrm{\char 37\relax}$ testing set. The training data
is further divided into a training set and a validation set. The training set
is used to train the regression models and the validation set is used to tune
the model hyper-parameters. The testing set is used for unbiased evaluation of
the model. As a part of data preparation, we performed input feature scaling
and output feature transformation to tailor the data for the machine learning
models. During the exploratory data analysis, we observed that the output data
had a skewed distribution with chunks of repeated values and outliers. Hence
we transformed the distribution to follow a normal distribution that spread
out the most frequent values and also reduces the impact of outliers.

###  0.3.5 Evaluation Metrics

In this section, we list the metrics that are used to evaluate the
classification and regression models. The metrics are selected based on the
type of algorithms, scale, and distribution of the input and output features
that reveal the fairness, robustness, and reliability of ML models.

#### Binary Classification

For the equations used below, we define TP, FP, TN, and FN as true positives,
false positives, true negatives, and false negatives respectively.

  1. 1.

Accuracy: It measures how many observations were correctly classified overall.
This is expressed as a percentage. Accuracy is calculated as below.

$\displaystyle P$ $\displaystyle=\frac{TP+TN}{TP+FP+TN+FN}$

  2. 2.

Precision (P): It measures the proportion of positive classifications by the
model that were actually correct. Precision is calculated as below.

$\displaystyle P$ $\displaystyle=\frac{TP}{TP+FP}$

  3. 3.

Recall (R): It measures the proportion of actual positives classifications
that were identified correctly by the model. Recall is calculated as below.

$\displaystyle R$ $\displaystyle=\frac{TP}{TP+FN}$

#### Multi-output Regression

For the equations used below, we define $n$ as the number of data points,
$y_{i}$, and $\hat{y}_{i}$ as the actual and predicted values respectively for
the data point $i$ and $\bar{y}$ as the mean value of the data points.

  1. 1.

Mean absolute percentage error (MAPE) is a statistical measure to evaluate the
accuracy of a regression model. The error is independent of the scale of the
output as it measures the accuracy as a percentage. MAPE is calculated as
below.

$\displaystyle MAPE$
$\displaystyle=\frac{100}{n}\sum_{i=1}^{n}\left|\frac{y_{i}-\hat{y}_{i}}{y_{i}}\right|$

  2. 2.

Root Mean squared error (RMSE) measures the square root of the average squared
errors. The RMSE is a good estimate for ensuring that the ML model has no
outlier predictions with large errors since it puts a higher weight on these
errors due to the squaring. RMSE is calculated as below.

$\displaystyle RMSE$
$\displaystyle=\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{n}}$

  3. 3.

Coefficient of determination ($R^{2}$ Score) is the measure of how close the
data points are to the fitted regression line. It explains how much of the
variance of actual data is explained by the predicted values. $R^{2}$ Score is
calculated as below.

$\displaystyle R^{2}(y,\hat{y})$
$\displaystyle=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y}_{i})^{2}}$

Table 2: Functions provided by the bvpTune library Function |  Description
---|---
getSolvabiltyStatus(TestCaseType, Setting) |  returns whether the numerical setting is solvable or not for the specified TestCaseType of the boundary value problem.
getSolverPerformance(TestCaseType, Settings) |  returns the predicted solver statistics of the numerical setting for the specified TestCaseType of the boundary value problem.
getOptimalODEevals(TestCaseType) getOptimalGridPoints(TestCaseType) getOptimalResiduum(TestCaseType) |  returns the optimal number of ODE evaluations, grid points, and residuum for the specified TestCaseType of the boundary value problem.
getOptimalODEevalsAndGridPoints(TestCaseType) getOptimalGridPointsAndResiduum(TestCaseType) getOptimalResiduumAndODEevals(TestCaseType) |  return pareto front numerical settings as a data frame that optimizes the two criteria for the specified TestCaseType of the boundary value problem.
getOptimizedSettings(TestCaseType) |  returns pareto front numerical settings as a data frame that optimizes the ODE evaluations, grid points, and the residuum for the specified TestCaseType of the boundary value problem.
visualize(dataframe) |  visualize the 2d/3d pareto front plots for the provided data frame returned by the optimization functions.

###  0.3.6 Multi-objective Optimization

After determining the best machine learning models based on the evaluation
metrics, we use them for optimization. The objective in our problem setting is
to minimize the residual error obtained from the BVP solver by optimizing the
computational resources such as the number of ODE evaluations and the grid
points. Hence, we define the problem as a multi-objective optimization as
shown in equation (1).

$\min_{x\in X}(\phi_{eval}(x),\phi_{gp}(x),\phi_{res}(x))$ (1)

subject to the constraints,

$\displaystyle f_{lb}(i)\leq x_{i}\leq f_{ub}(i),\mbox{ for }i=1,...,n$

and

$\displaystyle\phi_{cl}(x)=1$

where $\phi_{eval}$, $\phi_{gp}$, and $\phi_{res}$ are the selected ML
regression algorithms for predicting ODE evaluations, grid points, and maximum
residuum respectively. $\phi_{cl}$ is the selected ML classification
algorithm. $X$ is the feasible set of input vectors and
$x=[x_{1},x_{2},....,x_{n}]$. The functions $f_{lb}(i)$ and $f_{ub}(i)$
provide the upper and lower bounds of the $i^{th}$ input parameter search
space respectively. Along with equation (1) which optimizes all three
objectives, we also examine one and two objective optimizations to cater to
the requirements based on the different problem setups where one or two
objectives can be relaxed.

We use the Optuna library optuna_2019  for solving our multi-objective
optimization. This library provides a framework to automatically search for
the optimal input vectors and efficiently search large spaces and prune
unpromising results using state-of-the-art algorithms such as Nondominated
Sorting Genetic Algorithm II deb2002fast , the Tree-structured Parzen
Estimator Algorithm ozaki2020multiobjective ; ozaki2022multiobjective , and
the Quasi Monte Carlo Sampling Algorithm bergstra2012random , etc. We provide
the objective function defined in the equation (1) and the best machine
learning model to the Optuna framework for optimization.

Figure 3: ROC curves for different classification algorithms

##  0.4 Experiments and Results

In this section, we present the performance of the proposed workflow on the
BVP solver. The results reveal the necessity and effectiveness of our workflow
in fine-tuning the numerical settings of the BVP solvers. We further discuss
the different quality aspects related to the proposed workflow to show its
scalability, reliability, and stability. The scalability aspects are computed
on a workstation with a 40-core Intel® Xeon® E5-2680 v2 (2.80 GHz) CPU.

Table 3: Performance of classification models on the dataset ML Algorithm |  Accuracy |  Precision |  Recall
---|---|---|---
Support Vector Machine |  0.8458 |  0.8722 |  0.9355
Light Gradient Boosting Machine |  0.8593 |  0.9048 |  0.9120
Random Forest |  0.8123 |  0.8696 |  0.8874
K Nearest Neighbor |  0.8206 |  0.875 |  0.893
Artificial Neural Net |  0.8562 |  0.9023 |  0.9106

###  0.4.1 Binary Classification

Figure 4: Precision-Recall curves for different classification algorithms

The machine learning algorithms to achieve binary classification were trained
independently of the types of boundary value problems discussed in Section
0.2. We were able to generalize the models for all the considered types of
problems given to the BVP solver. Table 3 shows the performance of the binary
classification models on the test dataset. While the LGBM algorithm provided
the best results in terms of accuracy and precision, the SVM algorithm
provided the best recall. For further investigation, we plotted ROC curves of
the ML models, as seen in Figure 3. These curves show the performance of
classification models at different classification thresholds. The area under
the ROC curve (AUC) measures how well predictions are ranked, rather than
their absolute values. As observed from the figure, LGBM and ANN have better
AUCs.

Table 4: Performance of regression models on the dataset ML Algorithm |  RMSE |  MAPE |  $R^{2}$ Score
---|---|---|---
Support Vector Machine |  10658.80 |  229.40 |  0.094
Light Gradient Boosting Machine |  2760.38 |  4.30 |  0.95
Random Forest |  1870.41 |  1.31 |  0.97
K Nearest Neighbor |  8315.08 |  69.65 |  0.48
Artificial Neural Net |  3954.85 |  7.90 |  0.90

The dataset generated by the BVP solver consisted of more positive cases
(solvable settings) compared to the negative cases (non-solvable settings). To
evaluate the classification models with this imbalanced dataset, we plotted
precision-recall curves as depicted in Figure 4. Also like ROC curves,
precision-recall curves provide a graphical representation of a classifier’s
performance across many thresholds rather than a single value. The blue dotted
line represents a baseline classifier that simply predicts that all instances
belong to the positive class. As we can see from the figure, LGBM, RF, and ANN
have the largest area under the precision-recall curve (AUC-PR).

Figure 5: Violin plot comparing the distribution of the actual and predicted
values of ODE evaluations, grid points and maximum residuum.

###  0.4.2 Multi-output Regression

As discussed in the proposed workflow, we train regression models to predict
the solver’s performance after estimating the solvability status of the
boundary value problems using a binary classification technique. During the
exploratory data analysis, we discovered that the solver-specific numerical
settings have little (problem 7) to no effect (problems 1, 3, and 4) on the
performance of linear problems. The effect is significant for non-linear
problems. In contrast to classification, regression techniques strongly depend
on the problem type to forecast its performance. Hence, the regression models
are trained based on the type of problem that is solved. Table 4 shows the
problem-specific performance of the regression models. The table shows that
the RF algorithm outperformed other algorithms in terms of evaluation metrics.
The comparison of the output distribution of the actual and predicted values
using the RF algorithm is shown in Figure 5.

Figure 6: Pareto Front plot showing the results of two-objective optimization
trials along the number of ODE evaluations and the maximum residuum.

###  0.4.3 Multi-objective Optimization

The ML models, which are chosen based on their performance, are used to
optimize the problem (1). The equation can be customized based on the
requirements to perform single-, two-, and three-objective optimizations.
Figure 6 shows the Pareto front for optimizing the number of ODE evaluations
and the maximum residuum. The 3D Pareto front for all three criteria
optimization is shown in Figure 7. In these plots, the blue dots represent the
regular optimization trials, and the red dots represent the best trials
chosen. Since predictions from the ML models are calculated in real-time, the
time required for optimization is reduced. To examine the reliability of the
results produced from our workflow, the best trials of the optimization were
tested against the actual output from the BVP solver.

Figure 7: Pareto Front plot showing the results of three-objective
optimization trials.

#### Fine-tuning library for the BVP solver

Based on the proposed workflow in the paper, we implemented a Python library,
bvpTune, to fine-tune the numerical settings of the BVP solver. Using the
library, users can retrieve the solvability of the desired problem based on
the solver settings. The library also facilitates the optimization of the
quality metrics based on the custom queries provided through the function
interfaces. The library can be installed as a Python package. The source code
and documentation are provided on
GitHub.111https://github.com/VictorVinySaajan/bvpTune The functionalities
provided by the library are described in Table 2.

##  0.5 Workflow Evaluation

In this section, we evaluate the scalability, stability, and reliability of
our workflow for fine-tuning the numerical settings of the BVP solvers. The
scalability of any workflow is important for real-world applications. The
scalability of our workflow mainly depends on the scalability of the ML model.
Hence, in our context, ”scalability” refers to the ML model’s capacity for
handling large volumes of data and performing numerous computations rapidly
and effectively. Since the training is only performed once, we primarily
concentrate on the prediction time. We evaluated the computing time of the
predictions for a significant data set in order to ascertain whether the
models scale well with the data. Figures 8 and 9 compare the time taken by the
classification and regression models for predicting 10,000 and 1,000,000
outputs, respectively. We did not include the SVM algorithms in the analysis
as they did not scale well with the data. The figures show that, with the
exception of the KNN, all the models scale reasonably well with the data for
both classification and regression. We observe that the LGBM algorithm is the
fastest for classification and the RF algorithm is the fastest for regression.
To examine the reliability of the workflow, the outputs from the bvpTune
library are tested against the actual outputs from the BVP solver. We found
that both observations were in agreement with each other.

As mentioned previously, the regression models are trained specifically to the
type of reference boundary value problems. We do this because the effect of
the numerical settings on the solver’s performance varies significantly as the
problem changes. To confirm this, we compared the output distributions of our
reference problems for the same set of inputs. We quantified the drift using
the population stability index ($PSI$) yurdakul2018statistical  that measures
population stability between two distribution samples. For samples divided
into $B$ bins, the $PSI$ is calculated as below:

Figure 8: Bar plot showing the scalability of the classification models.

$PSI=\sum_{i=1}^{B}(y_{i}-y_{ri})\ln(y_{i}/y_{ri})$ (2)

Table 5: Pairwise combination of PSI values for the reference boundary value problems defined in soetaert2010package . The values of $i^{\text{th}}$ row and $j^{\text{th}}$ column are calculated by keeping output distributions of Reference Problem Number(i) and Reference Problem Number(j) as initial populations respectively. Reference Problem Number |  7 |  19 |  20 |  22 |  23 |  24 |  33
---|---|---|---|---|---|---|---
7 |  0.0000 |  2.5384 |  4.5560 |  2.4123 |  1.5956 |  2.5570 |  1.6833
19 |  6.0619 |  0.0000 |  1.9525 |  4.2577 |  3.9291 |  6.2365 |  2.0148
20 |  5.2940 |  0.6268 |  0.0000 |  4.6000 |  2.8879 |  6.6869 |  1.4494
22 |  5.7805 |  3.8481 |  3.5217 |  0.0000 |  2.1835 |  1.3254 |  3.0429
23 |  3.3124 |  2.2905 |  2.5908 |  2.1386 |  0.0000 |  0.7574 |  3.6245
24 |  4.2097 |  3.6440 |  3.6107 |  1.5194 |  0.7938 |  0.0000 |  4.3400
33 |  5.7000 |  4.9902 |  5.8062 |  6.6005 |  7.5409 |  7.0110 |  0.0000

The $y_{i}$ is the proportion of the new output distribution that falls in the
$i^{th}$ bin and $y_{ri}$ is the proportion of the initial output distribution
that falls in the $i^{th}$ bin. According to the widely used rule of thumb
proposed in siddiqi2012credit , the drift between the populations is stable if
the $PSI$ value is less than 0.1. A value between 0.1 to 0.25 indicates a
small shift between the populations and a higher than 0.25 value represents
two substantially different populations. The PSI values observed from Table 5,
suggest that all the reference problems have significantly different behaviors
corresponding to the same inputs and cannot be generalized.

Figure 9: Bar plot showing the scalability of the regression models.

##  0.6 Conclusion

In this chapter, an ML-based optimization workflow for tuning the numerical
settings of boundary value problem solvers is proposed, and models based on
binary classification and multi-output regression are established. We
demonstrated the optimization of quality metrics exemplarily for a BVP solver
using these models. Furthermore, several machine learning algorithms for
parameter tuning are explored, and the scalability and reliability of these
models are discussed. Experimental results show that LGBM has good accuracy
and computational performance for classification, while Random Forests have
good accuracy and computational performance for regression. Additionally, the
stability of the models for unseen problem types is discussed and can be
continued for different boundary value problem solvers in future work.

##  0.7 Acknowledgement

This work was developed in Transport Processes department of Fraunhofer ITWM.
Furthermore, we acknowledge the support of the High Performance Center for
Simulation and Software-Based Innovation.

##  0.8 Authors Contribution

The authors confirm their contribution to the paper as follows: study
conception and design: Viny Saajan Victor, Manuel Ettmüller, Andre Schmeißer,
Heike Leitte and Simone Gramsch; data collection: Viny Saajan Victor and
Manuel Ettmüller; analysis and interpretation of results: Viny Saajan Victor;
draft manuscript preparation: Viny Saajan Victor. All authors reviewed the
results and approved the final version of the manuscript.

##  0.9 Conflict of Interest

All the named authors of this chapter declare that they have no conflicts of
interest, financial or otherwise.

## References

  * (1) Victor, Viny Saajan and Schmeißer, Andre and Leitte, Heike and Gramsch, Simone, “Visual Parameter Space Analysis for Optimizing the Quality of Industrial Nonwovens“ IEEE Computer Graphics and Applications, IEEE, vol. 42, no. 2, pp. 56–67, 2022.
  * (2) Ettmüller, Manuel and Arne, Walter and Marheineke, Nicole and Wegener, Raimund, “On flow-enhanced crystallization in fiber spinning: asymptotically justified boundary conditions for numerics of a stiff viscoelastic two-phase model“ Journal of Non-Newtonian Fluid Mechanics, Elsevier, vol. 296, pp. 104636, 2022.
  * (3) Arne, Walter and Marheineke, Nicole and Schnebele, Johannes and Wegener, Raimund, “Fluid-fiber-interactions in rotational spinning process of glass wool production“ Journal of Mathematics in Industry, SpringerOpen, vol. 1, no. 1, pp. 1–26, 2011.
  * (4) Kierzenka, Jacek and Shampine, Lawrence F, “A BVP solver based on residual control and the Matlab PSE“ ACM Transactions on Mathematical Software (TOMS), ACM New York, NY, USA, vol. 27, no. 3, pp. 299–316, 2001.
  * (5) Soetaert, Karline and Cash, Jeff and Mazzia, Francesca, “Package bvpSolve, solving testproblems“ 2010.
  * (6) Noble, William S, “What is a support vector machine¿‘ Nature biotechnology, Nature Publishing Group, vol. 24, no. 12, pp. 1565–1567, 2006.
  * (7) Ke, Guolin and Meng, Qi and Finley, Thomas and Wang, Taifeng and Chen, Wei and Ma, Weidong and Ye, Qiwei and Liu, Tie-Yan, “Lightgbm: A highly efficient gradient boosting decision tree“ Advances in neural information processing systems, vol. 30, 2007.
  * (8) Breiman, Leo, “Random forests“ Machine learning, Springer, vol. 45, no. 1, pp. 5–32, 2001.
  * (9) Fix, Evelyn and Hodges, Joseph Lawson, “Discriminatory analysis. Nonparametric discrimination: Consistency properties“ International Statistical Review/Revue Internationale de Statistique, JSTOR, vol. 57, no. 3, pp. 238–247, 1989.
  * (10) Cover, Thomas and Hart, Peter, “Nearest neighbor pattern classification“ IEEE transactions on information theory, vol. 13, no. 1, pp. 21–27, 1967.
  * (11) Akiba, Takuya and Sano, Shotaro and Yanase, Toshihiko and Ohta, Takeru and Koyama, Masanori, “Optuna: A Next-generation Hyperparameter Optimization Framework“ Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 2623–2631, 2019.
  * (12) Deb, Kalyanmoy and Pratap, Amrit and Agarwal, Sameer and Meyarivan, TAMT, “A fast and elitist multiobjective genetic algorithm: NSGA-II“ IEEE transactions on evolutionary computation, vol. 6, no. 2, pp. 182–197, 2002.
  * (13) Ozaki, Yoshihiko and Tanigaki, Yuki and Watanabe, Shuhei and Onishi, Masaki, “Multiobjective tree-structured parzen estimator for computationally expensive optimization problems“ Proceedings of the 2020 genetic and evolutionary computation conference, pp. 533–541, 2020.
  * (14) Ozaki, Yoshihiko and Tanigaki, Yuki and Watanabe, Shuhei and Nomura, Masahiro and Onishi, Masaki, “Multiobjective Tree-Structured Parzen Estimator“ Journal of Artificial Intelligence Research, vol. 73, pp. 1209–1250, 2022.
  * (15) Bergstra, James and Bengio, Yoshua, “Random search for hyper-parameter optimization.“ Journal of machine learning research, vol. 13, no. 2, 2012.
  * (16) Yurdakul, Bilal, “Statistical properties of population stability index“ Western Michigan University, 2018.
  * (17) Siddiqi, Naeem, “Credit risk scorecards: developing and implementing intelligent credit scoring“ John Wiley & Sons, vol. 3, 2012.