A detailed analysis of the origin of deep-decoupling oscillations

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.

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

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

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

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

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

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$.

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

which are separeted by the switching line

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.

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

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

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

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].

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.

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

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 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).

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.

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].