# A detailed analysis of the origin of deep-decoupling oscillations

John Bailie, Henk A. Dijkstra and Bernd Krauskopf

(May 2024)

###### Abstract

The variability of the strength of the Atlantic Meridional Overturning
Circulation is influenced substantially by the formation of deep water in the
North Atlantic. In many ocean models, so-called deep-decoupling oscillations
have been found, whose timescale depends on the characteristics of convective
vertical mixing processes. Their precise origin and sensitivity to the
representation of mixing have remained unclear so far. To study this problem,
we revisit a conceptual Welander model for the evolution of temperature and
salinity in two vertically stacked boxes for surface and deep water, which
interact through diffusion and/or convective adjustment. The model is known to
exhibit several types of deep-decoupling oscillations, with phases of weak
diffusive mixing interspersed with strong convective mixing, when the
switching between them is assumed to be instantaneous. We present a
comprehensive study of oscillations in Welander’s model with non-instantaneous
switching between mixing phases, as described by a smooth switching function.
A dynamical systems approach allows us to distinguish four types of
oscillations, in terms of their phases of diffusive versus convective mixing,
and to identify the regions in the relevant parameter plane where they exist.
The characteristic deep-decoupling oscillations still exist for non-
instantaneous switching, but they require switching that is considerably
faster than needed for sustaining oscillatory behaviour. Furthermore, we
demonstrate how a gradual freshwater influx can lead to transitions between
different vertical mixing oscillations. Notably, the convective mixing phase
becomes shorter and even disappears, resulting in long periods of much reduced
deep water formation. The results are relevant for the interpretation of
ocean-climate variability in models and (proxy) observations.

##  1 Introduction

Decades of studies of the North Atlantic Ocean circulation with Ocean General
Circulation Models [1, 2, 3, 4, 5] have shown that convective mixing, in
particular, in the Greenland Sea and the Labrador Sea, is involved in the
variability of the Atlantic Meridional Overturning Circulation (AMOC) on
timescales of decades to millennia. This type of strong mixing between the top
layer and deeper parts of the ocean modifies density gradients locally in the
respective geographic regions, which lead to a large scale AMOC response via
associated pressure gradients.

A particular example are the so-called Dansgaard-Oeschger (D-O) events, which
display variability [6] on a millennial timescale. Their modelling has shown
that the variability arises through the interaction of ocean-atmosphere, sea-
ice, and convective mixing, which couples changes in atmospheric buoyancy
fluxes to the AMOC variations. It has been speculated that so-called deep-
decoupling oscillations may describe the D-O events. These oscillations are
characterized by a short phase of convective mixing in the North Atlantic
followed by a very long deep-decoupling phase. During the deep-decoupling
phase, mixing is diffusive and the AMOC re-adjusts, leading to overall low-
frequency variations [7]. For this variability, the AMOC’s response to surface
buoyancy flux anomalies is sensitive to the representation of the vertical
mixing component [8, 9].

Conceptual models have played an important role in investigating the
sensitivity of the AMOC to representations of vertical mixing [10, 11]. The
model by [12] for vertical mixing in the North Atlantic takes a central place
here, because it exhibits periodic solutions with the relevant deep-decoupling
properties [13, 7, 14]. More precisely, it takes the form of a differential
equation for the evolution of temperature and salinity, and its oscillations
are characterized by two phases: a deep-coupling phase, where convective
mixing between the surface and deep ocean in the North Atlantic strongly
‘couples’ the surface and deep ocean, and a deep-decoupling phase when a polar
halocline prevents convective mixing. The [12] model has been analysed
extensively by [10] for the case of instantaneous switching between convective
and diffusive mixing phases, as represented by a step function. [15] studied
the mathematical mechanism of how the periodic solution (dis)appears in this
piecewise-smooth model. A comprehensive bifurcation analyses of all dynamical
regimes of this limiting case of the Welander model has been presented in
[16].

While instantaneous switching between different phases is an interesting
limiting case, most models will feature a continuous and smooth switching
function to represent switching behaviour. In the Welander model, which is
introduced in detail in Sec. 2, switching is realised by a $\tanh$ function
with a switching timescale parameter $\varepsilon$. This switching function
becomes a step function for $\varepsilon=0$, but otherwise the model is
smooth. The existence of periodic solutions for $\varepsilon>0$ was also
studied in [16]; the arguably surprising result is that the Welander model
supports oscillations only when the switching between convective and
diffucsive mixing phases is sufficiently fast.

We present here an investigation into the nature of the oscillations of the
Welander model for positive $\varepsilon$ in terms of the existence of deep-
coupling phases where mixing is convective, and deep-decoupling phases
characterised by weaker diffusive mixing. We take a mathematical approach
grounded in dynamical system theory for systems of smooth differential
equations [17, 18] in combination with the continuation software AUTO [19]. We
define the boundaries of the switching zone for any $\varepsilon\geq 0$, and
this enables us to determine when the periodic solution has a component in the
zones of deep coupling and/or of deep decoupling. For the limiting case
$\varepsilon=0$, we find only oscillations with both a deep-coupling and a
deep-decoupling phase [10, 15, 16], which we refer to as _Welander
oscillations_. For $\varepsilon>0$, however, these are no longer the only type
of oscillation. We present a bifurcation analysis that distinguishes different
sub-regions in the relevant parameter plane (of virtual salinity flux and
density threshold), where the corresponding oscillation has qualitatively
different features. In particular, we find that the existence of Welander
deep-(de)coupling oscillations requires even faster switching than needed to
sustain periodic solutions.

The paper is organised as follows. Section 2 introduces Welander’s model, and
this is followed in Sect. 3 by a brief presentation of the oscillatory regime
for the piecewise-smooth limiting case $\varepsilon=0$. In Sec. 4, we define
the boundary of the switching zone and the zones of deep-(de)coupling for
$\varepsilon>0$. Section 5 then discusses the different types of oscillations
and where they can be found for the representative value $\varepsilon=0.009$.
We start by showing that Welander deep-(de)coupling oscillations still exist,
and show in Sec. 55.1 that they are lost when the underlying periodic orbit
becomes tangent to the boundary of the switching zone. In Sec. 55.2, we then
present bifurcation diagrams that show for which combinations of virtual
salinity flux and density threshold Welander and other types of oscillations
exist; their properties in terms of different oscillation phases are detailed
in Sec. 55.3. By way of an interpretation of these results, Sec. 6 discusses
how the nature of the observed oscillation changes gradually when the virtual
salinity flux slowly decreases, as would be the case with increasing influx of
meltwater from the Greenland ice sheet. Section 7 then shows how the
bifurcation diagram in the parameter plane changes with increasing
$\varepsilon$, as the region of periodic solutions shrinks and then
disappears. The final Sec. 8 discusses our findings and points out directions
for future investigations.

##  2 Conceptual mixing model

Figure 1: Deep water formation in the North Atlantic and its modelling. Panel
(a) presents the Greenland Ice Sheets (GIS), and the Labrador and Nordic seas
with the average mixing depth in meters during January 1958 using ORAS5 global
ocean reanalysis data [20]; here, dark blue regions illustrate deeper mixing
depths, indicating areas of active deep water formation. Panel (b) is a
schematic of the two-box Welander model as described by Eqs. (1).

The [12] model comprises of two vertically stacked boxes representing the
surface and deep ocean near deep water formation sites in the subpolar North
Atlantic, as shown in Fig. 1(a). It describes the evolution of salinity $S$
and temperature $T$ in the surface box of height $H$ that is coupled to a deep
water box with prescribed and fixed temperature $T_{0}$ and salinity $S_{0}$.
Figure 1(b) presents a schematic of this model. Temperature and salinity
between the surface and bottom boxes interact according to a mixing process
with the exchange function $K_{\varepsilon}(\rho)$ that depends on the density
$\rho$ at the surface relative to that of the bottom box. Moreover, the
salinity $S$ in the surface box experiences the virtual salinity flux
$S_{0}F_{0}$, while the temperature $T$ is relaxing to the atmospheric
temperature $T_{a}$ at rate $\gamma$. The overall model takes the form

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

Here, we parameterize the vertical mixing process between the two boxes with
enhanced diffusion, as also done in earlier studies [21, 22]; specifically,
the convective exchange function

$\displaystyle
K_{\varepsilon}(\rho)=k_{1}+\frac{1}{2}(k_{2}-k_{1})\left(1+\tanh{\left(\frac{\rho-\rho_{0}-g^{*}}{\varepsilon}\right)}\right),$
(2)

smoothly transitions between diffusive and convective vertical mixing phases,
given by vertical mixing coefficients $k_{1}$ and $k_{2}$, where
$0<k_{1}<k_{2}$. Further, $g^{*}$ is the density threshold that controls the
transition between the two types of mixing. When the density difference
$\rho-\rho_{0}$ between the surface and bottom box exceeds $g^{*}$, the water
column becomes unstably stratified. Consequently, denser surface water
overlies a lighter deep water, promoting strong, vigorous convective mixing
between the boxes, represented by $k_{2}$. In contrast, when $\rho-\rho_{0}$
falls below $g^{*}$, the water column is stably stratified, resulting in weak,
diffusive mixing represented by $k_{1}$. Note that the slope of
$K_{\varepsilon}(\rho)$ at $\rho-\rho_{0}=g^{*}$ is $1/\varepsilon$. The
smaller $\varepsilon$, the steeper is this slope and, hence, the faster the
transition, which is why we refer to $\varepsilon$ as the switching timescale
parameter. Indeed, for $\varepsilon=0$ the transition is ‘infinitely fast’,
that is, instantaneous.

To complete the model and express $K_{\varepsilon}(\rho)$ as a function of $S$
and $T$, we follow standard practice in many conceptual box models [12, 10,
23] and employ a linear equation of state

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

Here, $\alpha_{S}$ and $\alpha_{T}$ are the saline expansion and thermal
compression coefficients, respectively, and the bottom box is taken as a
prescribed reference state to measure the relative dynamics of the surface
box.

The analysis of system (1) with $K_{\varepsilon}(\rho)=K_{\varepsilon}(S,T)$
as given by (2) and (3) is greatly simplified by non-dimensionalising [10,
16], which reduces the number of parameters. The resulting model is given by

$\displaystyle\dot{x}=1-x-\mathcal{K}_{\varepsilon}(x,y)x,$ (4)
$\displaystyle\dot{y}=\mu-\mathcal{K}_{\varepsilon}(x,y)y,$

for the non-dimensional temperature $x$ and salinity $y$; here, the dot
represents the derivative with respect to the rescaled time and $\mu$ is the
rescaled virtual salinity flux. Additionally, the convective exchange function
now takes the form

$\displaystyle\mathcal{K}_{\varepsilon}(x,y)=\kappa_{1}+\frac{1}{2}(\kappa_{2}-\kappa_{1})\left(1+\tanh{\left(\frac{y-x-\eta}{\varepsilon}\right)}\right),$
(5)

where $\kappa_{1}$ and $\kappa_{2}$ are the non-dimensional vertical mixing
coefficients and $\eta$ is the non-dimensional density threshold.
Conveniently, the non-dimensional density takes the simple form $\rho=y-x$, so
that we have $\mathcal{K}_{\varepsilon}(\rho)=\mathcal{K}_{\varepsilon}(x,y)$.

We refer to system (4) with (5) as the adjusted Welander model, and it is our
main object of study. Specifically, we fix the vertical mixing coefficients at
$\kappa_{1}=0.1$ and $\kappa_{2}=1.0$, as in [16], and investigate the
location in the $(\mu,\eta)$-plane of oscillations with both deep-decoupling
and deep-coupling phases, and how their existence depends on the switching
timescale parameter $\varepsilon$. To this end, we make use of the fact that
the adjusted Welander model is an example of a ‘switched’ slow-fast system
[24, 25]. More specifically, the convective exchange function
$\mathcal{K}_{\varepsilon}(\rho)$ in (5) is a switch from $\kappa_{1}$ for
small $\rho$ to $\kappa_{2}$ for large $\rho$; the switching happens over a
$\rho$-interval around $\rho=\eta$, whose width shrinks to $0$ with
$\varepsilon$.

##  3 Welander oscillations of the piecewise-smooth limit

Figure 2:  Bifurcation diagrams and representations of the stable periodic
solution of system (4) with $\varepsilon=0$. Panel (a) is the bifurcation
diagram in the $(\mu,\eta)$-plane with curves $\mathrm{BE_{1}}$,
$\mathrm{BE_{2}}$, $\mathrm{BB}$, and $\mathrm{PS}$ of discontinuity-induced
bifurcations (coloured as in [16]) that bound the region $W$ of Welander
oscillations (shaded purple). Panel (b) shows the one-parameter bifurcation
diagram in $\mu$ for fixed $\eta=-0.17$, with branches of equilibria (black;
solid when stable, dashed when unstable) and the branch of periodic solutions
represented by their minimum and maximum values $y_{\rm min}$ and $y_{\rm
max}$ (green curves). Panel (c) shows the time series of temperature $x$
(blue) and salinity $y$ (green) for $(\mu,\eta)=(0.219,-0.17)$, at the marker
in panel (a1). Panel (d1) shows the corresponding stable periodic solution
$\Gamma$ (green) on the piecewise constant graph of
$\mathcal{K}_{\varepsilon}$ over the $(x,y)$-plane; here, the graph is
coloured blue where $\mathcal{K}_{\varepsilon}=\kappa_{1}$ over $R_{1}$ and
red where $\mathcal{K}_{\varepsilon}=\kappa_{2}$ over $R_{2}$, which are
connected by the instantaneous vertical switching plane $\Sigma$ (grey).
Panels (d2) and (d3) are the correspondingly coloured time series of
$\mathcal{K}_{\varepsilon}$ and $\rho$ along $\Gamma$, respectively.

In the limit $\varepsilon=0$, the switching between $\kappa_{1}$ and
$\kappa_{2}$ is instantaneous and system (4) is no longer smooth; rather, it
is now a piecewise linear system with a discontinuity along the line where
$\rho=y-x=\eta$. As the starting point of our investigation, we recall and
present where in the $(\mu,\eta)$-plane one finds a piecewise smooth periodic
orbits, corresponding to Welander oscillations (see [16] for a complete
bifurcation analysis of this limiting case).

For $\varepsilon=0$, system (4) with (5) is known as a planar Filipov system
[26, 27, 28] with the two adjacent regions or zones

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

which are separeted by the switching line

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

Zone $R_{1}$ represents a stably stratified water column with diffusive mixing
and, conversely, zone $R_{2}$ represents an unstably stratified water column
with associated vigorous convective mixing that couples the surface and deep
water boxes. The dynamics in $R_{1}$ is determined by setting
$\mathcal{K}_{\varepsilon}(\rho)\equiv\kappa_{1}$ in (4), and that in $R_{2}$
by setting $\mathcal{K}_{\varepsilon}(\rho)\equiv\kappa_{2}$. Both of these
systems are linear, and this implies that any periodic solutions of the planar
Filipov system (4) with $\varepsilon=0$ must consist of a segment in zone
$R_{1}$ and a segment in zone $R_{2}$ and, thus, cross the switching line
$\Sigma$ twice.

Figure 2 illustrates that this type of oscillating solution exists, namely for
parameter combinations of $\mu$ and $\eta$ in the shaded parameter region
labeled W of the $(\mu,c)$-plane shown in panel (a). Along the curves
$\mathrm{FF}$, $\mathrm{BE}_{1}$, $\mathrm{BB}$ and $\mathrm{BE_{2}}$ that
bound this region one finds discontinuity-induced bifurcations [27, 29, 26]
that lead to the emergence or disappearance of the periodic solution [16].
When $\eta$ is fixed at a correpsonding value, a stable periodic solution
$\Gamma$ exist for $\mu$-values between the discontinuity-induced bifurcations
$\mathrm{BE_{1}}$ and $\mathrm{BE_{2}}$. This is illustrated in Fig. 2(b) for
$\eta=-0.17$, where $\Gamma$ is represented by its maximal and minimal
$y$-values; note that system (4) with $\varepsilon=0$ converges to stable
equilibria outside the range of the stable periodic solution. Panel (c) shows
the stable oscillating solution for $(\mu,\eta)=(0.219,-0.17)$ as the time
series of temperature $x$ and salinity $y$. Notice the piecewise smooth nature
of these Welander oscillations, with corners (discontinuities of their
derivative) at the maxima and minima; they are indeed consistent with
oscillations found in earlier studies [10, 12].

To clarify which parts of the stable periodic orbit are in $R_{1}$ and in
$R_{2}$, respectively, we show $\Gamma$ in Fig. 2(d1) on the piecewise-
constant graph of the convective exchange function $\mathcal{K}_{\varepsilon}$
with $\varepsilon=0$; the switching line $\Sigma$ is represented here by a
vertical plane connecting the plateaux of constant $\kappa_{1}$ over $R_{1}$
and $\kappa_{2}$ over $R_{2}$. As is illustrated also by the time series of
$\mathcal{K}_{\varepsilon}$ in Fig. 2(d2), the periodic solution $\Gamma$
indeed features an alternation of deep-decoupling and deep-coupling phases of
considerable lengths, with instantaneous switching at $\Sigma$ between them.
This oscillating behaviour is represented in panel (d3) by the time series of
the density $\rho$ in the surface box. Notice the non-smooth nature of the
oscillation whenever the threshold $\rho=\eta$ defining $\Sigma$ is crossed.
The density is above this threshold during the deep-coupling phase, where the
mixing is convective and vigorous. This phase ends abruptly when $\Gamma$
crosses $\Sigma$, which results in an immediate switch to a stable
stratification of the water column. This is the beginning of the deep-
decoupling phase with only weak, diffusive mixing, which ends when $\Gamma$
crosses $\Sigma$ again and there is, hence, an immediate switch back to the
deep-coupling phase. We remark that the deep-decoupling phase of $\Gamma$ is
longer when $\mu$ is closer to the boundary $\mathrm{BE}_{1}$ in Fig. 2(a) and
(b). Conversely, higher values of $\mu$ near $\mathrm{BE_{2}}$ lead to a
longer deep-coupling phase.

##  4 Switching zone of the smooth model

Figure 2 shows that any periodic orbit of system (4) with $\varepsilon=0$
necessarily has both a deep-coupling phase and a deep-decoupling phase, which
is the characterising property of Welander oscillations. We now address the
question whether and where Welander oscillations can be found when
$\varepsilon>0$. Our starting point is the analysis in [16], which showed that
there exists a region of the $(\mu,\eta)$-plane where a stable periodic orbit
$\Gamma$ exists, provided $\varepsilon$ is sufficiently small, namely below
$\varepsilon^{*}\approx 0.147$. However, no statement has been made previously
as to when this periodic orbit is actually a Welander oscillation.

For non-zero $\varepsilon>0$ the transition between the deep-decoupling phase
near $\kappa_{1}$ and the deep-decoupling phase near $\kappa_{2}$ is no longer
instantaneous. Instead, it occurs within an intermediate (narrow) switching
zone of the $(x,y)$-plane that corresponds to the steep part of the now smooth
convective exchange function $\mathcal{K}_{\varepsilon}$ from (5). To make
this idea precise, one needs to decide when the switching begins and ends.
This can be done in different ways, such as by setting threshold values for
$\mathcal{K}_{\varepsilon}$ or by considering a percentage deviation from
$\kappa_{1}$ and $\kappa_{2}$.

We adopt here a geometric approach and define the switching part to lie in
between the two points of maximal curvature of the graph of
$\mathcal{K}_{\varepsilon}$ as a function of the density $\rho$. These points
can be found from (5) as the zeros of the derivative of the curvature along
$\mathcal{K}_{\varepsilon}(\rho)$; their $\rho$-values are given as the roots
$\rho^{\pm}$ of the function

$\displaystyle G(\rho)$
$\displaystyle=8\varepsilon^{2}-4(\kappa_{1}-\kappa_{2})$ (6)
$\displaystyle+\left(9\varepsilon^{2}+8(\kappa_{1}-\kappa_{2})^{2}\right)\cosh\left(\frac{2(\rho-\eta)}{\varepsilon}\right)$
$\displaystyle-\varepsilon^{2}\cosh\left(\frac{6(\rho-\eta)}{\varepsilon}\right),$

and the corresponding $\kappa$-values are
$L^{\pm}=\mathcal{K}_{\varepsilon}(\rho^{\pm})$ according to (5) with
$\rho=y-x$.

Figure 3:  The maximum curvature points $L^{\pm}$ of the convective exchange
function $\mathcal{K}_{\varepsilon}(\rho)$. Panel (a) shows them for
$\varepsilon=0.02$ and $\eta=-0.17$ on the graph of
$\mathcal{K}_{\varepsilon}(\rho)$, which is coloured blue up to $L^{-}$ and
red from $L^{+}$, with transitional, fainter colouring in between. Panels (b)
and (c) show, respectively, the $\rho$-values $\rho^{\pm}$ and corresponding
$\kappa$-values of $L^{\pm}$ as functions of $\varepsilon$ over the range
$[0,0.3]$.

Figure 3 illustrates the location of the maximum curvature points, which we
also denote $L^{\pm}$ for notational convenience. Panel (a) displays
$\mathcal{K}_{\varepsilon}(\rho)$ for $\varepsilon=0.02$ with $\eta=-0.17$,
with $L^{-}$ and $L^{+}$ clearly lying where the curvature of this graph is
largest. The color scheme indicates that $\mathcal{K}_{\varepsilon}(\rho)$ is
near $\kappa_{1}$ to the left of $L^{-}$ and near $\kappa_{2}$ to the right of
$L^{+}$, with a rapid but smooth switch between the two over the $\rho$-range
in between $L^{-}$ and $L^{+}$. Panels (b) and (c) of Fig. 3 show how the
points $L^{\pm}$ depend on $\varepsilon$ over a considerable range, which
includes the value $\varepsilon^{*}\approx 0.147$ beyond which system (4) is
no longer able to exhibit oscillations [16]. In the limit $\varepsilon=0$ the
associated $\rho$-values coincide at $\rho^{\pm}=\eta$ and the $\kappa$-values
are $\kappa_{1}$ and $\kappa_{2}$, respectively. When $\varepsilon$ is
increased from $0$, the distance in $\rho$ between $L^{-}$ and $L^{+}$
increases quite rapidly, while the corresponding $\kappa$-values remain very
close to $\kappa_{1}$ and $\kappa_{2}$.

Overall, Fig. 3 shows that it quite a natural choice, for any $\varepsilon\geq
0$, to consider the maximal curvature points $L^{\pm}$ as the boundary of the
switching part of the convective exchange function
$\mathcal{K}_{\varepsilon}(\rho)$. It results in the division of the
$(x,y)$-plane of system (4) into the three distinct open regions

$\displaystyle R_{1}$ $\displaystyle=\\{(x,y)\in\mathbb{R}^{2},\
\mathcal{K}_{\varepsilon}(x,y)<L^{-}\\},$ (7) $\displaystyle S$
$\displaystyle=\\{(x,y,)\in\mathbb{R}^{2},\
L^{-}<\mathcal{K}_{\varepsilon}(x,y)<L^{+}\\},$ (8) $\displaystyle R_{2}$
$\displaystyle=\\{(x,y)\in\mathbb{R}^{2},\
L^{+}<\mathcal{K}_{\varepsilon}(x,y)\\},$ (9)

by the two lines along which $\mathcal{K}_{\varepsilon}(x,y)=L^{\pm}$. To
avoid confusion with regions in the parameter $(\mu,\eta)$-plane we discuss in
Sec. 5, we again rather speak of zones: the deep-decoupling zone $R_{1}$, the
switching zone $S$, and the deep-coupling zone $R_{2}$.

For $\varepsilon=0$ this division of the phase plane reduces to that into only
the two zones $R_{1}$ and $R_{2}$, since then $S=\emptyset$ with
$L^{\pm}=\Sigma$. For $\varepsilon>0$, zone $R_{1}$ indeed still represents
the deep-decoupling regime with $\mathcal{K}_{\varepsilon}(x,y)$ very near
$\kappa_{1}$, and zone $R_{2}$ the deep-coupling regime with
$\mathcal{K}_{\varepsilon}(x,y)$ very near $\kappa_{2}$. However, there is now
a local timescale separation: in zones $R_{1}$ and $R_{2}$ the dynamics is
much slower, compared to the fast switching between $\kappa_{1}$ and
$\kappa_{2}$ in zone $S$. We remark that the dynamics in the intermediate
region $S$ represents a slow-fast desingularization or “blow up" with
$\varepsilon$ of the switching line $\Sigma$ [30, 31].

##  5 Welander oscillations with smooth switching

Figure 4: Welander oscillation at $(\mu,\eta)=(0.14,-0.3)$ of system (4) with
$\varepsilon=0.009$. Panel (a1) shows the corresponding stable periodic
solution $\Gamma$ (green) on the graph of $\mathcal{K}_{\varepsilon}$ (colored
as in Fig. 3(a)) with $L^{\pm}$, and panels (d2) and (d3) are the
corresponding time series of $\mathcal{K}_{\varepsilon}$ and $\rho$. Compare
with Fig. 2(d1)–(d3).

Welander oscillations can be found for sufficiently small $\varepsilon>0$, and
an example with $\varepsilon=0.009$ is shown in Fig. 4. Panel (a1) shows the
underlying stable periodic orbit $\Gamma$ on the graph of
$\mathcal{K}_{\varepsilon}$. It indeed features a deep-decoupling phase in
$R_{1}$ and a deep-coupling phase in $R_{2}$, now with transitions through the
switching zone $S$. The time series of $\mathcal{K}_{\varepsilon}$ in panels
(a2) illustrate that this oscillation indeed spends considerable amounts of
time in $R_{1}$ and in $R_{2}$, respectively, with rapid transitions within
$S$ on a faster timescale when $L^{\pm}$ are crossed. The density $\rho$ in
panel (a3) shows the same properties but is overall a more gradual
oscillation: $\rho$ changes quite slowly during the deep-decoupled phase in
$R_{1}$, where it has a minimum, and in the deep-coupling phase in $R_{2}$,
where it has a maximum; the faster transtions through the switching zone $S$
occur in the range between $\rho^{-}$ and $\rho^{+}$. Comparison with Fig.
2(d1)–(d3) shows that the Welander oscillation shown in Fig. 4 is the natural
counterpart for $\varepsilon>0$ of those of the piecewise smooth limit for
$\varepsilon=0$. Indeed, its characterizing feature is that the coexistence of
these deep-(de)coupling phases generates extended periods of inactive and
active deep water production; notice that the actual switching time is short
here, namely about 8% of the period of the oscillation.

###  5.1 Tangencies of periodic solution with $L^{\pm}$

Figure 5: Tangencies of the periodic orbit $\Gamma$ of system (4) with
$\varepsilon=0.009$, shown as in Fig. 4. Panels (a1)–(a3) show the tangency
$T^{+}$ of $\Gamma$ with $L^{+}$ at $(\mu,\eta)=(0.1616,-0.17)$, and panels
(b1)–(b3) show its tangency $T^{-}$ with $L^{-}$ at
$(\mu,\eta)=(0.3092,-0.17)$.

The stable periodic orbit $\Gamma$ in Fig. 4 can be followed or continued as
the virtual salinity $\mu$ and the density threshold $\eta$ are varied. In the
process, the lengths of its deep-decoupling phase in $R_{1}$ and its deep-
coupling phase in $R_{2}$ change. In fact, either phase may shrink down to
zero and disappear, and this happens when $\Gamma$ becomes tangent to $L^{-}$
or $L^{+}$. Figure 5 illustrates two cases for fixed $\eta=-0.17$ in the style
of Fig. 4. Panels (a1)–(a3) of Fig. 5 show the situation when $\Gamma$ is
tangent to the upper threshold $L^{+}$, and we refer to this tangency as
$T^{+}$. While the periodic orbit still crosses the entire switching zone $S$
to just reach $\rho^{+}$ at the maximum of $\rho$, it no longer enters the
zone $R_{2}$ and, hence, $\Gamma$ no longer has a deep-coupling phase.
Physically, the water column is almost sufficiently destabilized: changing the
virtual salinity $\mu$ slightly will result in one of two changes: either a
deep-coupling phase will form in zone $R_{2}$, evolving on a slower timescale
and with a density sufficient for convective mixing to occur, or the solution
will be completely contained within zones $S$ and $R_{1}$. Fig. 5(ab)–(b3)
presents the converse scenario of the tangency $T^{-}$, where $\Gamma$ is
tangent to the lower threshold $L^{-}$. Here the density has a minimum at
$\rho^{+}$ and changing $\mu$ slightly will result in either the formation of
a deep-decoupling phase, or the oscillation being restricted to zones $S$ and
$R_{2}$.

###  5.2 Existence region in the $(\mu,\eta)$-plane

Figure 6: Bifurcation diagrams of system (4) with $\varepsilon=0.009$. Panel
(a) shows the region in the $(\mu,\eta)$-plane where there exists a stable
periodic orbit $\Gamma$; it is bounded by branches $\mathrm{H}$ of Hopf
bifurcation (green) and $\mathrm{S}$ of saddle-node bifurcation (gray, and
black when a SNIC), which meet at the shown points $\mathrm{BT}$ of Bogdanov-
Takens and $\mathrm{N}$ of non-central saddle-node bifurcation. This region is
divided by the tangency curves $T^{+}$ and $T^{-}$ into the four sub-regions
$P_{S}$ (shaded gray), $P_{1}$ (shaded blue), $P_{2}$ (shaded red), and $W$
(shaded purple) of different types of oscillations. One-parameter bifurcation
diagrams in $\mu$ along the accordingly labeled horizontal lines at
$\eta=0.05$, $\eta=-0.17$ and $\eta=-0.3$ are presented in panels (b), (c) and
(d), respectively. Shown are the values of $\rho$ of branches of equilibria
(black; solid when stable, dashed when unstable), and the minimum and maximum
values $\rho_{\rm min}$ and $\rho_{\rm max}$ of $\Gamma$ between the two
points $\mathrm{H}$, which are coloured by the value of
$\mathcal{K}_{\varepsilon}(\rho)$ as in Fig. 3(a); the tangency points $T^{+}$
and $T^{-}$ delimit the $\mu$-ranges of the deep-decoupling phase (blue) and
the deep-coupling phase (red), respectively.

The stable periodic orbit $\Gamma$ of system (4) exists, for sufficiently
small $\varepsilon$, in a region of the $(\mu,\eta)$-plane that is the
continuation of the triangular region $W$ in Fig. 2(a). Figure 6(a) shows this
region for $\varepsilon=0.009$. It is bounded almost entirely by a curve
$\mathrm{H}$ of Hopf bifurcation, all the way from the point $\mathrm{N}$ of
non-central saddle-node bifurcation to the point $\mathrm{BT}$ of Bogdanov-
Takens bifurcation (from where $\mathrm{H}$ emerges). Between $\mathrm{BT}$
and $\mathrm{N}$ the boundary is formed by the part of a curve $\mathrm{S}$ of
saddle-node bifurcations, where the saddle-node lies on the periodic orbit
$\Gamma$ (this is also known as a SNIC or SNIPER bifurcation [17, 18]).
Additional bifurcation curves emerge from the points $\mathrm{N}$ and
$\mathrm{BT}$, but they bound such narrow regions [16] that they are not of
interest for the discussion here.

A crucial difference with the limiting case $\varepsilon=0$ is that the stable
periodic orbit $\Gamma$ for $\varepsilon>0$ is not necessarily a Welander
oscillation. In fact, when it is born at the Hopf bifurcation $\mathrm{H}$ the
periodic orbit is initially very small and surrounds an equilibrium in the
switching zone $S$ [16]; therefore, it is not a Welander oscillation near the
curve $\mathrm{H}$. However, as $\Gamma$ grows, it has tangencies $T^{\pm}$
with $L^{\pm}$ and then enters the respective zones $R_{1}$ and/or $R_{2}$. In
the $(\mu,\eta)$-plane these two types of tangencies from Fig. 5 occur along
the curves $T^{+}$ and $T^{-}$. These two curves can be found by continuing
$\Gamma$ together with the condition that it is tangent to $L^{+}$ and
$L^{-}$, respectively; the boundary value problem setup required for this
computation is described briefly in Appendix A.

As Fig. 6(a) shows, the curves $T^{\pm}$ divide the region of existence of
$\Gamma$ into the four parameter sub-regions $P_{S}$, $P_{1}$, $P_{2}$ and
$W$. As we discussed in Sec. 5.1, when crossing $T^{+}$ and $T^{-}$, the
periodic orbit $\Gamma$ gains or loses its deep-coupling or deep-decoupling
phase. Hence, the bifurcation diagram in the $(\mu,\eta)$-plane provides a
complete characterization of the nature of oscillations in system (4), where
each sub-region corresponds to a different combination in terms of the
existence of deep-(de)coupling phases. In particular, only in region $W$ one
finds Welander oscillations with both a deep-coupling and a deep-decoupling
phase. Note from the shape of region $W$ in Fig. 6(a) that Welander
oscillations require lower values of the density threshold $\eta$, with lower
values of the virtual salinity flux $\mu$ for lower $\eta$. To illustrate
further how the different types of oscillations arise, we present in panels
(b)–(d) of Fig. 6 the one-parameter bifurcation diagrams in $\mu$ for the
three fixed values of $\eta$ that are indicated by the corresponding
horizontal lines in panel (a).

For $\eta=0.05$ as in Fig. 6(b), the periodic orbit $\Gamma$ appears and
disappears for increasing $\mu$ at the Hopf bifurcations $\mathrm{H}$. In
between these two points, the tangency curves $T^{\pm}$ are not encounterd.
Therefore, $\Gamma$ remains entirely in the switching zone $S$ along the
entire shown branch of periodic solutions, which is represented again by the
maximum and minimum values of $\rho$ along $\Gamma$. Figure 6(c) shows the
bifurcation diagram in $\mu$ for $\eta=-0.17$, which is a value for which the
tangency curves $T^{\pm}$ in panel (a) are encountered, namely twice each.
Immediately after the Hopf bifurcation that creates $\Gamma$ in Fig. 6(c), the
oscillation is confinced to the switching zone $S$. However, the tangency
$T^{-}$ with $L^{-}$ takes place already for a slightly larger value of $\mu$.
Subsequently, $\Gamma$ has a deep-decoupling phase, over the $\mu$-range up to
the second encounter with the curve $T^{-}$ at $\mu\approx 0.3092$. The
tangency $T^{+}$ with $L^{+}$ takes place only at $\mu\approx 0.1616$, and
$\Gamma$ has a deep-coupling phase until $T^{+}$ is encountered again, which
happens very close to the Hopf bifurcation $\mathrm{H}$ where $\Gamma$
disappears. Hence, there are Welander oscillations in the range
$\mu\in[0.1616,0.3092]$; the two tangencies shown in Fig. 5 are exactly those
at the two end points $T^{+}$ and $T^{-}$ of this range for $\eta=-0.17$.

In Fig. 6(d) for the even lower value of $\eta=-0.3$, the order in which the
curves $T^{\pm}$ are encountered remains unchanged; however, the two left-hand
instances of $T^{-}$ and $T^{+}$ are now extremely close to the left-most Hopf
bifurcation point $\mathrm{H}$ and, similarly, the right-hand points $T^{-}$
and $T^{+}$ are extremely close to the right-most point $\mathrm{H}$. This
means that, for $\eta=-0.3$, Welander oscillations with definite deep-coupling
and deep-decoupling phases exist, for all practical purposes, throughout the
entire $\mu$-range of existence of the stable periodic orbit $\Gamma$. As the
two-parameter bifurcation diagram in panel (a) shows, this is also the case
along horozontal lines for even lower $\eta$, as long as they are above the
point $\mathrm{N}$. Notice also that the size of $\Gamma$, which is
represented by $\rho_{\min}$ and $\rho_{\max}$ in Fig. 6(d), grows (or
shrinks) extremely rapidly over tiny $\mu$-intervals near the two points
$\mathrm{H}$ of Hopf bifurcation. This phenomenon is typical for slow-fast
systems and known as a canard explosion [32, 31].

###  5.3 The four types of mixing oscillations

Figure 7: The periodic orbit $\Gamma$ of system (4) with $\varepsilon=0.009$,
when it is of type $P_{S}$ at $(\mu,\eta)=(0.32,-0.05)$ in panels (a1)–(a3),
of type $P_{1}$ at $(\mu,\eta)=(0.11,-0.17)$ in panels (b1)–(b3), and of type
$P_{2}$ at $(\mu,\eta)=(0.32,-0.17)$ in panels (c1)–(c3), shown as in Figs. 4
and 5.

The example in Fig. 4 of a Welander oscillation with $(\mu,\eta)=(0.14,-0.3)$
actually shows $\Gamma$ at the star in region $W$ of Fig. 6(a). The
oscillations in the parameter sub-regions $P_{S}$, $P_{1}$ and $P_{2}$ are
shown in the same way in Fig. 7, as represented by $\Gamma$ at the stars in
these corresponding sub-regions of the $(\mu,\eta)$-plane in Fig. 6(a). These
illustrations of $\Gamma$ on the graph of $\mathcal{K}_{\varepsilon}$, with
the times series of $\mathcal{K}_{\varepsilon}$ and $\rho$, provide
comprehensive insight into the nature of the different phases that define each
type of oscillation.

Figure 7(a1) shows that the periodic orbit in sub-region $P_{S}$ is confined
to the switching zone $S$; its time series of $\mathcal{K}_{\varepsilon}$ in
panel (a2) and of $\rho$ in panel (a3) are almost sinosoidal and have a short
period, which is due to $\Gamma$ evolving on the faster timescale associated
with $S$. Indeed, the oscillation in $P_{S}$ lacks distinct deep-(de)coupling
phases, meaning that there is no discernible deep water formation. Figure7(b1)
shows $\Gamma$ in $P_{1}$, which is characterized by a deep-decoupled phase in
zone $R_{1}$ and fast excursions into zone $S$, with the absence of a deep-
coupling phase. The time series of $\mathcal{K}_{\varepsilon}$ in panel (b2)
shows a clear timescale separation: the deep-decoupled phase continues for
quite some time, and is periodically interrupted by short bursts of the period
of the oscillation in panel (a2), with no distinct coupling phase. Panel (b3)
shows that, during the deep-decoupled phase, the density $\rho$ has a value
below $\rho^{-}$ for long times. Hence, there is predominantly diffusive
mixing in a stably stratified water column, with only brief intervals during
which deep water production attempts to recover but fails. In sub-region
$P_{2}$ one finds the converse: $\Gamma$ now has a deep-coupled phase in zone
$R_{2}$ with fast excursions into zone $S$, but no deep-decoupled phase. This
is illustrated in Figure7(c1) on the graph of $\mathcal{K}_{\varepsilon}$ with
the associated time series in panel (c2). The quite long deep-coupling phase
with strong convective mixing is interrupted periodically by the water column
attempting to re-adjust to a stable stratified state. However, this is
unsuccessful and mixing returns to being convective as the periodic solution
re-enters the deep-coupled zone $R_{1}$. Panel (c3) illustrates that the
density is generally sufficient to maintain an active deep water formation
with brief lower-density episodes that are not causing a meaningful shutdown
of convective mixing.

##  6 Drifting virtual salinity flux

We now study the influence of a steadily increasing freshwater influx in the
North Atlantic, which corresponds to slowly decreasing virtual salinity flux,
represented by $\mu$. The focus here is on an intermediate range of the
density threshold $\eta$, where one encounters the sub-regions $P_{1}$, $W$
and $P_{2}$. Our analysis so far reveals that decreasing $\mu$ leads to
progressively shorter deep-coupling phases and longer deep-decoupling phases.
Hence, the deep water formation weakens and the AMOCs heat transport to the
Northern Hemisphere diminishes in a negative convective feedback loop [33].
Furthermore, the deep water re-adjusts for short periods when the surface and
deep-water boxes are decoupled, which can lead to low-frequency variations in
the AMOC.

To model the deepwater formation process in this open system, we allow the
virtual salinity $\mu$ to drift with time $t$ according to the linear ramp
function

$\displaystyle\mu(t)=(1-rt)\mu_{\rm start}+rt\mu_{\rm end}$ (10)

with rate (or slope of the ramp) $r$, which is chosen to be sufficiently small
to ensure that system (4) with (10) adiabatically tracks the respective stable
solutions between $\mu_{\rm start}$ and $\mu_{\rm end}$. In other words,
$\mu(t)$ changes slowly enough so that only gradual transitions between
different types of deep-(de)coupling oscillations are observed. This approach
enables us to identify qualitatively different density oscillations in system
(4) under a real-time gradual freshwater influx — as could be observed in the
modern climate due to an increased meltwater influx from the Greenland Ice
Sheet into the Labrador and Nordic seas [34, 35, 36].

Figure 8: Drifting dynamics of system (4) with $\varepsilon=0.009$ and for
fixed $\eta=-0.17$ when the virtual salinity decreases linearly from $\mu_{\rm
start}=0.32$ to $\mu_{\rm end}=0.11$ with rate of $r=0.001$ in 10. Panel (a1)
shows the ramping $\mu(t)$, and panel (a2) the corresponding time series of
$\rho$ for $t\in[0,1000]$. Also shown are the thresholds $T^{\pm}$ and three
chosen time intervals of width $20$ (shaded in grey). The time series of
$\mathcal{K}_{\varepsilon}$ in each of these time intervals are shown in
panels (b1), (c1) and (d1), and those of $\rho$ in panels (b2), (c2) and (d2).

Figure 8 illustrates the transition as the virtual salinity flux $\mu(t)$
gradually decreases. It begins at the star in sub-region $P_{2}$ and ends at
the star in sub-region $P_{1}$ from the bifurcation diagram in Figure 6(a); in
particular, $\varepsilon=0.009$ and $\eta=-0.17$ here. Specifically, we
consider $\mu(t)$ as given by (10), with the start point $\mu_{\rm
start}=0.32$, as in Figure 7(c1)–(c3), and end point $\mu_{\rm end}=0.11$, as
Figure 7(b1)–(b3), and drift rate $r=0.001$. This downward ramp is illustrated
in Fig. 8(a1) over the $1000$ time units it takes to complete the drift from
start to finish. The corresponding time series of the density $\rho$ of system
(4) is shown in panel (a2), where the initial condition was a point on the
periodic orbit $\Gamma$ in Figure 7(c1). Observe in Fig. 8(a2) the overall
decreasing trend of $\rho$, which is expected to accompany a freshening of the
North Atlantic. This goes along with a shortening of the deep-coupling phase
above $\rho^{+}$, and a considerable lengthening of the deep-coupling phase
below $\rho^{-}$. Indeed, the former phase only exists until $\mu(t)$ drifts
past the tangency $T^{+}$, and the latter only exists past $T^{-}$.

The drift rate $r$ is so small that the time series is effectively periodic
over short time intervals; hence, the properties of the oscillation can be
studied locally near a ‘frozen’ value of $\mu$. We now examine the time series
of $\mathcal{K}_{\varepsilon}$ and $\rho$ in the three windows, of 20 time
units each, that are highlighted in Fig. 8(a2). Panels (b1) and (b2) present
these short time series in the first window with $t\in[5,25]$, where we
observe the characteristic properties of oscillations from sub-region $P_{2}$:
long deep-decoupling phases interrupted by shorter dips in
$\mathcal{K}_{\varepsilon}$ and $\rho$; compare with Figure 7(c2) and (c3).
The second time window with $t\in[300,320]$ is well in between $T^{-}$ and
$T^{+}$, and the corresponding short time series in Fig. 8(c1) and (c2) are
indeed typical for Welander oscillations in sub-region $W$. They exhibit fast
switching between prolonged deep-coupling and deep-decoupling phases of
comparable length, which are here of about the same length here; compare with
Figure 4(a2) and (a3). As $\mu(t)$ monotonically decreases, the length of the
deep-decoupling phases increases, and the periodic convective flushes during
the deep-coupling phases completely vanish when $T^{-}$ is passed. As is
illustrated with Fig. 8(d1) and (d2) for $t\in[950,970]$, the short time
series then clearly consists of long deep-coupling phases with periodic
‘blips’ where the water column briefly looses its stratification, which is the
hallmark of sub-region $P_{1}$; compare with Figure 7(b2) and (b3).

##  7 Influence of the switching timescale parameter

Figure 9: Bifurcation diagram in the $(\mu,\eta)$-plane of system (4) with
$\varepsilon=0.02$ in panel (a), $\varepsilon=0.041$ in panel (b),
$\varepsilon=0.085$ in panel (c), and $\varepsilon=0.11$ in panel (d). Shown
are the sub-regions $P_{S}$ (grey), $P_{1}$ (light blue), $P_{2}$ (light red)
and $W$ (light purple), their bounding curves $T^{-}$ (blue), $T^{+}$ (red),
$\mathrm{H}$ (green) and $\mathrm{S}$ (gray, and black when a SNIC), and the
points $\mathrm{BT}$ and $\mathrm{N}$. Compare with Fig. 6(a) and Fig. 2(a).

The parameter $\varepsilon$ directly influences the switching time between the
convective and non-convective mixing phases, and it plays a crucial role in
determining the configuration of deep-(de)coupling oscillations observed in
system (4). Since any oscillation completely disappears for
$\varepsilon>0.147$ [16], we now address what this means for the fate of the
sub-regions $P_{S}$, $P_{1}$, $P_{2}$, and $W$.

Figure 9 presents the bifurcation diagrams of system (4) in the
$(\mu,\eta)$-plane for four increasingly larger values of $\varepsilon$. As
before, the relevant parts of the curves $\mathrm{H}$ of Hopf bifurcation and
$\mathrm{S}$ of saddle-node bifurcation bound the region of existence of the
stable periodic solution $\Gamma$. Comparison of Fig. 6(a) for
$\varepsilon=0.009$ with Fig. 9(a) for $\varepsilon=0.02$ and Fig. 9(b) for
$\varepsilon=0.041$ reveals that these three bifurcation diagrams have the
same qualitative features in terms of existence and locations of the sub-
regions $P_{S}$, $P_{1}$, $P_{2}$ and $W$. In particular, they all feature
Welander oscillations for lower values of $\eta$ above the point $\mathrm{N}$.
We conclude that the situation in Fig. 6 is representative of any positive
$\varepsilon$ over this initial range. In particular, as $\varepsilon$ is
decreased to $0$ the Hopf bifucation curve $\mathrm{H}$, and the tangency
curves $T^{-}$ and $T^{+}$ converge to the discontinuity indced bifurcations
$\mathrm{BE}_{1}$, $\mathrm{FF}$ and $\mathrm{BE}_{2}$ in Figure 2(a). In the
process, the sub-regions $P_{S}$, $P_{1}$, and $P_{2}$ disappear, and only
region $W$ of Welander oscillations remains in the limit $\varepsilon=0$.

While the bifurcation diagrams in Fig. 6(a) and panels (a) and (b) of Fig. 9
share their qualitative features, there are important quantitative
differences. The existence region exhibiting oscillations decreases in size
with increasing $\varepsilon$ (notice the different scales of the panels of
Fig. 9). At the same time, the relative sizes of the sub-regions $P_{S}$ and
$P_{1}$ increase, while $P_{2}$ and $W$ decrease significantly; they are quite
small in Fig. 9(b) for $\varepsilon=0.041$. As is shown in Fig. 9(c), when the
switching timescale parameter is increased to $\varepsilon=0.085$, the sub-
regions $P_{2}$ and $W$ and their bounding curve $T^{+}$ have actually
disappeared. This menas that the periodic orbit $\Gamma$ is always restricted
to the switching zone $S$ and the deep-decoupling zone $R_{1}$. As a result,
there are no longer any Welander oscillations, and deep water formation is
suppressed. Note that the disappearance of Welander oscillations occurs
substantially below $\varepsilon^{*}\approx 0.147$, where the region of
oscillations itself vanishes. Figure 9(d) presents the bifurcation diagram for
$\varepsilon=0.11$, where now also sub-region $P_{1}$ and its boundary curve
$T^{-}$ have disappeared. Hence, the periodic solution always lies in the (now
larger) switching zone $S$.

Overall, we observe from Fig. 9 that the deep-coupling phase becomes
increasingly limited with increasing $\varepsilon$. Already for ‘intermediate’
values of $\varepsilon$, system (4) is biased toward exhibiting deep-
decoupling oscillations. Moreover, Welander oscillations, and deep-coupling
phases more generally, only exist when the switching between zones $R_{1}$ and
$R_{2}$ is sufficiently fast — considerably faster than required for the
existence of oscillations in the first place.

##  8 Conclusions and outlook

We have provided a detailed analysis of deep-decoupling oscillations in the
Welander model (4) with the smooth convective exchange function given by (5),
which describes the transition between diffusive and convective mixing, the
timescale of which is determined by the additional parameter $\varepsilon$. In
the context of the deepwater formation, this model exhibits self-sustained
oscillations of temperature and salinity at the ocean surface. We conducted a
bifurcation analysis that identifies sub-regions in the $(\mu,\eta)$-plane of
virtual salinity influx and density threshold where one finds different types
of deep-(de)coupling oscillations. The Welander oscillations with coexisting
deep-decoupling and deep-coupling phases, which were shown to be the only type
of oscillation for $\varepsilon=0$, also exist for positive $\varepsilon$, but
only up to a value considerable below the value $\varepsilon^{*}$ where
oscillations disappear altogether. Indeed, this means that other types of
oscillations exist for $\varepsilon>0$, but which exhibit only a deep-
decoupling phase, only a deep-coupling phase, or neither of the two. We
delimited their sub-regions of existence in the $(\mu,\eta)$-plane for any
given $\varepsilon$ by computing curves of tangencies of the underlying stable
periodic orbit with the boundary of the switching zone. Our bifurcation
diagrams ‘extend’ the known results for $\varepsilon=0$ of instantaneous
switching to the more general case where the switching is more gradual and
described by a smooth convective exchange function with $\varepsilon>0$.

Overall, our results complete the analysis of the possible oscillatory
behaviour exhibited by the Welander model. As we have shown, they may provide
insights into the potential qualitative changes in the deepwater formation
under an influx of freshwater. Specifically, we found that large influxes
result in a shortening of the convective phase, leading to prolonged
adjustment during deep-decoupled phases. Conversely, with lower levels of
freshening, the deepwater formation strengthened due to longer phases of
convective mixing and reduced weaker mixing.

We distinguished the deep-decoupling, switching and deep-coupling zones by the
lines where the curvature of the convective exchange function is maximal. This
is convenient mathematical choice that works well for all $\varepsilon\geq 0$,
but other choices for the boundaries of the switching zone could be made, such
as using a percentage deviation from the minimum and maximum values
$\kappa_{1}$ and $\kappa_{2}$. Moreover, one might consider a smooth switching
function different from the hyperbolic tangent in (5). Whatever the choice,
the technique of finding the curves of tangency with these zone boundaries can
still be used to identify where the different types of oscillations exist.
Moreover, if the alternative boundaries are close to the lines $L^{\pm}$ we
considered here, then the respective sub-regions will only differ a little.
Hence, what we presented is broadly representative for any reasonable choice
of smooth switching function and switching thresholds. Since a form of
convective adjustment is a common way to parameterize convective vertical
mixing in (particular relatively low resolution) ocean models, our general
approach can also be applied to the underlying convective adjustment scheme in
more complex models. This would allow one to identify thresholds between deep-
coupling and deep-decoupling regimes with realistic representations of the
ocean and atmosphere dynamics. In particular, it will be interesting to
determine the timescales of these two phases, in comparison with that of the
(generally faster) switching between them; indeed, their ratio would give a
‘physical interpretation’ to the switching timescale parameter $\varepsilon$
of the Welander model. In this way, it may also be possible to provide
insights into the role of deep-decoupling oscillations in modelled D-O events
[9].

A next step in the vein of conceptual modelling would be to consider the
competition between convective and salt-advective feedback mechanisms in the
AMOC [23]. This can be achieved by adding additional boxes to the current
schematic in Figure 1. Alternatively, one could incorporate an advective term
with a time delay to represent the AMOCs adjustment period. However, the
latter yields a model in the form of a delay differential equation with an
infinite dimensional state space, the analysis of which comes with its own
challenges [37].

##  Appendix A Computing tangency curves

Equilibria, periodic orbits and their bifurcations can be found and followed
in parameters with well-established numerical continuation techniques; see,
for example, [38] as an entry point to the literature. We use here the
continuation software package AUTO [19] to find the stable periodic orbit
$\Gamma$, as well as the bifurcation curves that bound the region where it
exists; see [16] for more information on the exact nature of the bifurcations
involved.

The important new aspect here is the computation of the curves $T^{-}$ and
$T^{+}$ by continuation, which we achieve as follows. The periodic orbit
$\Gamma$ is continued from a point of Hopf bifurcation into the region $W$,
where it has two intersetion points with each of the switching lines $L^{-}$
and $L^{+}$. The periodic solution $\Gamma$ is found as an orbit segment
$\bm{u}$ that satisfies a boundary value problem (BVP); more specifically,
$\bm{u}(t)=(\bm{u}_{x}(t),\bm{u}_{y}(t))$ satisfies the differential equation
(4) (after rescaling with the period $T_{\Gamma}$ of $\Gamma$) subject to the
periodicity conditions

$\displaystyle\begin{cases}\bm{u}_{x}(0)=\bm{u}_{x}(1),\\\
\bm{u}_{y}(0)=\bm{u}_{y}(1).\end{cases}$ (11)

Note that the point $\bm{u}(0)$ could be any point along $\Gamma$ for
conditions (11) to be satisfied. This is why, to ensure this BVP is well
defined and has a unique solution, one needs to impose a so-called phase
condition [38, Chapter 1 by E. J. Doedel]. AUTO uses an integral phase
condition [19], but we impose now the alternative condition

$G(\bm{u}_{y}(0)-\bm{u}_{x}(0))=0,$ (12)

where the function $G$ is as defined in (6). Hence, the point $\bm{u}(0)$ is
now a point on $L^{-}$ or $L^{+}$, and this also defines the solution $\bm{u}$
uniquely.

We remark that an initial solution $\bm{u}$ satisfying (11) and (12) can be
found in AUTO by ‘rotating’ the data representing $\Gamma$ until (12) is
satisfied. This orbit segement $\bm{u}$ can then be continued in a parameter,
say, in $\mu$. When $\Gamma$ becomes tangent to the corresponding line $L^{-}$
or $L^{+}$, the continuation detects a fold point in the parameter (which is a
double root of $G$). This fold can then be continued in an additional
parameter, say, in $\eta$, yielding the respective curve $T^{-}$ and $T^{+}$
in the $(\mu,\eta)$-plane.

## References

  * [1] U. Mikolajewicz, E. Maier-Reimer, Internal secular variability in an ocean general circulation model, Climate Dynamics 4 (1990) 145–156.
  * [2] A. J. Weaver, J. Marotzke, P. F. Cummins, E. S. Sarachik, Stability and variability of the thermohaline circulaion, Journal of Physical Oceanography (1993) 39–60.
  * [3] T. Delworth, M. Mann, Observed and simulated multidecadal variability in the Northern Hemisphere, Climate Dynamics 16 (9) (2000) 661–676.
  * [4] D. L. Bars, J. P. Viebahn, H. A. Dijkstra, A Southern Ocean mode of multidecadal variability, Geophysical Research Letters 43 (2016) 1–9.
  * [5] J. Lohmann, H. A. Dijkstra, M. Jochum, V. Lucarini, P. D. Ditlevsen, Multistability and Intermediate tipping of the Atlantic Ocean Circulation, Science Advances 10 (2024) eadi4253.
  * [6] L. C. Menviel, L. C. Skinner, L. Tarasov, P. C. Tzedakis, An ice–climate oscillatory framework for Dansgaard–Oeschger cycles, Nature Reviews Earth & Environment 1 (12) (2020) 677–693.
  * [7] M. Winton, E. S. Sarachik, Thermohaline oscillations induced by strong steady salinity forcing of ocean general circulation models, Journal of Physical Oceanography 23 (7) (1993) 1389–1410.
  * [8] W. Weijer, H. A. Dijkstra, Multiple oscillatory modes of the global ocean circulation, Journal of Physical Oceanography 33 (2003) 2197–2213.
  * [9] G. Vettoretti, P. Ditlevsen, M. Jochum, S. O. Rasmussen, Atmospheric CO2 control of spontaneous millennial-scale ice age climate oscillations, Nature Geoscience 15 (4) (2022) 300–306.
  * [10] P. Cessi, Convective adjustment and thermohaline excitability, Journal of Physical Oceanography 26 (4) (1996) 481–491.
  * [11] A. Colin de Verdière, A simple model of millennial oscillations of the thermohaline circulation, Journal of Physical Oceanography 37 (5) (2007) 1142–1155.
  * [12] P. Welander, Thermohaline effects in the ocean circulation and related simple models, in: Large-scale transport processes in oceans and atmosphere, Vol. 190 of NATO ASI Series, Springer, 1986, pp. 163–200.
  * [13] M. Winton, Energetics of deep-decoupling oscillations, Journal of Physical Oceanography 25 (3) (1995) 420–427.
  * [14] A. Sima, A. Paul, M. Schulz, The younger dryas—an intrinsic feature of late pleistocene climate change at millennial timescales, Earth and Planetary Science Letters 222 (3-4) (2004) 741–750.
  * [15] J. Leifeld, Non-smooth homoclinic bifurcation in a conceptual climate model, European Journal of Applied Mathematics 29 (5) (2018) 891–904.
  * [16] J. Bailie, B. Krauskopf, Bifurcation analysis of a conceptual model for vertical mixing in the North Atlantic, Physica D: Nonlinear Phenomena 460 (2024) 134077.
  * [17] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Vol. 112 of Applied Mathematical Sciences, Springer, 1998.
  * [18] S. H. Strogatz, Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering, CRC press, 2018.
  * [19] E. J. Doedel, A. R. Champneys, F. Dercole, T. F. Fairgrieve, Y. A. Kuznetsov, B. Oldeman, R. C. Paffenroth, B. Sandstede, X. J. Wang, C. H. Zhang, Auto-07p: Continuation and bifurcation software for ordinary differential equations, Tech. rep. (2007).
  * [20] Copernicus Climate Change Service (C3S), ORAS5 global ocean reanalysis monthly data from 1958 to present, Climate Data Store (CDS), 2020.
  * [21] M. Vellinga, Multiple equilibria in ocean models as a side effect of convective adjustment, Journal of Physical Oceanography 28 (4) (1998) 621–633.
  * [22] M. den Toom, H. A. Dijkstra, F. W. Wubs, Spurious multiple equilibria introduced by convective adjustment, Ocean Modelling 38 (1-2) (2011) 126–137.
  * [23] W. Weijer, W. Cheng, S. S. Drijfhout, A. V. Fedorov, A. Hu, L. C. Jackson, W. Liu, E. L. McDonagh, J. V. Mecking, J. Zhang, Stability of the atlantic meridional overturning circulation: A review and synthesis, Journal of Geophysical Research: Oceans 124 (8) (2019) 5336–5375.
  * [24] M. Wechselberger, Geometric Singular Perturbation Theory Beyond the Standard Form, Vol. 6 of Frontiers in Applied Dynamical Systems: Reviews and Tutorials, Springer, 2021.
  * [25] K. U. Kristiansen, P. Szmolyan, Relaxation oscillations in substrate-depletion oscillators close to the nonsmooth limit, Nonlinearity 34 (2) (2021) 1030.
  * [26] M. di Bernardo, C. J. Budd, A. R. Champneys, P. Kowalczyk, Piecewise-Smooth Dynamical Systems: Theory and Applications, Vol. 163 of Applied Mathematical Sciences, Springer, 2008.
  * [27] A. F. Filippov, Differential equations with discontinuous righthand sides: control systems, Vol. 18 of Mathematics and its Applications, Springer Science & Business Media., 2013.
  * [28] M. Guardia, T. M. Seara, M. A. Teixeira, Generic bifurcations of low codimension of planar Filippov systems, Journal of Differential Equations 250 (4) (2011) 1967–2023.
  * [29] 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.
  * [30] M. A. Teixeira, P. R. da Silva, Regularization and singular perturbation techniques for non-smooth systems, Physica D: Nonlinear Phenomena 241 (22) (2012) 1948–1955.
  * [31] C. Kuehn, Multiple Time Scale Dynamics, Vol. 191 of Applied Mathematical Sciences, Switzerland: Springer International Publishing., 2015.
  * [32] 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.
  * [33] H. A. Dijkstra, Nonlinear Physical Oceanography: a Dynamical Systems Approach to the Large Scale Ocean circulation and El Niño, Vol. 28, Springer, 2005\.
  * [34] J. Bamber, M. Van Den Broeke, J. Ettema, J. Lenaerts, E. Rignot, Recent large increases in freshwater fluxes from Greenland into the North Atlantic, Geophysical Research Letters 39 (19) (2012).
  * [35] L. D. Trusel, S. B. Das, M. B. Osman, M. J. Evans, B. E. Smith, X. Fettweis, J. R. McConnell, B. P. Noël, M. R. van den Broeke, Nonlinear rise in Greenland runoff in response to post-industrial arctic warming, Nature 564 (7734) (2018) 104–108.
  * [36] I. Sasgen, M. van den Broeke, J. L. Bamber, E. Rignot, L. S. Sørensen, B. Wouters, Z. Martinec, I. Velicogna, S. B. Simonsen, Timing and origin of recent regional ice-mass loss in Greenland, Earth and Planetary Science Letters 333 (2012) 293–303.
  * [37] A. Keane, B. Krauskopf, C. M. Postlethwaite, Climate models with delay differential equations, Chaos 27 (11) (2017) 114309.
  * [38] B. Krauskopf, H. M. Osinga, J. Galán-Vioque (Eds.), Numerical Continuation Methods for Dynamical Systems: Path following and boundary value problems, Understanding Complex Systems, Springer, 2007.