InTeGrate Modules and Courses >Modeling Earth Systems > Student Materials > Unit 8 Reading
InTeGrate's Earth-focused Modules and Courses for the Undergraduate Classroom
showLearn More
These materials are part of a collection of classroom-tested modules and courses developed by InTeGrate. The materials engage students in understanding the earth system as it intertwines with key societal issues. The collection is freely available and ready to be adapted by undergraduate educators across a range of courses including: general education or majors courses in Earth-focused disciplines such as geoscience or environmental science, social science, engineering, and other sciences, as well as courses for interdisciplinary programs.
Explore the Collection »
show Download
The student materials are available for offline viewing below. Downloadable versions of the instructor materials are available from this location on the instructor materials pages. Learn more about using the different versions of InTeGrate materials »

Download a PDF of all web pages for the student materials

Download a zip file that includes all the web pages and downloadable files from the student materials

For the Instructor

These student materials complement the Modeling Earth Systems Instructor Materials. If you would like your students to have access to the student materials, we suggest you either point them at the Student Version which omits the framing pages with information designed for faculty (and this box). Or you can download these pages in several formats that you can include in your course website or local Learning Managment System. Learn more about using, modifying, and sharing InTeGrate teaching materials.
Initial Publication Date: January 7, 2016

Unit 8 Reading: Modeling Thermohaline Circulation

by David Bice, The Pennsylvania State University

In this exercise, we will build a STELLA model of the thermohaline circulation in the North Atlantic Ocean based on Stommel's (1961) classic paper on this topic. This is a fascinating model in the sense that it seems to be very simple and yet its behavior is surprisingly complex. By complex, we mean that the system has more than one steady state and the evolution of the model can be very sensitive to the initial conditions.

We will see that there are thresholds in this system; crossing one of these thresholds leads to rapid, inexorable change in the state of the system. We will explore this notion of thresholds by "perturbing" the system with changes that represent realistic variability in the system.

In the real world, this system plays a critical role in the global and regional climate, so changes in its state can have important implications for climate change. In particular, we will consider the case of the Younger Dryas abrupt climate change event of about 11 ka, and use the model to evaluate some hypotheses for what may have caused this event. First, however, we will cover some background on seawater density and thermohaline circulation.


Surface ocean currents are driven mainly by the winds, but the deeper circulation of the oceans is driven primarily by density differences caused by changes in salinity and temperature. This kind of temperature- and salinity-driven flow is often called thermohaline circulation, and it is a very important feature of the oceans, exchanging water (along with all of its dissolved constituents) between the surface and the vast deepwater reservoir.

From the relationship between temperature, salinity, and density (Fig. 1) we see that colder waters are denser than warmer waters, and saltier water is denser than fresher water.

In order to initiate thermohaline circulation, the water at the surface of the ocean must become denser than the water beneath it, and then gravity will pull it down until it reaches the bottom of the ocean or until it reaches a depth where the water has the same density (Fig. 2). If water sinks in one area, it must rise in another, completing the circulation.

So, how can the surface waters become denser than the underlying water? One way is by strong cooling, and the other way is by increasing the salinity. Cooling is pretty easy to understand, but there are a couple of ways to increase the salinity. One way is to have very strong evaporation that is not balanced by precipitation; this is generally something that happens in warmer settings. Another way is to form ice from the seawater — ice formation removes pure H2O from seawater, leaving a higher concentration of salts in the remaining seawater — this is what happens at the edge of Antarctica.

In the North Atlantic, the salinity of the surface water is significantly higher in the south than in the north (Fig. 3), a fact largely due to differences in evaporation and precipitation. Evaporation increases salinity, while precipitation and river input decrease salinity. On the whole, there is a net export of evaporated water from the Atlantic into the Pacific, and in fact, the Atlantic is more saline than any other ocean, making it a prime location for the formation of deep water currents.

The Gulf Stream (Fig. 4) delivers some of this warm, salty water into the North Atlantic, where it cools and becomes dense enough to sink to the bottom of the ocean, initiating one of the main deep currents in the global ocean, the North Atlantic Deep Water (NADW). As this dense water sinks, it tends to draw in additional warm salty water from the Gulf Stream, which also cools, becomes dense and sinks.

The NADW is an enormous flow, with a rate of about 20 Sverdrups (one Sverdrup = 1e6 m3/s; Sverdrup was a famous oceanographer). For comparison, the flow rate of the Amazon is about 0.2 Sverdrups, and all of the major rivers of the world combined have a flow rate of about 1.2 Sverdrups. The actual velocity of the NADW flow is quite low, however — around 1 to 1.5 cm/s.

We can map out this flow through the world's oceans because it has a distinctive set of physical and chemical characteristics that do not change very quickly — it retains a sort of fingerprint — enabling us to see distinct water masses in the deep ocean. A series of transects (Figs. 5 – 6) through the Atlantic Ocean illustrate how the chemistry allows us to distinguish different water masses.

Using these kinds of fingerprints, the pattern of deep thermohaline circulation has been mapped out, as shown in Figure 7.

This NADW eventually returns to the surface and makes its way back into the North Atlantic, where it rejoins the Gulf Stream, creating what Broecker (1991) calls the Great Ocean Conveyor, which is of great importance to the global climate system.

To get a sense for how important this system is, consider the following rough calculation of how much heat is delivered to the polar region.

Volumetric flux or discharge of NADW formation = 2e7 m3/s

Density of seawater = 1e3 kg/m3

Mass flux (Fm) = (2e7 m3/s) x (1e3 kg/m3) = 2e10 kg/s

Temperature change during cooling (` DeltaT`) = 25°C

Heat capacity of water (Cp) = 4184 Joules .kg-1 .°C-1

Energy flux = Cp x ` DeltaT` x Fm = 2.09e15 Joules/s (same units as a Watt)

This is a big number — a lot of energy — but what does it mean? If you take this energy and spread it out over the whole polar region above 60°N latitude, you come up with 61 W/m2. For the sake of comparison, the sunlight up there provides an average of 190 W/m2, so this process provides about 30% of the total energy reaching that region, which is enough to cause ice sheets to grow or shrink, and as they grow or shrink, they can change the albedo to the point where it can influence the climate in the polar region. This Arctic region ice-albedo feedback mechanism can then influence the thermohaline circulation and eventually the greenhouse effect, thus producing a global climate response. This chain of feedbacks is the main reason why the Pleistocene glacial-interglacial cycles are synchronous in both the northern and southern hemispheres.

Clearly, this is a critically important part of the climate system, so we should be curious about how it works. Some questions we might ask include: does it always operate at the same rate, or does the conveyor belt circulation speed up and slow down? If the latter, why? And, what might the consequences be if the circulation were to speed up or slow down? There are two main approaches we might take in pursuing these questions. One would be to study sea floor sediments, looking for some kind of geochemical or sedimentological proxy for the operation of this current. Another approach would be to create numerical models of this system and then experiment with those models — this is the approach we will take here.

Back in 1961, Henry Stommel, a famous American oceanographer, studied thermohaline circulation in the form of a simple system, one in which there are two boxes or reservoirs of water — one warm and salty, the other cold and fresher. Water is exchanged between these boxes by means of both a surface flow and a deep flow. By exploring the behavior of this system, Stommel discovered some unexpected complexity that is now widely believed to be the cause of some dramatic and abrupt climate shifts that have occurred in the last 100,000 years. These abrupt climate changes are revealed by the oxygen isotope record of temperature supplied by ice cores from Greenland (Fig. 8).

Let us examine this further. The oxygen in water can occur as either 16O (by far the most common), 17O, or 18O (all of which are stable), so there are slight mass differences in water molecules, depending on which oxygen isotope they contain. As we learned in Unit 5, because the processes involved in the movement of water (evaporation and precipitation) are sensitive to these mass differences, the relative concentrations of lighter and heavier water can change over space and time. Heavier water is precipitated more easily, so by the time an air mass moves into the interior of an ice sheet, it has lost most of its heavy water vapor (with 18O) and is relatively enriched in the lighter water (with 16O). The snow falling on the ice sheet will therefore be much lighter than the ocean water from which the atmospheric vapor first evaporated. The relative abundance of the heavy and light isotopes in glacial ice are reported in terms of a ratio that looks like this:

` del{::}^{18}"O" = (( ({::}^{18}"O") / ({::}^{16}"O"))_"sample" - ( ({::}^{18}"O") / ({::}^{16}"O"))_"standard")/(({::}^{18}"O") / ({::}^{16}"O"))_"standard" * 1000 `   

The standard is often normal ocean water, so the precipitation falling on the ice sheet will always be negative (since the sample isotopic ratio will be smaller than the ocean water isotopic ratio), and it turns out that the colder the atmospheric temperature is at the time of precipitation, the more negative the isotopic value will be. In fact, there is such a good relationship between these two variables that we can use the oxygen isotope values to determine past temperatures.

Teams of scientists have done just this for ice cores taken from around the world and have found something very interesting in Greenland. The oxygen isotopic record of the past 20 kyrs or so indicates that as Earth began warming at the end of the last ice age, a sudden cooling occurred, followed by abrupt warming (Fig. 8). The cooling, known as the Younger Dryas event (Fig. 8) lasted about a thousand years, and the transitions between cold and warm states appear to have occurred in a matter of years to decades. The magnitudes of the temperature changes (~ 10 °C) are shocking, and the Younger Dryas is not the only such event to have been recorded in Greenlandic ice cores; there are at least four other events like it in the past 50 kyr. These events lead us to some important questions regarding the causes of this climatic flipping, namely: What triggers the switch? How does it work? The general scientific consensus (e.g., Rahmstorf, 2002) is that the abrupt and dramatic northern hemisphere climate changes shown in this Greenland oxygen isotope record are most likely related to abrupt changes in the thermohaline circulation system. We can make some progress toward understanding how this works with a simple model based on Stommel's work.

The Model

Like any good modeler, Stommel began with a conceptual model and then figured out how to express the processes in the model with some simple, yet realistic, equations. Simplicity was especially important considering that Stommel did not have a computer! The conceptual diagram below (Fig. 9) is adapted from Stommel's paper.

It helps to study Fig. 9 carefully before starting to build the model. The two boxes or regions of the ocean (equatorial and polar) have inherently different salinities and temperatures due to the global climate — the equatorial box is hot and salty, while the polar box is cold and less saline (fresher) — and the global climate will always maintain these differences. These differences are ` DeltaT_(eq) ` and ` DeltaS_(eq) ` in the equations below. These inherent differences in salinity and temperature create a density difference between the two boxes (`Deltarho `; `rho`, pronounced rho, is the Greek letter used for density), which then drives a flow between the two boxes (Q). The flow between the boxes mixes them and thus tends to reduce the differences between the two. The mixing flow is divided equally into a surface flow (i.e., the Gulf Stream) and a deep flow (i.e., the NADW) that are closely connected — you cannot have one without the other.

It is important to understand an important negative feedback in this system involving Q and the pole-equator temperature difference ` DeltaT `. When the temperature difference, ` DeltaT `, is great, then Q will be strong, but this transports more heat to the polar region, thus lessening ` DeltaT `, which in turn weakens Q. On the other hand, if Q is weak, then there is little heat transport to the polar region, allowing this region to grow colder, thus increasing ` DeltaT `, which can then in turn leads to a stronger Q and thus more heat to the polar region.

You might look at the conceptual model (Fig. 9) and think that we will need to keep track of four different quantities (i.e., four reservoirs) — the salinity and temperature of the polar box, and the salinity and temperature of the equatorial box. But Stommel realized that what really matters is the salinity and temperature differences between the two regions, because it is these differences that lead to a density difference, which in turn controls the strength of the mixing flow. It is worth pausing for a bit and appreciating that these differences are unusual kinds of reservoir quantities.

Stommel described this system in terms of two basic differential equations for `DeltaT` and `DeltaS` — the temperature and salinity differences between the polar and tropical surface ocean:

` (dDeltaT)/(dt) = (DeltaT_(eq) - DeltaT) - QDeltaT `      Eqn. 1 ` (dDeltaS)/(dt) = ((DeltaS_(eq) - DeltaS))/delta - QDeltaS `      Eqn. 2 Here, `DeltaT_(eq) ` and `DeltaS_(eq) ` are the equilibrium temperature and salinity differences between the two parts of the ocean maintained by the difference in heating of the ocean at low versus high latitudes and by differences in the relative importance of evaporation and precipitation associated with the atmospheric Hadley cell circulation; they are normalized (meaning they vary from 0 to 1) and dimensionless numbers. These equilibrium differences determine where `DeltaT ` and `DeltaS ` would end up if you turned off the flux of ocean water, Q, between the two reservoirs. In other words, given the climatological differences between the polar region and the tropics, there are natural differences in ocean temperature and salinity, and the system will tend to move toward those equilibrium differences. Equation 2 above includes a parameter ` delta`, which is a factor (set as 6 here) that expresses the fact that salinity changes more slowly than temperature — it diffuses more slowly than heat. If we ignore the flux (Q) in equations 1 and 2, we can see that for `DeltaT `, the response time should be one time unit and for `DeltaS `, the response time should be about 6 time units. The fact that the two reservoirs here have different response times is one of the main reasons behind some of the surprisingly complex behaviors that you will see in this system.

Q is another dimensionless number representing the flow of water between the two parts of the ocean — the thermohaline circulation; it is defined as:

` Q = (Deltarho)/lambda `     Eqn. 3
where ` lambda ` is essentially the resistance to flow (set to 0.2), and `Deltarho` is the density difference between the two parts of the ocean:
` Deltarho = |RDeltaS - DeltaT| `     Eqn. 4

where R is a ratio (with a value of 2) reflecting the fact that density is more sensitive to salinity changes than to temperature changes. In Stommel's model, one unit of time represents about 200 years.

The flow between the polar and equatorial regions, Q, is important in this model since it helps to determine the temperature and salinity differences between the two regions. If Q is large, the model vigorously mixes water from the two regions, bringing their temperatures and salinities closer together, so that `DeltaT ` and `DeltaS ` become smaller. On the other hand, if Q is small, then mixing is reduced and `DeltaT ` and `DeltaS ` become greater and approach the values of `DeltaT_(eq) ` and `DeltaS_(eq) ` .

Complexity in Systems Dynamics

Before we embark on building and experimenting with a model of this system, it might help to spend a bit of time considering some important characteristics of complex systems. As mentioned earlier, one characteristic of complex systems is that they have more than one steady state. This idea is illustrated in the simple cartoon in Fig. 10.

Along the x-axis, we have some system variable — like the concentration of some chemical in water — and the blue curve represents something like the potential energy of the system. You can also conceptualize this figure as though it were a topographic profile. There are four places where the slope of the blue curve equals 0 — these are equilibrium points. Points 1 and 3 are stable in the sense that if you "perturb" or alter the variable when it is in this position, its value will fall back down to the original position. Furthermore, there is effectively a "capture area" for each equilibrium point — a range of values the system variable can have that will still return the system to the equilibrium point. We call these stable equilibrium points steady states. You can see that steady state 1 is more stable than steady state 3 — 1 has a big capture area, and you need a bigger boost of energy to get it over the hump that separates 1 and 3. Points 2 and 4 are also equilibrium points, but they are unstable in the sense that the tiniest change or perturbation in the system variable will cause it to move into one of the stable points (think of a ball sitting at the crest of a hill). These unstable points are like topographic divides separating the capture areas of the steady states; in systems dynamics, these conditions are called thresholds or tipping points. Once a threshold or tipping point is crossed, the system may rapidly shift into a new steady state. If the system is near a threshold, it is obviously very sensitive to perturbations. Lenton et al. (2007) provide a nice summary of some important tipping points relevant to the climate system.

Another characteristic of complex systems is a sensitivity to initial conditions. In Fig. 10, the initial conditions can be thought of as a point along the x-axis where the system begins. If the initial conditions place the system near one of the divides, then a tiny shift one way or another will make a big difference on what steady state the system evolves to. On the other hand, if the initial conditions fall well inside one of the capture areas, the end result is not really sensitive to changes in the initial conditions.

There are a couple of important points related to this notion of thresholds or tipping points. One is that if you assume that the Younger Dryas event (Fig. 8) is related to a change in the thermohaline circulation, the fact that the transition is so fast makes us suspect that it is related to crossing a threshold. If such a change were to happen today, it would have grave consequences for our climate and this in turn would result in very serious economic consequences (i.e., climate-related damages). This makes the business of risk assessment related to climate change very difficult — the fact that there are these kinds of thresholds with severe consequences means that we cannot count on there being a simple relationship between the amount of global warming and the amount of climate damages that the global economy must absorb.

Obviously, knowing the details about these thresholds is of critical importance, and modeling provides some insight, but real-time monitoring of the state of the thermohaline circulation is equally important. At the current time, the Gulf Stream (i.e., Q in our model) is relatively strong and this delivers a lot of heat to the polar region, but this current state of the system could conceivably cross a threshold if there is a polar warming and/or a decrease in the salinity of the surface waters of the North Atlantic (Rahmstorf, 2002). In fact, the polar region is warming and the salinity is decreasing (Curry et al., 2005), both of which are moving us closer to the tipping point. If the salinity is lowered enough, and/or if the polar region is warmed enough (precisely how much is an unresolved question), then Q could be greatly diminished and the result would be a very strong, sudden cooling throughout most of the northern hemisphere, and a bit of warming in the southern hemisphere (Rahmstorf, 2002); the economic consequences of such a change would be staggering, which makes the thermohaline circulation system so important to understand.


Curry, R., Dickson, B, and Yashayaev, I., 2003, "A change in the freshwater balance of the Atlantic Ocean over the past four decades", Nature, v. 426, p. 826-829. doi:10.1038/nature02206.

Lenton, T., and others, 2007, "Tipping elements in the Earth system", PNAS, p.1786–1793, doi: 10.1073/pnas.0705414105.

Rahmstorf, S., 2002, "Ocean circulation and climate during the past 120,000 years", Nature , v. 419., p. 207-214.

Stommel, H., 1961, "Thermohaline Convection with Two Stable Regimes of Flow," Tellus, v. 13, p. 224–230. DOI: 10.1111/j.2153-3490.1961.tb00079.x

These materials are part of a collection of classroom-tested modules and courses developed by InTeGrate. The materials engage students in understanding the earth system as it intertwines with key societal issues. The collection is freely available and ready to be adapted by undergraduate educators across a range of courses including: general education or majors courses in Earth-focused disciplines such as geoscience or environmental science, social science, engineering, and other sciences, as well as courses for interdisciplinary programs.
Explore the Collection »