In his work (Stommel 1961) did not account for evaporation or precipitation, so there is no virtual salt flux, \(H = 0\), but each box does exchange salinity with its surrounding environment, so \(d > 0\). The set of four differential equations presented in Chapter 23 (Equation 23.6) can be reduced two coupled differential equations for the mean temperature and salinity. We notice that if we subtract the pairs of equations, then we can write \[ \begin{aligned} \frac{\mathrm{d}\Delta T}{\mathrm{d}t} = c(\Delta T^* - \Delta T) - 2|q|\Delta T, \\ \frac{\mathrm{d}\Delta S}{\mathrm{d}t} = d(\Delta S^* - \Delta S) - 2|q|\Delta S, \end{aligned} \tag{24.1}\]
where \(\Delta T^* = 2T^*\), \(\Delta S^* = 2S^*\), and we have defined the two new unknowns, \[ \begin{aligned} \Delta T &= T_2 - T_1, \\ \Delta S &= S_2 - S_1, \end{aligned} \] and now, the flow \(q\) is given by \[ q = k(\alpha \Delta T - \beta \Delta S). \]
We can then nondimensionalise the system by setting \[ \Delta T = {\Delta T^*}y, \quad \Delta S = \Delta S^* x, \qquad t = [t]t', \] where \(t'\) is nondimensional time. Annoyingly, notice that we have flipped the order of the equations, i.e. above we presented temperature and salinity in that order, and now this is associated with \(y\) and \(x\).
For example, the equation in \(\Delta T\) now becomes \[ \frac{\mathrm{d} y}{\mathrm{d}t'} = c[t](1 - y) - 2|q|[t] y \] and we see that we should select \[ [t] = \frac{1}{c}. \] If we now examine the remaining factor: \[ 2|q|[t] = \frac{2k \alpha \Delta T^*}{c} \left\lvert -y + \frac{\beta \Delta S^*}{\alpha \Delta T^*}\right\rvert = \frac{1}{\lambda}(Rx - y), \] we can see where the constants \(\lambda\) and \(R\) are defined from. Notice that there has been a negative sign flip within the absolute values.
Dropping the primes henceforth, we have the following set of non-dimensional equations to study for the unknowns, \(x = x(t)\), and \(y = y(t)\): \[ \begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t} &= \delta(1 - x) - |f(x, y)|x, \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= 1 - y - |f(x, y)|y, \end{aligned} \tag{24.2}\] where we have introduced the function, \[ f(x, y; R, \lambda) = \frac{1}{\lambda}(Rx - y), \tag{24.3}\] where there are now three non-dimensional parameters given by \[ \begin{aligned} \delta &= d/c, \\ \lambda &= c/(2\alpha k \Delta T^*), \\ R &= \beta \Delta S^*/(\alpha \Delta T^*). \end{aligned} \]
Together Equation 24.2 and Equation 24.3 form a system of equations for \((x(t), y(t))\), with parameters \(\delta\), \(\lambda\), and \(R\).
Although there are many parameters, it is important to note that the quantity \(f\) is crucial since it corresponds to the direction of flow in the bottom pipe of the Ocean model. By convention, the system was set up so that \(f > 0\) corresponds to flow through the bottom pipe goes in the direction of the equator (Box 1) as a result of higher densities at the high latitudes. We are subsequently interested in whether it is possible for \(f\) to switch sign, which would correspond to a bottom flow going from equator to poles, and the entire system switching direction.
Let the equilibrium states be given by \((x^*, y^*)\) with \(f(x^*, y^*) = f^*\). Then setting the right hand sides of (Equation 24.2) to zero and solving, we find \[ x^* = \frac{\delta}{\delta + |f^*|} \quad \text{and} \quad y^* = \frac{1}{1 + |f^*|}. \tag{24.4}\] We can solve for \(f^*\) in both equations and equate the result to one another. This gives \[ \delta \frac{1-x^*}{x^*} = \frac{1 - y^*}{y^*} = |f^*|. \] Therefore the equilibrium points lie along the above curve. Let us generate different values of \(f^*\) and plot the combination.
import numpy as np
import matplotlib.pyplot as plt
delt = 1/6
fmat = np.linspace(1e-3, 10, 101)
xs = delt/(delt + np.abs(fmat))
ys = 1/(1 + np.abs(fmat))
plt.plot(xs, ys, '.-')
plt.xlim((0,1))
plt.ylim((0,1))
plt.xlabel("$x^*$"); plt.ylabel("$y^*$");
plt.grid(1)
The above plot shows the steady-state values at \(\delta = 1/6\). Each point refers to a different value of \(|f^*|\). Note that as \(|f^*| \to \infty\), then \((x^*, y^*) \to (0, 0)\) while as \(|f^*| \to 0\), then \((x^*, y^*) \to (1,1)\). Now the above does not tell us what equilibrium states will exist at a specific value of \(\delta\) since it requires information about \(f^*\). We take (Equation 24.3) and combine with (Equation 24.4). We conclude that the equilibrium points must satisfy \[ \lambda f^* = \phi(f^*), \] where we have defined the function \(\phi\) according to \[ \phi(f^*) = \frac{\delta R}{\delta + |f^*|} - \frac{1}{1 + |f^*|}. \tag{24.5}\] Therefore, for each combination of the parameters \((\lambda, \delta, R)\), we must solve \(\delta f^* = \phi(f^*)\) to determine \(f^*\). Once the values of \(f^*\) are known, then the steady-states \((x^*, y^*)\) are also known. Here is a typical graph showing the potential intersections at the values of \(\lambda = 1/5\), \(\delta = 1/6\) and \(R = 2\).
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as sciopt
lam = 1/5;
delt = 1/6;
R = 2;
f = np.linspace(-2.5, 2.5, 101)
phi = delt*R/(delt + np.abs(f)) - 1/(1 + np.abs(f))
plt.plot(f, phi, 'k', label='$\phi$')
plt.plot(f, lam*f, 'b', label='$\lambda f$' )
plt.ylim((-1.5,2))
plt.xlabel("$f$")
plt.ylabel("$\phi$")
plt.legend()
plt.grid(1)
Although there are many parameters that can alter the shape of the overall graphs, note that if \(\lambda\) is increased from its current value of \(\lambda = 1/6\), then two of the left roots will disappear, leaving only a single root.
For our purposes, we will mainly proceed using numerics, and bypass the need to study the above results analytically.
In the accompanying lecture script, lecture29-StommelPhasePlane.ipynb, we will obtain the eigenvalues numerically, and generate the phase plane. At the test point \(R = 2\), \(\delta = 1/6\), and \(\lambda = 1/5\), we generate the following phase plane picture.
Once the above stability analysis, it is possible to sketch the final bifurcation diagram.
First, we can note that as \(\lambda\) increases, the graph in Section 24.2 indicates that there may either be three intersections (as shown) or if \(\lambda\) increases then two of the intersections may disappear. We can note that:
| \(f_i^*\) | \(\lambda \to 0^+\) | \(\lambda \to \infty\) |
|---|---|---|
| \(f_1^*\) | \(-\infty\) | none |
| \(f_2^*\) | fixed point < 0 | none |
| \(f_3^*\) | fixed point > 0 | 0 |
We can also ascertain that there should be a point, say \(\lambda = \lambda^*\) where the two roots coalesce, \(f_1^* = f_2^*\), and thereafter disappear. This allows us to now plot the bifurcation \((\lambda, f^*)\)-plane:
What is the final takeaway message(s)?
It is possible to pose a conceptual model of the ocean and the thermohaline circulation, which governs the exchange of hot-and-cold waters and salty waters from equatorial-to-pole zones.
This conceptual model suggests that tipping points and hysteresis are possible in such systems. In particular, it is possible, according to these models, for the flow directions to entirely change!