# Effects of the Asymmetry between Surface and Interior Flow on the Dynamics of a Thermohaline Loop

Friederike Pollmann, Fabien Roquet, Gurvan Madec

► To cite this version:

Friederike Pollmann, Fabien Roquet, Gurvan Madec. Effects of the Asymmetry between Surface and Interior Flow on the Dynamics of a Thermohaline Loop. Journal of Physical Oceanography, 2015, 45 (10), pp.2544-2563. 10.1175/JPO-D-15-0022.1. hal-01234077

HAL Id: hal-01234077

https://hal.science/hal-01234077

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Effects of the Asymmetry between Surface and Interior Flow on the
Dynamics of a Thermohaline Loop

FRIEDERIKE POLLMANN

Institut für Meereskunde, Universität Hamburg, Hamburg, Germany

FABIEN ROQUET

Department of Meteorology (MISU), Stockholm University, Stockholm, Sweden

GURVAN MADEC

LOCEAN Sorbonne Universités (UPMC, University of Paris 6)/CNRS/IRD/MNHN, Paris, France

(Manuscript received 31 January 2015, in final form 23 June 2015)

ABSTRACT

Large-scale overturning cells in the ocean typically combine an essentially horizontal surface branch and an interior branch below, where the circulation spans both horizontal and vertical scales. The aim of this study is to analyze the impact of this asymmetry between the two branches by “folding” a one-dimensional thermohaline loop, such that its lower part remains vertical while its upper part is folded down into the horizontal plane. It is found that both the transitory response and the distribution of thermohaline properties are modified significantly when the loop is folded. In some cases, velocity oscillations are induced during the spinup that were not seen in the unfolded case. This is because a circular loop allows for compensations between the density torques produced above and below the heat forcing level, while such compensations are not possible in the folded loop because of the horizontal direction of the surface circulation. Furthermore, the dynamical effects associated with nonlinearities of the equation of state are significantly altered by the folding. Cabbeling tends to decelerate the flow in the folded loop, instead of accelerating it as in the circular case, and can also act to dampen velocity oscillations. Thermobaricity also alters the loop circulation, although comparatively less.

**1. Introduction**

Overturning cells of various shapes and sizes are a pervasive feature of the global ocean, ranging from the small eddy-induced circulations to the global-scale meridional overturning circulation (MOC). Overturning circulations are often analyzed in simplified frameworks like the Stommel box model, which consists of two connected, well-mixed reservoirs of different temperature and salinity and features different circulation modes depending on the relative strength of the thermal and haline forcings (Stommel 1961).

The thermohaline loop is a variation of the Stommel box model in which properties can vary continuously along a one-dimensional loop, allowing for a more realistic representation of the advective-diffusive balance (e.g., Welander 1967; Huang 1999; Wunsch 2005). It does not only find application in physical oceanography but in the form of the so-called thermosyphons also in nuclear and solar energy engineering (e.g., Zvirin 1982; Miljkovic and Wang 2011). Despite their simple setup, thermohaline loops exhibit a variety of dynamical behaviors including instability and chaos (Malkus 1972; Welander 1986; Yuan and Wunsch 2005).

Although it is unclear to what extent analogies between these conceptual representations and the real ocean can be drawn, it has been argued that by capturing the basic physics these idealized models can indeed provide an intuitive understanding of more complex, higher-dimensional processes (Bryan 1986; Ruddick and Zhang 1996; Manabe and Stouffer 1999; Held 2005).

Denotes Open Access content.

Corresponding author address: Friederike Pollmann, Institut für Meereskunde, Universität Hamburg, Bundesstr. 53, 20146 Hamburg, Germany.
E-mail: friederike.pollmann@uni-hamburg.de

DOI: 10.1175/JPO-D-15-0022.1
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Schematic description of the fluid loop as proposed by Wunsch (2005) and (b) the modified folded version. The loop coordinate is represented as $l = a\Phi$, where the angle $\Phi$ is measured clockwise from the top and $a$ denotes the radius. Sources and sinks of temperature and salinity are applied at $l_+$ and $l_-$, respectively, which are located at the same height $Z_f$. The azimuthal velocity $w$ is constant around the loop for reasons of continuity. In the modified version, the loop is folded at the level of source and sink until the part above is purely horizontal. The steady-state temperature distribution is also shown; blue and red colors denote cold and warm fluid, respectively, with dashed lines representing the strong temperature gradients within the diffusive boundary layers.</figcaption></figure>

Wunsch (2005) devised a circular loop model with point sources and sinks of salinity and temperature applied on either side of the loop, typically at an angle of 45° from the loop's top. The case of applying sources and sinks at the same height is of particular interest for ocean studies, as it represents the case of horizontal convection best (Hughes and Griffiths 2008). This means that in the Wunsch model part of the circulation occurs at a level above the position of sources and sinks. In the real ocean, however, the thermohaline forcing (leaving aside geothermal heating) is applied directly at the ocean surface, introducing an asymmetry between an essentially horizontal surface flow and a deeper, both horizontal and vertical, one in the interior. We therefore propose a simple modification of the Wunsch model where the upper portion of the loop above the level of source and sink is folded down until it becomes purely horizontal (cf. Fig. 1).<fnref n="1" />

In this study, we investigate the dynamical modifications introduced when folding the loop, considering a simplified nonlinear EOS similar to the one in Vallis (2006, p. 34). In most idealized studies, including the one by Wunsch (2005), a linear EOS is assumed. Yet, in the real ocean, the EOS is nonlinear mainly through a strong dependence of the seawater thermal expansion on temperature and pressure, giving rise to cabbeling and thermobaricity, respectively. Both effects induce vertical advection through neutral surfaces (albeit only downward in the case of cabbeling) and have been shown to considerably affect the circulation and water mass properties (e.g., McDougall 1987; Klocker and McDougall 2010; Schanze and Schmitt 2013; Hieronymus and Nycander 2013). One key aspect of this study is to explore how the impact of these two nonlinear effects changes when the loop is folded.

This paper is organized as follows: In section 2, we describe the loop model, discuss it analytically, and present the numerical implementation. We then describe the dynamics of the linear EOS scenario (section 3) before addressing the nonlinear case in section 4. An extensive discussion of the results and their applicability is given in section 5, and a brief summary is given in section 6.

<footnote n="1">¹ Note that also the lowest portion of the loop could be folded into the horizontal plane since many overturning cells feature both a horizontal surface and a horizontal bottom branch. However, as the bulk of the buoyancy forcing is confined to the ocean's surface, little dynamical effect is expected from a folded lower branch, which is why we restrict ourselves to folding the loop's upper part only.</footnote>
## 2. The loop model

### a. Governing equations

The model consists of a closed fluid loop with an infinitesimal section area, so that radial velocities can be neglected. The loop natural coordinate $l$ is defined as the distance along the circle from the top of the loop, increasing clockwise. The total length of the loop is given by $L = 2\pi a$, where $a$ denotes the radius. The focus of this study is on heat-induced horizontal convection, where heating and cooling occur at the same height and no salinity forcing is applied. Two geometries of the loop will be considered (cf. Fig. 1):

1) The circular loop: Essentially the same model as in Wunsch (2005). The loop coordinate is $l = a\Phi$, where $\Phi$ is simply the angle from the top, which is defined positive for clockwise turning. Height is thus given as $z(l) = a\cos(l/a)$. Sinks and sources of temperature and salinity are applied at positions $l_− = a\Phi_− \le L/2$ and $l_+ = a\Phi_+ = L - l_−$, respectively, that is, at the same height $Z_f$.

2) The folded loop: The part above source and sink ($z \ge Z_f$ in the circular loop) is now folded down horizontally. The loop coordinate is kept unchanged, only the height function is modified, with $z(l) = \min[Z_f, a\cos(l/a)]$. By analogy, an equivalent angle $\Phi = l/a$ can still be defined, although it must be kept in mind that it is now associated with two different circular paths, one vertical for the interior flow and one horizontal for the surface flow.

Loop equations are derived based on the seawater Boussinesq approximation (e.g., Young 2010; Roquet 2013), which differs from the standard Boussinesq approximation in its ability to use a nonlinear EOS. The continuity equation simplifies to the constraint that the loop velocity $w(l) = Dl/Dt$ is constant around the loop:

$$ \frac{\partial w}{\partial l} = 0, \qquad (1) $$

with an immediate consequence that the momentum advection term is always null in this model.

The momentum equation is given by

$$ \frac{\partial w}{\partial t} = -\frac{1}{\rho_o} \frac{\partial p}{\partial l} + g\sigma\Pi(l) - \varepsilon w + \tau, \qquad (2) $$

where $\sigma = (\rho - \rho_o)/\rho_o$ represents the density anomaly with respect to a reference value $\rho_o$, $g$ is the gravitational acceleration, and $p$ is the pressure. The quantity $\Pi(l)$ is the so-called curvature term defined as

1) $\Pi(l) = \sin(l/a)$ for the circular loop, and

2) $\Pi(l) = m(l)\sin(l/a)$ for the folded loop,

with $m(l)$ denoting a mask function, equal to 1 in the lower branch of the loop (i.e., where $z \le Z_f$) and 0 elsewhere (i.e., in the horizontal, folded part). The parameter $\varepsilon$ describes a constant Rayleigh friction coefficient, and $\tau$ is an applied stress averaged along the loop, representing the integrated effect of the wind.

Since velocity is constant, integrating Eq. (2) around the loop produces a simplified form without the pressure gradient term:

$$ \frac{\partial w}{\partial t} = \overline{g\sigma\Pi} - \varepsilon w + \tau, \qquad (3) $$

where the overbar represents the loop mean: $\overline{T} = (1/L) \int T dl$. Whereas in the circular loop each point around the loop contributes to the buoyancy term $\overline{g\sigma\Pi}$, the surface branch of the folded loop becomes irrelevant to dynamics as it is completely horizontal.

Incidentally, the loop model's pressure is a diagnostic quantity and follows a nonhydrostatic balance:

$$ \frac{\partial p}{\partial z} = -\rho_o g \left( \sigma - m \frac{\overline{\sigma\Pi}}{\sin(l/a)} \right), \qquad (4) $$

although the nonhydrostatic contribution is typically negligible. Indeed, the ratio of hydrostatic to non-hydrostatic vertical pressure variations, given by $\sigma\Pi/\overline{\sigma\Pi}$, is on the order of $10^{11}$ for a typical length scale of $10^7$ m, a time scale of 600 yr, and a scale of $10^{-3}$ for the density anomaly $\sigma$ (a discussion of suitable scaling parameters for the ocean is provided in section 5 and shown below in Table 3).

The EOS describes the density anomaly $\sigma$ as a function of temperature $\Theta$, salinity $S$, and, in the seawater Boussinesq approximation, height $z$:

$$ \sigma = \tilde{\sigma}(\Theta, S, z). \qquad (5) $$

The temperature variable $\Theta$ represents Conservative Temperature (McDougall 2003) and the salinity variable $S$ is Absolute Salinity (McDougall et al. 2009). Trends in buoyancy are given by

$$ \frac{D\tilde{\sigma}}{Dt} = -\alpha \frac{D\Theta}{Dt} + \beta \frac{DS}{Dt} - \gamma \frac{Dz}{Dt}, \qquad (6) $$

where $D/DT = \partial/\partial t + w\partial/\partial l$ is the Lagrangian derivative. The parameters $\alpha$, $\beta$, and $\gamma$ describe the coefficients of thermal expansion, haline contraction, and compressibility, respectively. These are defined as

$$ \alpha = -\left.\frac{\partial \tilde{\sigma}}{\partial \Theta}\right|_{S,z}, \quad \beta = \left.\frac{\partial \tilde{\sigma}}{\partial S}\right|_{\Theta,z} \quad \text{and} \quad \gamma = -\left.\frac{\partial \tilde{\sigma}}{\partial z}\right|_{\Theta,S}. \qquad (7) $$
Here, we assume the following nonlinear EOS based on the formulation by Vallis (2006, p. 34), which captures the main nonlinear EOS effects:

$$ \tilde{\sigma} = -\alpha_o \left[ 1 + \frac{\lambda}{2}(\Theta - \Theta_o) - \mu z \right] (\Theta - \Theta_o) + \beta_o (S - S_o), \quad (8) $$

where $\Theta_o$, $S_o$, $\alpha_o$, $\beta_o$, $\lambda$, and $\mu$ are constant values. The cabbeling effect involves a quadratic term in temperature, whose strength is set by the parameter $\lambda$. The thermobaric effect, on the other hand, is related to the dependence of the thermal expansion coefficient on pressure, here represented as a depth-temperature product term of magnitude $\mu$. Setting $\lambda = \mu = 0$ consequently produces a linear EOS.

The tracer equations are given by

$$ \frac{D\Theta}{Dt} = \kappa \frac{\partial^2\Theta}{\partial l^2} + \eta_\Theta \Delta_f, \quad \text{and} \quad (9) $$

$$ \frac{DS}{Dt} = \kappa \frac{\partial^2 S}{\partial l^2} + \eta_S \Delta_f, \quad (10) $$

where $\kappa$ represents eddy diffusivity assumed constant around the loop and equal for salinity and temperature. The terms $\eta_\Theta$ and $\eta_S$ control the magnitude of the heat and salt forcing, respectively, whose distribution is given by

$$ \Delta_f = 2\pi[\delta(l - l_+) - \delta(l - l_-)], \quad (11) $$

where $\delta$ denotes the Dirac delta function. Only fixed flux forcings will be considered in this study; the case of relaxation forcing will be addressed in future work.

## b. Nondimensionalization

Equations are nondimensionalized assuming that the magnitudes of density anomalies are scaled by the thermal source strength. Scaling constants are denoted by a superscript $s$ and the nondimensional quantities are denoted by a prime with $T = T_o + T' T^s$, where reference values are equal to zero except in the case of salinity and temperature. The scaling constants are set to

$$ l^s = a, \quad w^s = \frac{\sigma^s g}{\varepsilon} \Rightarrow t^s = \frac{l^s}{w^s} = \frac{a\varepsilon}{g\sigma^s}. \quad (12) $$

From the equation of density [Eq. (6)], the scaling constants for temperature and salinity can be obtained:

$$ \Theta^s = \frac{\sigma^s}{\alpha_o}, \quad S^s = \frac{\sigma^s}{\beta_o}. \quad (13) $$

As a result, the following scaling expression for density is found:

$$ \eta_{\Theta}^{s} = \frac{\Theta^{s}}{t^{s}} \Rightarrow \sigma^{s} = \sqrt{\frac{a\varepsilon\alpha_{o}\eta_{\Theta}^{s}}{g}}. \quad (14) $$

The nondimensional equations of motion then take the form

$$ F \frac{\partial w'}{\partial t} = \overline{\sigma' \Pi'} - w' + \tau', \quad (15) $$

$$ \sigma' = - \left( 1 + \frac{\lambda'}{2} \Theta' - \mu' z' \right) \Theta' + S', \quad (16) $$

$$ \frac{D\Theta'}{Dt'} = R \frac{\partial^2\Theta'}{\partial l'^2} + \Delta_f, \quad \text{and} \quad (17) $$

$$ \frac{DS'}{Dt'} = R \frac{\partial^2 S'}{\partial l'^2} + \eta \Delta_f, \quad (18) $$

where $D/Dt' = \partial/\partial t' + w'/\partial l'$.

Apart from geometrical parameters (i.e., circular vs folded cases and the vertical position of forcings $Z_f$), several nondimensional parameters are defined in the loop model:

1) the inverse Rayleigh number $R = \kappa \sqrt{\varepsilon/(ga^3\alpha_0\eta_\Theta^s)}$, which corresponds to the ratio of diffusive to advective time scales and controls the size of the diffusive boundary layers near the source and sink (as will be shown later);

2) the momentum stress forcing $\tau'$ as a representation of frictional forcing;

3) the Grashof number $F = \sqrt{(g\alpha_0\eta_\Theta^s)/(a\varepsilon^3)}$, which represents the ratio of frictional to advective time scales. Assuming inertialess flow, that is, $F \to 0$, is equivalent to neglecting frictional oscillations. Equation (15) then simplifies to

$$ w' = \overline{\sigma' \Pi'} + \tau'; \quad (19) $$

4) the salt forcing strength parameter $\eta = (\beta_o \eta_S^s)/(\alpha_o \eta_\Theta^s)$, which scales the magnitude of the salinity forcing relative to that of the temperature forcing; when it is set to one, the temperature and salinity forcings cancel exactly when a linear EOS is assumed. Note that a circulation can be generated even if $\eta = 1$ if the EOS is nonlinear. No salt forcing will be considered in this paper, that is, $\eta$ will be kept to zero; and

5) finally, the cabbeling and thermobaricity parameters $\lambda'$ and $\mu'$.

The only formal difference between the circular and folded loop model concerns the buoyancy term, which is masked in the folded loop case. Another, somewhat minor difference is that the definition of height is changed with consequences for the thermobaric term in the EOS.
In the rest of the paper, we will refer to non-dimensionalized quantities only, and the primes will be dropped for clarity.

### c. The numerical model

Following the approach by Wunsch [2005, cf. his Eqs. (11)–(14)], the system of equations—Eqs. (16)–(19), using the inertialess assumption—can be solved analytically for the equilibrium velocity $w$ if a linear EOS is assumed. This procedure can be extended to nonlinear scenarios including only thermobaricity; when adding the cabbeling term or folding the loop, however, the ability to compute a simple analytical solution is lost. Therefore, a numerical investigation of loop dynamics is required.

The loop model is implemented in Fortran90, using a forward differencing scheme for the diffusive part and a leapfrog scheme for the advection and forcing terms. With a time increment of $\Delta t = 10^{-4}$ and 360 grid points around the loop (i.e., $\Delta\Phi = 1^\circ$), the diffusion scheme is stable for $R < 0.4$. Since the velocity overshoot is typically below $w = 2$, the Courant–Friedrich–Levy criterion $w \le \Delta\Phi/\Delta t \approx 175$ is easily fulfilled.

Our standard parameter settings are $R = 0.1$ (weak diffusion), $\tau = 0$ (no wind stress forcing), and $\eta = 0$ (no salinity forcing). The simulations are started at rest ($w = 0$) with temperatures set to zero everywhere. The thermal forcing, applied at the same height (i.e., $\Phi_+ + \Phi_- = 2\pi$), is switched on at $t = 0$. Our main focus is on the scenario $Z_f = 0.5$ (i.e., $\Phi_- = \pi/3$; cf. Fig. 1). With these settings, equilibrium is typically reached at $t \approx 15$ in the circular loop model and at $t \approx 20$ in the folded loop model. Unless otherwise noted, the simulations involve $10^6$ iterations (i.e., $t_{end} = 100$), which take less than a minute to perform on a standard computer.

The numerical model is validated based on scenarios that allow for the computation of an analytical solution for $w$. We find that for our standard parameter settings, the difference between numerical and analytical solutions is several orders of magnitude smaller than the respective analytical velocity. The remaining inaccuracies are therefore considered too small to influence the conclusions drawn here.

## 3. Loop dynamics with a linear EOS

Considering a linear EOS, density<fnref n="2" /> is given by $\sigma = -\Theta$, and velocity is given by $w = -\Theta\Pi$ (again, we

focus on a thermally driven circulation, i.e., $\eta = \tau = 0$). The curvature term $\Pi$ weights various parts of the loop differently as to their influence on velocity; changes in density around $z = 0$, where $\sin(\Phi)$ is maximum, have the greatest impact, while events at the loop's top or bottom are irrelevant for velocity. Since the sine function is antisymmetric with respect to the vertical axis, its product with any symmetric function integrates to zero around the loop—a property that is preserved when folding the loop. This has important consequences for the buoyancy term $\sigma\Pi$, as density can be decomposed into the symmetric horizontal mean $\sigma^m$ and the antisymmetric horizontal anomaly $\sigma^a$:

$$ 
\begin{aligned}
\sigma^m &= \frac{\sigma(\Phi) + \sigma(-\Phi)}{2}, \\
\sigma^a &= \frac{\sigma(\Phi) - \sigma(-\Phi)}{2},
\end{aligned}
\qquad (20) $$

with $\sigma^m(\Phi) = \sigma^m(-\Phi)$ and $\sigma^a(-\Phi) = -\sigma^a(\Phi)$. As a consequence, only the distribution of the horizontal density anomaly $\sigma^a$ can affect the loop velocity:

$$ \overline{\sigma\Pi} = \overline{\sigma^a\Pi}. \qquad (21) $$

Note that the vertical integral of the density anomaly as a function of height is simply equal to half of the instantaneous velocity:

$$ w = 2 \int_{-1}^{1} \sigma^a(z) dz. \qquad (22) $$

This means that positive values tend to accelerate the flow (and vice versa), and compensations can arise when both negative and positive density anomalies are observed simultaneously.

### a. The transitory regime

The simulations are started from a state of rest with temperatures set to zero everywhere. When the thermal forcing is turned on at $t = 0$, temperature anomalies are generated at the source and sink (cf. Fig. 2) that induce horizontal density anomalies and thus a nonzero velocity. This will in turn advect the temperature anomalies away from the forcing region and further modify velocity. Because of the influence of the curvature term, the highest velocities can be obtained when source and sink are applied at $Z_f = 0$, where $\sin(\Phi) = 1$. The circulation reaches its maximum speed when the horizontal density anomaly, integrated around the loop, is largest in accordance with Eq. (19). For example, when source and sink are applied at $Z_f = 0$ in the circular loop model, this is given when the warm front reaches the top and the cold front reaches the bottom of the loop. As velocity

<footnote n="2">² Note that from here on we refer to the nondimensional density and mass anomalies simply as density and mass for the sake of brevity. This is the reason why “negative densities” are sometimes obtained.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Initial evolution of temperature as a function of time and position for the circular and the folded loop model. (a), (c) When source and sink are applied at $Z_f = 0$, the temperature amplitudes are the same at source and sink. (b), (d) Shifting the source and sink up or down, for example, with $Z_f = 0.5$, leads to an asymmetric temperature distribution. (e), (f) When the loop is folded, the temperature distribution is similar but exhibits higher amplitudes because of the significant velocity decrease, which is underlined by the temperature difference $\Delta\Theta = \Theta(\text{full}) - \Theta(\text{folded})$. Contour lines are placed at intervals of 1 in (a)–(d) or 0.5 in (e) and (f).</figcaption></figure>

decreases, the fluid is advected more slowly past the source and sink, which increases their efficiency and leads to higher temperature amplitudes (and vice versa). There are further oscillations in temperature (cf. Figs. 2 and 3) that are associated with the arrival of warm (cold) temperature anomalies at the sink (source) and whose amplitudes are steadily decreased by diffusion. Yet, those temperature oscillations do not yield velocity

oscillations because they tend to compensate almost exactly upon integration around the loop.

When the loop is folded, significant differences in magnitude are observed. The velocity is generally reduced as expected, and it takes longer for the system to reach a steady state (cf. Fig. 3). Also, the temporal variations of velocity and temperature are modified substantially. For example, in the circular loop the
<figure><img src="image_4.png" /><figcaption>FIG. 3. Variation of velocity and temperature at the source and sink with time in the linear EOS scenario for the (a) circular loop and the (b) folded loop model. The highest equilibrium velocities are obtained when source and sink are applied at $Z_f = 0$, where the curvature term $\Pi$ is maximum. The corresponding temperature variation demonstrates that faster advection results in lower temperature amplitudes. Significant changes in the variability of the transitory regime are observed when the loop is folded down.</figcaption></figure>

maximum velocity overshoot is obtained for the case $Z_f = 0$ (Fig. 3a, red curve), while in the folded loop the velocity maximum is highest for the setting $Z_f = 0.5$ (Fig. 3b, green curve). Moreover, in the latter case oscillations in velocity are observed that do not exist in the circular loop case. This is because once the loop is folded, buoyancy anomalies produced in the upper part of the loop can no longer compensate anomalies in its lower part (see Fig. 3).

b. Steady-state properties

The steady state is reached when a diffusive-advection temperature balance is achieved everywhere around the loop. Then, temperatures are negative in most of the lower loop between the sink and source and positive in most of the upper loop between source and sink (cf.

Figs. 1 and 2). Because of diffusion, however, the temperature changes sign slightly above the sink and slightly below the source, which gives rise to a small region around $Z_f$ where $\Theta$ and thus $\sigma$ change sign horizontally. This region is pivotal for maintaining the circulation in equilibrium; if the temperature and hence the density distribution were exactly symmetrical around the vertical axis, the horizontal density anomaly would be zero everywhere and there could be no circulation.

Simple scaling laws can be derived for the size of these diffusive boundary layers and the observed amplitude of temperature anomalies. Except at the point of source and sink, the source balance in equilibrium is given by

$$w \frac{\partial \Theta}{\partial l} = R \frac{\partial^2 \Theta}{\partial l^2}, \qquad (23)$$
which describes an exponential variation of temperature above the sink and below the source. The thickness of the boundary layers $\delta_{BL}$ thus scales as $\delta_{BL} \propto R/w$, which implies that their size is increased when diffusion is stronger and decreased when the circulation is faster.

Assuming thin boundary layers allows for another approximation: Based on the constraint of heat conservation, temperature anomalies have to cancel around the loop in equilibrium. When diffusion is weak, the temperature is nearly constant in both the warm and cold branch, and that constraint can be expressed as

$$L^+ \Delta T^+ + L^- \Delta T^- \approx 0, \quad (24)$$

where $L^+$ denotes the length of the warm upper path, $L^-$ is that of the cold path below $Z_f(L^+ + L^- = L)$, and $\Delta T^+$ and $\Delta T^-$ are their respective temperatures. Note that if diffusion is strong (i.e., $R \approx 1$ or more), the boundary layers cannot be neglected and the temperature distribution within each path can no longer be approximated by a constant.

Equation (24) indicates that the relative length of the warm path determines its temperature amplitude, that is, the temperature amplitude at the source, relative to that at the sink. In the symmetric case with $Z_f = 0$, the warm path is as long as the cold path ($L^+ = L^-$), and consequently, $\Delta T^+ = -\Delta T^-$. Positioning source and sink at $Z_f = 0.5$, on the other hand, means that the warm branch is half as long as the cold branch, which requires $\Delta T^+ = -2\Delta T^-$ in equilibrium.

The heat budget of a section of the loop including the source and the associated boundary layer, but not the sink, is given by

$$2\pi + w\Delta T^- - w\Delta T^+ = 0. \quad (25)$$

Combining Eqs. (24) and (25), we obtain

$$\Delta T^+ = \frac{L^-}{w}. \quad (26)$$

This relation is confirmed in numerical experiments for both model versions when the standard parameter settings are applied (cf. Table 1).

The diffusive boundary layer of the loop's right branch is situated above the sink, which in the folded loop model lies within the folded part that does not contribute to the buoyancy torque. In the dynamically relevant part of the folded loop, the section where horizontal temperature anomalies are high is then approximately half as large as in the circular model (cf. Fig. 1). This partly explains the observed reduction in velocity compared to the circular model ($\Delta w = -18.3\%$ for $Z_f = 0.5$ and $\Delta w = -22.4\%$ for $Z_f = 0$). However, the velocity

<table><caption>TABLE 1. Steady-state velocity and temperature values for different positions of source and sink in the circular and the folded loop model.</caption><thead><tr><th>Model</th><th>Z<sub>f</sub></th><th>Φ<sub>-</sub></th><th>w</th><th>Θ(Φ<sub>+</sub>)</th><th>Θ(Φ<sub>-</sub>)</th></tr></thead><tbody><tr><td rowspan="5">Circular loop</td><td>0.87</td><td>π/6</td><td>0.46</td><td>11.34</td><td>-2.29</td></tr><tr><td>0.5</td><td>π/3</td><td>0.55</td><td>7.60</td><td>-3.80</td></tr><tr><td>0</td><td>π/2</td><td>0.58</td><td>5.42</td><td>-5.42</td></tr><tr><td>-0.5</td><td>2/3 π</td><td>0.55</td><td>3.80</td><td>-7.60</td></tr><tr><td>-0.87</td><td>5/6 π</td><td>0.46</td><td>2.29</td><td>-11.34</td></tr><tr><td rowspan="5">Folded loop</td><td>0.87</td><td>π/6</td><td>0.40</td><td>12.81</td><td>-2.61</td></tr><tr><td>0.5</td><td>π/3</td><td>0.45</td><td>9.36</td><td>-4.68</td></tr><tr><td>0</td><td>π/2</td><td>0.45</td><td>6.97</td><td>-6.97</td></tr><tr><td>-0.5</td><td>2/3 π</td><td>0.41</td><td>5.16</td><td>-10.31</td></tr><tr><td>-0.87</td><td>5/6 π</td><td>0.26</td><td>3.95</td><td>-18.28</td></tr></tbody></table>

decrease is somewhat smaller than expected considering the length of loop that is folded (30% and 50%, respectively). This is because the lower velocity induces, in turn, a distinct rise in temperature amplitudes as predicted by Eq. (26), leading to a compensation through slightly increased horizontal density gradients.

### c. Sensitivity to the forcing position $Z_f$

In the circular loop, the effect of vertically shifting the level of source and sink is symmetric about $z = 0$, that is, at any time, the same loop velocity is obtained for $Z_f$ and $-Z_f$ (cf. Table 1). Choosing $Z_f = -0.5$ instead of $Z_f = 0.5$, for example, means that the cold path is half instead of twice as long as the warm path. The relative length of upper and lower paths, however, is not affected; in consequence, the temperature distribution is the same but of opposite sign, and the velocity remains unaffected. This symmetry is broken in the folded loop model; since the loop is always folded at the level of source and sink, the size of the dynamically important part of the loop is varied when source and sink are shifted vertically.

An asymmetric positioning of source and sink (i.e., $Z_f \neq 0$) has important consequences for dynamics: First, the timing of the various dynamical events becomes asynchronous, as the warm front reaches the top or the sink earlier/later than the cold front reaches the bottom or the source (cf. Fig. 2b). Second, because of their different amplitudes, the warm or cold temperature anomalies no longer compensate (the asymmetric temperature distribution begins to be established at $t \approx 0.6$ for $Z_f = 0.5$). This causes oscillations and even sign changes in the horizontal density anomaly (cf. Fig. 4), which approaches its equilibrium distribution smoothly in both model setups when $Z_f = 0$ (not shown).

## 4. Cabbeling and thermobaricity

Introducing cabbeling only, the EOS assumes the form $\sigma = \sigma_{\text{lin}} + \sigma_{\text{cabb}}$. Unlike the linear term $\sigma_{\text{lin}} = -\Theta$,
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time evolution of the horizontal density anomaly $\sigma^a$ as a function of height in the (a) full circular loop and the (b) folded loop model for $Z_f = 0.5$ and $\Phi \in (0, \pi)$. Contour lines are placed at intervals of 1. Note that the vertical integral of $\sigma^a$ is equal to half the loop velocity.</figcaption></figure>

the nonlinear contribution $\sigma_{cabb} = -0.5\lambda\Theta^2$ has the same sign everywhere, leading to a buoyancy gain throughout the dynamically relevant part of the loop. Including only thermobaricity in the EOS leads to the form $\sigma = \sigma_{\text{lin}} + \sigma_{\text{therm}}$, where the nonlinear contribution $\sigma_{\text{therm}} = \mu z \Theta$ changes sign both with height and temperature. Consequently, cabbeling and thermobaricity are expected to induce different dynamical modifications compared to the linear scenario.

### a. Effect on the steady state

As illustrated in Fig. 4, the region around source and sink is pivotal for loop dynamics; the impact of the nonlinear effects at that level will therefore determine how they affect the steady-state velocity of the circular loop. We will first discuss the case of introducing cabbeling only.

At the source, density is negative and reduced even further due to cabbeling (assuming $\lambda > 0$). This leads to an increased tendency to rise and hence a faster circulation (represented by the longer red arrow in the schematic of the circular loop model in Fig. 5a). At the sink, on the other hand, density is positive, so that the nonlinear buoyancy gain due to cabbeling diminishes the fluid's tendency to sink (represented by the shorter blue arrow in Fig. 5a). This causes a deceleration of the flow. With a completely symmetric geometry, where source and sink are applied at $z = 0$, the buoyancy increase by cabbeling is equally strong in the warm and the cold branch. In consequence, the effects at source and sink compensate exactly, and cabbeling does not influence the steady-state velocity in the circular loop

model, which is indeed observed and can even be shown analytically.

For $Z_f \neq 0$, the temperature distribution determines which of the opposing effects of cabbeling at source and sink ultimately dominates; with an asymmetric position of source and sink, temperature amplitudes in the warm branch differ from those in the cold branch. Because cabbeling is a quadratic function of temperature, it has the highest effect at the source and within the warm branch for the standard setting of $Z_f = 0.5$, for which it induces a velocity increase (cf. the schematic illustration in Fig. 5a). This is underlined by the variation of steady-state velocity with $\lambda$ (cf. Fig. 6a): a positive cabbeling parameter leads to an acceleration of the flow, while $\lambda < 0$ decelerates the circulation. In the former case, cabbeling induces a decrease in mass $\bar{\sigma}$, and in the latter, cabbeling induces an increase (cf. Fig. 6a); with a linear EOS, that is, $\sigma = -\Theta$, the total mass is zero for reasons of heat conservation.

This compensation between these opposite effects is also illustrated in the difference in the horizontal density anomaly with respect to the linear scenario [cf. Fig. 7a for the circular loop model with $\lambda = 0.1$ and $\Phi \in (0, \pi)$]: in the upper loop, where temperatures are positive, cabbeling acts to increase the horizontal density anomaly of the right branch, while in the lower right branch, $\Delta\sigma^a$ assumes negative values, indicating a reduced density anomaly due to cabbeling. The higher amplitudes in the warm section then lead to a stronger effect of cabbeling in the upper than in the lower loop, resulting in the observed velocity increase.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Schematic description of the effect of (a) cabbeling and (b) thermobaricity in both model versions with $Z_f = 0.5$. Full arrows represent the linear scenario, while dashed arrows denote the nonlinear contribution. Blue and red colors describe cold and warm fluid, respectively. The effect on steady-state velocity is illustrated in black.</figcaption></figure>

When the loop is folded, cabbeling (and thermobaricity) can influence the steady-state velocity even when $Z_f = 0$, since a perfect cancellation between warm and cold loop parts is no longer possible. For better comparison, however, source and sink are applied at $z = 0.5$ like in the circular loop. Moreover, we again focus on the scenario $\lambda = 0.1$, which allows for linear effects to dominate loop dynamics, while the influence of cabbeling is expected to be sufficiently strong to be discernible. In that case, cabbeling induces a 5% deceleration of the circulation in contrast to the 4.6% increase in velocity in the circular loop model (cf. Table 2 and Fig. 5a for a schematic description of how cabbeling affects dynamics in the folded loop model). The variation of velocity with $\lambda$ (cf. Fig. 6b) underlines that the influence of cabbeling is reversed compared to the circular loop model. This is because the upper, warm branch is no longer relevant for velocity. The horizontal density anomaly difference to the linear EOS case (cf. Fig. 7b) illustrates that in most of the dynamically important part of the loop a decrease in $\sigma^a$ is induced when cabbeling is included in

the EOS. There, temperatures are mostly negative, and the nonlinear buoyancy gain acts to retard the flow, which causes the observed velocity decrease. However, the downward diffusion of heat from the source leads to a small section where temperatures are positive and cabbeling acts to accelerate the flow ($\Delta\sigma^a > 0$); hence even in the folded loop a compensation between the retarding and the accelerating effect of cabbeling is observed.

The velocity decrease due to cabbeling results in higher temperature amplitudes (cf. Table 2), which is also why mass is affected more strongly than in the circular loop (cf. Fig. 6b). In consequence, the amplitude of the horizontal density anomaly difference between the nonlinear and the linear scenario is also higher than in the circular loop (cf. Fig. 7).

In both model setups, an oscillatory instability can be provoked when the cabbeling parameter is sufficiently large or in the circular loop for any choice of $\lambda$ when the linear density forcings by salinity and temperature cancel ($\eta = 1$). Critical values are $\lambda < -0.385$ and $\lambda > 0.951$.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Velocity and mass as a function of $\mu$ and $\lambda$. The linear reference value is shown in red. When the amplitude of the cabbeling parameter $\lambda$ is too large, an oscillatory instability is induced, which is why the full range between $-1 \le \lambda \le 1$ is not depicted for the circular loop model. In the folded loop model, higher values of $\lambda$ lead to ever lower velocities; the different behavior for $\lambda \ge 0.3$ is related to the shift from an advection- to a diffusion-dominated regime. In that regime, the circulation is very slow, which results in very high, negative values for mass.</figcaption></figure>

for $\eta = 0$ in the circular loop model, which might vary slightly when the simulations are run arbitrarily long (here, $t_{end} = 500$). Nevertheless, they suggest that stability is affected differently when the nonlinear and linear contributions to density do not act in the same way ($\lambda < 0$) and when they do ($\lambda > 0$). In the folded loop, different parameter settings induce this instability. A detailed analysis of this phenomenon will be deferred to future work.

Considering only thermobaricity in the EOS, the nonlinear contribution $\sigma_{therm} = \mu z \Theta$ varies both with temperature and height. Throughout most of the loop, $\Theta$ and $z$ have the same sign, so that the thermobaric term is positive for $\mu < 0$. In that case, thermobaricity induces an increase in mass, which is illustrated in Fig. 6a. At the source, where $\Theta > 0$, the thermobaric term is positive if $\mu > 0$ in the standard scenario with $Z_f = 0.5$, which leads to a denser fluid and thus a reduced tendency to rise. At the sink, where $\Theta < 0$, the thermobaric term is negative and acts to increase buoyancy. Both processes lead to a

velocity decrease, which is shown schematically in Fig. 5b. When $\mu < 0$, the opposite effect is obtained, as the variations of mass and velocity with $\mu$ (cf. Fig. 6a) show.

When source and sink are not applied at the loop's equator, the thermobaric term changes sign twice in the circular loop model. With $\mu > 0$ and $Z_f = 0.5$, it causes a deceleration in the upper part of the loop ($z > 0$) and an acceleration for $z < 0$. The overall impact of thermobaricity on steady-state velocity is governed by its effect at source and sink and in that part of the loop with the highest temperature amplitudes; for $\mu = 0.1$ and the standard parameter settings, velocity is decreased by 1.6% in the circular loop model. When the loop is folded, there still is a compensation between the counteracting effects of thermobaricity in the upper and lower part of the loop, but the total impact on steady-state velocity is reduced: the deceleration observed in the circular loop for $z > Z_f = 0.5$ is irrelevant for dynamics in the folded loop model. In consequence, the
<figure><img src="image_8.png" /><figcaption>FIG. 7. Difference in the horizontal density anomaly $\sigma^a$ between the case with $\lambda = 0.1$ and the one based on a linear EOS ($\lambda = 0$, depicted in Fig. 4) for the (a) circular loop and the (b) folded loop with $Z_f = 0.5$. In the circular loop, the largest differences are observed in the upper part, where temperature amplitudes are highest, especially at the level $Z_f = 0.5$. When the loop is folded, the largest differences are found at and slightly below $Z_f = 0.5$, that is, where the diffusive boundary layer of the left branch is located. Contour lines are placed at intervals of 0.4, and again only $\Phi \in (0, \pi)$ is considered in the computation of $\sigma^a$.</figcaption></figure>

accelerating effect of thermobaricity below $z = 0$ becomes more important, leading to a velocity change of $\Delta w = -0.9\%$ compared to the linear case (illustrated in Fig. 5b). Contrary to the cabbeling scenario, the region around $Z_f$ still determines the effect of thermobaricity on equilibrium velocity in the folded loop—the section where $\Delta\sigma^a < 0$ extends from $z = 0.5$ to $z = 0.2$ (cf. Fig. 8b), while for cabbeling $\Delta\sigma^a$ changes sign already at $z = 0.4$ (cf. Fig. 7b).

When source and sink are applied below $z = 0$, the effect of thermobaricity on velocity is reversed in both model versions (cf. the positive horizontal density difference below $z = 0$ shown in Fig. 8): the dynamically important region, where temperatures change sign horizontally (i.e., within the boundary layers around the level of source and sink), as well as that part of the loop where temperature amplitudes are highest are then associated with a negative value of $z$.

The final results for velocity, temperature, and mass enlisted in Table 2 show that in both model geometries the effects of cabbeling and thermobaricity are superimposed almost linearly in a combined scenario with $\lambda = \mu = 0.1$. For example, velocity is increased by 4.6% due to cabbeling in the circular loop, decreased by 1.6% due to thermobaricity, and increased by 3.2% in a combined scenario. In the folded loop, on the other hand, both effects add up with respect to their influence in velocity, leading to a deceleration by 6.0% for $\lambda = \mu = 0.1$. Mass on the other hand is still altered oppositely by the two nonlinear processes, which have a

stronger effect than in the circular loop model because of the higher temperature amplitudes caused by the folding and the resultant velocity decrease.

In the parameter ranges discussed here, cabbeling has a somewhat stronger effect on dynamics and steady-state properties in both model versions. This is because the thermobaric term is only a linear, not a quadratic, function of temperature and because $z$ is maximum where the curvature term $\Pi$ is minimum and vice versa. Yet, for other parameter choices the influence of thermobaricity and cabbeling could be comparable, or their relative importance might even be reversed. Note, too, that with respect to mass both nonlinear effects have a significant influence—when a linear EOS is used, mass is zero in equilibrium for reasons of heat conservation.

<table><caption>TABLE 2. Final results for velocity, mass, temperature, and density at the source in the folded and the circular loop model for various forms of the EOS.</caption><thead><tr><th>Model</th><th>λ</th><th>μ</th><th>w</th><th>Mass</th><th>Θ(Φ<sub>+</sub>)</th><th>σ(Φ<sub>+</sub>)</th></tr></thead><tbody><tr><td rowspan="4">Circular loop</td><td>0</td><td>0</td><td>0.55</td><td>0</td><td>7.60</td><td>-7.60</td></tr><tr><td>0.1</td><td>0</td><td>0.58</td><td>-1.15</td><td>7.26</td><td>-9.90</td></tr><tr><td>0</td><td>0.1</td><td>0.54</td><td>0.31</td><td>7.72</td><td>-7.33</td></tr><tr><td>0.1</td><td>0.1</td><td>0.57</td><td>-0.89</td><td>7.36</td><td>-9.70</td></tr><tr><td rowspan="4">Folded loop</td><td>0</td><td>0</td><td>0.45</td><td>0</td><td>9.36</td><td>-9.36</td></tr><tr><td>0.1</td><td>0</td><td>0.43</td><td>-2.02</td><td>9.84</td><td>-14.69</td></tr><tr><td>0</td><td>0.1</td><td>0.44</td><td>0.28</td><td>9.44</td><td>-8.97</td></tr><tr><td>0.1</td><td>0.1</td><td>0.42</td><td>-1.77</td><td>9.95</td><td>-14.41</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>FIG. 8. As in Fig. 7, but varying the thermobaricity parameter $\mu$ instead of the cabbeling parameter $\lambda$. Differences are shown between the case where $\mu = 0.1$ and the linear case $\mu = 0$ (depicted in Fig. 4) and contour lines are placed at intervals of 0.05.</figcaption></figure>

b. Effects on the mass distribution

As stated in section 1, cabbeling is associated with a densification of water masses. In both model setups, however, we observe a decrease in mass when cabbeling is included (cf. Table 2). This apparent contradiction can be explained in the following manner: Adding the cabbeling term to the EOS implies that the thermal expansion coefficient is no longer constant:

$$ \alpha = -\left.\frac{\partial\sigma}{\partial\Theta}\right|_{S,z} = 1 + \lambda\Theta. \qquad (27) $$

Buoyancy trends [cf. Eqs. (6) and (9)] are then given by

$$ \frac{D\sigma}{Dt} = -(1 + \lambda\Theta)\left(R\frac{\partial^2\Theta}{\partial l^2} + \Delta_f\right). \qquad (28) $$

The term related to temperature diffusion can be rewritten as

$$ -\alpha R \frac{\partial^2 \Theta}{\partial l^2} = -\frac{\partial}{\partial l} \left( \alpha R \frac{\partial \Theta}{\partial l} \right) + \lambda \left( \frac{\partial \Theta}{\partial l} \right)^2, \qquad (29) $$

where the first term can be interpreted as a diffusive density flux. The second term only depends on the sign of $\lambda$ and hence constitutes a density source in our standard case of $\lambda = 0.1$. It is proportional to the temperature gradient and thus strongest where mixing is significant (i.e., in the boundary layers). This mixing-related density source is balanced by the density sink associated with the nonlinear density forcing at source and sink (density diffusion and advection integrate to zero in equilibrium). Integrating

the second term in Eq. (28) around the loop gives the total density forcing:

$$ \begin{aligned} \overline{-\alpha\Delta_f} &= -\int_0^{2\pi} (1 + \lambda\Theta)[\delta(l - l_+) - \delta(l - l_-)] dl \\ &= -\lambda[\Theta(l_+) - \Theta(l_-)]. \end{aligned} \qquad (30) $$

This shows that the linear contribution to the density forcing integrates to zero, while the nonlinear one due to cabbeling is negative for $\lambda > 0$ and thus constitutes a net sink of density. The variation of mass with time (cf. Fig. 9) underlines that this buoyancy source is the dominant mechanism during the first time steps, when the final temperature distribution with the strong gradients in the diffusive boundary layers is not yet established and the mixing-induced buoyancy sink is weaker. That nonlinear buoyancy source also determines the final mass of the fluid since the densification due to cabbeling is confined to a small section of the loop. The same holds true for the folded loop, with the only difference being that mass oscillates more strongly and is affected to a larger degree because of the higher temperature amplitudes and the asymmetry introduced by the folding.

With a nonlinear EOS including only thermobaricity, where $\alpha = 1 - \mu \min[Z_f, \cos(l)]$ and $\gamma = -\mu\Theta$ [cf. Eq. (7)], trends in buoyancy [cf. Eqs. (6) and (9)] are given by

$$ \frac{D\sigma}{Dt} = -\alpha \left( R \frac{\partial^2\Theta}{\partial l^2} + \Delta_f \right) + \Pi w \gamma. \qquad (31) $$

In the same manner as for cabbeling, the temperature diffusion term can be rewritten as a diffusive density flux and a nonlinear source or sink of density:
<figure><img src="image_10.png" /><figcaption>FIG. 9. Variation of velocity, temperature, and mass with time in the (a) circular loop and the (b) folded loop model for different versions of the EOS with $Z_f = 0.5$. For clarity, we show the scenario $\mu = 0.5$; the behavior is qualitatively the same as in the case with $\mu = 0.1$, which can hardly be distinguished from the scenario with a linear EOS. The amplitudes of the steady-state velocity (temperature) are decreased (increased) by less than 10%, and the final mass is increased by a factor of 5 when considering $\mu = 0.5$ instead of $\mu = 0.1$.</figcaption></figure>

$$ -\alpha R \frac{\partial^2 \Theta}{\partial l^2} = -\frac{\partial}{\partial l} \left( \alpha R \frac{\partial \Theta}{\partial l} \right) + \mu \Pi R \frac{\partial \Theta}{\partial l}. \quad (32) $$

The second term depends not only on the sign of the thermobaric parameter but also on the temperature gradient and the sine function. It can thus constitute either a source or sink of density, which is also why both upward and downward motion through neutral surfaces are associated with thermobaricity as described in the introduction.

The nonlinear density forcing due to thermobaricity, however, vanishes when source and sink are applied at the same height:

$$ \begin{aligned} \overline{-\alpha \Delta_f} &= - \int_0^{2\pi} \alpha[\delta(l - l_+) - \delta(l - l_-)] dl \\ &= \mu[\cos(l_+) - \cos(l_-)] = 0. \end{aligned} \quad (33) $$

Consequently, the mixing-related buoyancy source or sink described in Eq. (32) is in equilibrium balanced only by the compressibility term [the last term in Eq. (31)]. Through most of the loop, the fluid is denser than in the linear case when $\mu = 0.1$, concordant with the mass increase observed for $t > 0$ (cf. Fig. 9); locally, however, lower density values can be obtained since height and temperature change sign at different levels in the scenario with $Z_f = 0.5$.

### c. The transitory regime

To assess the influence of the nonlinear effects during the transient regime, the difference in the horizontal density anomaly with respect to the linear reference case is analyzed. As in the previous subsection, we first analyze the scenario with $\lambda = 0.1$ (cf. Fig. 7 for the difference
in $\sigma^a$ using $Z_f = 0.5)$. In both model versions, oscillations between positive and negative $\Delta\sigma^a$ can be observed for $t < 10$, indicating that cabbeling intensifies or weakens certain dynamical events. For example, when the cold front propagates down the right branch right after the onset of the forcing, the buoyancy gain due to cabbeling leads to a lighter right side and hence a reduced horizontal density difference compared to the linear case. When the warm front arrives at the sink and travels down the right branch for $t > 1.1$, the right side becomes lighter; this is intensified by cabbeling so that the horizontal density anomaly is higher than when a linear EOS is considered. These modifications of the horizontal density anomaly are the main reason for the different variation of velocity with time shown in Fig. 9. The magnitude of the changes brought about by cabbeling differs in the full and folded loops because the latter exhibits higher temperature amplitudes as a consequence of the folding.

This effect of cabbeling on buoyancy has interesting consequences in the folded loop model, where velocity oscillations can be observed in the linear case when an asymmetric positioning of source and sink is chosen (cf. Fig. 9 for the scenario with $Z_f = 0.5)$. These oscillations, however, are dampened when cabbeling is accounted for in the EOS: the velocity increase for $t \ge 3.8$ is brought about by the arrival of the warm minimum at the sink, causing anomalously cold and thus dense fluid, which increases the horizontal density anomaly and hence accelerates the flow. In the presence of cabbeling, this densification is counteracted by the nonlinear buoyancy gain. In consequence, $\sigma^a$ is not increased as much as when a linear EOS is considered, and the acceleration is weaker. This is confirmed by the negative density anomaly difference observed in the lower loop for $t \ge 3.8$ (cf. Fig. 7).

Considering now the variation of the horizontal density anomaly difference between the nonlinear scenario with $\mu = 0.1$ and the linear one (cf. Fig. 8), we see that, like cabbeling, the thermobaric effect intensifies and weakens the various dynamical events. Especially during the first time steps, the additional dependence on height is imminent; once the cold front propagating down the right branch for $t > 0$ passes below $z = 0$, the thermobaric term becomes positive, resulting in a densification of the fluid and therefore an increase in $\Delta\sigma^a$ and velocity compared to the linear case.

### d. Effect of salinity and wind stress forcing

We now briefly discuss the influence of introducing salinity or wind stress forcing. Figure 10 shows the variation of velocity, temperature, and mass for $\eta = 0.5$, that is, with a linear temperature forcing twice as strong

as the linear salinity forcing, which effectively halves the net linear density forcing compared to the standard case ($\eta = 0$). Significant changes in the transient behavior can be observed compared to the purely thermally driven circulation (cf. Fig. 9). There also is a noticeable effect on the equilibrium properties; for the settings $\lambda = \mu = 0.1$, the final velocity is increased by 7.6% relative to the linear scenario in the full loop (compared to 3.2% for $\eta = 0$) and decreased by 16.8% in the folded loop (compared to 6.0% for $\eta = 0$). This stronger influence of cabbeling and thermobaricity is expected because when the linear contributions of salinity and temperature to the buoyancy torque are opposed, the relative importance of nonlinear effects is naturally increased.

The influence of the wind stress forcing $\tau$ on the system's behavior and stability properties is discussed in detail by Wunsch (2005) and Yuan and Wunsch (2005) for a linear equation of state. Taking $\lambda = \mu = 0.1$ and setting for example $\tau = \eta = 0.5$ induces a circulation that is 60.8% faster than in the linear scenario without salinity or wind stress forcing in the full loop and 65.7% faster in the folded loop. A detailed investigation of the impact of wind and salinity forcing or their interaction, however, is beyond the scope of the present paper and will be deferred to future work.

## 5. Discussion

Either folded or not, the loop model is without doubt a simplified description of large-scale oceanic overturning cells such as the three-dimensional, temporally varying meridional overturning circulation. It features the simplest dynamics one could think of, consisting in a single time-dependent scalar for velocity and consequently no net momentum advection. Yet, it can potentially embed the full richness and complexity of ocean thermodynamics, through the ability to use any equation of state, and also allows the study of more exotic fluid thermodynamics. Apart from being interesting in its own right, the loop model thus has great educational value when it comes to the analysis of overturning flow behavior or the demonstration of the Sandström theorem (Wunsch 2005).

A natural classification of overturning models can be obtained based on the number of dimensions used to describe the tracer fields. The Stommel box model, involving only a scalar to represent a tracer field in each box, is essentially a zero-dimensional model, while the thermohaline loop is the one-dimensional equivalent. Two-dimensional overturning models are often used in the study of horizontal convection, either in numerical simulations or laboratory experiments (Hughes and Griffiths 2008), and, finally, three-dimensional models
<figure><img src="image_11.png" /><figcaption>FIG. 10. As in Fig. 9, but for η = 0.5; that is, the linear forcing by salinity is half as strong as that by temperature. Since temperature and salinity influence density in a counteracting way, the relative importance of nonlinear effects is enhanced, which is clearly visible in the larger differences between linear and nonlinear scenarios compared to the purely thermally driven case shown in Fig. 9, especially during the transient response.</figcaption></figure>

simply correspond to the real ocean simulated with general circulation models. In this perspective, it comes as evidence that the thermohaline loop features a wider (and somewhat more realistic) range of dynamical behaviors than the Stommel box model. In particular, it allows for a more transparent investigation of loop dynamics insofar as the loop does not feature discontinuous, infinitely large, or zero Rayleigh numbers (Dewar and Huang 1996; Wunsch 2005). Importantly, this also implies that the Stommel box model is a limit case of the thermohaline loop (in practice, when the inverse Rayleigh number tends toward zero), indicating that steady states observed in a box model can be reproduced with a thermohaline loop but not necessarily vice versa.

The possibility of three-dimensional analogs of the thermohaline loop is discussed by Wunsch (2005, cf. his

section 7 and the references given therein), who underlines that the great parameter sensitivity and the significantly fewer degrees of freedom of the one-dimensional system might make its stability characteristics, time histories, and responses to external disturbances very different from higher-dimensional models. Without a careful investigation of the connection between three-dimensional circulations and the flow in one-dimensional loops (which is beyond the scope of the present paper), quantitative conclusions about the real ocean cannot be drawn based on the findings of this study. However, interpreting the thermohaline loop as a “metaphor for the circulation instead of the circulation itself” (Wunsch 2005, p. 97), some basic understanding about the circulation’s fundamental properties can indeed be gained. In the present case, the asymmetry between the interior and
the surface flow was shown to considerably affect loop dynamics and properties, especially with respect to the influence of the nonlinear effects of cabbeling and thermobaricity.

The relevance of this result can be tested by comparing our standard parameter settings to observations. Because of the highly idealized nature of the loop model, it could symbolize various kinds of overturning cells. We will here focus on the Atlantic meridional overturning circulation but propose similar scalings for Labrador Sea convective plumes and the wind-driven shallow overturning circulation of the subtropical Atlantic in Table 3.

In the North Atlantic, the mean thermocline temperature is typically about $Θ_o = 10°C$, and we estimate 20°C for the temperature in the upper thermocline (i.e., $Θ = Θ_o + 10°C$) and 5°C below the thermocline (i.e., $Θ = Θ_o - 5°C$). This roughly corresponds to the temperature distribution in the loop model when the warm section is half as long as the cold section, a case achieved when $Z_f = 0.5$, using the nondimensionalization $Θ = Θ_o + Θ' \times 1°C$. With a circulation time scale of 600 yr (Toggweiler and Key 2003) and 100–200 days as commonly accepted values for the frictional decay time scale of the deep ocean, the Grashof number is on the order of $F \sim 10^{-3}$. The inertialess assumption thus appears reasonable.

Using a circulation time scale of 600 yr and a vertical length scale of $a_v \sim 10^3$ m, we obtain a vertical diffusivity of $\kappa_v \sim 10^{-5} m^2 s^{-1}$ for $R = 0.1$, a number in general agreement with observations (Toole et al. 1994). Choosing instead a length scale that corresponds to the global overturning's horizontal length scale $a_h \sim 10^7$ m produces an equivalent horizontal diffusivity of $\kappa_h \sim 10^3 m^2 s^{-1}$ for $R = 0.1$. Hence, it is arguably a reasonable zero-order approximation to use a single non-dimensional number for representing horizontal and vertical diffusivities. We hypothesize that it may be possible to slightly modify the loop model using an elliptic rather than a circular shape to represent the large aspect ratio characterizing large-scale ocean overturning cells; however, we do not attempt to develop such a model extension here. Our hypothesis is comforted by the fact that we find similar inverse Rayleigh numbers for the horizontal and vertical components of the overturning circulations under consideration (cf. Table 3). In a generalized loop model, it would be useful to be able to apply different Rayleigh numbers for horizontal and vertical circulations, similar to the approach by Hazewinkel et al. (2012) in the context of stressed horizontal convection. On the other hand, even the discrimination between (fixed) horizontal and vertical diffusivities would only yield a poor description of

<table><caption>TABLE 3. Suitable parameter settings for application of the thermohaline loop to Labrador Sea convective plumes (Phenomenon 1) and the shallow overturning circulation of the subtropical Atlantic (Phenomenon 2). Refer, for example, to Klinger et al. (1996) and Marshall and Schott (1999) for characteristics of convective plumes and to Schott et al. (2004) and Boccaletti (2005) for those of the North Atlantic subtropical shallow overturning cells. The notation is as in the main text; we furthermore introduce $\tau_o$ for the overturning time scale, $Θ^+$ for the temperature estimate of the upper, and $Θ^-$ for that of the lower branch of the cell.</caption><thead><tr><th></th><th colspan="2">Phenomenon</th></tr><tr><th></th><th>1</th><th>2</th></tr></thead><tbody><tr><td>S (g kg<sup>-1</sup>)</td><td>35.00</td><td>37.17</td></tr><tr><td>Θ<sub>0</sub> (°C)</td><td>3.4</td><td>23.0</td></tr><tr><td>Θ<sup>+</sup> (°C)</td><td>6.0</td><td>26.0</td></tr><tr><td>Θ<sup>-</sup> (°C)</td><td>2.7</td><td>20.0</td></tr><tr><td>τ<sub>o</sub></td><td>4 weeks</td><td>10 yr</td></tr><tr><td>a<sub>h</sub> (km)</td><td>1.5</td><td>2000.0</td></tr><tr><td>a<sub>v</sub> (km)</td><td>1.5</td><td>0.2</td></tr><tr><td>Z<sub>f</sub></td><td>0.8</td><td>0.5</td></tr><tr><td>λ'</td><td>0.14</td><td>0.03</td></tr><tr><td>μ'</td><td>0.30</td><td>0.01</td></tr><tr><td>κ<sub>h</sub> (m<sup>2</sup> s<sup>-1</sup>)</td><td>10</td><td>10<sup>3</sup></td></tr><tr><td>κ<sub>v</sub> (m<sup>2</sup> s<sup>-1</sup>)</td><td>10</td><td>10<sup>-5</sup></td></tr><tr><td>R<sub>h</sub></td><td>11</td><td>0.1</td></tr><tr><td>R<sub>v</sub></td><td>11</td><td>0.1</td></tr></tbody></table>

oceanic turbulence, which is inhomogeneous in both space and time (Polzin et al. 1997; Whalen et al. 2012). By mixing density in the vertical, this turbulence is an important contributor to driving the global overturning circulation (e.g., Munk and Wunsch 1998; Visbeck 2007; Talley 2013), and in the past years, efforts have been made to develop more elaborate parameterizations that, for example, involve the sources of this mixing [e.g., breaking internal gravity waves, as in Müller and Natarov (2003) and Olbers and Eden (2013)].

The cabbeling parameter $\lambda = (1/\alpha)(\partial\alpha/\partial\Theta)$ and the thermobaricity parameter $\mu = (1/\alpha)(\partial\alpha/\partial z)$ are here approximated by their particular values at $S = 35 \text{ g kg}^{-1}$, $\Theta = 10°C$, and $p = 0 \text{ dbar}$, based on the International Thermodynamic Equation Of Seawater—2010 (TEOS-10) values (IOC et al. 2010) and using the approximate equivalence $z = p \times 1 \text{ m dbar}^{-1}$. In this case, the cabbeling parameter is approximately $\lambda = 6 \times 10^{-2} \text{ K}^{-1}$, and the thermobaric parameter $\mu = 1.2 \times 10^{-4} \text{ m}^{-1}$. For a temperature scaling factor of $\Theta^s = 1°C$, the non-dimensional cabbeling parameter can then be specified as $\lambda' = \lambda \Theta^s \approx 0.1$, and the nondimensional thermobaric parameter, setting the length scale to $a \sim 10^3 \text{ m}$, can be specified as $\mu \approx 0.1$. Note that both the cabbeling and thermobaricity parameter are to the first order inversely proportional to the thermal expansion, which itself increases with temperature. Hence, nonlinear EOS effects are expected to be more important for colder
circulations; for example, when studying convective plumes in the Labrador Sea, where the average temperature is a few degrees only, the cabbeling parameter would have to be increased by a factor of 1.5 and the thermobaric parameter by a factor of 3 compared to our standard values (cf. Table 3). As depicted in Fig. 10, they are also more important when a salinity forcing is applied, too, because the linear contributions by temperature and salinity to the density forcing counteract. The implications of neglecting salinity for the cabbeling term itself, however, are minor for our standard parameter settings: the “densification upon mixing” (cabbeling) is here solely attributed to the variance of temperature [cf. McDougall and Garrett 1992], for the discussion of cabbeling in a system where density is a function of temperature alone], which is a reasonable approximation of conditions in the real ocean, where the nonlinear contraction due to salinity mixing is approximately an order of magnitude smaller than that due to the mixing of heat (McDougall and Garrett 1992; Schanze and Schmitt 2013).

**6. Summary and outlook**

The asymmetry between the essentially horizontal surface flow and the interior circulation, which spans both horizontal and vertical scales, is a fundamental characteristic of many oceanic overturning cells. In this study, we have investigated the effect of this asymmetry on dynamics and fluid properties in the context of a thermohaline loop. To that end, we have introduced a simple modification of the loop model devised by Wunsch (2005): by folding the loop at the level where the point sources and sinks of temperature and salinity are applied, until the upper part is purely horizontal, this asymmetry can be accounted for in a simple way. In such a setting, the forcing is concentrated at the uppermost level of the loop, which is in good agreement with the observation that in reality the buoyancy fluxes mainly occur at the ocean's surface.

In our analysis, we focused on weak diffusion (small values of the inverse Rayleigh number $R$) and a weakly nonlinear equation of state (relatively small nonlinear parameters). The circulation was purely thermally driven (no salinity or wind stress forcing applied), with the source and sink of heat located at the same level to best represent the case of horizontal convection. Folding the loop led to significant changes in the transitory behavior and steady-state properties, including velocity oscillations and considerably higher temperature amplitudes associated with a slower circulation.

The influence of the two nonlinear EOS effects of cabbeling and thermobaricity on loop dynamics and steady-state properties differed significantly in the two

model versions; for example, cabbeling was found to dampen the velocity oscillations observed in the folded loop, and its effect on the equilibrium velocity was reversed compared to the circular loop. This is because in the folded loop there is no compensation between density torques arising above and below the level of the heat forcing—an unrealistic property of the symmetric, circular loop model since the surface flow is basically horizontal in the real ocean. Folding the loop also induced a modification of the system’s stability properties, with dramatically different dynamic responses when large cabbeling parameters were used (cf. Fig. 6b). A detailed investigation of these differences will be the subject of future work, as the focus here was on weakly nonlinear forms of the EOS.

From the study of the one-dimensional fluid loop, only a basic, qualitative understanding of real ocean dynamics can be gained. Because of the wide range of parameter choices, different interpretations are possible; a useful analogy with the MOC, however, is encouraged by the observation of the following features in the folded loop model:

1) a thermocline just below the warm source, with a typical thickness of a few hundred of meters, below which cold waters are sitting;

2) a deep-reaching, cold vertical column, homogeneous almost to the bottom, below the cold source, reminiscent of the cold weakly stratified polar regions; and

3) a surface region where the temperature is gradually decreasing poleward in a transition zone between a large body of warm surface waters and the region of deep convection.

Nevertheless, several essential aspects of the global overturning circulation are missing from this picture. We have only briefly addressed the influence of salinity and wind stress forcing, which were seen to strongly affect the transient behavior and steady-state properties. In particular, the relative importance of nonlinear effects was enhanced when a salinity forcing was accounted for because the linear temperature and salinity forcings counteract. Moreover, the real MOC can be regarded as a superposition of various overturning cells with different aspect ratios rather than a single cell (Talley 2003), suggesting the implementation of a more complex geometry, for example, with an additional loop representing the Antarctic Bottom Water cell. Furthermore, boundary conditions that combine temperature relaxation with a fixed flux of salinity are generally considered more realistic than the fixed flux conditions for both temperature and salinity applied in this study (Dewar and Huang 1996; Arzel et al. 2006). The different forcing types for salinity and temperature are
regarded as pivotal for the emergence of multiple equilibria in the Stommel box model (Marotzke 2000), and it would be interesting to test whether such multiple equilibrium states can be found with the thermohaline loop, and if so, how these are affected by nonlinear EOS effects or the folding of the loop.

*Acknowledgments*. We thank Jonas Nycander, Carl Wunsch, and two anonymous reviewers for their helpful comments. This work has been funded by the Bolin Centre for Climate Research (RA1).

REFERENCES

Arzel, O., T. Huck, and A. Colin de Verdière, 2006: The different nature of the interdecadal variability of the thermohaline circulation under mixed and flux boundary conditions. *J. Phys. Oceanogr.*, **36**, 1703–1718, doi:10.1175/JPO2938.1.

Boccaletti, G., 2005: Timescales and dynamics of the formation of a thermocline. *Dyn. Atmos. Oceans*, **39**, 21–40, doi:10.1016/j.dynatmoce.2004.10.010.

Bryan, F., 1986: High-latitude salinity effects and interhemispheric thermohaline circulations. *Nature*, **323**, 301–304, doi:10.1038/323301a0.

Dewar, W., and R. Huang, 1996: On the forced flow of salty water in a loop. *Phys. Fluids*, **8**, 954–970, doi:10.1063/1.868874.

Hazewinkel, J., F. Paparella, and W. Young, 2012: Stressed horizontal convection. *J. Fluid Mech.*, **692**, 317–331, doi:10.1017/jfm.2011.514.

Held, I. M., 2005: The gap between simulation and understanding in climate modeling. *Bull. Amer. Meteor. Soc.*, **86**, 1609–1614, doi:10.1175/BAMS-86-11-1609.

Hieronymus, M., and J. Nycander, 2013: The buoyancy budget with a nonlinear equation of state. *J. Phys. Oceanogr.*, **43**, 176–186, doi:10.1175/JPO-D-12-063.1.

Huang, R. X., 1999: Mixing and energetics of the oceanic thermohaline circulation. *J. Phys. Oceanogr.*, **29**, 727–746, doi:10.1175/1520-0485(1999)029<0727:MAOTO>2.0.CO;2.

Hughes, G. O., and R. W. Griffiths, 2008: Horizontal convection. *Annu. Rev. Fluid Mech.*, **40**, 185–208, doi:10.1146/annurev.fluid.40.111406.102148.

IOC, SCOR, and IAPSO, 2010: The International Thermodynamic Equation of Seawater—2010: Calculation and use of thermodynamic properties. Intergovernmental Oceanographic Commission, Manuals and Guides 56, 220 pp. [Available online at http://www.teos-10.org/pubs/TEOS-10_Manual.pdf.]

Klinger, B. A., and U. Send, 1996: Representation of convective plumes by vertical adjustment. *J. Geophys. Res.*, **101**, 18 175–18 182, doi:10.1029/96JC00861.

Klocker, A., and T. J. McDougall, 2010: Influence of the nonlinear equation of state on global estimates of dianeutral advection and diffusion. *J. Phys. Oceanogr.*, **40**, 1690–1709, doi:10.1175/2010JPO4303.1.

Malkus, W. V., 1972: Non-periodic convection at high and low Prandtl number. *Mem. Soc. Roy. Sci. Liege*, **4**, 125–128.

Manabe, B. S., and R. J. Stouffer, 1999: Are two modes of thermohaline circulation stable? *Tellus*, **51A**, 400–411, doi:10.1034/j.1600-0870.1999.t01-3-00005.x.

Marotzke, J., 2000: Abrupt climate change and thermohaline circulation: Mechanisms and predictability. *Proc. Natl. Acad. Sci. USA*, **97**, 1347–1350, doi:10.1073/pnas.97.4.1347.

Marshall, J., and F. Schott, 1999: Open-ocean convection: Observations, theory, and models. *Rev. Geophys.*, **37**, 1–64, doi:10.1029/98RG02739.

McDougall, T. J., 1987: Thermobaricity, cabbeling, and water-mass conversion. *J. Geophys. Res.*, **92**, 5448–5464, doi:10.1029/JC092iC05p05448.

—, 2003: Potential enthalpy: A conservative oceanic variable for evaluating heat content and heat fluxes. *J. Phys. Oceanogr.*, **33**, 945–963, doi:10.1175/1520-0485(2003)033<0945:PEACOV>2.0.CO;2.

—, and C. J. Garrett, 1992: Scalar conservation equations in a turbulent ocean. *Deep-Sea Res.*, **39A**, 1953–1966, doi:10.1016/0198-0149(92)90007-G.

—, D. R. Jackett, and F. J. Millero, 2009: An algorithm for estimating Absolute Salinity in the global ocean. *Ocean Sci. Discuss.*, **6**, 215–242, doi:10.5194/osd-6-215-2009.

Miljkovic, N., and E. N. Wang, 2011: Modeling and optimization of hybrid solar thermoelectric systems with thermosyphons. *Sol. Energy*, **85**, 2843–2855, doi:10.1016/j.solener.2011.08.021.

Müller, P., and A. Natarov, 2003: The internal wave action model (IWAM). *Near-Boundary Processes and Their Parameterization: Proc. 13th 'Aha Huliko'a Winter Workshop*, Honolulu, HI, University of Hawai'i at Manoa, 95–105.

Munk, W., and C. Wunsch, 1998: Abyssal recipes II: Energetics of tidal and wind mixing. *Deep-Sea Res. I*, **45**, 1977–2010, doi:10.1016/S0967-0637(98)00070-3.

Olbers, D., and C. Eden, 2013: A global model for the diapycnal diffusivity induced by internal gravity waves. *J. Phys. Oceanogr.*, **43**, 1759–1779, doi:10.1175/JPO-D-12-0207.1.

Polzin, K., J. Toole, J. Ledwell, and R. Schmitt, 1997: Spatial variability of turbulent mixing in the abyssal ocean. *Science*, **276**, 93–96, doi:10.1126/science.276.5309.93.

Roquet, F., 2013: Dynamical potential energy: A new approach to ocean energetics. *J. Phys. Oceanogr.*, **43**, 457–476, doi:10.1175/JPO-D-12-098.1.

Ruddick, B., and L. Zhang, 1996: Qualitative behavior and nonoscillation of Stommel's thermohaline box model. *J. Climate*, **9**, 2768–2777, doi:10.1175/1520-0442(1996)009<2768:QBANOS>2.0.CO;2.

Schanze, J. J., and R. W. Schmitt, 2013: Estimates of cabbeling in the global ocean. *J. Phys. Oceanogr.*, **43**, 698–705, doi:10.1175/JPO-D-12-0119.1.

Schott, F. A., J. P. Mccreary, and G. C. Johnson, 2004: Shallow overturning circulations of the tropical-subtropical oceans. *Earth's Climate, Geophys. Monogr.*, Vol. 147, Amer. Geophys. Union, 261–304.

Stommel, H., 1961: Thermohaline convection with two stable regimes of flow. *Tellus*, **13A**, 224–230, doi:10.1111/j.2153-3490.1961.tb00079.x.

Talley, L. D., 2003: Shallow, intermediate, and deep overturning components of the global heat budget. *J. Phys. Oceanogr.*, **33**, 530–560, doi:10.1175/1520-0485(2003)033<0530:SIADOC>2.0.CO;2.

—, 2013: Closure of the global overturning circulation through the Indian, Pacific, and Southern Oceans: Schematics and transports. *Oceanography*, **26**, 80–97, doi:10.5670/oceanog.2013.07.

Toggweiler, J., and R. Key, 2003: Thermohaline circulation. *Encyclopedia of Atmospheric Sciences*, A. C. Judith and J. A. Pyle, Eds., Elsevier, 1549–1555.

Toole, J. M., R. W. Schmitt, and K. L. Polzin, 1994: Estimates of diapycnal mixing in the abyssal ocean. *Science*, **264**, 1120–1123, doi:10.1126/science.264.5162.1120.
Vallis, G. K., 2006: *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 745 pp.

Visbeck, M., 2007: Oceanography: Power of pull. *Nature*, **447**, 383–383, doi:10.1038/447383a.

Welander, P., 1967: On the oscillatory instability of a differentially heated fluid loop. *J. Fluid Mech.*, **29**, 17–30, doi:10.1017/S0022112067000606.

—, 1986: Thermohaline effects in the ocean circulation and related simple models. *Large-Scale Transport Processes in Oceans and Atmosphere*, J. Willebrand and D. L. T. Anderson, Eds., Springer, 163–200.

Whalen, C., L. Talley, and J. MacKinnon, 2012: Spatial and temporal variability of global ocean mixing inferred from Argo

profiles. *Geophys. Res. Lett.*, **39**, L18612, doi:10.1029/2012GL053196.

Wunsch, C., 2005: Thermohaline loops, Stommel box models, and the Sandström theorem. *Tellus*, **57A**, 84–99, doi:10.1111/j.1600-0870.2005.00093.x.

Young, W. R., 2010: Dynamic enthalpy, Conservative Temperature, and the seawater Boussinesq approximation. *J. Phys. Oceanogr.*, **40**, 394–400, doi:10.1175/2009JPO4294.1.

Yuan, S., and C. Wunsch, 2005: Stress-driven thermohaline loops. *Phys. Fluids*, **17**, 066601, doi:10.1063/1.1927887.

Zvirin, Y., 1982: A review of natural circulation loops in pressurized water reactors and other systems. *Nucl. Eng. Des.*, **67**, 203–225, doi:10.1016/0029-5493(82)90142-X.