# Bifurcation analysis of a conceptual model for the Atlantic Meridional
Overturning Circulation

John Bailie and Bernd Krauskopf

###### Abstract

The Atlantic Meridional Overturning Circulation (AMOC) distributes heat and
salt into the Northern Hemisphere via a warm surface current toward the
subpolar North Atlantic, where water sinks and returns southwards as a deep
cold current. There is substantial evidence that the AMOC has slowed down over
the last century. We introduce a conceptual box model for the evolution of
salinity and temperature on the surface of the North Atlantic Ocean, subject
to the influx of meltwater from the Greenland ice sheets. Our model, which
extends a model due to Welander, describes the interaction between a surface
box and a deep-water box of constant temperature and salinity, which may be
convective or non-convective, depending on the density difference. Its two
main parameters $\mu$ and $\eta$ describe the influx of freshwater and the
threshold density between the two boxes, respectively.

We use tools from bifurcation theory to analyse two cases of the model: the
limiting case of instantaneous switching between convective or non-convective
interaction, where the system is piecewise-smooth (PWS), and the full smooth
model with more gradual switching. For the PWS model we perform a complete
bifurcation analysis by deriving analytical expressions for all bifurcations.
The resulting bifurcation diagram in the $(\mu,\eta)$-plane identifies all
regions of possible dynamics, which we show as phase portraits — both at
typical parameter points, as well as at the different transitions between
them. We also present the bifurcation diagram for the case of smooth switching
and show how it arises from that of the PWS case. In this way, we determine
exactly where one finds bistability and self-sustained oscillations of the
AMOC in both versions of the model. In particular, our results show that
oscillations between temperature and salinity on the surface North Atlantic
Ocean disappear completely when the transition between the convective and non-
convective regimes is too slow.

##  1 Introduction

The Atlantic Meridional Overturning Circulation (AMOC) is a large conveyor
belt of water that spans the entire Atlantic Ocean. Light surface currents
transport relatively warm and saline waters northward to high latitudes. Here,
the water becomes denser, leading to downward convection and mixing with the
deep ocean, and subsequent formation of deepwater masses. A deep current then
transports this water back to lower latitudes, where it upwells to the
surface, thus closing the circulation loop [1]. The strength of the AMOC is
governed by the interplay between two proposed upwelling mechanisms [2, 3].
The first perspective is that turbulent mixing across surfaces of equal
density results in the upwelling of deepwater to the surface ocean in low
latitudes [4, 5]. The second perspective suggests that strong circumpolar
winds induce upwelling in the South Atlantic Ocean [6]. Regardless of the
mechanism, the process of deepwater formation is crucial in determining the
shape and strength of the associated return current — making it a critical
factor for the stability of the AMOC.

This paper focuses on the deepwater formation sites in the North Atlantic.
Specifically, as illustrated in Figure 1(a), the convection of highly saline
water from the surface to the deep ocean in the Labrador and Nordic seas forms
the North Atlantic Deep Water (NADW); it has an associated return current
referred to as the NADW overturning cell. Several climate processes, such as
salt rejection and atmospheric cooling, facilitate this convection by
preconditioning the subpolar North Atlantic to have relatively high salinity
[7]. Furthermore, an advective process transports saline water to the North
Atlantic and, thus, stimulates the convection and formation of the NADW [8].
An inherent negative feedback loop forms: weaker convection results in a
smaller NADW and, consequently, a weaker overturning cell. This weaker cell
then advects less salt to the North Atlantic, further weakening the
convection.

Figure 1: Panel (a) shows a simplified sketch of the North Atlantic component
of the AMOC, inspired by [9]. The surface flow is displayed in red, and the
NADW overturning cell in blue. Deepwater formation sites are denoted
$\mathrm{L}$ and $\mathrm{N}$ in the Labrador sea and the Nordic seas,
respectively. Panel (b) shows the two-box model setup for the interaction
between surface water and cold deep water at the sites $\mathrm{L}$ and
$\mathrm{N}$.

Evidence from proxy and sea surface temperature measurements indicate that the
AMOC has weakened over the twentieth century [10]. There was also a
particularly abrupt change in overturning strength during the 1970s [11]
attributed to a large-scale influx of fresh water into the North Atlantic;
this is known as the Great Salinity Anomaly and is linked to Arctic sea-ice
export [12]. A weakened AMOC has significant consequences on the Earth’s
climate system since it leads to reduced northern heat transport, which lowers
the oceanic and atmospheric temperature in the Northern hemisphere [13] via a
weakening or even shutdown of the northern deepwater formation. Some
significant implications drawn from simulations are a widespread cooling in
Europe [10], the possible collapse of the North Atlantic plankton stocks [14],
and a rise in the sea level [15]. As a result of external environmental
factors, the AMOC is likely to weaken further, and a complete shutdown of the
deep water formation in the Labrador Sea is a possibility [10]. In particular,
meltwater from the melting Greenland ice sheets contributes to a large influx
of freshwater into the subpolar North Atlantic [16]. As freshwater is strictly
non-saline, it dilutes the ocean surface water by lowering its salinity, thus,
inhibiting the deep-water formation and, hence, the NADW overturning cell
strength.

Of particular interest in this context is the modelling of the underlying deep
water formation itself — with the aim of understanding the possible long-term
behaviour of the AMOC in response to freshwater influx. Climate models form a
hierarchy of complexity, and the choice of model depends on the nature of the
question that is being asked. We study here a conceptual model of low
complexity for investigating the AMOC in regard to deep water formation —
specifically, from the class of box models that consider only a few variables
in a relatively small number of interacting boxes, each representing a body of
water of concern. While they are not designed to be used for prediction, box
models are simple enough to be amenable to mathematical analysis, including
with tools from dynamical systems theory [17].

The stability of the AMOC was first investigated by Stommel with a two-box
model [18]; it considers the circulation between a subtropical box and a
subpolar box, where a capillary flow represents the advection of water between
the two boxes. Stommel’s model features three qualitatively different regimes.
In the first regime, the AMOC is driven by salinity differences between the
boxes, and surface currents move water toward the equator. Temperature
differences are the main driver in the second regime, and the surface currents
move water toward the poles. The final regime features bistability, where the
AMOC may tip to either of the described equilibrium states. Stommel laid the
foundation for several advective models, which add more boxes and physical
processes; see, for example, [19, 20, 21].

A two-box model presented by Welander [22] attempts to describe self-sustained
oscillations of temperature and salinity on the ocean surface in the presence
of external forcing. The boxes interact by exchanging heat and salt via a
mixing process. When the water in the surface box is sufficiently dense, the
mixing is convective (strong). When the water in the boxes has comparable
density, on the other hand, the mixing is non-convective (weak) and may happen
via several climate processes, such as double-diffusion [23]. In this setup,
an atmospheric basin with fixed properties interacts with the surface box,
which is modelled by Newton’s transfer law. The model by Welander is described
for two cases: when the transition between convective and non-convective
mixing is modelled as a continuous change and, alternatively, when it is
instantaneous and discontinuous. In both cases, self-sustained oscillations
are observed, which are characterised by a convective and a non-convective
phase. Welander’s model was re-examined by Leifeld [24] with the aim of
formalising the previous analysis by using a modern approach of piecewise-
smooth (PWS) dynamical systems. They undertook a preliminary stability
analysis and made a first comparison between the smooth and non-smooth models;
however, this work falls short of describing the full bifurcation picture and,
to the best of our knowledge, there is as yet no complete analysis of the
Welander model, nor any closely related models.

###  1.1 The adjusted Welander model

We take this as the starting point of our study of an _adjusted Welander
model_ that also considers the impact of a freshwater influx into the North
Atlantic ocean. Following on from work in [25], where the external forcing
enters in the form of Newton’s transfer law, we consider here a direct
freshwater flux that dilutes the salinity in a surface ocean box at the North
Atlantic, which is coupled to a box of deep water of constant lower
temperature and salinity. As is illustrated by the schematic in Figure 1(b),
the model takes the form of a planar system of ordinary differential equations
for temperature $T$ and salinity $S$ in the surface ocean box, which is given
by

$\displaystyle\frac{dT}{dt}=-\gamma(T-T_{a})-k_{\varepsilon}(\rho)(T-T_{0}),$
(1)
$\displaystyle\frac{dS}{dt}=\frac{F_{0}}{H}S_{0}-k_{\varepsilon}(\rho)(S-S_{0}).$

The atmosphere externally drives the surface ocean box to a thermal
equilibrium $T_{a}$ at rate $\gamma$, which is Newton’s transfer law. The
salinity, on the other hand, is directly forced by the freshwater flux $F_{0}$
at the rate $\frac{F_{0}}{H}S_{0}$, where $H$ is the depth of the surface
ocean box. Moreover, $T_{0}$ and $S_{0}$ are the (fixed) temperature and
salinity of the deep-ocean box that drive $T$ and $S$, respectively, as given
by the convective exchange function $k_{\varepsilon}$. This function
determines the coefficient for Newton’s transfer law and takes as its argument
the density $\rho$ of the surface ocean box given (in linear approximation) by

$\displaystyle\frac{\rho}{\rho_{0}}=1+\alpha_{S}(S-S_{0})-\alpha_{T}(T-T_{0}).$
(2)

Here, the constant $\rho_{0}$ is the density of the bottom box, and the
coefficients $\alpha_{S}$ and $\alpha_{T}$ are, respectively, the saline
expansion and thermal compression constants [17].

Figure 2: The switching functions
$\mathcal{H}_{\varepsilon}(\rho-\rho_{0}-g^{*})$ with $\varepsilon=0.1$ in
panel (a) and $\varepsilon=0$ in panel (b). Blue shading indicates mixing is
mainly non-convective, and red shading that mixing is mainly convective. In
the smooth transition region in panel (a), the colour changes from blue, via
white, to red.

The convective exchange function is a key ingredient in (1) and describes the
transition between the two regimes when the vertical mixing between the two
boxes is non-convective at rate $k_{1}>0$ and when it is convective at rate
$k_{2}>k_{1}$. This transition is modeled in its general form as

$\displaystyle
k_{\varepsilon}(\rho)=k_{1}+\mathcal{H}_{\varepsilon}(\rho-\rho_{0}-g^{*})(k_{2}-k_{1}),$
(3)

where $\mathcal{H}_{\varepsilon}$ is a suitable switching function from zero
to one, whose switching time depends on the switching-time parameter
$\varepsilon$. Hence, when the density difference $\rho-\rho_{0}$ is
(sufficiently) greater than the density threshold $g^{*}$, mixing between the
boxes is mainly convective; on the other hand, it is mainly non-convective
when $\rho-\rho_{0}$ is (sufficiently) smaller than $g^{*}$. Different
switching functions have been used in the literature [24, 20], including those
based on the $\arctan$ function. In this paper, we define
$\mathcal{H}_{\varepsilon}$ as

$\displaystyle\mathcal{H}_{\varepsilon}(u)=\frac{1}{2}\left(1+\tanh\left(\frac{u}{\varepsilon}\right)\right).$
(4)

Note that $\mathcal{H}_{\varepsilon}$ has a switching time of order
$\varepsilon>0$ and maximal rate of switching $\frac{1}{\varepsilon}$ given by
the derivative of $\mathcal{H}_{\varepsilon}$ at zero. Moreover, the limiting
case for $\varepsilon\xrightarrow{}0$ is an instantaneous switch, represented
by the Heaviside function $\mathcal{H}_{0}$. Figure 2 shows the resulting
convective exchange functions $k_{\varepsilon}(\rho)$ from (3) for
$\varepsilon=0.1$ and $\varepsilon=0$.

The adjusted Welander model in the form (1) has ten parameters, making a
direct analysis impractical. The first step in our analysis is to non-
dimensionalise the system by introducing rescaled temperature, salinity and
time

$x=\frac{T-T_{0}}{T_{a}-T_{0}},\quad
y=\frac{\alpha_{S}(S-S_{0})}{\alpha_{T}(T_{a}-T_{0})},\quad\tau=\gamma t,$ (5)

and parameters

$\kappa_{i}=\frac{k_{i}}{\gamma},\quad\mu=\frac{F_{0}S_{0}\alpha_{S}}{\gamma\alpha_{T}(T_{a}-T_{0})H}\quad\eta=\frac{g^{*}(\kappa_{2}-\kappa_{1})}{\gamma\alpha_{T}(T_{a}-T_{0})\rho_{0}}.$
(6)

This transforms (1) into

$\displaystyle\dot{x}=1-(1+\kappa_{1}+\mathcal{H}_{\varepsilon}(y-x-\eta)(\kappa_{2}-\kappa_{1}))x,$
(7)
$\displaystyle\dot{y}=\mu-(\kappa_{1}+\mathcal{H}_{\varepsilon}(y-x-\eta)(\kappa_{2}-\kappa_{1}))y,$

where the dot represents the derivative with respect to the rescaled time.

###  1.2 Outline of the work

The adjusted Welander model in the form (7) is our central object of study. We
first perform in Section 2 a (non-smooth) bifurcation analysis for the
limiting case of system (7) with the Heaviside switching function
$\mathcal{H}_{0}$. Specifically, we determine and catalogue all of the
possible dynamics, by presenting analytical expressions for all codimension-
one and codimension-two bifurcations; the corresponding proofs and derivations
can be found in Appendix A. The rescaled freshwater flux $\mu$ and density
threshold $\eta$ are the bifurcation parameters, and we show the complete
bifurcation diagram in the $(\mu,\eta)$-plane for a reasonable choice of the
vertical mixing coefficients $0<\kappa_{1}<\kappa_{2}$. Moreover, we present
representative phase portraits in the $(x,y)$-plane for all open regions, of
which there are eight, for the different types of transitions of codimension
one between them, as well as at the five codimension-two points that organise
the bifurcation diagram. In particular, we identify the parameter regime where
the system exhibits bistability between states, where deep-water convection is
either substantial or shut down, which is characteristic behaviour of several
models of different complexity [8]. Moreover, we determine the parameter
regime with self-sustained relaxation-type oscillations that have been
observed in [25] and also in Welander’s original work [20]. In addition to
this earlier work, we determine this region analytically and clarify the
nature of the different possible transitions to/from this oscillatory regime.
Our results also show that the bifurcation diagram in the $(\mu,\eta)$-plane
is toplogically the same for any fixed $0<\kappa_{1}<\kappa_{2}$.

Section 3 is then concerned with the smooth case of system (7) with
$\mathcal{H}_{\varepsilon}$ for $\varepsilon>0$. Here we first present the
bifurcation diagram in the $(\mu,\eta)$-plane for $\varepsilon=0.1$ (for the
same choice of $0<\kappa_{1}<\kappa_{2}$). This requires computing the
relevant bifurcation curves by making use of established bifurcation theory
[26] in conjunction with the continuation software package AUTO-07p [27]. We
focus here on the main parameter regimes, especially those that feature
bistability and self-sustained oscillations, for which we show representative
phase portraits. We then present a partial bifurcation analysis in
$(\mu,\eta,\varepsilon)$-space that clarifies the convergence of the
bifurcation diagram in the $(\mu,\eta)$-plane as $\varepsilon$ approaches $0$.
Moreover, we show that there is a codimension-three bifurcation at a quite low
value of the switching-time parameter $\varepsilon$, at which the region with
self-sustained oscillations completely disappears from the $(\mu,\eta)$-plane.
In other words, the switching between the regimes with strong convective
mixing and with weak non-convective mixing needs to be sufficiently fast for
relaxation-type oscillations to occur in the adjusted Welander model (7). In
particular, this shows the relevance of the non-smooth limiting system with
$\mathcal{H}_{0}$ for explaining this oscillatory behaviour.

In the final Section 4 we summarise our findings, briefly discuss their
significance for the dynamics of AMOC, and point out some direction for future
work.

##  2 Bifurcation analysis of the PWS model for $\varepsilon=0$

In the limiting case of an instantaneous transition with the transition
function $\mathcal{H}_{0}$ from (4), system (7) reduces to the piecewise-
smooth linear Filippov system

$\displaystyle\binom{\dot{x}}{\dot{y}}=\left\\{\begin{matrix}&f_{1}(x,y),\ \ \
y<x+\eta\\\ &f_{2}(x,y),\ \ \ y>x+\eta\end{matrix}\right.$ (8)

with

$\displaystyle f_{i}(x,y)=\binom{1-(1+\kappa_{i})x}{\mu-\kappa_{i}y}.$ (9)

The switching manifold

$\displaystyle\Sigma=\\{(x,y)\in\mathbb{R}^{2},\ \ y=x+\eta\\}$ (10)

is a straight line that partitions the phase space of (8) into the open
regions

$\displaystyle R_{1}=\\{(x,y)\in\mathbb{R}^{2},\ \ y<x+\eta\\},$ (11)
$\displaystyle R_{2}=\\{(x,y)\in\mathbb{R}^{2},\ \ y>x+\eta\\},$ (12)

where $f_{1}$ and $f_{2}$ apply, respectively.

We now perform a bifurcation analysis of the piecewise-smooth AMOC model (8).
To this end, we use tools from the bifurcation theory for this class of non-
smooth systems from the relevant literature [28, 29, 30], which we largely
follow also in terms of notation and where more details can be found. More
specifically, we determine analytic expressions for all (non-smooth)
bifurcations, which is possible because of the simple expression for the
switching manifold, and the fact that $f_{1}$ and $f_{2}$ are linear. We
present these results in the form of propositions, whose proofs can be found
in Appendix A. The associated curves of codimension-one bifurcations divide
the $(\mu,\eta)$-plane into eight open regions, denoted $\mathrm{I-VIII}$. We
also present the corresponding phase portraits in the $(\mu,\eta)$-plane, as
well as those at the different types of bifurcations. The vertical mixing
coefficients are fixed here to $\kappa_{1}=0.1$ and $\kappa_{2}=1.0$. This
choice is suitable for our purposes and in the realistic range [25], yet
slightly different from the values found in the literature [24, 22]. Moreover,
as can be seen from the expressions in Section 2.3, the bifurcation diagram in
the $(\mu,\eta)$-plane is qualitatively the same for any
$0<\kappa_{1}<\kappa_{2}$.

###  2.1 Sliding properties and pseudo-equilibria

We start by introducing some relevant notions from the theory of PWS systems.
An equilibrium $p_{i}$ of the vector field $f_{i}$ that lies in region $R_{i}$
is an equilibrium of the overall system and called _admissible_. An important
part of the bifurcation theory of planar Filippov systems is the interaction
of equilibria and other invariant objects of $f_{1}$ and $f_{2}$ with the
switching manifold $\Sigma$ [28]. First of all, orbits may cross the switching
manifold at the _crossing segment_ $\Sigma_{c}\subset\Sigma$, along which the
vector fields $f_{1}$ and $f_{2}$ are both transverse and have the same sign.
The set of points where $f_{1}$ and $f_{2}$ are transversal but have opposite
signs is the _sliding segment_ $\Sigma_{s}\subset\Sigma$, which we also refer
to as $\Sigma_{s}^{a}$ when it is attracting and as $\Sigma_{s}^{r}$ when it
is repelling. These different segments of the switching manifold are bounded
by _tangency points_ $F_{1}$ and $F_{2}$, where either $f_{1}$ or $f_{2}$ is
tangent to $\Sigma$, respectively. Generically, such a tangency of $f_{i}$ is
quadratic and isolated, and it is called visible if nearby parabolic orbits
lie in $R_{i}$, and invisible otherwise. For system (8) we have the following.

######  Proposition 1 (Tangency points and sliding segments).

System (8) has a single sliding segment $\Sigma_{s}$ that is delimited by two
tangency points $F_{1}$ and $F_{2}$ at

$\displaystyle
F_{i}=\binom{1-\mu+\eta\kappa_{i}}{1-\mu+(1+\kappa_{i})\eta}\in\Sigma,$ (13)

which are quadratic when

$\displaystyle\mu+(\mu-\eta-1)\kappa_{i}-\eta\kappa_{i}^{2}\neq 0.$ (14)

The (quadratic) tangency point $F_{1}$ is visible for

$\displaystyle\mu+(\mu-\eta-1)\kappa_{i}-\eta\kappa_{i}^{2}<0,$ (15)

and the (quadratic) tangency point $F_{2}$ is visible for

$\displaystyle\mu+(\mu-\eta-1)\kappa_{i}-\eta\kappa_{i}^{2}>0.$ (16)

Otherwise, the (quadratic) tangency at $F_{i}$ is invisible. For $\eta\neq 0$,
system (8) has a sliding segment $\Sigma_{s}$. When $\eta>0$ the sliding
segment is attracting, denoted $\Sigma_{s}^{a}$ and given by

$\displaystyle\Sigma_{s}^{a}=\\{s\in\Sigma,F_{1}<s<F_{2}\\},$ (17)

and when $\eta<0$ it is repelling, denoted $\Sigma_{s}^{r}$ and given by

$\displaystyle\Sigma_{s}^{r}=\\{s\in\Sigma,F_{2}<s<F_{1}\\}.$ (18)

Here, in a slight abuse of notation, we mean the ordering on the line
$\Sigma_{s}$, as given by the $x$-component.

A crucial ingredient of the theory is the extension of the flow to the sliding
segment $\Sigma_{s}$ by defining the sliding vector field $f_{s}$. This is
achieved with Filippov’s convex method by forming a weighted sum of the
adjoining vector fields $f_{1}$ and $f_{2}$ such that $f_{s}$ is in the
direction of (the tangent to) $\Sigma_{s}$ [28, 29, 30]. With this definition,
a PWS orbit is the union of orbit segments induced by the vector fields
$f_{1}$ on $R_{1}$, $f_{2}$ on $R_{2}$, and $f_{s}$ on $\Sigma_{s}$. Moreover,
every point of the phase plane lies on a unique PWS orbit of the planar
Filippov system; see [30] for details. Orbits that remain in $R_{1}\cup
R_{2}\cup\Sigma_{c}$ are called _regular_ , and orbits with segments on
$\Sigma_{s}$ are called _sliding orbits_. Here we use a common convention that
sliding orbits continue into $R_{1}$ or $R_{2}$ when the end of the sliding
segment $\Sigma_{s}$ is reached (in forward or backward time, respectively, by
following the trajectory from the respective tangency point) [30, 31].
However, the end of $\Sigma_{s}$ may not be reached because the sliding vector
field may have equilibria, called _pseudo-equilibria_ , which are referred to
as _admissible_ when they lie $\Sigma_{s}$.

The properties of all equilibria of system (8) can be stated as follows.

######  Proposition 2 (Equilibria, sliding vector field and pseudo-
equilibria).

System (8) has the following equilibria and pseudo-equilibria for
$0<\kappa_{1}<\kappa_{2}$.

  1. 1.

The vector field $f_{i}$ has the stable nodal equilibrium

$\displaystyle
p_{i}=\left(\frac{1}{1+\kappa_{i}},\frac{\mu}{\kappa_{i}}\right).$ (19)

The equilibrium $p_{1}$ is admissible when

$\displaystyle\eta>\frac{\mu}{\kappa_{1}}-\frac{1}{1+\kappa_{1}},$ (20)

and $p_{2}$ is admissible when

$\displaystyle\eta<\frac{\mu}{\kappa_{2}}-\frac{1}{1+\kappa_{2}}.$ (21)

The admissible equilibrium $p_{i}$ has a strong stable manifold
$W^{ss}(p_{i})$ defined by the piecewise-smooth orbit along the linear strong
stable direction $W^{ss}_{loc}=\textnormal{span}\binom{1}{0}$.

  2. 2.

The _sliding vector field_ defined on the sliding segment $\Sigma_{s}$ is
given by

$\displaystyle f_{s}(x)$
$\displaystyle=\frac{1}{\eta}\left(\mu+(\mu-\kappa_{2}\eta-1)x+x^{2}\right)\binom{1}{1},$
(22)

where $\Sigma$ is parametrised by $x$.

  3. 3.

There are two pseudo-equilibria, that is, equilibria of $f_{s}$, given by

$\displaystyle
q^{\pm}=\frac{1}{2}\left(1-\mu-\eta\pm\sqrt{(\eta+\mu+1)^{2}-4\mu}\right)\binom{1}{1}+\binom{0}{\eta}.$
(23)

When $\eta>0$, the pseudo-equilibrium $q^{-}$ is asymptotically unstable and
$q^{+}$ is asymptotically stable on $\Sigma_{s}$. On the other hand, when
$\eta<0$, the pseudo-equilibrium $q^{+}$ is asymptotically unstable and
$q^{-}$ is asymptotically stable on $\Sigma_{s}$. The admissibility of these
pseudo-equilibria is presented and described in Section 2.2.

Global invariant manifolds of admisible equilibria are defined in complete
analogy to those of smooth systems, but with regard to the piecewise-smooth
flow $\varphi^{t}$ constructed in [30]. Each admissible equilibrium $p_{i}\in
R_{i}$ of system (8) is attracting with real eigenvalues and, hence, has a
strong stable manifold $W^{ss}(p_{i})$ consisting of the two orbits that
approach $p_{i}$ tangent to the strong eigenspace. Since system (8) is
piecewise linear, $W^{ss}(p_{i})$ is actually a straight line locally near
$p_{i}$; however, this is not the case globally since the strong stable
manifold typically crosses the switching manifold $\Sigma$.

We also consider here global invariant manifolds of admissible pseudo-
equilibria, which we define as follows. If $q\in\Sigma_{s}$ is a saddle
pseudo-equilibrium then its stable manifold $W^{s}(p)$ or unstable manifold
$W^{u}(p)$ is the union of the two _arriving or departing orbits_ in $R_{1}$
and $R_{2}$, consisting of points that reach $q$ under the piecewise-smooth
flow $\varphi^{t}$ in finite forward or backward time, respectively. The
saddle pseudo-equilibrium $q$ then also has associated _generalised (un)stable
manifolds_ $W^{u}_{g}(q)$ or $W^{s}_{g}(q)$. These generalised manifolds
consist of segments on $\Sigma_{s}$ of points that converge to $q$ under the
sliding flow (in backward and forward time, respectively), together with their
globalisation under $\varphi^{t}$, which generally consists of departing and
arriving orbits to tangency points that bound $\Sigma_{s}$. When an admissible
pseudo-equilibrium $q$ is a nodal attractor, its arriving orbits form the
strong stable manifold $W^{ss}(p)$; similarly, a nodal repellor
$q\in\Sigma_{s}$ has the strong unstable manifold $W^{uu}(p)$ consisting of
its pair of departing orbits.

###  2.2 Bifurcation diagram and structurally stable phase portraits

The bifurcation diagram of system (8) consists of curves of (piecewise-smooth)
bifurcations that divide the $(\mu,\eta)$-plane into eight open regions
$\mathrm{I}$ to $\mathrm{VIII}$, which are equivalence classes of topological
equivalence where the phase portraits are structurally stable. This
classification is based on the following common notion [30, 31]: two planar
Filippov systems $f$ and $\tilde{f}$ with switching manifolds $\Sigma$ and
$\widetilde{\Sigma}$, respectively, are _topologically equivalent_ if there
exists an orientation preserving homeomorphism
$h:\mathbb{R}^{2}\xrightarrow{}\mathbb{R}^{2}$ that maps $\Sigma$ to
$\widetilde{\Sigma}$ and orbits of $f$ to orbits of $\tilde{f}$. Note that
this definition is a direct and natural extension of that for smooth dynamical
systems. In particular, a bifurcation of a planar Filippov system concerns a
topological change, and its codimension is given (colloquially speaking) by
the number of parameters one needs to find it generically at an isolated
point. More information and formal definitions can be found as part of the
broad classification in [31] of discontinuity-induced bifurcations in planar
Filippov systems.

Figure 3: Two-parameter bifurcation diagram of system (8) in the
$(\mu,\eta)$-plane with $\kappa_{1}=0.1$ and $\kappa_{2}=1.0$. The curves of
boundary equilibrium bifurcation $\mathrm{BE}_{1}$ and $\mathrm{BE}_{2}$,
fold-fold bifurcation $\mathrm{FF}$, and pseudo-saddle-node bifurcation
$\mathrm{PS}$ from Proposition 3 bound regions $\mathrm{I}$ to $\mathrm{VIII}$
with structurally stable phase portraits shown in Figures 4–6. Grey shading
indicates the existence of a (crossing) periodic orbit, and blue shading
bistability between equilibria. These curves intersect at the codimension-two
points $\mathrm{FB_{1}}$, $\mathrm{FB_{1}}$, $\mathrm{BB}$, $\mathrm{GB_{1}}$,
and $\mathrm{GB_{2}}$ from Proposition 4, which generates segments of
different bifurcation types shown in Figures 7–11. Specifically, the curves
$\mathrm{BE}_{1}$ and $\mathrm{BE}_{2}$ consist of segments
$\mathrm{BE}_{1}^{P}$, $\mathrm{BE}_{2}^{P}$, $\widehat{\mathrm{BE}}_{1}^{P}$,
$\widehat{\mathrm{BE}}_{2}^{P}$ and $\widetilde{\mathrm{BE}}_{1}^{P}$ of
persistence boundary equilibrium bifurcation, and $\mathrm{BE}_{1}^{F}$,
$\mathrm{BE}_{2}^{F}$ and $\widehat{\mathrm{BE}}_{2}^{F}$ of non-smooth fold
boundary equilibrium bifurcation. The curve $\mathrm{FF}$ consists of segments
$\mathrm{FF_{1}}$ and $\mathrm{FF_{2}}$ of fold-fold bifurcation and
$\mathrm{FU}$ of fused-focus bifurcation.

Figure 3 shows the bifurcation diagram of system (8) in the
$(\mu,\eta)$-plane, for the fixed values $\kappa_{1}=0.1$ and $\kappa_{2}=1.0$
of the vertical mixing rates, with the regions $\mathrm{I}$ to
$\mathrm{VIII}$. Their boundaries are formed by bifurcation curves
$\mathrm{BE}_{1}$ and $\mathrm{BE_{2}}$ of boundary equilibrium bifurcation,
$\mathrm{FF}$ of fold-fold bifurcation and $\mathrm{PS}$ of pseudo-saddle-node
bifurcation that are formally presented and determined in Proposition 3. More
precisely, these curves cross or meet at codimension-two bifurcation points
$\mathrm{FB_{1}}$, $\mathrm{FB_{1}}$, $\mathrm{BB}$, $\mathrm{GB_{1}}$, and
$\mathrm{GB_{2}}$. As is spelled out in Proposition 4, these points divide the
curves of codimension-one bifurcations into the segments of different
bifurcation types that are shown and labeled in Figure 3.

We first present and discuss the structurally stable phase portraits of system
(8) in regions $\mathrm{I}$ to $\mathrm{VIII}$; the different bifurcations
between them are analysed and illustrated in subsequent sections. The eight
cases of phase portraits are shown in Figures 4–6 in a suitable part of the
$(x,y)$-plane. In every phase portrait, the switching manifold $\Sigma$
appears as a straight grey line that partitions phase space into the open
regions ${R_{1}}$ and ${R_{2}}$. Admissible equilibria of system (8) are shown
in black, and non-admissible equilibria in grey. Non-admissible equilibria
outside the frame of interest (far away from the switching manifold) are not
shown. There exists a sliding segment in each region: attracting sliding
segments $\Sigma_{s}^{a}$ are coloured blue, and repelling sliding segments
$\Sigma_{s}^{r}$ are coloured orange. In either case, the sliding segment is
bounded by the quadratic tangency points $F_{1}$ and $F_{2}$, which are
coloured cyan when visible and grey when invisible. Admissible pseudo-
equilibria $q^{-}$ and $q^{+}$ are coloured by their stability: stable pseudo-
equilibria are green, and unstable pseudo-equilibria are red. Admissible
equilibria and pseudo-equilibria may have (strong) invariant manifolds that
are coloured blue when stable and red when unstable. Some representative
trajectories are shown in black, and they were obtained numerically with an
integrator based on event-detection, as described in [32].

Figure 4: Representative phase portraits in regions $\mathrm{I}$ to
$\mathrm{III}$ along the horizontal slice $\eta=0.35$ of the
$(\mu,\eta)$-plane, each with an attracting sliding segment $\Sigma_{s}^{a}$
bounded by quadratic tangency points $F_{1}$ and $F_{2}$. Panel (a) for
$\mu=0.0225$ shows the admissible equilibrium $p_{1}$ with its strong stable
manifold $W^{ss}(p_{1})$, as well as the non-admissible equilibrium $p_{2}\in
R_{1}$. Panel (b1) for $\mu=0.385$ and magnification (b2) near the sliding
segment show $p_{2}\in R_{1}$ and the attracting pseudo-node $q^{+}$ with
strong stable manifold $W^{ss}(q^{+})$. Panel (c) for $\mu=0.975$ shows the
admissible equilibrium $p_{2}\in R_{2}$ with $W^{ss}(p_{2})$.

Figure 4 presents phase portraits of system (8) in regions $\mathrm{I}$ to
$\mathrm{III}$, which all feature and attracting sliding segment
$\Sigma_{s}^{a}$. In the phase portrait in region $\mathrm{I}$, shown in panel
(a), $\Sigma_{s}^{a}$ is bounded by a visible quadratic tangency point $F_{1}$
on the left and an invisible quadratic tangency point $F_{2}$ on the right.
Neither of the pseudo-equilibria $q^{-}$ and $q^{+}$ are on $\Sigma_{s}^{a}$
and, hence, they are non-admissible (and not shown). The equilibrium $p_{2}$
lies in region $R_{1}$ and is non-admissible, while $p_{1}\in R_{1}$ is
admissible and a global attractor. Orbits in $R_{1}$ and $R_{2}$ are either
regular and converge to $p_{1}$ or hit the attracting sliding segment
$\Sigma_{s}^{a}$, along which sliding orbits approach $F_{1}$ and then depart
into $R_{1}$ to converge to $p_{1}$. Note, that the strong stable manifold
$W^{ss}(p_{1})$ of $p_{1}$ is composed of a horizontal component in ${R_{1}}$
and the corresponding arriving orbit in $R_{2}$. Crossing the segment
$\mathrm{BE}_{1}^{P}$ of boundary equilibrium bifurcation results in $p_{1}$
becoming non-admissible by moving into $R_{2}$ through $F_{1}$; at the same
time, a pseudo-equilibrium $q^{+}\in\Sigma_{s}^{a}$ emerges from $F_{1}$,
where the tangency is now invisible. The resulting phase portrait in region
$\mathrm{II}$ is shown in Figure 4(b1) with a magnification near the sliding
segment in panel (b2). The pseudo-equilibrium $q^{+}$ is a global attractor:
all orbits in $R_{1}$ and $R_{2}$ hit the sliding segment $\Sigma_{s}^{a}$,
along which the sliding orbits converge to $q^{+}$. Moreover, $q^{+}$ has the
strong stable manifold $W^{ss}(q^{+})$, consisting of the two arriving orbits
to $q^{+}$ from within $R_{1}$ and $R_{2}$, respectively; see panel (b2). When
the segment $\mathrm{BE}_{2}^{P}$ is crossed there is again a boundary
equilibrium bifurcation, but now of $p_{2}$ at $F_{2}$: as Figure 4(c) shows,
in region $\mathrm{III}$ the pseudo-equilibrium $q^{+}$ moved off
$\Sigma_{s}^{a}$ through $F_{2}$, and $p_{2}$ with strong stable manifold
$W^{ss}(p_{2})$ is now admissible and the global attractor.

Figure 5: Representative phase portraits in regions $\mathrm{IV}$ to
$\mathrm{VI}$ along a horizontal slice $\eta=-0.2$ of the $(\mu,\eta)$-plane,
each with a repelling sliding segment $\Sigma_{s}^{r}$ bounded by quadratic
tangency points $F_{1}$ and $F_{2}$. Panel (a) for $\mu=0.0225$ is similar to
Figure 4(a), but now the sliding segment is repelling. Panel (b1) for
$\mu=0.25$ and the magnification near the sliding segment (b2) features a
repelling pseudo-node $q^{+}$ with a strong unstable manifold $W^{uu}(q^{+})$
and a (crossing) periodic orbit $\Gamma$ that encircles $\Sigma_{s}^{r}$.
Panel (c) for $\mu=0.525$ is similar to Figure 4(c), but now the sliding
segment is repelling.

Phase portraits in regions $\mathrm{IV}$ to $\mathrm{VI}$ are presented Figure
5; as was the case for regions $\mathrm{I}$ to $\mathrm{III}$, this also
concerns the the transition from $p_{1}$ to $p_{2}$ being the global
attractor, with the difference that there is now a repelling sliding segment
$\Sigma^{r}_{s}$. In the phase portrait in region $\mathrm{IV}$, shown in
panel (a), $\Sigma^{r}_{s}$ is bounded by an invisible quadratic tangency
point $F_{2}$ on the left and by a visible quadratic tangency point $F_{1}$ on
the right; the pseudo-equilibria $q^{-}$ and $q^{+}$ are on $\Sigma_{c}$ and
non-admissible (and not shown). The only admissible equilibrium is $p_{1}\in
R_{1}$, and it is a global attractor. Sliding orbits on $\Sigma_{s}^{r}$
approach $F_{2}$, where they depart into $R_{1}$ and converge to $p_{1}$;
however, in contrast to region $\mathrm{III}$, no forward orbits hit
$\Sigma_{s}^{r}$ as this sliding segment is repelling. When crossing segment
$\widehat{\mathrm{BE}}_{1}^{P}$, we find again a (persistence) boundary
equilibrium bifurcation where $p_{1}$ moves through $F_{1}$ and becomes non-
admissible. As Figure 5(b1) and the magnification in panel (b2) show, in
region $\mathrm{V}$ this results again in the pseudo-equilibrium $q^{+}$ being
admissible. However, $q^{+}\in\Sigma_{s}^{r}$ is now a repelling node with
strong unstable manifold $W^{uu}(q^{+})$, consisting of the departing orbits
from $q^{+}$ in $R_{1}$ and $R_{2}$, respectively. Importantly, in region
$\mathrm{V}$ there is a stable (crossing) periodic orbit $\Gamma$, which is
composed of orbit segments of $f_{1}$ in $R_{1}$ and $f_{2}$ in $R_{2}$ that
join on the crossing segment $\Sigma_{c}$. All points except $q^{+}$ converge
to this periodic orbit; in particular, $W^{uu}(q^{+})$ accumulates on
$\Gamma$, while initial conditions on $\Sigma_{s}^{r}\setminus\\{q^{+}\\}$
move to an end point $F_{2}$ or $F_{1}$ of $\Sigma_{s}^{r}$, where they depart
into $R_{1}$ or $R_{2}$, respectively, to converge to $\Gamma$; see panel
(b2). Crossing segment $\widehat{\mathrm{BE}}_{2}^{P}$ concerns a second
(persistence) boundary equilibrium bifurcation, but now of $p_{2}$ at the
tangent point $F_{2}$. As a result, the now admissible equilibrium $p_{2}$ is
indeed the global attractor in region $\mathrm{VI}$, as is shown in Figure
5(c).

Figure 6: Representative phase portraits in regions $\mathrm{VII}$ and
$\mathrm{VIII}$, both with a repelling sliding segment $\Sigma_{s}^{r}$
bounded by quadratic tangency points $F_{1}$ and $F_{2}$. Panel (a1) for
$(\mu,\eta)=(0.0987,-0.463)$ and the magnification near the sliding segment
(a2) show $p_{2}$ with $W^{ss}(p_{2})$, pseudo-saddle-equilibrium $q^{-}$ with
a generalised stable manifold $W^{s}_{g}(q^{-})$ and unstable manifold
$W^{u}(q^{-})$, and the repelling pseudo-node $q^{+}$ with strong unstable
manifold $W^{uu}(q^{+})$. Panel (b) for $(\mu,\eta)=(-0.115,-0.95)$ shows the
simultaneously admissible equilibria $p_{1}$ with $W^{ss}(p_{1})$ and $p_{2}$
with $W^{ss}(p_{2})$, and the pseudo-saddle-equilibrium $q^{-}$ with
$W^{s}_{g}(q^{-})$ and $W^{u}(q^{-})$.

Phase portraits for regions $\mathrm{VII}$ and $\mathrm{VIII}$ are presented
in Figure 6; they both still feature the repelling sliding segment
$\Sigma^{r}_{s}$, bounded by the quadratic tangency points $F_{2}$ on the left
and $F_{1}$ on the right. In region $\mathrm{VII}$, as in panel (a1) with a
magnification in panel (a2), we find the repelling pseudo-equilibrium
$q^{+}\in\Sigma_{s}^{r}$ as in region $\mathrm{V}$. However, due to the
transition through the bounding segment $\widehat{\mathrm{BE}}_{2}^{F}$ of
(non-smooth fold) boundary equilibrium bifurcation, the equilibrium $p_{2}$ is
now in open region $R_{2}$ and admissible, and the second pseudo-equilibrium
$q^{-}$ now also lies on $\Sigma_{s}^{r}$. The point $p_{2}$ attracts all
points, apart from those on the generalised stable manifold $W^{s}_{g}(q^{-})$
of $q^{-}$, which is composed of sliding orbits on $\Sigma^{r}_{s}$
approaching $q^{-}$ and the arriving orbit to $F_{2}$. The unstable manifold
$W^{u}(q^{-})$ of $q^{-}$ and strong unstable manifold $W^{uu}(q^{+})$ both
converge to the attractor $p_{2}$. Note that $W^{ss}(p_{2})$ is composed of a
horizontal component in $R_{1}$, the corresponding sliding orbit in
$\Sigma^{r}_{s}$ and the arriving orbit to $F_{1}$. When crossing segment
$\widetilde{\mathrm{BE}}_{1}^{P}$ into region $\mathrm{VIII}$, there is a
(persistence) boundary equilibrium bifurcation, at which $p_{1}$ becomes
admissible and the pseudo-equilibrium $q^{+}$ becomes non-admissible by moving
through $F_{1}$ onto the crossing segment $\Sigma_{c}$. As the phase portrait
in Figure 6(b) shows, both $p_{1}\in R_{1}$ and equilibrium $p_{2}\in R_{2}$
are now attractors in region $\mathrm{VIII}$; hence, this is the region of
bistability. The generalised stable manifold $W_{g}^{s}(q^{-})$ of the saddle
pseudo-equilibrium $q^{-}\in\Sigma^{r}_{s}$ is now composed of
$\Sigma^{r}_{s}$ and the arriving orbits to both $F_{1}$ and $F_{2}$, and it
forms the boundary between the basins of attraction of the attractors $p_{1}$
and $p_{2}$. Indeed, the lower branch of $W^{u}(q^{-})$ converges to $p_{1}$,
and its upper branch to $p_{2}$.

###  2.3 Codimension-one and codimension-two bifurcations

We now present analytical expressions for all (non-smooth) bifurcations of
system (8) of codimension one and two in Propositions 3 and 4, respectively.
The respective proofs can be found in Appendix A.

######  Proposition 3 (Codimension-one bifurcations).

System (8) has the following codimension-one bifurcations, further information
on which can be found in [31].

  1. 1.

_Boundary equilibrium bifurcations_ of $p_{1}$ and $p_{2}$ occur,
respectively, along the straight lines

$\displaystyle\mathrm{BE}_{1}:\ \ (\mu,\eta)=\left(\mu,\
\frac{\mu}{\kappa_{1}}-\frac{1}{1+\kappa_{1}}\right),$
$\displaystyle\mathrm{BE}_{2}:\ \ (\mu,\eta)=\left(\mu,\
\frac{\mu}{\kappa_{2}}-\frac{1}{1+\kappa_{2}}\right).$

At $\mathrm{BE}_{i}$ the equilibrium $p_{i}$ collides with the tangency point
$F_{i}$, changing its visibility. The tangency point $F_{1}$ is visible in
regions $\mathrm{I,IV}$ and $\mathrm{VIII}$. Similarly, the tangency point
$F_{2}$ is visible in regions $\mathrm{III,VI,VII}$ and $\mathrm{VIII}$. The
pseudo-equilibrium $q^{+}$ is admissible in regions $\mathrm{II}$,
$\mathrm{V}$ and $\mathrm{VII}$. Similarly, the pseudo-equilibrium $q^{-}$ is
admissible in regions $\mathrm{VII}$ and $\mathrm{VIII}$.

  2. 2.

_Fold-fold bifurcations_ occur along the horizontal line

$\displaystyle\mathrm{FF}:\ \ (\mu,\eta)=(\mu,\ 0).$

At $\mathrm{FF}$ the tangency points $F_{1}$ and $F_{2}$ coincide at a
singular tangency point $F^{*}$ and switch places on the sliding segment
boundary, resulting in the sliding segments changing between being attracting
and repelling [31]; see also Proposition 1 for a description of the tangency
points.

  3. 3.

A _pseudo-saddle-node bifurcation_ occurs along the curve segment

$\displaystyle\mathrm{PS}:\ \ (\mu,\eta)=(\mu,\ -(\mu+1)+2\sqrt{\mu}),\ \ \
\frac{\kappa_{1}^{2}}{(\kappa_{1}+1)^{2}}<\mu<\frac{\kappa_{2}^{2}}{(\kappa_{2}^{2}+1)^{2}}.$

Along $\mathrm{PS}$ the pseudo-equilibria $q^{-}$ and $q^{+}$ form a saddle-
node at

$\displaystyle q^{*}=(1-\sqrt{\mu})\binom{1}{1}+\binom{0}{\eta}$

on the repelling sliding segment $\Sigma_{s}^{r}$.

The curves $\mathrm{BE}_{1}$, $\mathrm{BE}_{2}$, $\mathrm{FF}$ and
$\mathrm{PS}$ from Proposition 3 intersect or meet at codimension-two
bifurcation points. These points divide $\mathrm{BE}_{1}$, $\mathrm{BE}_{2}$,
$\mathrm{FF}$ into the segments shown in Figure 3, along which the respective
codimension-one bifurcation manifests itself in a topologically different way,
as follows.

######  Proposition 4 (Codimension-two bifurcations).

System (8) has the following codimension-two bifurcations for
$0<\kappa_{1}<\kappa_{2}$.

  1. 1.

_Fold-boundary equilibrium bifurcations_

$\displaystyle\mathrm{FB_{1}}:\ \
(\mu,\eta)=\left(\frac{\kappa_{1}}{1+\kappa_{1}},\ 0\right),$ (24)
$\displaystyle\mathrm{FB_{2}}:\ \
(\mu,\eta)=\left(\frac{\kappa_{2}}{1+\kappa_{2}},\ 0\right),$ (25)

occur at the intersection point of the curve $\mathrm{FF}$ with the curves
$\mathrm{BE_{1}}$ and $\mathrm{BE_{2}}$, respectively. At the point
$\mathrm{FB}_{i}$, the equilibrium $p_{i}$ collides with the singular tangency
point $F^{*}$.

The point $\mathrm{FB}_{i}$ divides the curve $\mathrm{FF}$ locally into
segments $\mathrm{FF}_{1}$ and $\mathrm{FF_{2}}$, which is the case of fold-
fold bifurcation of type VI1 as presented in [31], and the segment
$\mathrm{FU}$ of fused-focus bifurcation [31, 33] along which the (crossing)
periodic orbit $\Gamma$ (dis)appears. Both of these fold-fold bifurcations
result in the sliding segment changing between being repelling and attracting,
and the quadratic tangency points $F_{i}$ switching places as the sliding
segment boundaries.

The point $\mathrm{FB}_{i}$ also divides the curve $\mathrm{BE}_{i}$ locally
into segment $\mathrm{BE}^{P}_{i}$, where there is a standard persistence
boundary equilibrium bifurcation with a nodal equilibrium as presented in
[31], and a segment $\widehat{\mathrm{BE}}_{i}^{P}$ along which a stable
(crossing) periodic orbit $\Gamma$ (dis)appears in a homoclinic-like
persistence boundary equilibrium bifurcation.

  2. 2.

A _double-boundary equilibrium bifurcation_

$\displaystyle\mathrm{BB}:\ \
(\mu,\eta)=\left(\frac{\kappa_{1}\kappa_{2}}{(\kappa_{1}+1)(\kappa_{2}+1)},\ \
\frac{-1}{(\kappa_{1}+1)(\kappa_{2}+1)}\right),$ (26)

occurs at the intersection of the curves $\mathrm{BE_{1}}$ and
$\mathrm{BE_{2}}$. At point $\mathrm{BB}$ the equilibria $p_{1}$, and $p_{2}$
simultaneously collide at the two different quadratic tangency points $F_{1}$
and $F_{2}$, respectively.

The point $\mathrm{BB}$ divides the curve $\mathrm{BE_{1}}$ locally into
segment $\widehat{\mathrm{BE}}_{1}^{P}$, along which a stable (crossing)
periodic orbit $\Gamma$ (dis)appears, and a segment
$\widetilde{\mathrm{BE}}_{1}^{P}$, where there is a standard persistence
boundary equilibrium bifurcation with a nodal equilibrium as presented in
[31]. Similarly, the curve $\mathrm{BE_{2}}$ is divided locally by
$\mathrm{BB}$ into the segment $\widehat{\mathrm{BE}}_{2}^{F}$, along which
the (crossing) periodic orbit $\Gamma$ (dis)appears, and a segment
$\mathrm{BE}^{F}_{2}$, where there is the standard non-smooth fold boundary
equilibrium bifurcation with a nodal equilibrium [31].

  3. 3.

_Generalized boundary equilibrium bifurcations_ [30]

$\displaystyle\mathrm{GB_{1}}:\ \
(\mu,\eta)=\left(\frac{\kappa_{1}^{2}}{(\kappa_{1}+1)^{2}},\ \
-\frac{1}{(\kappa_{1}+1)^{2}}\right),$ (27) $\displaystyle\mathrm{GB_{2}}:\ \
(\mu,\eta)=\left(\frac{\kappa_{2}^{2}}{(\kappa_{2}+1)^{2}},\ \
-\frac{1}{(\kappa_{2}+1)^{2}}\right),$ (28)

occur at end points of the curve $\mathrm{PS}$, respectively, on the curves
$\mathrm{BE}_{1}$ and $\mathrm{BE}_{2}$. At the point $\mathrm{GB}_{i}$,
equilibrium $p_{i}$ collides with the quadratic tangency point $F_{i}$. At the
same time, a pseudo-saddle-node bifurcation takes place at $F_{i}$, resulting
in a generalised boundary equilibrium bifurcation with respect to $f_{i}$.

The point $\mathrm{GB_{1}}$ separates the curve $\mathrm{BE_{1}}$ locally into
the segment $\widetilde{\mathrm{BE}}^{P}_{1}$ and the segment
${\mathrm{BE}}^{F}_{1}$ of boundary equilibrium bifurcations. Similarly, the
point $\mathrm{GB_{2}}$ separates the curve $\mathrm{BE_{2}}$ locally into the
segment $\widehat{\mathrm{BE}}^{P}_{2}$ and the segment
$\widehat{\mathrm{BE}}^{F}_{2}$ of boundary equilibrium bifurcations.

###  2.4 Phase portraits at codimension-one bifurcations

We now present in Figures 7–11 phase portraits in the $(x,y)$-plane for each
segment of codimension-one bifurcation introduced in Proposition 4 and shown
and labeled accordingly in Figure 3. Here, we take a global view of each such
transition to allow for comparison with the respective neighbouring
structurally stable phase portraits in Figures 4–6.

####  2.4.1 Fold-fold and pseudo-Hopf bifurcations

Figure 7: Representative phase portraits of segments
$\mathrm{FF_{1}},\mathrm{FU}$ and $\mathrm{FF_{2}}$ along $\eta=0$, each with
a singular tangency point $F^{*}$. Panel (a) for $\mu=0.0225$ on
$\mathrm{FF}_{1}$ shows the global attractor $p_{1}$ with $W^{ss}(p_{1})$ and
the equilibrium $p_{2}$. Panel (b1) for $\mu=0.25$ on $\mathrm{FU}$ and the
magnification (b2) shows the weakly attracting fold-fold point $F^{*}$; a
representative orbit is highlighted in purple. Panel (c) for $\mu=0.7$ on
$\mathrm{FF}_{2}$ shows the global attractor $p_{2}$ with $W^{ss}(p_{2})$.

Figure 7 shows the phase portraits along segments $\mathrm{FF_{1}}$,
$\mathrm{FU}$, and $\mathrm{FF_{2}}$, each of which with a singular tangency
point $F^{*}$. The phase portrait along segment $\mathrm{FF_{1}}$, which
separates regions $\mathrm{I}$ and $\mathrm{IV}$, is shown in panel (a). The
admissible equilibrium $p_{1}\in R_{1}$ is a global attractor with strong
stable manifold $W^{ss}(p_{1})$. The singular tangencyat the point $F^{*}$ is
invisible to $f_{2}$ and visible to $f_{1}$, and orbits of system (8) are
collinear at $F^{*}$. The phase portrait along segment $\mathrm{FU}$,
separating regions $\mathrm{II}$ and $\mathrm{V}$, is shown in panel (b1) with
a magnification near $F^{*}$ in panel (b2). The singular tangency at $F^{*}$
is now invisible to both vector fields $f_{1}$ and $f_{2}$, and orbits are
anti-collinear at $F^{*}$. Therefore, orbits spiral inward toward $F^{*}$ (at
a very slow rate), and this point is a global attractor. This situation is
reminiscent of a (supercritical) Hopf bifurcation for smooth dynamical
systems, which is why this bifurcation is also known as a pseudo-Hopf
bifurcation [33]. The phase portrait along segment $\mathrm{FF_{2}}$, which
separates regions $\mathrm{III}$ and $\mathrm{VI}$, is presented in panel (c).
The singular tangency at $F^{*}$ is now visible to $f_{2}$ and invisible to
$f_{1}$, and $p_{2}\in R_{2}$ is admissible and the global attractor.

####  2.4.2 Boundary equilibrium and pseudo-saddle-node bifurcations

The phase portraits along the segments of the curves $\mathrm{BE}_{1}$ and
$\mathrm{BE}_{2}$ from Proposition 3 are characterised by an equilibrium of
system (8) being on the switching manifold, but they have different global
manifestations.

Figure 8: Representative phase portraits along the segments
$\mathrm{BE}_{1}^{P}$ and $\mathrm{BE}_{2}^{P}$ on the horizontal slice
$\eta=0.35$ of the $(\mu,\eta)$-plane, each with an attracting sliding segment
$\Sigma_{s}^{a}$. Panel (a) for $\mu=0.1259$ on $\mathrm{BE}_{1}^{P}$ shows
boundary equilibrium $p_{1}$ with $W^{ss}(p_{1})$. Panel (b) for $\mu=-0.25$
on $\mathrm{BE}_{2}^{P}$ shows the boundary equilibrium $p_{2}$ with
$W^{ss}(p_{2})$.

The phase portrait along segment $\mathrm{BE}_{1}^{P}$, which separates
regions $\mathrm{I}$ via $\mathrm{II}$, is shown in Figure 8(a). It has the
attracting sliding segment $\Sigma_{s}^{a}$ bounded by the boundary-node
$p_{1}$ on the left, and by the invisible quadratic tangency point $F_{2}$ on
the right; both pseudo-equilibria $q^{-}$ and $q^{+}$ are non-admissible (and
not shown). The equilibrium $p_{2}$ is not admissible and the boundary-node
$p_{1}$ is a global attractor; note that its strong stable manifold
$W^{ss}(p_{1})$ consists only of the horizontal arriving orbit to $p_{1}$ in
$R_{1}$. The phase portrait in panel (b) along segment $\mathrm{BE}_{2}^{P}$,
separating regions $\mathrm{II}$ and $\mathrm{III}$, is the corresponding
situation but for the boundary-node $p_{2}$: this point is now the global
attractor with strong stable manifold $W^{ss}(p_{2})$ in $R_{2}$, and it
bounds $\Sigma_{s}^{a}$ together with the invisible quadratic tangency point
$F_{1}$.

Figure 9: Representative phase portraits along segments
$\widehat{\mathrm{BE}}_{1}^{P}$ and $\widehat{\mathrm{BE}}_{2}^{P}$ on the
horizontal slice $\eta=-0.2$ in the $(\mu,\eta)$-plane, each with a repelling
sliding segment $\Sigma_{s}^{r}$. Panel (a1) for $\mu=0.0709$ on
$\widehat{\mathrm{BE}}_{1}^{P}$ and the magnification (a2) show the boundary
equilibrium $p_{1}$ with $W^{ss}(p_{1})$, and homoclinic connections
$\gamma^{*}_{1}$ and $\Gamma^{*}_{1}$. Panel (b1) for $\mu=0.3$ on
$\widehat{\mathrm{BE}}_{2}^{P}$ and the magnification (b2) show boundary
equilibrium $p_{2}$ with $W^{ss}(p_{2})$, and homoclinic connections
$\gamma_{2}^{*}$ and $\Gamma_{2}^{*}$.

Figure 9(a1) shows the phase portrait along segment
$\widehat{\mathrm{BE}}_{1}^{P}$, separating regions $\mathrm{IV}$ and
$\mathrm{V}$, with a magnification in panel (a2) near the sliding segment,
which is now repelling. Here $\Sigma_{s}^{r}$ is bounded by the invisible
quadratic tangency point $F_{2}$ on the left and by the boundary-node $p_{1}$
on the right, with both pseudo-equilibria non-admissible (and not shown). The
point $p_{1}$ is globally attracting with strong stable manifold
$W^{ss}(p_{1})$ in $R_{1}$. However, the vector field $f_{2}$ is transverse to
$\Sigma$ at $p_{1}$, and this departing orbit forms a non-sliding homoclinic
connection $\Gamma_{1}^{*}$ back to the boundary-node $p_{1}$. Observe in
panel (a2) that $\Gamma_{1}^{*}$ bounds a region of a family of homoclinic
orbits that involve sliding (in backward time) along the repelling sliding
segment $\Sigma_{s}^{r}$. The orbit labeled $\gamma^{*}_{1}$, consisting of
$\Sigma_{s}^{r}$ and the departing orbit from $F_{2}$, is the maximal sliding
homoclinic orbit: it divides this region inside $\Gamma_{1}^{*}$ into
homoclinic orbit that remain in $R_{1}$ from those that have segments in both
$R_{1}$ and $R_{2}$. The phase portrait along segment
$\widehat{\mathrm{BE}}_{2}^{P}$, which separates regions $\mathrm{V}$ and
$\mathrm{VI}$, is shown similarly in Figure 9(b1) and (b2). The overall
picture is effectively that same, but now $p_{2}$ is the globally attracting
boundary equilibrium on $\Sigma_{s}^{r}$, with analogous non-sliding and
maximal sliding homoclinic orbits $\Gamma_{2}^{*}$ and $\gamma^{*}_{2}$,
respectively. The characterising feature of this bifurcation is the existence
of a non-sliding and crossing homoclinic connection $\Gamma_{i}^{*}$, from
which the stable (crossing) periodic orbit $\Gamma$ in region $\mathrm{V}$
bifurcates; compare with Figure 5(b). This type of (persistence) boundary
equilibrium bifurcation is hardly discussed in the literature; to our
knowledge, it has only been observed in the related Welander’s box model in
[24], where it is referred to as a homoclinic-like boundary equilibrium
bifurcation.

Figure 10: Representative phase portraits along segments
$\widehat{\mathrm{BE}}_{2}^{F}$, $\widetilde{\mathrm{BE}}_{1}^{P}$ and along
the curve $\mathrm{PS}$, each with a repelling sliding $\Sigma_{s}^{a}$. Panel
(a1) for $(\mu,\eta)=(0.14,-0.36)$ on $\widehat{\mathrm{BE}}_{2}^{F}$ and the
magnification (a2) show the boundary equilibrium $p_{2}$ with $W^{ss}(p_{2})$
and homoclinic connection $\Gamma_{2}^{*}$, and the pseudo-equilibrium $q^{+}$
with $W^{uu}(q^{+})$. Panel (b) for $(\mu,\eta)=(-0.0141,-1.05)$ on
$\widetilde{\mathrm{BE}}_{1}^{P}$ shows the admissible equilibrium $p_{2}\in
R_{2}$ with $W^{ss}(p_{2})$, the boundary equilibrium $p_{1}$ with
$W^{ss}(p_{1})$, and the pseudo-equilibrium $q^{-}$ with $W^{u}(q^{-})$ and
$W^{s}_{g}(q^{-})$. Panel (c) for $(\mu,\eta)=(0.09,-0.7636)$ on $\mathrm{PS}$
shows the equilibrium $p_{2}$ with $W^{ss}(p_{2})$ and the singular pseudo-
equilibrium $q^{*}$ with $W^{uu}(q^{*})$.

The phase portrait along segment $\widehat{\mathrm{BE}}^{F}_{2}$, which
separates regions $\mathrm{V}$ and $\mathrm{VII}$, is shown in Figure 10(a1)
with a magnification near the repelling sliding segment in panel (a2). Here,
$\Sigma_{s}^{r}$ is bounded by the attracting boundary-node $p_{2}$ on the
left and by an invisible tangency $F_{1}$ on the right; moreover, it contains
the admissible and repelling pseudo-equilibrium $q^{+}\in\Sigma_{s}^{r}$
(while the pseudo-equilibrium $q^{-}$ is non-admissible and not shown). As was
the case along segment $\widehat{\mathrm{BE}}_{1}^{F}$, the phase portrait in
Figure 10(a) features a (crossing) homoclinic orbit $\Gamma^{*}_{2}$ of
$p_{2}$. However, due to the existence of $q^{+}$ on $\Sigma_{s}^{r}$, this
special orbit does now not bound a region with further (sliding) homoclinic
orbits. Regardless, $\Gamma^{*}_{2}$ is still the limit of the stable
(crossing) periodic orbit $\Gamma$ in region $\mathrm{V}$. Note that all
points inside the region bounded by $\Gamma^{*}_{2}$ converge in backward time
to the unstable pseudo-equilibrium $q^{+}$, whose strong unstable manifold
$W^{uu}(q^{+})$ converges to $p_{2}$; see panel (a2). Segment
$\widetilde{\mathrm{BE}}_{1}^{P}$ separates regions $\mathrm{VII}$ and
$\mathrm{VIII}$, and the phase portrait along it is shown in Figure 10(b).
Here, $p_{1}$ is the attracting boundary-node, and the quadratic tangency
$F_{2}$ is visible. The equilibrium $p_{2}\in R_{2}$ is admissible and also
attracting. Moreover, the pseudo-equilibrium $q^{-}$ lies on the repelling
sliding section $\Sigma_{s}^{r}$, and it is a saddle. Its generalised stable
manifold $W^{s}_{g}(q^{-})$ consists of $\Sigma_{s}^{r}$ and the arriving
orbit to $F_{2}$. The boundary between the basins of attraction of $p_{1}$ and
$p_{2}$ is formed by the union of $W^{s}_{g}(q^{-})$ and the strong stable
manifold $W^{ss}(p_{1})$ in $R_{1}$. Points below these curves and including
$W^{ss}(p_{1})$ converge to $p_{1}$, while points above these curves converge
to $p_{2}$.

The pseudo-saddle-node bifurcation along the curve $\mathrm{PS}$ separates
regions $\mathrm{VI}$ and $\mathrm{VII}$, and its phase portrait is shown in
Figure 10(c). As the name suggests, there is a saddle-node $q^{*}$ of pseudo-
equilibria on the repelling sliding segment $\Sigma_{s}^{r}$, which is the
limiting point where the admissible pseudo-equilibria
$q^{-},q^{+}\in\Sigma_{s}^{r}$ in region $\mathrm{VII}$ (dis)appear. Note that
$q^{*}$ is semi-stable on $\Sigma_{s}^{r}$ and has the strong unstable
manifold $W^{uu}(q^{*})$. The points on $\Sigma_{s}^{r}$ in between the
visible quadratic tangency point $F_{2}$ and $q^{*}$ end up at $q^{*}$ under
the sliding flow; all other points in the $(x,y)$-plane converge to the
admissible and stable equilibrium $p_{2}\in R_{2}$ with strong stable manifold
$W^{ss}(p_{2})$.

Figure 11: Representative phase portraits along the segments
${\mathrm{BE}}_{2}^{F}$ and $\mathrm{BE}_{1}^{F}$, each with a repelling
sliding segment $\Sigma_{s}^{r}$. Panel (a) for $(\mu,\eta)=(-0.05,-0.55)$ on
$\mathrm{BE}_{2}^{F}$ shows the boundary equilibrium $p_{2}$ with
$W^{ss}(p_{2})$ and the admissible equilibrium $p_{1}\in R_{1}$ with
$W^{ss}(p_{1})$. Panel (b) for $(\mu,\eta)=(0.0209,-0.7)$ on
$\mathrm{BE}_{1}^{F}$ similarly shows the boundary equilibrium $p_{1}$ with
$W^{ss}(p_{1})$ and the admissible equilibrium $p_{2}\in R_{2}$ with
$W^{ss}(p_{2})$.

Finally, the boundary equilibrium bifurcations along segments
$\mathrm{BE}_{2}^{F}$ and $\mathrm{BE}_{1}^{F}$ are encountered in the
transition from region $\mathrm{IV}$ via region $\mathrm{VIII}$ to region
$\mathrm{VI}$. In the phase portrait for $\mathrm{BE}_{2}^{F}$ in Figure
11(a), the attracting boundary-node $p_{2}$ bounds the repelling sliding
segment $\Sigma_{s}^{r}$ on the left, while a visible quadratic tangency point
$F_{1}$ bounds it on the right. There are no pseudo-equilibria on
$\Sigma_{s}^{r}$, and the equilibrium $p_{1}\in R_{1}$ is admissible and also
attracting. Trajectories above and including the union of $W^{ss}(p_{2})$ in
$R_{2}$, $\Sigma_{s}^{r}$ and the arriving orbit to $F_{1}$ in $R_{1}$
converge to $p_{2}$, and orbits below this union converge to $p_{1}$. The
phase portrait along segments $\mathrm{BE}_{1}^{F}$ in Figure 11(b) is
effectively the same with the roles of $p_{1}$ and $p_{2}$ exchanged. Here,
the orbits below and including the union of the arriving orbit to $F_{2}$ in
$R_{2}$, $\Sigma_{s}^{r}$ and $W^{ss}(p_{1})$ in $R_{1}$ converge to the
attracting boundary node $p_{1}\in\Sigma_{s}^{r}$, while orbits above this
union converge to the attracting equilibrium $p_{2}\in R_{2}$.

##  3 Bifurcation analysis of the smooth model

We now investigate the smooth model (7) for small $\varepsilon>0$. Here we
again fix the vertical mixing coefficients to $\kappa_{1}=0.1$ and
$\kappa_{2}=1.0$, to enable a direct comparison of the bifurcation diagram of
system (7) with that of the limiting case of system (8).

Figure 12: Two-parameter bifurcation diagram in the $(\mu,\eta)$-plane of
system (7) for $\varepsilon=0.1$ and with $\kappa_{1}=0.1$, $\kappa_{2}=1.0$.
Shown are curves of Hopf bifurcation $\mathrm{H}$ (red, solid when
supercritical, dashed when subcritical), saddle-node bifurcation $\mathrm{S}$
(black when on periodic orbit and grey otherwise) and homoclinic bifurcation
$\mathrm{h}_{1}$ (green), which are the main curves that divide the
$(\mu,\eta)$-plane into the large regions $A,B,C$ and $D$. Also shown are
codimension-two points $\mathrm{CP}$, $\mathrm{BT_{1}}$, $\mathrm{BT_{2}}$,
$\mathrm{GH_{1}}$, $\mathrm{GH_{2}}$, $\mathrm{N_{1}}$ and $\mathrm{N_{2}}$;
grey shading indicates the existence of a stable periodic orbit, and blue
shading bistability between equilibria.

We first consider the bifurcation diagram in the $(\mu,\eta)$-plane of system
(7) for the fixed value of $\varepsilon=0.1$. It is shown in Figure 12 and was
obtained by computing the shown bifurcation curves and codimension-two points
with the continuation package AUTO-07p [27], guided by established bifurcation
theory [26]. One clearly observes four main open regions, denoted $A,B,C$ and
$D$, on which we focus here; associated phase portraits are shown in Figures
13–15.

A main element of the bifurcation diagram in Figure 12 is a curve $\mathrm{S}$
of saddle-node bifurcation with two branches that meet at the cusp point
$\mathrm{CP}$. Along each branch of $\mathrm{S}$ there are points
$\mathrm{BT}_{1}$ and $\mathrm{BT}_{2}$ of Bogdanov-Takens bifurcation (one
close to $\mathrm{CP}$). From these points a curve $\mathrm{H}$ of Hopf
bifurcation emerges, which is the second main element of the bifurcation
diagram. Together, the curves $\mathrm{S}$ and $\mathrm{H}$ effectively form
the boundaries of the four main regions $A$, $B$, $C$ and $D$. Additional
ingredients are: the change of criticality of $\mathrm{H}$ at generalised Hopf
points $\mathrm{GH_{1}}$ and $\mathrm{GH_{2}}$; a curve $\mathrm{h}_{1}$ of
homoclinic bifurcation; and a segment of $\mathrm{S}$, bounded by points
$\mathrm{N_{1}}$ and $\mathrm{N_{2}}$ of non-central homoclinic bifurcation
[34], where the saddle-node bifurcation occurs on a periodic orbit (also known
as SNIC or SNIPER). We remark that the complete bifurcation diagram in the
$(\mu,\eta)$-plane involves subtle additional bifurcation phenomena near the
points $\mathrm{GH_{1}}$, $\mathrm{GH_{2}}$, $\mathrm{N}_{1}$ and
$\mathrm{N}_{2}$ that are indistinguishable on the scale of Figure 12; these
include very narrow regions bounded by additional curves of homclinic
bifurcation and of saddle-node bifurcation of periodic orbits, and their
discussion is beyond the scope of this paper.

###  3.1 Phase portraits in the main regions of the $(\mu,\eta)$-plane

Figure 13: Phase portraits at the points $A_{1}$ at $\mu=0.01$ and $A_{2}$ at
$\mu=0.1$ with $\eta=-0.4$ from region $A$ in Figure 12. Panels (a1) and (b1)
shows the phase portrait on the graph of $\mathcal{H}_{0.1}(x,y)$, and panels
(a2) and (b2) in the $(x,y)$-plane. Featured is the equilibrium $p$, its
strong stable manifold $W^{ss}(p)$ (blue curve) when it exists, and some
representative trajectories (purple curves). Figure 14: Phase portraits at the
point $A_{3}$ at $(\mu,\eta)=(0.1,-0.555)$ and from region $B$ at
$(\mu,\eta)=(0.07,-0.50)$, shown as in Figure 13 and featuring a periodic
orbit $\Gamma$ in panels (b1) and (b2).

Region $A$ of Figure 12 is bounded by the respective (supercritical) part of
the curves $\mathrm{S}$ and $\mathrm{H}$. Comparison with Figure 3 shows that
$A$ is the largest region and ‘covers’ the five regions $\mathrm{I}$,
$\mathrm{II}$, $\mathrm{III}$, $\mathrm{IV}$ and $\mathrm{VI}$ of the PWS
system (8). Throughout region $A$, there is a single attracting equilibrium,
denoted $p$, which may correspond to distinct mixing states: weak (non-
convective) mixing near $\kappa_{1}$, an intermediate state in between
convective and non-convective mixing, or strong (convective) mixing near
$\kappa_{2}$. This is illustrated in Figures 13 and 14(a) with phase portraits
at the parameter points labeled $A_{1}$, $A_{2}$ and $A_{3}$ in Figure 12. We
show all phase portrait of system (7) in two ways to indicate when the
dynamics corresponds to $\kappa_{1}$ or $\kappa_{2}$: on the graph of
$\mathcal{H}_{\varepsilon}(y-x-\eta)$ over the $(x,y)$-plane and on the
$(x,y)$-plane itself, where we use coloring as in Figure 2. At parameter point
$A_{1}$ as in Figure 13(a), the single stable equilibrium $p$ lies in the
region with $\mathcal{H}_{0.1}(y-x-\eta)$ near $0$ (that is, the dynamics of
system (7) is near $\kappa_{1}$), and it has real eigenvalues and a strong
stable manifold $W^{ss}(p)$; hence, $p$ corresponds here to the equilibrium
$p_{1}\in R_{1}$ from regions $\mathrm{I}$ and $\mathrm{IV}$ of the PWS system
(8). Moving to parameter point $A_{2}$ as in Figure 13(b), the equilibrium $p$
now lies in the transition region where the graph of
$\mathcal{H}_{0.1}(y-x-\eta)$ is steep; moreover, it is an attracting focus
with complex conjugate eigenvalues. Finally, at parameter point $A_{3}$ as in
Figure 14(a), the attracting point $p$ has again real eigenvalues and a strong
stable manifold $W^{ss}(p)$, and now lies in the region of the phase plane
with $\mathcal{H}_{0.1}(y-x-\eta)$ near $1$ (that is, the dynamics is now near
$\kappa_{2}$). Hence, $p$ now corresponds to the equilibrium $p_{2}\in R_{2}$
in either regions $\mathrm{III}$ and $\mathrm{VI}$ of the PWS limit. We
conclude that the gradual transition from $A_{1}$ to $A_{3}$ within region $A$
is very reminiscent of that from region $\mathrm{I}$, via region
$\mathrm{II}$, to region $\mathrm{III}$ of system (8); compare with Figure 4.

Region $B$ is bounded by the supercritical part of the curve $\mathrm{H}$ and
the SNIPER-part of $\mathrm{S}$, and it is the ‘smooth version’ of region
$\mathrm{V}$ of system (8). The phase portrait in Figure 14(b), at the marked
parameter point in Figure 12, shows that in region $B$ there is indeed a
stable periodic orbit $\Gamma$ surrounding the now unstable equilibrium $p$.
Observe that $\Gamma$ ‘lives’ in the switching region; that is, it lies on the
steep part of the graph of $\mathcal{H}_{0.1}(y-x-\eta)$. Note further that
the periodic orbit $\Gamma$ bifurcates at the supercritical part of the Hopf
bifurcation curve $\mathrm{H}$ from the attracting focus $p$ of the phase
portrait at $A_{2}$ in Figure 13(b). As $\eta$ is decreased within region $B$,
the periodic orbit $\Gamma$ grows and develops two segments that lie in the
region with $\mathcal{H}_{0.1}(y-x-\eta)$ near $0$ and near $1$, respectively;
these segments correspond to the two segments of the segments periodic orbit
$\Gamma$ in $\mathrm{V}$ of the limiting PWS system (8) in Figure 5(b).

Figure 15: Phase portraits in the region $C$ at $(\mu,\eta)=(0.04,-0.575)$ and
region $D$ at $(\mu,\eta)=(0.01,-0.65)$. Shown in the same manner as in Figure
13, now featuring equilibria $p_{1},p_{3}$, and $q$ with the stable manifold
$W^{s}(q)$ and unstable manifold $W^{u}(q)$.

Region $C$ of Figure 12 is bounded by segments of the two branches of
$\mathrm{S}$ and by the homoclinic bifurcation curve $\mathrm{h}_{1}$ (which
follows closely a subcritical part of the curve $\mathrm{H}$). Figure 15(a)
shows the representative phase portrait at the marked parameter point in
Figure 12. There is an attracting equilibrium, labeled $p_{2}$, with a high
value of the transition function $\mathcal{H}_{0.1}$, as well as a saddle-
equilibrium $q$ and a repelling equilibrium $p_{1}$ with an intermediate value
of $\mathcal{H}_{0.1}$. Note that $q$ has the stable manifold $W^{s}(q)$ and
unstable manifold $W^{u}(q)$, which converges to $p_{2}$. Region $C$ is the
‘smooth version’ of region $\mathrm{VII}$ of the limiting PWS system (8) in
the following way: $p_{2}$ of the smooth system (7) corresponds to $p_{2}\in
R_{2}$, and the equilibria $q$ and $p_{1}$ correspond to the pseudo-saddle-
equilibrium $q^{-}$ and pseudo-equilibrium $q^{+}$ on the repelling sliding
segment $\Sigma^{r}_{s}$, respectively; compare with Figure 6(a).

Finally, region $D$ is bounded by the other segments of the two branches of
$\mathrm{S}$ and the homoclinic bifurcation curve $\mathrm{h}_{1}$. As the
representative phase portrait in Figure 15(b) at the marked parameter point in
Figure 12 shows, it is the region of bistability and corresponds to region
$\mathrm{VIII}$ of the system (8). The attractor $p_{2}$ in Figure 15(b) is
still at a high value of $\mathcal{H}_{0.1}$, and the saddle-equilibrium $q$
is unchanged. However, in contrast to region $C$, the equilibrium $p_{1}$ is
at lower value of the transition function and, moreover, it is now an
attractor. The lower branch of $W^{u}(q)$ converges to $p_{1}$ and its upper
branch to $p_{2}$, meaning that the stable manifold $W^{s}(q)$ forms the
boundary between the basins of attraction of $p_{1}$ and $p_{2}$; compare with
Figure 6(b).

###  3.2 Partial bifurcation analysis in $(\mu,\eta,\varepsilon)$-space

Figure 16: Partial three-parameter bifurcation diagram in
$(\mu,\eta,\varepsilon)$-space of system (7) for $\varepsilon>0$ and system
(8) at $\varepsilon=0$, with $\kappa_{1}=0.1$ and $\kappa_{2}=0.1$.
Represented are curves $\mathrm{S}$ of saddle-node bifurcation (black) and
$\mathrm{H}$ of Hopf bifurcation (red) for fixed values of $\varepsilon>0$,
together with the bifurcation diagram for $\varepsilon=0$ from Figure 3. The
diagram also illustrates the curve $\mathrm{CP}$ of cusp bifurcation (grey),
along with branches $\mathrm{BT}_{1}$ and $\mathrm{BT}_{2}$ of Bogdanov-Takens
bifurcation (dark purple) that meet at the point $\mathrm{DBT}$. Included are
also curves $\mathrm{GH_{1}}$ and $\mathrm{GH_{2}}$ of generalised Hopf
bifurcation (pink).

We now describe how the main elements of the bifurcation diagram of system (7)
in the $(\mu,\eta)$-plane change with the switching-time parameter
$\varepsilon$. Our focus here is on the curves $\mathrm{S}$ and $\mathrm{H}$,
which meet at the Bogdanov-Takens points $\mathrm{BT}_{1}$ and
$\mathrm{BT}_{2}$ and effectively delimit the two main regions of interest,
namely region $B$ characterised by stable oscillations and region $D$
exhibiting bistability. Figure 16 presents the partial three-parameter
bifurcation diagram in $(\mu,\eta,\varepsilon)$-space for
$\varepsilon\in[0,0.16]$ and ranges of $\mu$ and $\eta$ as in Figure 3.
Specifically, Figure 16 shows the bifurcation diagram of the limiting system
(8) for $\varepsilon=0$ together with the curves $\mathrm{S}$ and $\mathrm{H}$
of system (7) as computed for 31 equidistant slices of fixed $\varepsilon>0$.
In this way, the corresponding surfaces $\mathrm{S}$ and $\mathrm{H}$ of
saddle-node and Hopf bifurcation are visualised in
$(\mu,\eta,\varepsilon)$-space with a ‘see-through effect’. Also shown in
Figure 16 is the curve $\mathrm{CP}$ of cusp bifurcation, and the curves
$\mathrm{BT}_{1}$ and $\mathrm{BT}_{2}$ of Bogdanov-Takens bifurcation, along
which the surface $\mathrm{H}$ ends on the surface $\mathrm{S}$. These
codimension-two bifurcation curves were computed directly by numerical
continuation in $(\mu,\eta,\varepsilon)$-space. In particular, this shows that
$\mathrm{BT}_{1}$ and $\mathrm{BT}_{2}$ form a single curve with a maximum at
the point $\mathrm{DBT}$ at $\varepsilon\approx 0.147$, which we identified as
a codimension-three degenerate Bogdanov-Takens point of focus type [35, 36].
Additionally, curves $\mathrm{GH_{1}}$ and $\mathrm{GH}_{2}$ of generalised
Hopf bifurcation are shown in Figure 16; they were found by identifying the
points the corresponding bifurcation points on the curves of Hopf bifurcation
in the individual slices for fixed $\varepsilon$. We observe for increasing
$\varepsilon$ that the curve $\mathrm{GH_{1}}$ ends at the point
$\mathrm{DBT}$. The curve $\mathrm{GH_{2}}$, on the other hand, terminates
where the curves $\mathrm{CP}$ and $\mathrm{BT_{2}}$ intersect at the
codimension-three point $\mathrm{GBC}$ at $\varepsilon\approx 0.1208$.
Similarly, the codimension-two non-central homoclinic bifurcation
$\mathrm{N_{1}}$ and $\mathrm{N_{2}}$ (not shown in Figure 16) are found to
vanish as $\varepsilon$ increases, prior to $\varepsilon$ reaching the value
$\varepsilon\approx 0.147$ of the point $\mathrm{DBT}$. The disappearance of
$\mathrm{N_{1}}$ and $\mathrm{N_{2}}$ involves a sequence of codimension-three
bifurcations, whose analysis is beyond the scope of this paper.

We first consider the relevance of the PWS limiting system (8) for the
bifurcation diagram of the smooth system (7). While the continuation of Hopf
and saddle-node bifurcation of system (7) becomes very challenging for small
values of $\varepsilon$ near $0$, we managed to compute the respective curves
$\mathrm{S}$ and $\mathrm{H}$ in the slice at $\varepsilon=0.005$. As
illustrated in Figure 16, this turns out to be sufficient for determining the
convergence of $\mathrm{S}$ and $\mathrm{H}$ to the corresponding non-smooth
bifurcation as $\varepsilon$ approaches $0$. Specifically, the lower boundary
of the surface $\mathrm{S}$ of saddle-node bifurcation in the
$(\mu,\eta)$-plane at $\varepsilon=0$ is the union of the non-smooth fold
boundary equilibrium bifurcation curve segments
$\mathrm{BE}_{1}^{F},\mathrm{BE}_{2}^{F}$ and $\widehat{\mathrm{BE}}^{F}_{2}$,
and the pseudo-saddle-node bifurcation $\mathrm{PS}$. The surface $\mathrm{H}$
of Hopf bifurcation has as its boundary at $\varepsilon=0$ the union of the
curve segments $\widehat{\mathrm{BE}}_{1}^{P}$,
$\widehat{\mathrm{BE}}_{2}^{P}$ and $\widetilde{\mathrm{BE}}_{1}^{P}$ of
persistence boundary equilibrium bifurcation, and $\mathrm{FU}$ of fused-focus
bifurcation. Moreover, the curves $\mathrm{BT}_{1}$ and $\mathrm{BT}_{2}$ of
Bogdanov-Takens bifurcation converge to the points $\mathrm{GB}_{1}$ and
$\mathrm{GB}_{2}$ of generalised boundary equilibrium bifurcation,
respectively; the curve $\mathrm{CP}$ of cusp bifurcation also converges to
the point $\mathrm{GB}_{2}$. Similarly, the curves $\mathrm{GH}_{1}$ and
$\mathrm{GH}_{2}$ of generalised Hopf bifurcation converge to the points
$\mathrm{FB}_{1}$ and $\mathrm{FB}_{2}$ of fold-boundary equilibrium
bifurcation, respectively.

We now consider the influence of increasing the switching-time parameter
$\varepsilon$. Observe from Figure 16 that the surface $\mathrm{H}$ of Hopf
bifurcation ‘ends’ at the point $\mathrm{DBT}$ at $\varepsilon\approx 0.147$.
Specifically, the curve $\mathrm{H}$ in the $(\mu,\eta)$-plane for fixed
$\varepsilon<0.147$ shrinks to a point at $\mathrm{DBT}$ and disappears. Since
all other curves of codimension-two bifurcations have also disappeared, above
$\mathrm{DBT}$ one only finds the surface $\mathrm{S}$ of saddle-node
bifurcation with the curve $\mathrm{CP}$ of cusp bifurcation. Therefore, in
any slice for fixed $\varepsilon>0.147$ the remaining regions are: region $A$
with a single attracting equilibrium that can take any value of
$\mathcal{H}_{\varepsilon}$, and the bistability region $D$, where two stable
equilibria coexist, one associated with $\mathcal{H}_{\varepsilon}$ near $0$
and the other with $\mathcal{H}_{\varepsilon}$ near $1$. In particular, both
region $B$ with stable oscillations, and region $C$, no longer exists for
$\varepsilon>0.147$. Hence, we conclude that the existence of self-sustained
oscillations in the (adjusted) Welander model requires sufficiently fast
switching between convective and non-convective mixing of surface water with
the deep ocean.

##  4 Discussion and outlook

We studied the adjusted Welander model (7) with transition function
$\mathcal{H}_{\varepsilon}$ between weak and strong mixing between the warm
surface and cold deep ocean as given by (4). This conceptual model in the
context of the AMOC describes the evolution of temperature and salinity on the
ocean surface in the Labrador and Nordic seas. We performed a bifurcation
analysis with advanced tools from (non-smooth) dynamical systems theory, first
for the piecewise-smooth limiting case $\varepsilon=0$ when $\mathcal{H}_{0}$
is the Heaviside function, and then for the smooth case of
$\mathcal{H}_{\varepsilon}$ with small $\varepsilon>0$. Specifically, we
presented bifurcation diagrams in the $(\mu,\eta)$-plane of salinity versus
temperature flux ratio $\mu$ and density threshold $\eta$, where the rates
$\kappa_{1}$ of weak (non-convective) and $\kappa_{2}$ of strong (convective)
mixing were fixed at suitable values. For the PWS model with $\varepsilon=0$,
all curves of codimension-one bifurcations and points of codimension-two
bifurcations were determined analytically — resulting in a complete
description of all possible dynamics and the transitions between them. In this
way, we identified the respective discontinuity-induced bifurcations,
including the continuum of homoclinic orbits investigated in [24], and showed
how these are generated or lost as $\mu$ and $\eta$ change along different
paths. In fact, the bifurcation diagram in the $(\mu,\eta)$-plane we presented
for this case is complete and representative: it does not change in a
qualitative way when a different choice is made for $0<\kappa_{1}<\kappa_{2}$,
as the expressions we derived show. For the smooth case, we computed the
corresponding bifurcation diagram in the $(\mu,\eta)$-plane for
$\varepsilon=0.1$ by means of numerical continuation. While the bifurcation
diagram is complete, we concentrated here on four main regions of dynamics. In
particular, we identified the region with oscillations found in Welander’s
original model [20], as well as a region of bistability that resembles
previously described dynamics in a hierarchy of AMOC models [8].

We also performed a partial bifurcation analysis in
$(\mu,\eta,\varepsilon)$-space for small values of $\varepsilon$, which
focused on surfaces of Hopf and saddle-node bifurcations that (effectively)
bound the main regions. In this way, we showed how the bifurcation diagram for
$\varepsilon>0$ is ‘connected’ to that of the PWS limit. Here, the switching
time $\varepsilon$ plays the role of a parameter that desingularises the
limiting Heaviside swtiching function for $\varepsilon=0$. A direction for
future mathematical work would be to use tools from geometric singular
perturbation theory [37] to study via slow-fast regularisation [38] how
complicated smooth dynamics arises from the piecewise-linear limit. In this
context, we conjecture that the family of homoclinic orbits along the segment
$\widehat{\mathrm{BE}}_{1}^{P}$ will generate a singular Hopf bifurcation with
subsequent canard explosion to the Welander-type large periodic orbit — with
the maximally sliding orbit $\gamma^{*}_{1}$ being the limit of a maximal
canard [39].

Returning to the context of the AMOC, we found for large influxes of
freshwater (smaller $\mu$) that mixing is dominantly non-convective, with the
system approaching a stable equilibrium associated with $\kappa_{1}$.
Conversely, the mixing is dominated by convection for large influxes of
salinity (larger $\mu$), with convergence to a stable equilibrium associated
with $\kappa_{2}$. We found that the intermediate region of bistability in the
AMOC strength exists throughout and is rather independent of the switching
time parameter $\varepsilon$. In contrast, the region of oscillations, where
the AMOC strength changes periodically between strong and weak, does depend on
$\varepsilon$. In fact, oscillations are present only for sufficiently small
$\varepsilon$: when the switching between the two regimes of mixing regimes
becomes too slow, oscialltions are no longer observed.

More generally, the investigation of a conceptual model, such as system (7),
is a tool to uncover and highlight possible types of dynamics one may observe
in the the AMOC. Specifically, we considered here the issue of deep ocean
mixing in the North Atlantic in isolation from the larger climate system. Of
course, there are many other climate processes that influence the overall
state of the AMOC, and the analysis presented should be seen as forming a
basis for the investigation of possible extensions of the model. There are
several interesting directions for future research in this regard, all with
their own mathematical challenges. One option is to consider additonal boxes
in the model, such as an Equatorial box as in Stommel’s original setup [18],
or even to model the two deep-water convection sites in the Labrador sea and
the Nordic seas by separate boxed as in [21]; indeed, such models are of
higher dimensions, which makes their bifurcation analysis more involved.
Another direction is to incorporate seasonal changes, for example, by periodic
forcing the freshwater influx parameter $\mu$, which leads to a non-autonomous
model. Finally, the AMOC displays a number of feedback loops, such as the
salt-advection into the subpolar North Atlantic. Incorporating feedback loops
leads to the study of conceptual climate models in the form of delay
differential equations, the study of which is possible but challenging because
they have an infinite-dimensional phase space [40].

## Acknowledgements

This work was supported in part by Royal Society Te Apārangi Marsden Fund
grant #19-UOA-223. We thank Henk Dijkstra for many helpful discussion,
especially regarding the form of the adjusted Welander model we study here.

##  Appendix A Proofs of Propositions 1–4

We now state and then verify the required properties for the specific case of
system (8), with reference to the literature on planar Filippov systems where
applicable. For in-depth background on general Filippov system theory and the
associated formalism see [30, 31].

Proof of Proposition 1 (Sliding segments and tangency points).

The linear switching manifold $\Sigma$ is given as the zero set of the
switching function $g(x,y)=y-x-\eta$, and it has the constant normal vector
$\mathbf{n}=\binom{-1}{1}$. A tangency point $F_{i}$ occurs when
$(f_{i}\cdot\mathbf{n})(x,y)=0$ (more generally, when the first Lie derivative
of $g$ with respect to $f_{i}$ is zero [30]). With $y=x+\eta$ on $\Sigma$ we
obtain

$\displaystyle(f_{i}\cdot\mathbf{n})(x,x+\eta)=x-1+\mu-\eta\kappa_{i},$ (29)

which yields (13). The visibility of the tangency point $F_{i}$, when it is
quadratic, is determined by the curvature of the orbit of $f_{i}$ from $F_{i}$
relative to $\Sigma$. This is measured by the second Lie derivative of $g$
with respect to $f_{i}$ [30], which for system (8) is given by

$\displaystyle(f_{i}\cdot\nabla(f_{i}\cdot\mathbf{n}))(x,y)=(1+\kappa_{i})(1-(1+\kappa_{i})x)-\kappa_{i}\mu+\kappa_{i}^{2}y,$

where $\nabla$ is the gradient. Evaluating at $F_{i}$ gives

$\displaystyle(f_{i}\cdot\nabla(f_{i}\cdot\mathbf{n}))(F_{i})=\mu+\kappa_{i}(\mu-\eta-1)-\eta\kappa_{i}^{2},$

which yields the genericity condition (14), and the visibility conditions (15)
and (16).

The tangency points $F_{1}$ and $F_{2}$ bound $\Sigma_{s}$ and (29) implies

$\displaystyle(f_{1}\cdot\mathbf{n})(x,x+\eta)>0\ \mathrm{for}\
x>F_{1}\quad\mathrm{and}\quad(f_{2}\cdot\mathbf{n})(x,x+\eta)<0\ \mathrm{for}\
x<F_{2}.$

From (13) we know that $F_{1}<F_{2}$ for $\eta>0$, while $F_{2}<F_{1}$ for
$\eta<0$, which yields (17) and (18). ∎

Proof of Proposition 2 (Equilibria, sliding vector field and pseudo-
equilibria).

  1. 1.

Expression (19) immediately follow from setting $f_{i}(x,y)=0$, and conditions
(20) and (21) are immediate consequences from the definition of $R_{i}$ in
(11) and (12), respectively. The Jacobian

$\displaystyle J_{f_{i}}(x,y)=\left[\begin{matrix}-(1+\kappa_{i})&0&\\\
0&-\kappa_{i}&\end{matrix}\right]$ (30)

of $f_{i}$ has two negative real eigenvalues $\lambda_{ss}=-(1+\kappa_{i})$
and $\lambda_{s}=\kappa_{i}$, which implies that $p_{i}$ is a stable node.
Since $\lambda_{ss}$ has eigenvector $\binom{1}{0}$, the statement on
$W^{ss}_{loc}(p_{i})$ follows.

  2. 2.

The sliding vector field on the line $\Sigma_{s}$ is given by

$\displaystyle f_{s}(x,x+\eta)$
$\displaystyle=\left((1-\lambda(x))f_{1}+\lambda(x)f_{2}\right)(x,x+\eta),$
(31)

where $\lambda(x)\in[0,1]$ is chosen such that the vector $f_{s}$ is in the
(constant) direction $\binom{1}{1}$ of $\Sigma_{s}$. This means that both
components of the vector $f_{s}(x,x+\eta)$ are equal, which this is the case
for

$\displaystyle\lambda(x)=\frac{x+\mu-\kappa_{2}\eta-1}{\eta(\kappa_{1}-\kappa_{2})}.$

Insertion into (31) and simplification yields $f_{s}$ as given in (22).

  3. 3.

Setting $f_{s}(x,x+y)=0$ means solving the quadratic equation

$\displaystyle Q(x):=\mu+(\mu-\kappa_{2}\eta-1)x+x^{2}$

in (22), which gives the expressions for $q^{\pm}$ in (23). The stated
properties follow from evaluating $\frac{dQ(x)}{dx}$ at the $x$-values of
$q^{-}$ and $q^{+}$, respectively. ∎

Proof of Proposition 3 (Codimension-one bifurcations).

  1. 1.

The equilibrium $p_{i}$ from Proposition 2 collides with the switching
manifold $\Sigma$ when

$\displaystyle g(p_{i})=\frac{\mu}{\kappa_{i}}-\frac{1}{\kappa_{i}+1}-\eta=0.$

Solving this for $\eta$ gives the stated expression for $\mathrm{BE}_{1}$ and
$\mathrm{BE}_{2}$. According to (13) the respective boundary equilibrium
bifurcations happens at the tangency point $F_{i}=p_{i}$, and (23) shows that
this involves the (dis)appearance of an admissible pseudo-equilibria through
$F_{i}$. Simultaneously, there is a change in visibility of $F_{i}$ [41, 31],
as can be seen from (15) and (16). It follows that the visibility of the
tangency points $F_{i}$ and the presence of admissible pseudo-equilibria
$q^{\pm}$ in the different regions of the $(\mu,\eta)$-plane are as stated.
See Proposition 4 for details regarding genericity conditions and different
manifestations of the boundary equilibrium bifurcations along the curves
$\mathrm{BE}_{1}$ and $\mathrm{BE}_{2}$.

  2. 2.

For $\eta=0$ we have $F^{*}=F_{1}=F_{2}$ according to (13), which is the
defining property of the fold-fold bifurcation $\mathrm{FF}$; for genericity
conditions and resulting different manifestations see Proposition 4.

  3. 3.

A saddle-node bifurcation of pseudo-equilibria occurs when the square root in
(23) is zero, which gives

$\displaystyle\eta=-(\mu+1)+2\sqrt{\mu}$

and, hence, $\mathrm{PS}$ as stated, and also $q^{*}$ as in (3). The saddle-
node is generic since $\frac{d^{2}Q(x)}{dx^{2}}=2\neq 0$. Since
$(\mu+1)^{2}-4\mu=(\mu-1)^{2}>0$, we know that $\eta<0$ along the curve
$\mathrm{PS}$. Hence, $q^{*}$ lies on $\Sigma_{s}^{r}$ with
$F_{2}<q^{*}<F_{1}$, and the stated bounds for $\mu$ follow. ∎

Proof of Proposition 4 (Codimension-two bifurcations).

  1. 1.

The expressions for $\mathrm{FB}_{i}$ follow immediately from Proposition 3 by
requiring that the curve $\mathrm{FF}$ intersects the curves $\mathrm{BE}_{1}$
and $\mathrm{BE}_{2}$, respectively, yielding $p_{i}=F^{*}$. Note that these
curves intersect transversely at $\mathrm{FB}_{i}$, and the genericity
conditions for $\mathrm{FF}$ and $\mathrm{BE}_{i}$ are satisfied, which means
that the fold-boundary equilibrium bifurcations are generic; see [42].

The points $\mathrm{FB}_{i}$ divide the fold-fold curve $\mathrm{FF}$ into
segments $\mathrm{FF}_{i}$, where $f_{1}(F^{*})$ and $f_{2}(F^{*})$ are
colinear, and a segement $\mathrm{FU}$ where they are not. With

$\displaystyle
f_{i}\cdot\nabla(f_{i}\cdot\mathbf{n})(F^{*})=\mu+(\mu-1)\kappa_{i},$

we conclude that along $\mathrm{FF}_{i}$ the fold-fold bifurcation is for a
visible and and invisible quadratic tangency, which is exactly the case VI1
described in [31]. It also follows that along $\mathrm{FU}$ the fold-fold
bifurcation is for two invisible quadratic tangencies, and with nearby flows
in opposite directions; this identifies this case as a fused-focus bifurcation
according to [31]. The bifurcating (crossing) periodic orbit $\Gamma$ is
stable as demonstrated by the phase portraits presented in Section 2.2. We
remark that the stability of $\Gamma$ can be determined by considering the
(local) return map around $F^{*}$ [30, 31], but this is beyond the scope of
this paper.

The point $\mathrm{FB}_{i}$ also divide $\mathrm{BE}_{i}$ locally as stated;
this follows from the change of stability of the sliding segement and the
associated change from $F_{1}<F_{2}$ for $\eta>0$ to $F_{2}<F_{1}$ for
$\eta<0$; see Proposition 1 and the illustrated and discussed in depth in
Section 2.4.

  2. 2.

At the point $\mathrm{BB}$ of double-boundary equilibrium bifurcation there
are boundary equilibrium bifurcations similtaneously at $p_{1}\neq p_{2}$, and
its location is readily found by equating expressions in Proposition 3 for the
curves $\mathrm{BE_{1}}$ and $\mathrm{BE_{2}}$, which intersect transversally.
It follows that the division of the curves $\mathrm{BE}_{i}$ as are stated;
this is illustrated and discussed in Section 2.4.

  3. 3.

The point $\mathrm{GB}_{i}$ is found by equating the expressions for the
curves $\mathrm{BE}_{i}$ and $\mathrm{PS}$ from Proposition 3. Whether the
boundary equilibrium bifurcation $\mathrm{BE}_{i}$ is of non-smooth fold or
persistence type depends on the sign of the higher-order term [42]

$\displaystyle(\mathbf{n}\cdot(J_{f_{j}})^{-1}\cdot
f_{i})(p_{2})=\frac{1}{\kappa_{j}}\left(\frac{{\kappa_{i}}^{2}}{(\kappa_{i}+1)^{2}}-\mu\right)$

Here $j\neq i\in\\{1,2\\}$ is the respective other index and $J_{f_{j}}$ is
the Jacobian from (30). Hence, a sign change for the curve $\mathrm{BE}_{i}$
happens at the point $\mathrm{GB}_{i}$; specifically, $\mathrm{BE}_{i}$ is of
persistence type for $\mu>\frac{\kappa_{i}^{2}}{(\kappa_{i}+1)^{2}}$ and of
non-smooth fold type for $\mu<\frac{\kappa_{i}^{2}}{(\kappa_{i}+1)^{2}}$. ∎

##  Appendix B Phase portraits at codimension-two bifurcations

Figure 17: Representative phase portrait at codimension-two points
$\mathrm{FB}_{1}$ and $\mathrm{FB}_{2}$. Panel (a) for $(\mu,\eta)=(0.0909,0)$
at $\mathrm{FB}_{1}$ shows the globally stable boundary-node $p_{1}$ with
strong stable manifold $W^{ss}(p_{1})$. Panel (b) for $(\mu,\eta)=(0.5,0)$ at
$\mathrm{FB}_{2}$ shows the globally stable boundary-node $p_{2}$ with the
strong stable manifold $W^{ss}({p_{2}})$.

We present in Figures 17 and 18 phase portraits at the points
$\mathrm{FB_{1},FB_{2}}$, $\mathrm{BB,GB_{1}}$ and $\mathrm{GB}_{2}$ from
Proposition 4. This illustrates how these codimension-two bifurcation points
give the nearby codimension-one boundary equilibrium bifurcations
$\mathrm{BE}_{i}$ and fold-fold bifurcations $\mathrm{FF}$ their different
flavours.

Figure 17 presents phase portraits at the fold-boundary equilibrium
bifurcation points $\mathrm{FB_{1}}$ and $\mathrm{FB_{2}}$. The phase portrait
at $\mathrm{FB}_{1}$ is shown in panel (a). It features an attracting
boundary-node $p_{1}$ that is simultaneously a singular tangency point, which
is invisible for $f_{2}$. All orbits converge to $p_{1}$ along the weak
eigendirection in $R_{1}$. The equilibrium $p_{2}$ is in $R_{1}$ and non-
admissible. The phase portrait at $\mathrm{FB}_{2}$ is shown in panel (b). The
equilibrium $p_{2}$ is now the attracting boundary-node and an invisible
tangency point for $f_{1}$ that attracts all orbits along the weak
eigendirection in $R_{2}$. In both cases, the strong stable manifold
$W^{ss}(p_{i})$ of $p_{i}$ is the corresponding arriving orbit in $R_{i}$.

Figure 18: Representative phase portraits at codimension two points
$\mathrm{GB_{2}},\mathrm{BB}$ and $\mathrm{GB_{1}}$. Each phase portrait
features a repelling sliding segment $\Sigma_{s}^{r}$. Panel (a) for
$(\mu,\eta)=(0.0083,-0.8264)$ at $\mathrm{BB}$ shows the boundary-nodes
$p_{1}$ and $p_{2}$ with corresponding strong stable manifolds $W^{ss}(p_{1})$
and $W^{ss}(p_{2})$. A heteroclinic connection is shown in green. Panel (b)
for $(\mu,\eta)=(0.0455,-0.4545)$ at $\mathrm{GB}_{1}$ shows the quadratic
tangency point $F_{2}$ and a generalized boundary-node $p_{1}$. The
equilibrium $p_{2}$ is shown in $R_{2}$. Panel (c1) for
$(\mu,\eta)=(0.25,-0.25)$ at $\mathrm{GB_{2}}$ and a magnification (c2) near
the sliding segment shows an invisible quadratic tangency point $F_{1}$ and a
generalized boundary-node $p_{2}$ with a strong stable manifold
$W^{ss}(p_{2})$. There are homoclinic connections to $p_{2}$ denoted
$\Gamma^{*}_{2}$ and $\gamma^{*}_{2}$.

Figure 18 presents phase portraits at the remaining codimension two points.
The phase portrait at the double boundary equilibrium bifurcation
$\mathrm{BB}$ at the intersection of the $\mathrm{BE_{1}}$ and
$\mathrm{BE}_{2}$ curves is shown in panel (a). It features a repelling
sliding segment $\Sigma_{s}^{r}$ bounded on the left by the attracting
boundary node $p_{2}$ and on the right by the attracting boundary node
$p_{1}$. The pseudo-equilibria $q^{-}$ and $q^{+}$ are both on $\Sigma_{c}$
and non-admissible (and not shown): the pseudo-equilibrium $q^{-}$ is at the
left hand boundary $p_{2}$, and $q^{+}$ is at the right hand boundary $p_{1}$.
There is a heteroclinic connection between $p_{1}$ and $p_{2}$ composed of
orbit segments in $R_{1}$ and $R_{2}$, respectively. Moreover, a (sliding)
heteroclinic connection between the equilibria is composed of the sliding
orbit from $p_{1}$ to $p_{2}$. If we interpret the departing orbits from
$\Sigma_{s}^{r}$ as having a sliding component, then there is a continuum of
homoclinic connections to $p_{1}$ in $R_{1}$ composed of departing orbits from
$\Sigma_{s}^{r}$ and a corresponding sliding component. Note that the
boundary-node $p_{1}$ has a strong stable manifold $W^{ss}(p_{1})$ composed of
a horizontal component in $R_{1}$. Similarly, the boundary-node $p_{2}$ also
has a strong stable manifold $W^{ss}(p_{2})$, composed of a horizontal
component in $R_{2}$. Overall, the boundary equilibrium bifurcations occurring
simultaneously leads to a pseudo-equilibrium emerging on the sliding segment
along $\widetilde{\mathrm{BE}}_{1}^{P}$ and $\widehat{\mathrm{BE}}_{2}^{F}$;
see Section 2.4. The phase portrait at the generalised boundary equilibrium
bifurcation $\mathrm{GB}_{1}$ is shown in panel (b). It features a repelling
sliding segment $\Sigma_{s}^{r}$ bounded on the left by the visible quadratic
tangency point $F_{2}$ and on the right by the attracting generalised
boundary-node $p_{1}$ (shown in magenta). The pseudo-equilibria $q^{-}$ and
$q^{+}$ undergo a pseudo-saddle-node bifurcation at the right hand boundary-
node $p_{1}$. The phase portrait at the generalised boundary equilibrium
bifurcation $\mathrm{GB}_{2}$ is shown in panel (c1) with a magnification near
the sliding segment in panel (c2). The repelling sliding segment
$\Sigma_{s}^{r}$ is bounded on the left by generalised boundary node $p_{2}$
(shown in magenta) and on the right by the invisible quadratic tangency point
$F_{1}$. The homoclinic connections $\gamma^{*}_{2}$ and $\Gamma_{2}^{*}$ are
the same as described in Section 2.4; see also Figure 9(b). The departing
orbits from $\Sigma_{s}^{r}$ together with the respective sliding component
from $p_{2}$, form a continuum of homoclinic connections to $p_{2}$. In
particular, within $\gamma^{*}_{2}$ there are homoclinic connections composed
of a departing orbit from $\Sigma_{s}^{r}$ in $R_{2}$ and the corresponding
sliding orbit. There are also homoclinic connections inbetween
$\gamma_{2}^{*}$ and $\Gamma_{2}^{*}$, which feature departing orbits from
$\Sigma_{s}^{r}$ in $R_{1}$ that cross $\Sigma$ into $R_{2}$.

## References

  * [1] T. Kuhlbrodt, A. Griesel, M. Montoya, A. Levermann, M. Hofmann, S. Rahmstorf, On the driving processes of the Atlantic meridional overturning circulation, Reviews of Geophysics 45 (2) (2007).
  * [2] K. Speer, S. R. Rintoul, B. Sloyan, The diabatic deacon cell, Journal of Physical Oceanography 30 (12) (2000) 3212–3222.
  * [3] B. M. Sloyan, S. R. Rintoul, The southern ocean limb of the global deep overturning circulation, Journal of Physical Oceanography 31 (1) (2001) 143–173.
  * [4] H. Jeffreys, On fluid motions produced by differences of temperature and humidity, Quarterly Journal of the Royal Meteorological Society 51 (216) (1925) 347–356.
  * [5] J. W. Sandstrom, Meteorologische studien im schwedischen Hochgebirge, Goteborg Wettergren and Kerber 1916, 1916.
  * [6] J. Toggweiler, B. Samuels, On the ocean’s large-scale circulation near the limit of no vertical mixing, Journal of Physical Oceanography 28 (9) (1998) 1832–1852.
  * [7] R. Marsh, W. Hazeleger, A. Yool, E. J. Rohling, Stability of the thermohaline circulation under millennial CO2 forcing and two alternative controls on Atlantic salinity, Geophysical Research Letters 34 (3) (2007).
  * [8] W. Weijer, W. Cheng, S. S. Drijfhout, A. V. Fedorov, A. Hu, L. C. Jackson, W. Liu, E. McDonagh, J. Mecking, J. Zhang, Stability of the atlantic meridional overturning circulation: A review and synthesis, Journal of Geophysical Research: Oceans 124 (8) (2019) 5336–5375.
  * [9] M. Srokosz, M. Baringer, H. Bryden, S. Cunningham, T. Delworth, S. Lozier, J. Marotzke, R. Sutton, Past, present, and future changes in the Atlantic meridional overturning circulation, Bulletin of the American Meteorological Society 93 (11) (2012) 1663–1676.
  * [10] S. Rahmstorf, J. E. Box, G. Feulner, M. E. Mann, A. Robinson, S. Rutherford, E. J. Schaffernicht, Exceptional twentieth-century slowdown in atlantic ocean overturning circulation, Nature Climate Change 5 (5) (2015) 475–480.
  * [11] M. Dima, G. Lohmann, Evidence for two distinct modes of large-scale ocean circulation changes over the last century, Journal of Climate 23 (1) (2010) 5–16.
  * [12] R. R. Dickson, J. Meincke, S.-A. Malmberg, A. J. Lee, The “great salinity anomaly” in the northern North Atlantic 1968–1982, Progress in Oceanography 20 (2) (1988) 103–151.
  * [13] L. Jackson, R. Kahana, T. Graham, M. Ringer, T. Woollings, J. Mecking, R. Wood, Global and european climate impacts of a slowdown of the AMOC in a high resolution GCM, Climate dynamics 45 (11) (2015) 3299–3316.
  * [14] A. Schmittner, Decline of the marine ecosystem caused by a reduction in the Atlantic overturning circulation, Nature 434 (7033) (2005) 628–633.
  * [15] M. Gregory, K. Dixon, R. Stouffer, A. Weaver, E. Driesschaert, M. Eby, T. Fichefet, H. Hasumi, A. Hu, J. Jungclaus, et al., A model intercomparison of changes in the atlantic thermohaline circulation in response to increasing atmospheric CO2 concentration, Geophysical Research Letters 32 (12) (2005).
  * [16] S. Nghiem, D. Hall, T. Mote, M. Tedesco, M. Albert, K. Keegan, C. Shuman, N. DiGirolamo, G. Neumann, The extreme melt across the Greenland ice sheet in 2012, Geophysical Research Letters 39 (20) (2012).
  * [17] H. A. Dijkstra, Nonlinear Physical Oceanography: a Dynamical Systems Approach to the Large Scale Ocean circulation and El Niño, Vol. 28, Springer Science & Business Media, 2005.
  * [18] H. Stommel, Thermohaline convection with two stable regimes of flow, Tellus 13 (2) (1961) 224–230.
  * [19] C. Rooth, Hydrology and ocean circulation, Progress in Oceanography 11 (2) (1982) 131–149.
  * [20] P. Welander, Thermohaline effects in the ocean circulation and related simple models, in: Large-scale transport processes in oceans and atmosphere, Springer, 1986, pp. 163–200.
  * [21] A. Neff, A. Keane, H. A. Dijkstra, B. Krauskopf, Bifurcation analysis of a North Atlantic ocean box model with two deep-water formation sites (2023).  arXiv:2305.11975.
  * [22] P. Welander, A simple heat-salt oscillator, Dynamics of Atmospheres and Oceans 6 (4) (1982) 233–242.  doi:10.1016/0377-0265(82)90030-6.
  * [23] H. E. Huppert, J. S. Turner, Double-diffusive convection, Journal of Fluid Mechanics 106 (1981) 299–329.
  * [24] J. Leifeld, Nonsmooth homoclinic bifurcation in a conceptual climate model (2016).  arXiv:1601.07936.
  * [25] P. Cessi, Convective adjustment and thermohaline excitability, Journal of Physical Oceanography 26 (4) (1996) 481–491.
  * [26] Y. A. Kuznetsov, I. A. Kuznetsov, Y. Kuznetsov, Elements of Applied Bifurcation Theory, Vol. 112, Springer, 1998.
  * [27] E. J. Doedel, T. F. Fairgrieve, B. Sandstede, A. R. Champneys, Y. A. Kuznetsov, X. Wang, Auto-07p: Continuation and bifurcation software for ordinary differential equations, Tech. rep., Concordia University (2007).
  * [28] M. Bernardo, C. Budd, A. R. Champneys, P. Kowalczyk, Piecewise-Smooth Dynamical Systems: Theory and Applications, Vol. 163, Springer Science & Business Media, 2008.
  * [29] A. F. Filippov, Differential equations with discontinuous righthand sides: control systems, Vol. 18, Springer Science & Business Media, 2013.
  * [30] M. Guardia, T. Seara, M. A. Teixeira, Generic bifurcations of low codimension of planar Filippov systems, Journal of Differential Equations 250 (4) (2011) 1967–2023.
  * [31] Y. A. Kuznetsov, S. Rinaldi, A. Gragnani, One-parameter bifurcations in planar Filippov systems, International Journal of Bifurcation and Chaos 13 (08) (2003) 2157–2188.
  * [32] P. T. Piiroinen, Y. A. Kuznetsov, An event-driven method to simulate Filippov systems with accurate computing of sliding motions, ACM Transactions on Mathematical Software (TOMS) 34 (3) (2008) 1–24.
  * [33] J. Castillo, J. Llibre, F. Verduzco, The pseudo-Hopf bifurcation for planar discontinuous piecewise linear differential systems, Nonlinear Dynamics 90 (3) (2017) 1829–1840.
  * [34] S. H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering, Westview Press, 2000.
  * [35] B. Krauskopf, H. M. Osinga, A codimension-four singularity with potential for action, in: Mathematical Sciences with Multidisciplinary Applications, Springer, 2016, pp. 253–268.
  * [36] F. Dumortier, R. Roussarie, J. Sotomayor, H. Zoladek, Bifurcations of Planar Vector Fields: Nilpotent Singularities and Abelian Integrals, Springer, 2006.
  * [37] C. K. R. T. Jones, Geometric Singular Perturbation Theory, Dynamical systems (Montecatini Terme, 1994), Springer-Verlag, Berlin/New York, 1995.
  * [38] M. Jeffrey, D. Duarte Novaes, Regularization of hidden dynamics in piecewise smooth flows, Journal of Differential Equations 259 (9) (2015) 4615–4633.
  * [39] M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. Osinga, M. Wechselberger, Mixed-mode oscillations with multiple time scales, SIAM Review 54 (2) (2012) 211–288.
  * [40] A. Keane, B. Krauskopf, C. M. Postlethwaite, Climate models with delay differential equations, Chaos 27 (11) (2017) 114309.
  * [41] M. Di Bernardo, A. Nordmark, G. Olivar, Discontinuity-induced bifurcations of equilibria in piecewise-smooth and impacting dynamical systems, Physica D: Nonlinear Phenomena 237 (1) (2008) 119–136.
  * [42] F. Della Rossa, F. Dercole, Generalized boundary equilibria in n-dimensional Filippov systems: The transition between persistence and nonsmooth-fold scenarios, Physica D: Nonlinear Phenomena 241 (22) (2012) 1903–1910.