Regular and chaotic Welander oscillations in a four-dimensional conceptual model for the Atlantic Meridional Overturning Circulation

Abstract: The Atlantic Meridional Overturning Circulation (AMOC) is a key component of the Earth's climate. Evidence indicates a twentieth-century weakening, and enhanced freshwater input to the subpolar North Atlantic may further reduce overturning strength. We present and study a conceptual four-dimensional, single-hemisphere box model with three compartments: a tropical surface box, a subpolar surface box, and a large deep-water box. Advective exchange couples the surface boxes and vertical exchange with the deep ocean is represented by a smooth convective-adjustment scheme. A comprehensive bifurcation analysis reveals an equilibrium structure with up to four coexisting overturning states, together with regimes of bistability and tristability. We identify families of periodic solutions and chaotic attractors with a clear timescale separation: a millennial oscillation is modulated by faster decadal-to-centennial variability arising from episodic shutdowns of subpolar convection. As prescribed freshwater fluxes increase, shutdown events become more frequent and the background overturning weakens. Additionally, for certain values of freshwater influx, the dynamics become chaotic, producing an irregular on-off switching of convection.

• The paper develops a four-dimensional Stommel–Welander AMOC model and finds up to four coexisting stable equilibria through saddle-node, Hopf, and codimension-two bifurcations.
• The paper shows that freshwater forcing produces slow–fast Welander cycles with one to four convective shutdowns per period, as incomplete deep-water temperature recovery clusters events.
• The paper confirms a broad chaotic regime with positive Lyapunov exponents, periodic windows, and coexistence of regular and chaotic oscillations bounded partly by boundary crises.

Overview and motivation

The Atlantic Meridional Overturning Circulation (AMOC) is a recognized tipping element of the climate system, whose weakening under enhanced freshwater input to the subpolar North Atlantic could trigger abrupt changes in other climate subsystems. In this paper, Bailie, Subramanian, and Krauskopf develop a conceptual three-box, four-dimensional ordinary differential equation (ODE) model that couples Stommel-type density-driven overturning with Welander-style convective adjustment, largely following the setup of Zhang et al. The model comprises a tropical surface box (E), a subpolar North Atlantic surface box (N), and a deep-water box (D), with advective exchange between the surface boxes and smooth convective-adjustment coupling of each surface box to the deep ocean. The authors' stated aim is not to reproduce the present-day climate state but to provide a comprehensive dynamical systems analysis—equilibria, periodic orbits, and chaotic attractors—in a two-parameter plane of virtual salinity flux $\mu$ and density threshold $\eta$.

Model formulation

The full six-dimensional model tracks temperature and salinity in each box. Meridional exchange is parameterized by a Stommel-type closure $\Psi = \sigma(\rho_N - \rho_E)$, regularized near $\Psi = 0$ via $\Psi\tanh(\Psi/\vartheta)$ for smoothness. Vertical exchange uses a Welander-type convective adjustment function

$$\mathbf{K}_i = k_d + \tfrac{1}{2}(k_i - k_d)\left(1 + \tanh\!\left(\frac{\rho_i - \rho_D - \Delta\rho}{\varepsilon}\right)\right),$$

with background diffusivity $k_d$, convective strengths $k_i > k_d$ ordered as $0 < k_d < k_E < k_N$, and smoothing parameter $\varepsilon = 0.02$. Two timescale-motivated assumptions reduce the system to four dimensions: deep salinity $\eta$0 is held constant, and tropical temperature relaxes instantly to its atmospheric value $\eta$1. One unit of rescaled time corresponds to roughly 3.97 days. The primary bifurcation parameters are the dimensionless freshwater forcing $\eta$2 (negative values denote net freshwater input to the North Atlantic) and the density threshold $\eta$3, so that convection triggers before static instability. The authors note explicitly that both parameters are phenomenological aggregates of unresolved processes, and that the bifurcation structure shifts significantly as $\eta$4; they therefore explore a broader $\eta$5 range than direct physical interpretation would suggest.

Equilibrium structure

One-parameter continuations in $\eta$6 at fixed $\eta$7 reveal a consistent qualitative pattern: a single stable northward-overturning equilibrium ($\eta$8) exists for $\eta$9 near zero, while only a stable southward-overturning equilibrium ($\Psi = \sigma(\rho_N - \rho_E)$0, nearly zero in magnitude) persists for sufficiently negative $\Psi = \sigma(\rho_N - \rho_E)$1. Between these regimes, sequences of saddle–node (SN) and Hopf (H) bifurcations create intervals of bistability, tristability, and—for $\Psi = \sigma(\rho_N - \rho_E)$2—a narrow interval with up to four coexisting stable equilibria. This multistability exceeds what is found in classical two-box models and implies that freshwater forcing can produce overlapping basins of attraction among distinct steady overturning states.

The two-parameter analysis in the $\Psi = \sigma(\rho_N - \rho_E)$3-plane organizes these regions via curves SN and H, with codimension-two points—generalized Hopf (GH), Bogdanov–Takens (BT), and zero–Hopf (ZH)—acting as organizing centers. In particular, a ZH point of type IV governs a region where SN and H curves meet tangentially, carving out thin strips of additional stable equilibria. The overall structure is consistent with the salt–advection feedback of classical box models and with hosing experiments in OGCMs, which likewise predict weakened and multiple AMOC states under enhanced freshwater forcing.

Periodic Welander oscillations

Welander oscillations are self-sustained oscillations generated by the convective-adjustment scheme, in which North Atlantic mixing alternates between convective and diffusive phases. The authors classify mixing regimes geometrically using the curvature loci $\Psi = \sigma(\rho_N - \rho_E)$4 and $\Psi = \sigma(\rho_N - \rho_E)$5 of the surface defined by $\Psi = \sigma(\rho_N - \rho_E)$6 over the density-contrast plane; a convective shutdown event is an excursion from the convective sheet through the intermediate regime into the diffusive sheet.

Families of periodic solutions bifurcate from Hopf points $\Psi = \sigma(\rho_N - \rho_E)$7, $\Psi = \sigma(\rho_N - \rho_E)$8, and $\Psi = \sigma(\rho_N - \rho_E)$9 on equilibrium branches. Only those born at $\Psi = 0$0 and $\Psi = 0$1 admit stable segments; the $\Psi = 0$2 family remains unstable throughout the plotted range. Stability alternates at successive saddle–node (S) and period-doubling (PD) bifurcations, producing progressively narrower stable intervals as $\Psi = 0$3 decreases. In the two-parameter plane, five yellow-shaded regions of stable periodic solutions are bounded by S and PD curves; these regions shrink and accumulate toward lower $\Psi = 0$4, organized by homoclinic tangencies associated with a codimension-two Shilnikov–Hopf point at $\Psi = 0$5. This accumulation suggests a route to chaos, which the authors confirm numerically.

Representative oscillations exhibit a clear slow–fast structure: a millennial-scale preconditioning phase during which subpolar salinity, tropical salinity, and North Atlantic temperature decrease until the density contrast crosses the threshold $\Psi = 0$6, followed by fast shutdown–recovery episodes on decadal-to-centennial timescales. As $\Psi = 0$7 becomes more negative, successive yellow regions correspond to one additional shutdown event per period, from one up to four events; a single episode occupies about 5% of the period, while four episodes occupy roughly 25%. A key mechanism is the incomplete recovery of deep-water temperature after each shutdown: it acts as a slowly evolving memory variable that leaves the density contrast close to $\Psi = 0$8, promoting clustered subsequent shutdowns. This sawtooth pattern of slow preconditioning and rapid convective collapse is reminiscent of Dansgaard–Oeschger variability, though the model omits sea-ice feedbacks, teleconnections, and interhemispheric coupling.

Chaotic Welander oscillations

For certain parameter values, shutdown events occur irregularly along chaotic attractors. At $\Psi = 0$9, trajectories spiral onto a thin sheet from which they are intermittently ejected during shutdown events, producing a "barcode" time series of $\Psi\tanh(\Psi/\vartheta)$0 with irregular clusters of one, two, or three closely spaced events; larger local maxima of $\Psi\tanh(\Psi/\vartheta)$1 precede larger clusters.

A map of the maximal Lyapunov exponent $\Psi\tanh(\Psi/\vartheta)$2 over an $\Psi\tanh(\Psi/\vartheta)$3 grid identifies a substantial chaotic set in the $\Psi\tanh(\Psi/\vartheta)$4-plane, interspersed with periodic windows. Notably, this set overlaps parts of the second and third periodic-orbit regions, so that periodic and chaotic Welander oscillations coexist for some parameters—an implication being that deterministic freshwater forcing alone can select between regular and irregular convective shutdown without any change in forcing statistics. The chaotic set is bounded in part by two boundary-crisis curves $\Psi\tanh(\Psi/\vartheta)$5 and $\Psi\tanh(\Psi/\vartheta)$6, computed via continuation of homoclinic tangencies using Lin's method, meeting at a double-boundary-crisis point DBC. The remaining lower edge of the chaotic set appears organized by further boundary crises, which the authors do not resolve.

Limitations and open questions

The paper concedes several limitations. First, the non-dimensional parameter region corresponds to a physically large dimensional range, $\Psi\tanh(\Psi/\vartheta)$7, owing to basin-scale box volumes and the choice of $\Psi\tanh(\Psi/\vartheta)$8; precise bifurcation curve locations are sensitive to the smoothing parameter, so the results identify qualitative structure rather than physically definitive thresholds. Second, the reduction assumes constant deep salinity and instantaneous tropical temperature relaxation—assumptions whose relaxation could alter the dynamics. Third, the boundary crises bounding part of the chaotic set remain incompletely characterized, requiring identification of additional saddle periodic orbits and continuation of connecting orbits. Open questions raised by the paper include how subdividing the subpolar box into distinct convection sites (e.g., Irminger and Labrador seas) modifies the shutdown dynamics, how interhemispheric coupling affects the oscillation structure, and how stochastic atmospheric forcing induces transitions between the coexisting periodic and chaotic regimes.

Conclusion

This work provides a comprehensive bifurcation-theoretic account of a four-dimensional conceptual AMOC model combining Stommel advection with Welander convective adjustment. Its principal findings are: an equilibrium structure with up to four coexisting stable overturning states; families of periodic Welander oscillations whose number of convective shutdown events per cycle increases with freshwater forcing; a memory mechanism based on incomplete deep-temperature recovery that clusters shutdown events; and a substantial deterministic chaotic regime featuring irregular on–off switching of convection, partially bounded by computable boundary-crisis curves. Within its idealized setting, the model demonstrates that both regular and irregular convective shutdown dynamics arise in a purely deterministic framework, providing a mechanistic template against which higher-complexity model behavior can be interpreted.