## Abstract

A thorough analysis of the stability of the uncoupled Rooth interhemispheric three-box model of thermohaline circulation (THC) is presented. The model consists of a northern high-latitude box, a tropical box, and a southern high-latitude box, which correspond to the northern, tropical, and southern Atlantic Ocean, respectively. Restoring boundary conditions are adopted for the temperature variables, and flux boundary conditions are adopted for the salinity variables. This paper examines how the strength of THC changes when the system undergoes forcings that are analogous to those of global warming conditions by applying the equilibrium state perturbations to the moisture and heat fluxes into the three boxes. In each class of experiments, using suitably defined metrics, the authors determine the boundary dividing the set of forcing scenarios that lead the system to equilibria characterized by a THC pattern similar to the present one from those that drive the system to equilibria with a reversed THC. Fast increases in the moisture flux into the northern high-latitude box are more effective than slow increases in leading the THC to a breakdown, while the increases of moisture flux into the southern high-latitude box strongly inhibit the breakdown and can prevent it, as in the case of slow increases in the Northern Hemisphere. High rates of heat flux increase in the Northern Hemisphere destabilize the system more effectively than low ones; increases in the heat fluxes in the Southern Hemisphere tend to stabilize the system.

## 1. Introduction

Following the classical picture of Sandström’s (Sandström 1908; Huang 1999), ocean circulation can be divided, albeit somewhat ambiguously, into the wind-driven circulation, which is directly induced by the mechanical action of the wind and essentially affects the surface waters, and the thermohaline circulation (THC), which is characterized by relatively robust gradients in the buoyancy of the water masses (Weaver and Hughes 1992; Rahmstorf 2000, 2002; Stocker et al. 2001). Although some studies show that the THC is very sensitive to wind stresses (Toggweiler and Samuels 1995), more recent studies suggest that the buoyancy forces are dominant when more realistic boundary conditions are considered (Rahmstorf and England 1997; Bugnion and Hill 2003). The THC dominates the large-scale picture of global ocean circulation and is especially relevant for intermediate-to-deep water.

One of the main issues in the study of climate change is the fate of the THC in the context of global warming (Weaver and Hughes 1992; Manabe and Stouffer 1993; Stocker and Schmittner 1997; Rahmstorf 1997, 1999a,b, 2000; Wang et al. 1999a,b). The THC plays a major role in the global circulation of the oceans as pictured by the conveyor belt scheme (Weaver and Hughes 1992; Stocker 2001); therefore, the effect on the climate of a change of its pattern may be global (Broecker 1997; Manabe and Stouffer 1999b; Cubasch et al. 2001; Seager et al. 2002). The THC is sensitive to changes in the climate since the North Atlantic Deep Water (NADW) formation is sensitive to variations in air temperature and precipitation in the Atlantic basin (Rahmstorf and Willebrand 1995; Rahmstorf 1996). There are several paleoclimatic datasets indicating that changes in the patterns or collapses of the THC coincide with large variations in climate, especially in the North Atlantic region (Broecker et al. 1985; Boyle and Keigwin 1987; Keigwin et al. 1994; Rahmstorf 2002).

Box models have historically played a major role in understanding the fundamental dynamics of the THC (Weaver and Hughes 1992). The Stommel two-box oceanic model (Stommel 1961) parameterized the THC strength as proportional to the density gradient between the northern latitude and equatorial boxes. The oceanic model proposed by Rooth (1982) introduced the idea that the driver of the THC was the density difference between the high-latitude basins of the Northern and Southern Hemispheres, thus implying that the THC is an interhemispheric phenomenon. Both the Stommel and Rooth models allow two stable equilibria, one—the warm mode—characterized by the downwelling of water in the high latitudes of the North Atlantic, which resembles the present oceanic circulation, and the other—the cold mode—characterized by the upwelling of water in the high latitudes of the North Atlantic. Perturbations to the driving parameters of the system, that is, freshwater and heat fluxes in the oceanic boxes, can cause transitions between the two regimes.

Various GCM experiments have shown that multiple equilibria of the THC are possible (Bryan 1986; Manabe and Stouffer 1988; Stouffer et al. 1991; Stocker and Wright 1991; Manabe and Stouffer 1999a; Marotzke and Willebrand 1991; Hughes and Weaver 1994), and most GCMs have shown that the greenhouse gases (GHGs) radiative forcing could cause a weakening of the THC by the inhibition of the sinking of the water in the northern Atlantic. Large increases of the moisture flux and/or of the surface air temperature in the downwelling regions are the driving mechanisms of such a process. Nevertheless, in a recent global warming simulation run the THC remained almost constant since the destabilizing mechanisms have been offset by the increase of the salinity advected into the downwelling region due to net freshwater export from the whole North Atlantic (Latif et al. 2000).

Hemispheric box models along the lines of Stommel’s have been extensively studied, and the stability of their THC has been assessed in the context of various levels of complexity in the representation of the coupling between the ocean and the atmosphere (Weaver and Hughes 1992; Nakamura et al. 1994; Marotzke 1996; Krasovskiy and Stone 1998; Tziperman and Gildor 2002).

However, analysis of GCM data has indicated that the THC is an interhemispheric phenomenon, and in some cases it has been found that the THC strength is approximatively proportional to the density difference between the northern and southern Atlantic (Hughes and Weaver 1994; Rahmstorf 1996; Klinger and Marotzke 1999; Wang et al. 1999a), thus supporting Rooth’s approach.

In this paper we perform an analysis of the stability of the present pattern of the THC against changes to the heat and freshwater atmospheric fluxes using an uncoupled version of the Rooth oceanic model. We note that this model assumes that buoyancy is well mixed in the vertical, and thus avoids the need for parameterizing the mixing (Munk and Wunsch 1998). The model considered here has no real explicit coupling between the ocean and the atmosphere: the atmospheric freshwater fluxes are prescribed and the atmospheric heat fluxes relax the oceanic temperatures toward prescribed values (Rahmstorf 1996; Scott et al. 1999; Stone and Krasovskiy 1999; Titz et al. 2002a,b).

We explicitly analyze the role of the spatial patterns and the rates of increase of the forcings in determining the response of the system. We consider a suitably defined metric and determine what are the thresholds beyond which there is destabilization of the warm mode of the THC. We emphasize that our treatment goes beyond quasi-static analysis since the effect of the rate of forcing is explicitly addressed, so we analyze how the effects on the system of slow (Wang et al. 1999a) and rapid (Wiebe and Weaver 1999) changes contribute. In particular, in the limit of very slow changes, our results coincide with those that can be obtained with the analysis of the bifurcations of the system (Stone and Krasovskiy 1999; Scott et al. 1999). Relatively few studies with more complex models have explicitly addressed how the spatial (Rahmstorf 1995, 1996, 1997 or temporal (Stocker and Schmittner 1997) patterns of forcing determine the response of the system. However, these studies did not carry out a parametric study of the influence of both spatial and temporal patterns. We emphasize that our study is based on a process model that aims more at obtaining a qualitative picture of the properties of the THC than providing quantitative results.

Our paper is organized as follows: in section 2 we provide a description and the general mathematical formulation of the dynamics of the three-box model used in this study. In section 3 we describe the parameterization of the physical processes. In section 4 we analyze the feedbacks of the system. In section 5 we describe the stability properties of the system when it is forced with perturbations in the freshwater flux. In section 6 we describe the stability properties of the system when it is forced with perturbations in the heat flux. In section 7 we present our conclusions.

## 2. Model description

The three-box model consists of a northern high-latitude box (box 1), a tropical box (box 2), and a southern high-latitude box (box 3). The volume of the two high-latitude boxes is the same and is 1/*V* times the volume of the tropical box. We choose *V* = 2 so that box 1, box 2, and box 3 can be thought of as describing the portions of an ocean like the Atlantic north of 30°N, between 30°N and 30°S, and south of 30°S, respectively. At the eastern and western boundaries of the oceanic boxes there is land; our oceanic system spans 60° in longitude so that it covers *ε* = 1/6 of the total surface of the planet. The boxes are 5000 m deep, so that the total mass *M _{i}*

_{=}

_{1,3}=

*M*of the water contained in each of the high-latitude boxes is ≈1.1 × 10

^{20}kg, while

*M*

_{2}=

*V*×

*M*

_{i}_{=}

_{1,3}= 2 ×

*M*. The tropical box is connected to both high-latitude boxes; moreover, the two high-latitude boxes are connected by a deep current passage (which bypasses the tropical box) containing a negligible mass. The three boxes are assumed to be well mixed. The assumption of well-mixed boxes excludes the polar halocline disaster (Bryan 1986; Weaver and Hughes 1992; Zhang et al. 1993) from the range of phenomena that can be described by this model. In our model we have not included a deep box as done in Rahmstorf (1996). This choice allows an easier and more detailed control of the behavior of the system, but implies that distortions are introduced in the response of the system to fast perturbations (Weaver and Hughes 1992). The physical state of the box

*i*is described by its temperature

*T*and its salinity

_{i}*S*; the box

_{i}*i*receives from the atmosphere the net flux of heat

*H̃*

_{i}and the net flux of freshwater

*F̃*

_{i}; the freshwater fluxes

*F̃*

_{i}globally sum up to 0, since the freshwater is a globally conserved quantity. The box

*i*is subjected to the oceanic advection of heat and salt from the upstream box through the THC, whose strength is

*q̃*.

In Fig. 1 we present a scheme of our system in the northern sinking pattern: note that the arrows representative of the freshwater fluxes are arranged in such a way that the conservation law is automatically included in the graph. The dynamics of the system is described by the tendency equations for the heat and the salinity of each box. We divide the three heat tendency equations by *c _{p}* ×

*M*, where

_{i}*c*is the constant pressure specific heat of water per unit of mass, and in the salinity tendency equations we neglect the contribution of the freshwater fluxes in the mass balance, so that virtual salinity fluxes and freshwater fluxes are equivalent (Marotzke 1996; Rahmstorf 1996; Scott et al. 1999). We obtain the following final form for the temperature and salinity tendency equations for the three boxes (Scott et al. 1999):

_{p}where *q* = *ρ*_{0}*q̃*/*M*, *H*_{i} = *H̃*_{i}/(*c*_{p}*M*), and *F*_{i} = *ρ*_{0}*S*_{0}*F̃*_{i}/*M*.

We impose the condition that the average salinity is a conserved quantity; this means that *S̃*_{1} + *VṠ*_{2} + *S̃*_{3} = 0, which implies that

This conservation law holds at every time, and so rules out the possibility of including the effects of the melting of continental ice sheets in our study. We note that the system can (asymptotically) reach an equilibrium only if its feedbacks can drive it toward a state in which the following equation holds:

The strength of the specific THC is parameterized as proportional to the difference between the density of the water contained in the box 1 and the density of the water contained in the box 3, as in the Rooth (1982) model. Given that the water density can approximately be expressed as

where *α* and *β* are the thermal and haline expansion coefficients, respectively; set to 1.5 × 10^{−4} °C^{−1} and 8 × 10^{−4} psu^{−1}, respectively (where psu stands for practical salinity unit), and *ρ*_{0} is a standard reference density, we obtain for the normalized THC strength *q* the following relation:

where *k* is the hydraulic constant, chosen to obtain a reasonable value of the THC strength. We emphasize that the constant *k* roughly summarizes the wind-powered turbulent mixing taking place in the ocean and powering the THC, and we think that, in principle, it should be expressed as a function of the climatic state of the system rather than being a constant.

The northern (southern) sinking state of the circulation is therefore characterized by *q* > (<)0. Considering the case *q* > 0, we obtain a diagnostic relation for steady state *q* by setting to zero Eqs. (3) and (6), and using Eq. (9):

in the case *q* < 0, a diagnostic relation for the steady state *q* can be obtained using the same procedure as above but setting to zero Eqs. (1) and (4):

We conclude that the equilibrium value of the THC strength is determined by the equilibrium values of the heat and freshwater fluxes of the box where we have upwelling of water. These results generalize the expressions given in Rahmstorf (1996). The transient evolution of *q* is determined by its tendency equation:

We observe that the difference between the forcings applied to the freshwater and heat fluxes into boxes 1 and 3 directly affect the evolution of *q*; the presence of terms involving the gradient of temperature and salinity between box 2 and box 1 (box 3) if *q* > 0 (*q* < 0) breaks the symmetry of the role played by the two high-latitude boxes. In our experiments we perturb an initial equilibrium state by changing at various rates the parameters controlling the fluxes *H _{i}* and/or

*F*, and observe under which conditions we obtain a reversal of the THC. The reversal of the THC causes a sudden cooling and freshening in box 1 and a sudden warming and increase of the salinity in box 3, because the former loses and the latter receives the advection of the tropical warm and salty water.

_{i}## 3. Parameterization of the physical processes

We perform our study considering the original formulation of the Rooth (1982) model and thus consider mixed boundary conditions for the surface fluxes. For a more detailed description of the model see also Scott et al. (1999). The heat flux into the box *i* is described by a Newtonian relaxation law of the form *H _{i}* =

*λ*(

_{i}*τ*

_{i}−*T*), where the parameter

_{i}*τ*is the target temperature (Bretherton 1982), which represents a climatologically reasonable value toward which the box temperature is relaxed, and the parameter

_{i}*λ*quantifies the efficacy of the relaxation. Along the lines of Scott et al. (1999), and following the parameterization introduced by Marotzke (1996), we choose

_{i}*τ*

_{1}=

*τ*

_{3}= 0°C and

*τ*

_{2}= 30°C, and select

*λ*

_{1}=

*λ*

_{2}=

*λ*

_{3}=

*λ*= 1.29 × 10

^{−9}s

^{−1}, which mimics the combined effect of radiative heat transfer and of a meridional heat transport parameterized as being linearly proportional to the meridional temperature gradient (Marotzke and Stone 1995; Marotzke 1996; Scott et al. 1999). In physical terms,

*λ*corresponds to a restoring coefficient of ≈4.3 W m

^{−2}°C

^{−1}affecting the whole planetary surface, which translates into an effective

*λ̃*≈ 1/

*ε*× 4.3 Wm

^{−2 }

^{o}C

^{−1}= 25.8 Wm

^{−2 }

^{o}C

^{−1}relative to the oceanic surface fraction only.

Using Eqs. (1)–(3), we can obtain the following tendency equation for the globally averaged temperature *T _{M}* = (

*T*

_{1}+

*VT*

_{2}+

*T*

_{3})/(2 +

*V*):

where *τ _{M}* is the average target temperature; the average temperature relaxes toward its target value in about 25 yr, which is very close to the value considered by Titz et al. (2002b). We observe that in the case of an ocean with vertical resolution, our heat flux parameterization corresponds to a restoring time of ∼3 months for a 50-m-deep surface level (Marotzke and Stone 1995).

The virtual salinity fluxes *F _{i}* are not functions of any of the state variables of the system, and are given parameters. We choose

*F*

_{1}= 13.5 × 10

^{−11}psu s

^{−1}, and

*F*

_{3}= 9 × 10

^{−11}psu s

^{−1}, which correspond to net atmospheric freshwater fluxes of ∼0.40 S

*υ*(1 Sv ≡ 10

^{6}m

^{3}s

^{−1}) into box 1 and of ∼0.27 Sv into box 3, respectively. These values are based on the analysis of the hydrology of the two hemispheres presented by Baumgartner and Reichel (1975) and Peixoto and Oort (1992). Our model neglects moisture fluxes associated with the wind-driven circulation. If these were strong enough to reverse the sign of

*F*

_{3}, the model’s behavior would by very different (Rahmstorf 1996; Titz et al. 2002b), but would not be in accord with observations.

Choosing *k* = 1.5 × 10^{−6} s^{−1}, at equilibrium we have *q* ≈ 1.47 × 10^{−10} s^{−1} (Scott et al. 1999). In physical terms this value describes a THC in the northern sinking pattern with a strength of ≈15.5 Sv; this value agrees reasonably well with estimates coming from observations (Roemmich and Wunsch 1985; Macdonald and Wunsch 1996). We emphasize that this THC strength introduces an oceanic time scale *t _{s}* ≈ 250 yr, which is the ratio between the volumes of box 1 and 3 and the equilibrium state water transport. We refer to this time scale as the flushing or advective time of the ocean since it is the time it takes for the ocean circulation to empty out box 1 or box 3, as well as to propagate the information from one box to another. We emphasize that, since we have considered the approximation of well-mixed boxes, our model should not be considered very reliable for phenomena taking place in time scales ≪

*t*.

_{s}In this equilibrium state, boxes 1 and 3 receive a net poleward oceanic heat advection of ∼1.58 PW and ∼0.17 PW, respectively. The former agrees reasonably well with estimates. The latter agrees with the sign of the global transport in the Southern Hemisphere, but disagrees with the sign in just the southern Atlantic (Peixoto and Oort 1992; Macdonald and Wunsch 1996). The equilibrium temperature of box 1 is larger than the equilibrium temperature of box 3, and the same inequality holds for the salinities, the main reason being that box 1 receives directly the warm and salty tropical water. Taking into account Eq. (9) and considering that the actual equilibrium *q* is positive, we conclude that the THC is haline dominated; the ratio between the absolute value of the thermal and of the haline contribution on the right-hand side of Eq. (9) is ≈0.8. The occurrence of the haline dominated mode depends of course on the forcing parameters, which we have based on observations. The haline mode, with its positive value for q, is consistent with observations. We report in Table 1 the values of the main model constants, while in Table 2 we present the fundamental parameters characterizing this equilibrium state.

## 4. Feedbacks of the system

The Newtonian relaxation law for the temperature implies that the atmosphere has a negative thermic feedback: *T _{i}* increases ⇒

*H*decreases ⇒

_{i}*T*decreases.

_{i}The reaction of the ocean to a decrease in the strength of the THC can be described as follows: *q* decreases ⇒

*T*_{1}decreases more than*T*_{3}⇒*q*increases,the change is limited by the atmospheric feedback;

*S*_{1}decreases more than*S*_{3}⇒*q*decreases.

The heat advection feedback mechanism 1 is negative and dominates the positive salinity advection feedback mechanism 2. We emphasize that the sign of the salinity feedback depends on the asymmetry in the freshwater fluxes (Scott et al. 1999). We note that the strength of the second-order feedback mechanism 1b is relevant in establishing the overall stability of the THC (Tziperman et al. 1994); a very strong temperature-restoring atmospheric feedback like that in our model tends to decrease the stability of the system (Nakamura et al. 1994; Rahmstorf 2000). The feedback mechanisms 1 + 2 are governed by *q* so that their time scale is around the flushing time *t _{s}* ≈ 250 yr (Scott et al. 1999) of the oceanic boxes. We expect that a fast perturbation avoids those feedbacks, since it bypasses the ocean-controlled transfer of information needed to activate them.

In order to simulate global warming conditions, we perform two sets of experiments. In the first, we increase the freshwater fluxes into the two high-latitude boxes. The rationale for this forcing is that an increase in the global mean temperature produces an increase in the moisture capacity of the atmosphere, and this is likely to cause the enhancement of the hydrological cycle, especially in the transport of moisture from the Tropics to the high latitudes. In the second, we represent the purely thermic effects of global warming by setting the increase in the target temperatures of the two high-latitude boxes larger than that of the tropical boxes, since virtually all GCM simulations forecast a larger increase in temperatures in the high latitudes (the so-called polar amplification).

## 5. Freshwater flux forcing

We force the previously defined equilibrium state with a net freshwater flux into the two high-latitude boxes that increase linearly with time; this is obtained by prescribing

We impose the conservation of the globally averaged salinity of the oceanic system by requiring that Eq. (7) holds at all times. When the perturbation is over, we reach a final state that can be either qualitatively similar to the initial northern sinking equilibrium or radically different, that is, presenting a reversed THC. Each perturbation can be uniquely specified by a set of three parameters, such as the triplet 1 [*t*_{0}, Δ*F*_{3}/Δ*F*_{1}, Δ*F*_{1}], or the triplet 2 [*t*_{0}, Δ*F*_{3}/Δ*F*_{1}, *F*^{t}_{1}], where Δ*F*_{i} ≡ *F*^{t}_{i} × *t*_{0} is the total change of the freshwater flux into box *i*. In order to describe a global warming scenario, we explore the effect of positive Δ*F*_{1} and Δ*F*_{3}. Freshwater forcing into box 1 tends to destabilize the THC, since it induces a freshening, and so, a decrease of the density of the water. If the freshwater forcing is larger in box 3, we do not obtain a reversal of the THC, the main reason being that increases in the net freshwater flux into box 3 tend to stabilize the circulation, as can be seen in Eq. (11). Since we are interested in the destabilizing perturbations, we limit ourselves to the case Δ*F*_{3}/Δ*F*_{1} < 1. These perturbations cause an initial decrease in the THC strength *q* and so trigger the following feedbacks: *F*_{1}’s increase is larger than *F*_{3}’s ⇒ *S*_{1} decreases more than *S*_{3} ⇒ *q* decreases.

In Fig. 2 we present two trajectories of *q* describing the evolution of the system from the initial equilibrium state when two different forcings, one subcritical and one supercritical, are applied. The oscillations in the value of *q*, which are slowly decaying in the subcritical case, are caused by the action of the heat and salinity advection feedbacks, which act on time scales of the order of the flushing time of the oceanic boxes (*t _{s}* ≈ 250 yr). In the case of the supercritical perturbation, the oscillations of the value of

*q*grow in size until the system collapses to a southern-sinking equilibrium. The two trajectories separate ∼400 yr after the beginning of the perturbation, when the negative feedback due to the decrease of advection of warm, salty water in box 1 becomes of critical relevance. The threshold value for

*q*is of the order of 10 Sv; more complex models have shown that weak THCs are not stable (Rahmstorf 1995; Tziperman 2000). When the THC does not change sign, it recovers to its initial unperturbed value (if Δ

*F*

_{3}= 0) or overshoots it (Δ

*F*

_{3}> 0), as can be deduced from (11); this reminds us of results obtained with much more complex models (Wiebe and Weaver 1999).

### a. Hysteresis

The hysteresis describes the response of the THC against quasi-static freshwater flux perturbations, that is, with a very large *t*_{0}, in both boxes 1 and 3. In Fig. 3 we have plotted four curves, corresponding to Δ*F*_{3}/Δ*F*_{1} = [0, 0.1, 0.2, 0.3]. The obtained graphs essentially represent the stationary states for given values of the freshwater flux. For a given curve, each point (*x, y*) plotted represents a stable state of the system having *q* = *y* and *F*_{1} = *x*. We scale all the graphs so that the common initial equilibrium state is the point (1, 1). For a range of *F*_{1}, which depends on the value of Δ*F*_{3}/Δ*F*_{1}, the system is bistable; that is, it possesses two distinct stable states, one having *q* > 0, the other one having *q* < 0. The history of the system determines which of the two stable states is realized. The boundaries of the domain in *F*_{1} where the system is bistable correspond to subcritical Hopf bifurcations (Scott et al. 1999) that drive the system from the northern (southern)-sinking equilibrium to the southern (northern)-sinking equilibrium if *F*_{1} crosses the right (left) boundary of the domain of bistability. If Δ*F*_{3}/Δ*F*_{1} is large enough, slow increases of the freshwater fluxes do not destabilize the system, so that no bifurcations are present. It is interesting to note that if *F*_{3} is kept constant, the value of *q* does not depend on the value of *F*_{1} in the northern-sinking branch, while clearly the dependence is apparent in the southern-sinking branch. This can be explained by considering that at equilibrium, the value of *q* is mainly determined by the value of the freshwater flux into the upwelling box, as shown in Eq. (11). We note that since these hysteresis experiments involve simultaneous perturbations to the moisture fluxes in both hemispheres, they differ from those presented in Rahmstorf (1996), where the perturbations to the moisture fluxes occurred in one high-latitude box.

Albeit with differences in the experimental settings in some cases, the presence of a hysteresis for slow freshwater flux perturbations has been previously presented both for simple (Rahmstorf 1996; Scott et al. 1999; Titz et al. 2002b) and more complex interhemispheric models (Stocker and Wright 1991; Mikolajewicz and Maier-Reimer 1994; Rahmstorf 1995, 1996, while the inclusion of changes in the freshwater flux into both high-latitude boxes as presented here had not been previously considered.

### b. Critical forcings

Figures 4a and 4b present, using the coordinates described by the triplets 1 and 2, respectively, the contours of the height of the manifold dividing the subcritical from the supercritical forcings. We consider logarithmic scales instead of linear scales since the former allow us to capture more clearly the qualitative behavior of the system. In Fig. 4a we observe that, as expected, the lower the value of the ratio Δ*F*_{3}/Δ*F*_{1}, the lower the total change Δ*F*_{1} needed to obtain the reversal of the THC. For a given value of the ratio Δ*F*_{3}/Δ*F*_{1}, more rapidly increasing perturbations (larger *F*^{t}_{1}) are more effective in disrupting the circulation. For values of Δ*F*_{3}/Δ*F*_{1} ≤ 0.4 we see a changeover between a slow and a fast regime, geometrically described by a relatively high-gradient region dividing two portions of the manifold having little *x* dependence. Figure 4b shows more clearly that for large values of Δ*F*_{3}/Δ*F*_{1} (≥0.5), the collapse of the THC occurs only for fast increases because of the presence of a threshold in the rate of increase depicted by the little *t*_{0} dependence of the manifold. The changeover corresponds to *t*_{0} ≈ *t _{s}*, which essentially matches the flushing time of the oceanic boxes (Scott et al. 1999). For low values of Δ

*F*

_{3}/Δ

*F*

_{1}, the disruption of the initial THC pattern takes place even in the case of a very low

*F*

^{t}

_{1}, as essentially captured in the hysteresis study of the previous subsection. This occurs because, as shown in the analysis performed by Scott et al. (1999), there is a critical ratio

*F*

_{3}/

*F*

_{1}below which the northern-sinking equilibrium becomes unstable. Such a ratio can be reached after a long enough

*t*

_{0}if the slow perturbation is characterized by a low enough value of Δ

*F*

_{3}/Δ

*F*

_{1}.

### c. Relevance of the mixed boundary conditions

In order to explore the relevance of the mixed boundary conditions (BCs) for the stability of the THC pattern against freshwater flux perturbations (Tziperman et al. 1994; Nakamura et al. 1994), we have also considered a second version of the model where we replace the expression of the heat flux in the box *i* defined as *λ*(*τ** _{i}* −

*T*) with a free parameter

_{i}*H*, thus realizing a flux BC model. Choosing the value of the initial

_{i}*H*equal to the value of

_{i}*λ*(

*τ*

*−*

_{i}*T*) of the previous version of the model at equilibrium, we obtain the same equilibrium solution. We then perturb this equilibrium solution with freshwater flux forcing as in the previous case. Since the heat flux

_{i}*H*does not depend on

_{i}*T*, the positive feedback (condition 1b) described at the beginning of the section is off; therefore we expect this version of the system to be more stable than the previous one. Repeating the analysis of stability to freshwater flux perturbations for the flux BC model, it is indeed confirmed that temperature-restoring conditions decrease the stability of the model, in agreement with the results of Tziperman et al. (1994), Nakamura et al. (1994), and Rahmstorf (2000).

_{i}Comparing the hysteresis graphs obtained for the mixed and flux BC systems presented in Fig. 5, it is apparent that the mixed BC model is less stable with respect to freshwater flux perturbations. Excluding the temperature feedback mechanism, we do not obtain collapses of the THC for quasi-static perturbations for Δ*F*_{3}/Δ*F*_{1} = 0.3, while for each of the cases Δ*F*_{3}/Δ*F*_{1} = [0, 0.1, 0.2] the bifurcation point is much farther than in the mixed BC model from the point (1, 1), which describes the unperturbed system.

Figure 6 shows the ratio between the values of the critical values of Δ*F*_{1} presented in Fig. 4 and the corresponding values obtained with the flux BC model. The average value of this ratio over all the domain is ≈0.3. Observing Fig. 6 for Δ*F*_{3}/Δ*F*_{1} = [0.3, 0.4], we see that the ratio goes to zero with increasing *t*_{0} because in this subdomain for slow forcings we have destabilization for finite values of Δ*F*_{1} only in the mixed BC case. This occurs because for these values the threshold in the rate of increase is present only in the flux BC model; therefore in this model for large *t*_{0}, the quantity Δ*F*_{1} diverges, while in the mixed BC model it is finite.

## 6. Thermal forcing

We explore the stability of the northern-sinking equilibrium against perturbations to the initial heat fluxes expressed as changes in the target temperatures in the three boxes. We do not impose a constraint fixing the average target temperature *τ _{M}*, so that in this case not two but three parameters

*τ*

_{1},

*τ*

_{2},

*τ*

_{3}can be changed. It is possible to recast the system of Eqs. (1)–(6) (Scott et al. 1999) so that the heat forcing is described completely by the changes in

*τ*=

_{N}*τ*

_{1}−

*τ*

_{2}and in

*τ*=

_{S}*τ*

_{3}−

*τ*

_{2}. Therefore, we set

*τ*

_{2}=

*τ*

_{2}(0) at all times, so that

*τ*

^{t}

_{1}=

*τ*

^{t}

_{N}and

*τ*

^{t}

_{3}=

*τ*

^{t}

_{S}. We apply the following forcings:

We characterize each forcing experiment, which leads to a final state that is characterized by a northern- or a southern-sinking equilibrium, using the parameter set 1 defined as [*t*_{0}, Δ*τ*_{3}/Δ*τ*_{1}, Δ*τ*_{1}] and the parameter set 2 defined as [*t*_{0}, Δ*τ*_{3}/Δ*τ*_{1}, *τ*^{t}_{1}]; consistently with the previous case, we define Δ*τ*_{i} ≡ *τ*^{t}_{i}*t*_{0}, that is, the total change of the target temperature of the box *i*. In order to imitate global warming scenarios where polar amplification occurs, we consider Δ*τ _{i}* ≥ 0,

*i*= 1, 3. We observe that if Δ

*τ*

_{3}> Δ

*τ*

_{1}, no destabilization of the THC can occur because such a forcing reinforces the THC. The specified forcings trigger the previously described feedbacks of the system so that:

*τ*

_{1}’s increase is larger than

*τ*

_{3}’s ⇒

*Η*

_{1}’s increase is larger than

*H*

_{3}’s ⇒

*T*

_{1}increases more than

*T*

_{3}⇒

*q*decreases.

We present in Fig. 7 two trajectories of *q* starting from the initial equilibrium state, characterized by slightly different choices of the forcing applied to *τ*_{1} and *τ*_{3}. One of the forcings is supercritical, the other one is subcritical. The two trajectories are barely distinguishable up to the end of the perturbation, because the forcing dominates the internal feedbacks of the system thanks to the large value of *λ*. The oscillations in the value of *q* at the end of the perturbation are caused by the internal feedbacks of the system and have a period comparable with the time scale of the flushing of the oceanic boxes. In general, the mean flow dampens the increase in *T*_{1} and tends to keep the system in the northern-sinking state, although the decreasing *q* lowers the amount of warm, salty tropical water injected into box 1. When *q* gets too small (in this case about 4 Sv) the mean flow feedback becomes negligible and the system eventually collapses.

### a. Hysteresis

We present in Fig. 8 only the graphs relative to Δ*τ*_{3}/Δ*τ*_{1} = [0, 0.2, 0.5, 0.8] for the sake of simplicity; the other curves that can be obtained for other values of 0 ≤ Δ*τ*_{3}/Δ*τ*_{1} < 1 are not qualitatively different. The initial state is the point (0, 1). The subcritical Hopf bifurcation points are (∼5, ∼0.7), (∼6, ∼0.7) (∼9, ∼0.7), and (∼18, ∼0.6), respectively. This means that the bifurcation occurs when the THC has already declined down to ∼10 Sv. The bistable region is extremely large in all cases, because the *x* of the bifurcation points in the lower part of the circuits are below −50°C; they become farther and farther from the initial equilibrium value as Δ*τ*_{3}/Δ*τ*_{1} increases. This means that once the circulation has reversed, the newly established pattern is extremely stable and can hardly be changed again. To our knowledge, this is the first time such an analysis of hysteresis graphs for this class of forcing is presented in the literature.

### b. Critical forcings

Figure 9 presents the contours of the height of the manifold of the critical forcings. We can observe that for all the values of Δ*τ*_{3}/Δ*τ*_{1} there is a high-gradient region dividing two flat regions describing the fast and slow regimes. In the fast region, a smaller total increase Δ*τ*_{1} is required to achieve the reversal of the THC. For Δ*τ*_{3}/Δ*τ*_{1} = 0.9, which represents the case of highly symmetrical forcings, the critical Δ*τ*_{1} essentially does not depend on *t*_{0}. We emphasize that the time scale of such a changeover is larger than *t _{s}* ≈ 250 yr, the main reason being that because of the large value of

*λ*, which is about one order of magnitude larger than the initial

*q*, the forcing is stronger than any negative feedback even if it is relatively slow. For the same reason, in the case of perturbations to the target temperatures, we do not observe the presence of thresholds for the rate of increase of the perturbations as seen for freshwater forcings. If the increase in

*τ*

_{1}is slow enough, the negative feedbacks make the destabilization more difficult, but do not prevent it.

## 7. Conclusions

In this paper we have accomplished a thorough analysis of the stability of the THC as described by the Rooth model. We simulate global warming conditions by applying to the equilibrium state perturbations to the moisture and heat fluxes into the three boxes. In particular, the relevance of the role of changes in the freshwater flux into the Southern Ocean is determined, along the lines of Scott et al. (1999) for a similar box model and Wang et al. (1999a) for an OGCM.

We have presented the hysteresis graphs for freshwater flux perturbations that reproduce and extend the results previously given in the literature, because the effect of freshwater flux increases in the Southern Ocean are included, while for the first time hysteresis experiments for heat flux forcings, obtained by perturbation in the target temperature, are presented.

High rates of increase in the moisture flux into the northern high-latitude box lead to a THC breakdown at smaller final increases than low rates, while the presence of moisture flux increases into the southern high-latitude box strongly inhibit the breakdown. Similarly, fast heat flux increases in the Northern Hemisphere destabilize the system more effectively than slow ones, and again the enhancement of the heat fluxes in the Southern Hemisphere tend to drive the system toward stability. In all cases analyzed, slow forcings, if asymmetric enough, lead to the reversal of the THC.

In our model, the spatial pattern of forcing is critical in determining the response of the system: quasi-symmetric forcings, which seem more reasonable under realistic global warming conditions, are not very effective in destabilizing the THC. The simplicity of the feedbacks characterizing our model, and in particular the absence of stratification effects, may be responsible for this behavior.

We have shown that the restoring boundary conditions for temperature decrease the stability of the system because they introduce an additional positive feedback that tends to destabilize the system. Flux boundary conditions provide the system with much greater stability with respect to perturbations in the freshwater fluxes at all time scales.

Our study is based on a conceptual model; therefore it aims at obtaining a qualitative picture of the properties of the THC more than providing realistic results. Nevertheless, future increases in computer power availability could allow the adoption of this methodology for performing extensive parametric studies of the THC stability properties and evolution with more complex and realistic models in both paleoclimatic and global warming perspectives.

We conclude that the study of the THC response to perturbations performed with uncoupled models, of any level of complexity, should take into much greater account the spatial pattern of the freshwater and/or heat flux forcing. The THC is a highly nonlinear, nonsymmetric system, and we suggest that the effect of changing the rate of increase of the forcings should be explored in greater detail to provide a bridge between the instantaneous and quasi-static changes. The THC dynamics encompass very different time scales, which can be explored only if the timing of the forcings is varied.

We conclude by pointing out some possible improvements to the present model and possible extensions of this work. We have put emphasis on how the oceanic advection can reverse the symmetry of the pattern of the circulation in the context of a rigorously symmetric geometry. The system is likely to be very sensitive to changes in the volumes of the boxes and especially to the introduction of asymmetries between the two high-latitude boxes, which would also make the system more realistic since the southern mid- to high-latitude portion of the Atlantic is considerably larger than the northern mid- to high-latitude portion (Zickfeld et al. 2004). Asymmetries in the oceanic fractional areas would induce asymmetries in the values of *λ _{i}*, thus causing the presence of different restoring times for the various boxes.

The introduction of a noise component in the tendency equation would increase the model’s realism and make the model conceptually more satisfying since it would introduce a parameterization of the high-frequency variability. Recent studies have undertaken this strategy and shown that close to instability threshold the evolution of the THC has a very limited predictability (Wang et al. 1999a; Knutti and Stocker 2002), and that stochastic resonance could be responsible for glacial/interglacial climate shifts (Ganopolski and Rahmstorf 2002). It would be interesting to analyze how the intensity and color of the noise would influence the results obtained and the conclusions drawn in this work.

Finally, the model could be improved so that it would be able to represent the main features of the whole conveyor belt by adding other boxes, descriptive of other oceanic basins, and carefully setting the connections between the boxes, along the lines of some examples presented by Weaver and Hughes (1992). Moreover, the consideration of a deep ocean box would surely increase the model performance in terms of representation of the short time scales.

In an extension of this paper (Lucarini and Stone 2005), we perform a similarly structured stability analysis of the THC under global warming conditions using a more complex version of the Rooth model, where an explicit coupling between the oceanic and the atmospheric part is introduced and the main atmospheric physical processes responsible for freshwater and heat fluxes are formulated separately.

## Acknowledgments

The authors are grateful to J. Scott for interesting discussions and useful suggestions. One author (V.L.) wishes to thank R. Stouffer for having proposed improvements to an earlier version of the manuscript, and T. Stocker and E. Tziperman for having suggested a number of relevant references. This research was funded in part by the U.S. Department of Energy’s (DOE’s) Climate Change and Prediction Program and in part by the MIT Joint Program on the Science and Policy of Global Change (JPSPGC). Financial support does not constitute an endorsement by DOE or JPSPGC of the views expressed in this article.

## REFERENCES

**,**

_{2}upset the current balance?

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{2}on the ocean–atmosphere system.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{2}emission rates on the stability of the thermohaline circulation.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

* Current affiliation: Department of Mathematics and Computer Science, University of Camerino, Camerino, Italy

*Corresponding author address:* Dr. Valerio Lucarini, Dept. of Mathematics and Computer Science, University of Camerino, via Madonna delle Carceri 62032, Camerino, Italy. Email: lucarini@alum.mit.edu