Heaving modes in the world oceans

Part of climate changes on decadal time scales can be interpreted as the result of adiabatic motions associated with the adjustment of wind-driven circulation, i.e., the heaving of the isopycnal surfaces. Heat content changes in the ocean, including hiatus of global surface temperature and other phenomena, can be interpreted in terms of heaving associated with adjustment of wind-driven circulation induced by decadal variability of wind. A simple reduced gravity model is used to examine the consequence of adiabatic adjustment of the wind-driven circulation. Decadal changes in wind stress forcing can induce three-dimensional redistribution of warm water in the upper ocean. In particular, wind stress change can generate baroclinic modes of heat content anomaly in the vertical direction; in fact, changes in stratification observed in the ocean may be induced by wind stress change at local or in the remote parts of the world oceans. Intensification of the equatorial easterly can induce cooling in the upper layer and warming in the subsurface layer. The combination of this kind of heat content anomaly with the general trend of warming of the whole water column under the increasing greenhouse effect may offer an explanation for the hiatus of global surface temperature and the accelerating subsurface warming over the past 10–15 years. Furthermore, the meridional transport of warm water in the upper ocean can lead to sizeable transient meridional overturning circulation, poleward heat flux and vertical heat flux. Thus, heaving plays a key role in the oceanic circulation and climate.


Introduction
The oceanic general circulation consists of two major components: the wind-driven circulation in the upper kilometer and the thermohaline circulation over the entire depth of the ocean. Thermohaline circulation is intimately related to the thermohaline forcing on the air-sea interface; the mechanical energy requited for sustaining the thermohaline circulation comes from tidal dissipation, wind stress, and geothermal heat flux. However, wind-driven circulation is not directly linked to surface thermohaline forcing. In fact, wind-driven circulation in the ocean can be idealized as the adiabatic motions of sea water. Since motions in the upper ocean are mostly dominated by wind-driven circulation, many changes in the upper ocean can be understood in terms of the adiabatic motions associated with wind-driven circulation.
Climate variability is intimately related to many critically important aspects of our society; thus, it is desirable to identify the sources of climate variability. However, climate variability is complicated because it involves many dynamical processes. These include: the low-frequency variability in the atmosphere (such as the solar isolation), wind stress, air-sea heat flux and freshwater flux, plus internal variability of the climate system itself.
Abstract Part of climate changes on decadal time scales can be interpreted as the result of adiabatic motions associated with the adjustment of wind-driven circulation, i.e., the heaving of the isopycnal surfaces. Heat content changes in the ocean, including hiatus of global surface temperature and other phenomena, can be interpreted in terms of heaving associated with adjustment of winddriven circulation induced by decadal variability of wind. A simple reduced gravity model is used to examine the consequence of adiabatic adjustment of the wind-driven circulation. Decadal changes in wind stress forcing can induce three-dimensional redistribution of warm water in the upper ocean. In particular, wind stress change can generate baroclinic modes of heat content anomaly in the vertical direction; in fact, changes in stratification observed in the ocean may be induced by wind stress change at local or in the remote parts of the world oceans. Intensification of the equatorial easterly can induce cooling in the upper layer and warming in the subsurface layer. The combination of this kind of heat content anomaly with the general trend of warming of the whole water column under the increasing greenhouse effect may offer an explanation for the hiatus of global surface temperature and the accelerating subsurface warming over the past 10-15 years. Furthermore, the meridional transport of warm water in the upper ocean can lead to sizeable transient meridional

3
Many studies have been focused on the internal variability of the atmosphere-ocean coupled system, such as the basin modes (Cessi and Paparella 2001;Cessi and Primeau 2001), the oceanic response to stochastic wind (Cessi and Louazel 2001) and the oceanic teleconnections (Cessi and Otheguy 2003). Johnson and Marshall (2004) explored the global teleconnections of meridional overturning anomaly. They came to a conclusion that meridional overturning circulation (MOC hereafter) anomaly on decadal and shorter time scales are confined to the hemisphere basin where the perturbations originated. Thus, the simple ocean model used in our study may be able to catch the major part of the consequence induced by wind stress anomaly. The impact of thermohaline forcing anomaly on the thermohaline circulation was studied by many investigators, e.g., Johnson and Marshall (2002a, b), Zhai and Sheldon (2012), Zhai et al. (2013); however, this is not our focus in this study.
Due to the existence of the periodic channel and the relatively small radius of deformation, the dynamics of the Antarctic circumpolar current (ACC) is quite complicated. Up till now, a clear theory for the ACC remains a great challenge, e.g., Rintoul et al. (2001). This is still a research frontier; for example, recently the spin-up and adjustment of the ACC has been examined by Allison et al. (2011).
In general, climate anomalous signals in the world oceans can be classified into three basic categories: warming (cooling), freshening (salinification), and heaving. Isopycnal motions subjected to no heat and salinity exchange with the environment are called heaving, and such motions are induced by adjustment of wind-driven circulation, e.g., Bindoff and McDougall (1994); and they can be idealized as the movement of the main thermocline (called thermocline hereafter).
In particular, wind stress is one of the most variable components in the climate system. Due to the change of wind forcing wind-driven circulation also adjusts in response. Wind stress changes take place in many different parts of the world oceans, and the corresponding adjustments of wind driven circulation evolve with time in rather complicated ways. As a result, identifying individual adjustment process and its source from climate datasets for the world oceans can be a great challenge. Instead, we set a modest goal in this study to the examination of fundamental heaving modes in the world oceans.
Recently, hiatus in global warming has become a hotly debated issue related to global climate change. Over the past several decades, the global sea surface temperature (SST) kept increasing; in addition, there are decadal variability of global ocean heat content, e.g., Levitus et al. (2005Levitus et al. ( , 2009Levitus et al. ( , 2012, Easterling and Wehner (2009), Katsman and van Oldenborgh (2011), Kaufmann et al. (2011), Guemas et al. (2013), Watanabe et al. (2013). However, ocean warming is far from uniform. For example, in the North Atlantic, the tropics/subtropics have warmed, but that the subpolar ocean has cooled. These changes can be linked to NAO, and they are also directly linked to the gyre scale circulation changes and the MOC (Lozier et al. 2008(Lozier et al. , 2010. In particular, over the past 10-15 years, change in global SST is more or less leveled off, i.e., a hiatus in global SST record; on the other hand, subsurface heat content keeps increasing, with seemingly even high rate, e.g., Meehl et al. (2011Meehl et al. ( , 2013, Balmaseda et al. (2013), Kosaka and Xie (2013) and Chen and Tung (2014).
As discussed in previous studies, the hiatus of global SST may be due to many mechanisms. For example, stratospheric water vapor and aerosol layer may contribute, e.g., Solomon et al. (2010Solomon et al. ( , 2011. However, our focus in this study is to explore the possible mechanism directly linked to oceanic circulation, in particular the linkage to the variability in wind-driven circulation. In fact, early studies suggested that this may be an important connection, e.g., McGregor et al. (2012), Jones et al. (2011), England et al. (2014). As will be explained in this study, hiatus of global SST record in combination of accelerating subsurface heat content increasing may be explained in terms of the general trend of warming of the whole water column and stratification changes induced by the adjustment of winddriven circulation in response to global decadal wind stress perturbations.
In this study we will use idealized geometry and simple wind stress perturbations and focus on the dynamical consequence of large-scale adjustment of wind-driven gyres. One of the most important consequences of such adjustment is the basin scale quasi-horizontal transport of water mass. We assume that such movements take place within a relatively short time, on the order of inter-annual and decadal time scales. Neglecting the contributions due to the surface thermohaline forcing and internal diapycnal diffusion, such movements can be idealized as adiabatic, and they are commonly called heaving. Heaving can induce changes in the basin mean vertical stratification and the corresponding mean vertical heat content profile.
For example, change in wind induced Ekman pumping can lead to a shifting of warm water in the upper ocean and thus baroclinic modes of heat content anomaly, as illustrated in Fig. 1. When the Ekman pumping rate is enhanced, the slope of the thermocline increases, leading to a three dimensional redistribution of warm water in the basin (red curve and arrows in Fig. 1a). On the x-z plane, the anomalous circulation appears in the form of a zonal overturning cell rotating anticlockwise. In particular, this induces a westward shifting of heat content and a downward shifting of heat content, as depicted in Fig. 1b, c. It is readily seen that if the wind induced Ekman pumping rate is weakened, an opposite process should take place.
Similarly, if we look at the meridional plane, there is an anomalous MOC induced by wind stress anomaly, and this should give rise to a meridional shifting of heat content and the associated poleward heat flux.
Since our model is adiabatic, heat content in the whole basin has zero net gain or loss; thus, in the vertical direction heat content anomaly must appear in the form of baroclinic modes. These baroclinic modes of heat content anomaly are entirely due to the adiabatic adjustment of wind-driven circulation in the ocean, i.e., heaving.
The basic idea discussed above can be extended to a model ocean with multiple gyres. For an idealized twohemisphere basin, there are several wind-driven gyres, including the subpolar gyres, the subtropical gyres, and the equatorial gyres; the fundamental structure of the thermocline in a quasi-steady state is sketched by black curves in Fig. 2. Assuming the amount of warm water in the upper ocean remains unchanged, during the adjustment winddriven circulation warm water from one gyre is redistributed to other gyres in the model ocean.
For example, if the equatorial easterly relaxes, the equatorial thermocline moves up in response. The upward movement of the equatorial thermocline leads to the transportation of warm water in the upper ocean toward middle/ high latitudes, red arrows in Fig. 2a. Neglecting the small change in sea level, the total water column height at each station is nearly constant. Hence, in compensation, cold water in the lower layer moves toward the Equator (blue arrows) to fill up the space left behind by the poleward transport of warm water.
The movement of water in the upper and lower layers implies two very important physical processes. First, the poleward flow in the upper layer and the equatorward flow in the lower layer give rise to the anomalous MOC and poleward heat flux. Second, the stratification in the ocean is changed due to such exchange. Since the wind-driven circulation is considered as adiabatic, the total heat content for the whole ocean must be constant. Consequently,

Fig. 2
Symmetric and asymmetric modes of heaving induced motions in a two-hemisphere basin. Black arrows indicate the vertical movement of thermocline induced by wind stress changes; black curves depict the undisturbed thermocline and the red curves indicate the disturbed thermocline; red arrows indicate the movement of warm water above the thermocline and blue arrows indicate the movement of cold water below the thermocline basin mean vertical heat content anomaly must appear in the form of baroclinic modes. As shown in Fig. 1, in a single gyre basin heaving motions can induce a first baroclinic mode; for a two-hemisphere basin, the situation becomes more complex, and the second or even higher baroclinic modes can be generated as well.
The second example is for the case with wind stress forcing change in the subtropical basin of the Northern Hemisphere, Fig. 2b. If Ekman pumping in the subtropical basin is reduced, the thermocline in the subtropical basin shoals; this leads to the warm water transport toward both high and lower latitudes, depicted by the red arrows. In compensation, cold water in the lower layer flows toward the subtropical basin of the Northern Hemisphere, depicted by the blue arrows.
Similarly, if wind forcing in the subpolar basin in the Northern Hemisphere is enhanced, the cyclonic gyration is intensified, and the thermocline moves upward. The cyclonic gyre can no longer hold up the same amount of warm water; hence, warm water is transported from the northern subpolar basin southward to the rest of the basin. As a result, thermocline in other parts of the basin deepens in response, Fig. 2c.
In addition, wind stress in both hemispheres can change simultaneously. Such changes can be idealized in terms of the symmetric and asymmetric modes induced by symmetric and asymmetric wind stress perturbations, sketched in Fig. 2d, e.
From the point view of global water mass and heat distribution, one of the most important consequences associated with such basin-scale adjustment of wind-driven circulation is the meridional (vertical) redistribution of mass and heat, which can contribute to a substantial portion of the variability of the meridional (vertical) transport of mass and heat. Many relevant and interesting phenomena will be explored in this study.
The rest of this paper is organized as following. In Sect. 2, the formulation of a two-hemisphere reduced gravity model is presented, and the results from a series of numerical experiments forced by idealized wind stress perturbations will be discussed. In Sect. 3, the formulation of a Southern Hemisphere ocean model is presented and the results from a series of numerical experiments forced by idealized wind stress perturbations will be discussed. Finally, we conclude in Sect. 4.

Model set up
With the focus on the wind-driven circulation mostly confined to above the thermocline, the density structure in the ocean can be idealized as a step function in density coordinates. This two-hemisphere (labeled as 2H hereafter) model ocean consists of two layers: the upper (lower) layer has a constant density of ρ 0 − Δρ (ρ 0 ). The upper layer thickness is denoted as h, and the lower layer is infinitely deep and thus motionless.
The model is based on the rigid-lid approximation. As such, the effect of free surface ζ � = 0 is replaced by a nonconstant hydrostatic pressure p = p a at the flat surface z = 0. The pressure in layers beneath can be calculated using the hydrostatic relation.
The pressure gradient for a one-and-a-half layer model is where g ′ = g�ρ/ρ 0 is the reduced gravity. In approximation, the sea surface level is linked to the upper layer thickness Thus, the sea level can be inferred from the thermocline depth, with an unknown constant which can be determined by requiring the basin-integrated sea surface level is zero, i.e., Since the model has one active layer only, the time dependent momentum and continuity equations for the active upper layer are where (τ x , τ y ) are the zonal and meridional wind stress, (κ, A m ) are the parameters for the vertical and horizontal momentum dissipation, which are set to κ = 0.0015 m/s, A m = 1.5 × 10 4 m 2 /s. Note that an interfacial friction linearly proportional the velocity shear is used, which can be interpreted as a crude parameterization of baroclinic instability. There are many different ways of parameterizing the eddy's effects. For example, Greatbatch (1998) and Greatbatch and Lamb (1990) postulated a more sophisticated formulation for parameterizing the vertical momentum dissipation associated with eddies, which is essentially equivalent to the now commonly used lateral mixing induced by meso-scale eddies of Gent and McWilliams (1990).
The model ocean is formed on an equatorial beta-plane, with the Coriolis parameter defined as the linear function of latitude The 2H model is 150° wide in the zonal direction, and extends from 70°S to 70°N. Following the common practice of non-eddy resolving modeling, the model is based on the B-grid, with 152 × 142 grids and 1° × 1° resolution, which grid size is set to 110 km. In this model we chose β = 2.2367 × 10 −11 /m/s, so that the Coriolis parameter on a grid point at 35°N equals to f 35 • N = 8.36552 × 10 −5 /s, corresponding to the Coriolis parameter at this latitude in the spherical coordinates.
In this study, we will be focused on the role of zonal wind stress only. The undisturbed zonal wind stress profile (in N/m 2 ) applied for the 2H model is where θ is the latitude. The numerical model is based on the traditional leap-frog scheme with a time step of 3153.6 s for both the momentum equations and the continuity equation. The detail of numerical model can be found in Huang (1987).
To prevent the upper layer outcropping in the subpolar basin, a relatively large amount of warm water is specified in the initial state of rest which corresponds to a mean layer depth of 350 m. The reduced gravity model used in this study makes use of observations. According to WOA09 data (Antonov et al. 2010), the annual mean potential density referred to the sea surface, averaged over depth of 0-300 m for the global oceans, is estimated at σ 0 = 25.96 kg/m 3 , but the corresponding value over depth of 400-5500 m is estimated at σ 0 = 27.66 kg/m 3 ; thus, g ′ ≈ 1.7 cm/s 2 . Hence a round off value of g ′ = 1.5 cm/s 2 will be used in this model. The annual mean potential temperature averaged over depth of 0-300 m is 13.36 °C, and the corresponding value over depth of 400-5500 m is 2.555 °C; thus the temperature difference between the upper and lower layers is 10.8 °C.
Hence, in this model the temperature difference is set at 10 °C. The model was run for 300 years to reach a quasi-equilibrium reference state. This reference state is symmetric to the Equator, and it has three gyres in each hemisphere. The wind stress, the thermocline depth and the streamfunction of the reference state are shown in Fig. 3.
To explain the vertical profile of heat content anomaly, we also plot a meridional profile of the thermocline depth 7 grid points (770 km) east of the western boundary, which can represent the layer depth maximum or minimum as a function of latitude, Fig. 3b. The relatively wide western boundary layer is due to the low horizontal resolution and large frictional parameters used in our model.
In the subtropical basin, the maximum depth of the thermocline of the anticyclonic gyre is approximately 603 m, mimicking the situation in the North Pacific Ocean. The strength of the subtropical gyre is about 24 Sv (1 Sv = 10 6 m 3 /s), somewhat weaker than in the North Pacific Ocean. This relatively weak gyre is due to the idealized wind stress profile used in the model. In the subpolar basin, there is weak cyclonic gyre, giving rise to a domeshaped thermocline. Our main focus is to explore the fundamental structure induced by decadal variation of small amplitude wind stress perturbations; thus, the choice of the model parameters and wind stress profile should not qualitatively affect the main results from this study.

Zonal wind stress variability in the world oceans
Wind stress in the world oceans varies over a broad spectrum in space and time. As an example, zonal wind stress in the central Pacific Ocean is shown in Fig. 4, taken from the GODAS data provided by the NOAA/ OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at http://www.esrl.noaa.gov/psd/. As shown in Fig. 4b, c, on decadal time scales zonal wind stress in the central Pacific varies with amplitude on the order of 0.015-0.03 N/m 2 . Thus, it is reasonable to use decadal wind stress perturbations on the order to 0.015-0.02 N/ m 2 in numerical experiments exploring the dynamical consequences of heaving.

Numerical experiments
From the reference state, we carried out a series of numerical experiments. In each experiment, the model was restarted from the quasi-equilibrium reference state and forced by wind stress with small perturbations in the form of a Gaussian profile where �τ = ±0.015 N/m 2 is the amplitude, Δy = 1100 km. The wind stress profiles used in numerical experiments are shown in Fig. 5. Note that because the amplitude of wind stress anomaly is relatively small, perturbations to the wind-driven circulation are almost linearly proportional to the amplitude (including the sign) of wind stress anomaly; thus, the corresponding results for wind stress perturbations with opposite signs can be inferred from results presented in this study.
In each experiment wind stress perturbations were linearly increased from 0 at t = 0 to the specified strength at the end of 20 years. Afterward, the wind stress perturbations were kept constant and the model run for additional 20 years. Such numerical experiments may represent the typical cases for variability induced by decadal wind stress perturbations.

The pivotal case, Exp. 2H-A
In this case the equatorial easterly was enhanced. The adjustment of wind-driven circulation induces a three dimensional redistribution of warm water in the model ocean. The time evolution of the volume anomaly and heat content anomaly is shown in Fig. 6. Due to the intensification of the equatorial easterly, the slope of the equatorial thermocline enhances; the warm water above the thermocline is pushed toward the Equator from both hemispheres. Thus, the amount of warm water in the equatorial band increases with time, but it declines at middle/high latitudes in both hemispheres, Fig. 6a. At the end of the 40 year experiment, the meridional distribution of the volume anomaly is shown in Fig. 6b. A basic assumption made in the reduced gravity model is that the lower layer is infinitely deep; hence, the pressure gradient and velocity in this layer is negligible. In the ocean the poleward mass flux in the upper layer must be compensated by the equatorward return flow in the lower layer; thus, the adjustment of wind-driven gyre in the upper layer should induce an anomalous MOC. The equivalent MOC rate can be diagnosed as follows where x w and x e are the western and eastern boundaries, v is the meridional velocity. By convention, the MOC induced by a northward flow of warm water in the upper layer is defined as positive.
As shown in Fig. 6c, wind-driven circulation adjustment induces a pair of anomalous MOC asymmetric to the Equator, with a maximum rate of more than 0.3 Sv at year 20, when the wind stress perturbations reach the (10) M moc (y, t) = x e x w h(x, y, t)v(x, y, t)dx peak amplitude. Afterward, the anomalous MOC declines gradually. If the numerical experiment were run for much longer time, the model ocean would gradually reach a new quasi-equilibrium state, in which the anomalous MOC vanishes. Since MOC varies so much during the adjustment, it is more meaningful to use the MOC rate averaged over the entire 40 year of the numerical experiment. The corresponding meridional profile of the mean MOC rate, with maximum amplitude of 0.2 Sv, is shown in Fig. 6d.
The MOC in the ocean is primarily associated with surface thermal forcing, in particular with the thermohaline circulation; in a steady state the wind-driven circulation in combination of surface heating/cooling can also contribute to the MOC. Recent studies revealed a close link between the MOC and surface thermohaline and wind forcing. For example, the NAO cycle plays a vital role in generating the variability of the MOC in Atlantic Ocean, e.g., Lozier et al. (2010), Zhai et al. (2014).
However, a major point in our study is that during the adiabatic adjustment of wind-driven circulation anomalous MOC appears, which is not directly linked to the surface thermohaline forcing. Our numerical experiment indicated that, even taking the value averaged over the entire 40 years, the MOC associated with adiabatic adjustment of the wind-driven circulation may consist of a substantial portion of the variable MOC in the world oceans.
The anomalous MOC inferred from the model also gives rise to poleward heat flux, which is defined as where ρ o = 1035 kg/m 3 is the mean reference density, C p = 4186 J/kg/°C is the mean heat capacity under constant pressure, T = 10 °C is the temperature difference where T (x, y, z, t) is the instantaneous temperature and T ref (x, y, z) is the temperature in the reference state. The time evolution of heat content anomaly is shown in Fig. 6e. The heat content anomaly at the end of the 40 year experiment is shown in Fig. 6f.
In the reference state, thermocline near the western boundary at low latitudes is the deepest, approximately 603 m, Fig. 3b. Hence, positive thermocline perturbations at low latitudes induce positive anomaly of the basin-mean heat content at this depth range. On the other hand, thermocline along the eastern boundary and at high latitudes shoals, leading to negative heat content anomaly at the shallow depth, shown in the upper part of Fig. 6f. Since these motions are adiabatic, heat content anomaly must appear in the form of baroclinic modes.
Take the time derivative of the heat content anomaly gives rise to the mean heat content anomaly rate With T = 40 year, this gives rise to the mean rate (per meter) averaged over the 40 year experiment. The vertical integration of this time rate gives rise to the vertical heat flux As shown in Fig. 6g, for the present case the vertical heat flux is on the order to −8 TW (the negative sign indicates a downward shifting of heat content). Dividing by the total area of the model ocean leads to the vertical heat flux per unit area, on the order of 0.03 W/m 2 . It is to emphasize that such vertical heat flux is due to adiabatic adjustment of the water masses in the ocean; if we look through the one-dimensional potential temperature coordinate, there is no change at all.
At the end of the experiment, the horizontal structure of the perturbations to the wind-driven circulation is shown in Fig. 7. The equatorial thermocline deepens, in particular 10° off the Equator and near the western boundary, with maximum amplitude of 15 m, Fig. 7a; on the other hand, thermocline shoals at high latitudes, with a maximum value of −9.8 m. The streamfunction anomaly appears as a pair of gyres asymmetric to the Equator, with maximum amplitude of 2 Sv. The streamfunction of the anomalous circulation in the Northern Hemisphere is positive, i.e., it is an anomalous anticyclonic circulation. As a result, the original anticyclonic circulation, including its western boundary current, is intensified, Fig. 7b.
Within the framework of the reduced gravity model, the corresponding sea level anomaly in the final state can be inferred from the upper layer depth anomaly Δh, Eq. (2). The sea level anomaly at the end of the experiment has the same pattern as the thermocline thickness anomaly. For the present case, a negative equatorial wind anomaly induces a positive sea level anomaly in the western part of the basin at low latitudes, and a negative sea level anomaly in the eastern part of the basin. Observations, e.g., Merrifield (2011), Qiu and Chen (2012), indicate that the sea level anomaly in the Pacific Ocean over the past 10-20 years has the same signs as those shown in Fig. 7c. Hence, the strong positive (negative) sea level anomaly in the western (eastern) North Pacific may be linked to the stronger than normal easterly in the Equatorial Pacific. However, if this strong anomalous easterly is going to relax, the long-lived strong sea level anomaly in the Pacific Ocean may swing in the opposite direction.
It is well known that the adjustment of wind-driven circulation in a closed basin is carried out through wave motions, including the Rossby waves and Kelvin waves. In particular, the role of Kelvin waves and long baroclinic Rossby waves play vital roles in establishing the circulation, e.g., Anderson and Gill (1975), Hsieh et al. (1983), Wajsowicz and Gill (1986), Wajsowicz (1986), Hsieh and Bryan (1996), Marshall and Johnson (2013). Furthermore, for simplified geometry, the adjustment of the global ocean has been studied by Huang et al. (2000), Primeau (2002), Cessi et al. (2004).
In the present case, wind forcing anomaly along the Equator and low latitude band induces a decline of the eastwest slope of the thermocline at low latitudes. These signals reach the western boundary and form the coast-trapped Kelvin waves, which move toward the Equator and then propagate eastward along the Equator. After reaching the eastern boundary, these waves reflect and bifurcate into the poleward propagating waves. Although the waves moving along the eastern boundary with a speed close to the Kelvin waves, recent studies suggested that such waves should be interpreted in terms of long cyclonic Rossby waves, e.g., Marshall and Johnson (2013). On their poleward propagation path along the eastern boundary, these waves gradually shed their energy and mass, forming the westward baroclinic Rossby waves carrying the signals through the ocean interior. Since the dissipation along the eastern boundary is relatively low, the thermocline thickness perturbation is nearly constant along the entire eastern boundary, as shown in the right edge in Fig. 7a.
Because the total amount of warm water in the upper layer is conserved, thermocline perturbations along the eastern boundary must be negative in order to compensate the increase of warm water in the ocean interior at low latitudes; in the present case, layer thickness anomaly along the eastern boundary is about −4.8 m, Fig. 7a. Due to the negative wind stress anomaly applied to the Equator and low latitude band, thermocline slope anomaly in this region is positive. The sharp reduction of layer thickness perturbations near the western boundary indicates that the northward transport of the western boundary current at this latitude is enhanced. On the other hand, thermocline anomaly is negative at high latitudes, Fig. 7a.

Cases with wind stress perturbations at one latitudinal band
We continue to present the results from the first set of experiments, in which wind stress perturbations apply to a single latitudinal band, as shown in Fig. 5a-e. Exp. 2H-A has been discussed above; in Exp. 2H-B, a positive wind stress anomaly applied to the Equator. Since wind stress perturbations applied to these experiments have small amplitude, results from Exp. 2H-B are very close to those in Exp. 2H-A, but with opposite signs. For example, volume anomaly is now negative for the equatorial band, and it is positive at middle and high latitudes, black curve in Fig. 8a. In Exp. 2H-C, a positive wind stress anomaly applied to 20°N, reducing the Ekman pumping in the subtropical basin, and thus creating a negative volume anomaly at the latitudinal band from 15°N to 45°N; however, at other latitudes volume anomaly is positive, blue curve in Fig. 8a.
In Exp. 2H-D, a positive wind stress anomaly applied to 40°N band, creating a negative volume anomaly at a latitudinal band from 35°N to 60°N, red curve in Fig. 8a. Note that the volume anomaly created in this case is much larger than in the previous two cases. This will also lead to other stronger anomalous features. In compensation, at other latitudes the volume anomaly is positive.
In Exp. 2H-E, a positive wind stress anomaly applied to 60°N band, creating a negative volume anomaly north of 57°N. For latitudinal band from 30°N to 56°N the volume anomaly is positive; however, below 30°N, the volume anomaly is negative again, magenta curve in Fig. 8a.
The corresponding southward transport of warm water in the upper layer creates a southward MOC in Exp. 2H-C and 2H-D, with the maximum amplitude near the latitude of wind stress anomaly maximum, i.e., 20°N and 40°N respectively, Fig. 8b. In Exp. 2H-D, the mean MOC maximum reaches a large value of −0.64 Sv. The MOC associated with the thermohaline circulation in such a two-hemisphere model basin is likely to be on the order of 10 Sv; thus, the amplitude of perturbations in Exp. 2H-D is about a few percentages of the climatological mean MOC. In Exp. 2H-E, positive wind stress anomaly applied to 60°N, creating a relatively weak negative MOC, with its maximum value around 56°N; however, south of this negative MOC there is a large positive MOC occupying the rest of the basin, with a maximum value of 0.23 Sv. The anomalous MOC gives rise to a sizeable poleward heat flux in the basin, Fig. 8b. In particular, in Exp. 2H-D, the maximum southward heat flux is about 27.5 TW, which is a sizeable fraction of the heat flux variability for such a model basin.
Another critical important aspect of the wind-driven circulation adjustment is the redistribution of warm water in the vertical direction. The heat content anomaly in the vertical direction in these cases appears in the form of first baroclinic modes, as shown in Fig. 8c. In Exp. 2H-B, heat content anomaly is positive above the depth of 420 m, but and it is negative below this depth. In Exp. 2H-C, the zerocrossing of the first baroclinic mode is moved downward to the depth of 470 m.
The rate of warming in the upper ocean has been discussed in many recent publications, e.g., Lyman and Johnson (2008), Lyman et al. (2010), Abraham et al. (2013), Chen and Tung (2014). Due to the relatively spare data coverage, the rate of warming remains uncertain and for the following discussion we will use the mean rate reported in the comprehensive review by Abraham et al. (2013). According to their analysis, over the period of 1970-2012, the planetary heat storage in the upper 700 m is 0.27 ± 0.04 W/m 2 , equivalent to an increase of heat content in the global upper ocean of 1.9 × 10 23 J (or 27 × 10 19 J/m) or temperature change of 0.2 °C. This is a mean barotropic mode of warming, which can be used as a benchmark. It is readily seen that the magnitude of the baroclinic modes of heat content anomaly obtained from this set of experiments are about 1/3 of the mean global warming rate inferred from observations.
In Exp. 2H-D warm water at latitudes higher than 35°N is pushed toward lower latitudes, red curve in Fig. 8a. As a result, thermocline at lower latitudes deepens, so that heat content anomaly is positive below 390 m, but it is negative in the upper ocean, red curve in Fig. 8b. In Exp. 2H-E, the situation is opposite, i.e., heat content anomaly is negative below 380 m, but it is positive above this depth. Note that in terms of the vertical heat content anomaly, over most part of the depth the perturbations created in Exp. 2H-D are of opposite signs, compared with other three cases in this set of experiments. It is clear that the shape of heat content anomaly created by wind stress perturbations is the result of the delicate competition between the zonal mean thermocline depth and wind stress perturbations applied to the model ocean.
The strong heat content anomaly in the vertical direction corresponds to an equivalent vertical heat flux. As shown in Fig. 8d, the vertical heat flux averaged over 40 year in Exp. 2H-B, 2H-C, and 2H-E are all positive, indicating an upward heat content shifting; the vertical heat flux is on the order of 2-8 TW. On the other hand, the vertical heat flux in Exp. 2H-D is −11 TW (−0.04 W/m 2 ), and it is about 1/7 of mean warming rate of 0.27 ± 0.04 W/m 2 estimated by Abraham et al. (2013).
The slope of the thermocline is changed. When the wind stress perturbations apply at 20°N (Exp. 2H-C), the thermocline depth along the eastern boundary is increased in compensation with the shoaling of the thermocline at middle latitudes, Fig. 9a; the slope of the thermocline is reduced at the Equator (black curve in Fig. 9d), 20°N and 40°N (blue and red curves in Fig. 9a). In particular, the slope is greatly However, the slope of the thermocline at 60°N is slightly enhanced in the eastern basin, Fig. 9a. The meridional mean of layer depth anomaly is positive in the eastern basin, but it is negative at the western basin. This indicates that warm water is pushed from the western basin to the eastern basin. In Exp. 2H-D (wind perturbations apply at 40°N) the thermocline shoals greatly at 40°N. In compensation, thermocline depth along the eastern boundary is increased, and the thermocline depth at the Equator and 20°N is slightly increased, Fig. 9b, e. The thermocline at 60°N is slightly reduced at 60°N, Fig. 9b.
In Exp. 2H-E (wind perturbations apply at 60°N) the thermocline at 40°N deepens. In compensation, thermocline depth along the eastern boundary is reduced, and the thermocline depth at other latitudes is slightly reduced, except for a small region near the eastern boundary at 60°N, Fig. 9c.
Thermocline perturbations in the Southern Hemisphere are quite different from those in the Northern Hemisphere, lower panels in Fig. 9. First, perturbations in this hemisphere are much weaker because this region is remote from the forcing. The thermocline depth along the eastern boundary is globally constant, and this sets up the thermocline perturbations in the Southern Hemisphere. Thus, when the wind perturbations apply to 20°N and 40°N, thermocline depth perturbations in the Southern Hemisphere are positive; on the other hand, when wind stress perturbations apply to 60°N, thermocline perturbations in the Southern Hemisphere is negative and small.

Cases with symmetric wind stress perturbations
In this set of experiments, wind stress perturbations symmetric to the Equator apply to two latitudinal bands, labeled as Exp. 2H-F, 2H-G, and 2H-H shown in   Fig. 10a. When wind stress perturbations apply to 60°N/60°S, warm water volume anomaly is negative around 65°N/65°S latitude band, and it is positive around 45°N/45°S and negative between 35°N and 35°S, red curve in Fig. 10a. The meridional transport of warm water induces transient MOC and poleward heat flux, anti-symmetric to the Equator, Fig. 10b. The maximum anomalous MOC (0.48 Sv) and poleward heat flux (20 TW) appear in the case when wind stress perturbations apply to the 40°N/40°S latitude bands.
The combination of wind stress perturbations at two latitude bands induces a baroclinic mode of heat content anomaly in the vertical direction, which amplitude is larger than the cases with single wind perturbation band discussed above, Fig. 10c. These heat content anomalies correspond to relatively large vertical heat flux, Fig. 10d.
As discussed above, wind stress perturbations also induce changes of the thermocline depth, Fig. 11. Since wind stress anomaly is symmetric to the Equator, response in both hemispheres is the same; hence only changes in the Northern hemisphere are discussed here. In all these three experiments, changes of the thermocline at latitude 20°N and 40°N has feature similar to those produced in Exp. 2H-C, 2H-D and 2H-E. On the other hand, changes along the equatorial band and 60°N are much larger than in Exp. 2H-C, 2H-D and 2H-E. As a result, the basin mean zonal thermocline depth change is about double the size produced in Exp. 2H-C, 2H-D, and 2H-E.

Cases with asymmetric wind stress perturbations
In this set of experiments (Exp. 2H-I, 2H-J, and 2H-K in Fig. 5), wind stress perturbations asymmetric to the   Fig. 12b. These can be a sizeable contribution to the climate variability on decadal time scales. Furthermore, the asymmetric wind perturbations can induce heat content anomaly in the form of second baroclinic modes, Fig. 12c. The heat content anomaly corresponds to a vertical heat flux; but under the asymmetric wind perturbations, the heat content anomaly and vertical heat flux are much weaker than the cases with symmetric wind perturbations, Fig. 12d.
The reason of higher baroclinic modes is as follows. Due to the asymmetric wind perturbations, thermocline in one hemisphere deepens, but shoals in the other hemisphere. For example, in Exp. 2H-I, zonal winds stress is reduced along 20°S, leading to stronger Ekman pumping for the subtropical gyre in the Southern Hemisphere. With warm water transported from the Northern Hemisphere, thermocline in the whole Southern Hemisphere deepens, with the maximum gain at the lower part of the water column, the black curve in Fig. 13a.
On the other hand, around 20°N zonal wind is enhanced, leading to weaker Ekman pumping for the subtropical gyre in the Northern Hemisphere. As a result, warm water is removed from the subtropical gyre in the Northern Hemisphere, and there is negative heat content anomaly for the Northern Hemisphere, the blue curve in Fig. 13a. Due to the slight nonlinearity associated with the inverse of reduced gravity, the gain and loss of heat content profiles in two hemispheres are not exactly asymmetric, and thus resulting in a relatively small residual heat content profile, the red curve in Fig. 13a. When the wind stress perturbations apply to higher latitudes, the heat content profiles in each hemisphere gradually change their shape. For example, in Exp. 2H-J the heat content profile in each hemisphere appears in the form of a baroclinic mode, but in Exp. 2H-I and 2H-K heat content profile in each hemisphere has the same sign, i.e., there is no zero-crossing of the heat content profile in each hemisphere. The combination of heat content anomaly in two hemispheres leads to different shape of net heat content anomaly for the whole basin, as shown in Fig. 13.
The layer depth anomaly at different latitudes is shown in Fig. 14. Because wind stress perturbations are asymmetric to the Equator, layer thickness anomaly is asymmetric too, and the basin mean layer thickness anomaly along the Equator is zero; hence only layer thickness anomaly in the Northern Hemisphere is shown in Fig. 14. In Exp. 2H-I, layer thickness along 20°N and 40°N is greatly reduced because the weakened Ekman pumping due to the positive zonal wind stress perturbations at low latitudes. The corresponding layer thickness anomaly along 20°S and 40°S (not shown in this figure) should have large positive values. At high latitudes (60°N), there is virtually no change in layer thickness.
When zonal wind stress perturbations apply to 40°S/40°N (Exp. 2H-J), layer thickness anomaly along 40°N becomes more negative, while layer thickness anomaly along 20°N becomes slightly positive, depicted by the blue curve in Fig. 14b. There is now a small negative value of layer thickness anomaly along 60°N and near the eastern boundary. In Exp. 2H-K, layer thickness anomaly is mostly confined to the latitude band around 40°N (40°S) due to the decline of polar easterly and the weakening of the Ekman upwelling.

A Southern Hemisphere model ocean
Wind-driven circulation in the Southern Hemisphere has very special features because all sub-basins in the Southern Ocean are linked through the ACC. The existence of this periodic channel gives rise to some unique features of the heaving modes. As sketched in Fig. 15, in addition to the inter-gyre modes for a single basin discussed above, there are two new types of basic modes unique to the Southern Oceans, including the annular modes and the inter-basin modes.
First, we assume that zonal wind stress in the South Oceans changes along certain latitudinal bands. For example, if zonal wind stress over the ACC is intensified, the slope of the front in the ACC increases and the thermocline shoals. As a result, the mean depth of thermocline in the ACC declines, depicted as the change from the solid line to the dashed line, Fig. 15b. Warm water in the ACC band is pushed toward low latitudes, leading to slightly deeper thermocline in three sub-basins, depicted by the dashed lines.
A concrete example is as follows. We assume the zonal wind near the southern boundary of a model ocean is intensified, red curve in Fig. 16a. The front in the ACC moves northward and becomes steeper. The change of thermocline shape in the ACC pushes warm water in the upper ocean from the ACC toward lower latitudes. As a result, warm water volume in the ACC declines, but it increases in middle and lower latitudes, red curve in Fig. 16b.
As discussed above, the meridional movement of warm water in the upper layer implies that there must be  a compensating return flow below the thermocline, and thus an equivalent MOC in the reduced gravity model. Note that the MOC is zero at the beginning and end of the transition state associated with wind stress perturbations, and it is non-zero only during the transit state of the reduced gravity model. Therefore, we will use the rate averaged over the whole transition period to evaluate the equivalent MOC during the adjustment, the red curve in Fig. 16c. As discussed above, the wind stress change induced MOC also carries an equivalent poleward heat flux, which also contributes to the global climate change.
Another important consequence of wind stress induced adjustment is the change of the vertical stratification or the heat content profile. Due to the intensification of westerly over the ACC, warm water is transported from the ACC to lower latitudes. Since the thermocline is deep at lower latitudes and shallow at high latitudes, such movement of warm water in the upper ocean induces a negative volume anomaly at the shallow depth and a positive volume anomaly at the deep level, the red curve in Fig. 16d. The shifting of the heat content in the vertical direction implies a vertical heat flux, and the corresponding profile averaged over the adjustment is shown in Fig. 16e. If the zonal wind stress is weakened, an opposite process should take place, as depicted by blue curves in Fig. 16. Since the amplitude of wind stress perturbations are relatively small, the anomaly produced has nearly the same patterns, but with opposite signs. The details of these solutions will be explained in the following sections.
Second, we assume that zonal wind stress in individual sub-basins changes along certain latitudinal band. For example, if zonal wind stress in the Pacific basin is intensified, the thermocline there deepens, moving from the solid line to the dashed line in the upper middle part of Fig. 15c. Assuming the total volume of warm water in the upper ocean remains unchanged, the thermocline in both the Indian and Atlantic basins moves upward in compensation, the dashed lines in the upper left and right parts of Fig. 15c. The thermocline in the ACC band may move slightly, but such change is excluded in this sketch.

Model set up
Using the same basic equations in Sect. 2, a second model is formulated for an idealized Southern Hemisphere ocean (SH model hereafter). The model basin includes a 360° wide channel model subjected to a periodic condition, and it extends from 60°S to the Equator. The northern part of the model ocean is divided into three sub-basins, separated by three continents. The Indian and Atlantic basins are 60° wide, and the Pacific basin is 150° wide; each continent is 30° wide in longitude, Fig. 17. The southern part of the model ocean is occupied by a 15° wide periodic channel, corresponding to the ACC.
The model is also an equatorial beta-plane model, with β = 2.1 × 10 −11 /m/s this gives rise to a Coriolis parameter matches that in spherical coordinates at 45°S, the northern edge of the periodic channel. The other parameters of the model are set as follows: A m = 1.5 × 10 4 m 2 /s, κ = 0.005 m/s. For numerical stability, a rather high interfacial friction parameter is used for the SH model. The same time step of 3153.6 s was used in numerical experiments.
To capture the deep thermocline in the South Oceans a relatively large amount of warm water is specified in the initial state which corresponds to a mean layer depth of 750 m. The reduced gravity model used here makes use of observations. According to WOA09 data (Antonov et al. 2010), the annual mean potential density, averaged over depth of 0-700 m for the global oceans, is estimated at σ 0 = 26.61 kg/m 3 , but the corresponding value over depth of 800-5500 m is estimated at σ 0 = 27.72 kg/m 3 ; thus, g ′ ≈ 1.11 cm/s 2 . Hence a round off value of g ′ = 1.0 cm/ s 2 will be used in this model. The annual mean potential temperature averaged over depth 0-700 m is 9.55 °C, and the corresponding value over depth 400-5500 m is 1.87 °C; thus the temperature difference between the upper and lower layers is 7.68 °C. Hence, in this model the temperature difference is set at 7.5 °C. The zonal wind stress profile (Fig. 17a) is taken from the zonally mean zonal wind stress averaged over the 51 years of SODA 2.1.6 data (Carton and Giese 2008).
The model was run for 300 years to reach a reference state. As will be shown shortly, due to the selection of parameters, the model ocean can reach a quasi-equilibrium state within 100 years. The thermocline depth and the streamfunction of the reference state are shown in Fig. 17. This reference state has three subtropical gyres north of the periodic channel. The thermocline thickness is nearly constant along the all eastern boundaries of the model ocean. Along the southern edge of the model ocean thermocline shoals to the depth of less than 200 m. The thermocline is the deepest in the Pacific basin, reaching 956 m, red curve in Fig. 17b. The subtropical gyres in three sub-basins are quite weak, with maximum transport of 11 Sv (Indian basin), 16 Sv (Pacific basin) and 11 Sv (Atlantic basin); the transport of the modelled ACC is 26 Sv, Fig. 17d. It is clear that these values are much smaller than the corresponding values in the world oceans. For example, diagnosis based on climatological hydrographic data indicates that the thermocline depth is the deepest in the South Indian Ocean; on the other hand, the transport of ACC is on the order of 100 Sv. Such large differences between our model ocean and the world oceans are due to the fact that the model is formulated for a rather idealized geometry with uniform reduced gravity and forced by a simple zonal wind stress which is zonally constant over the whole Southern Oceans. Since the goal of this study is to explore the fundamental structure of the heaving modes, the difference in thermocline depth between the model and observations should not qualitatively affect the basic results of our study.
From this reference state, we carried out a series of numerical experiments. In the first set of experiments, the zonal wind stress was perturbed along certain latitudinal bands; such wind stress perturbations were chosen to explore the annular modes sketched in Fig. 15b. In the second set of experiments, the zonal wind stress perturbations were confined to individual basins; such choice of wind stress perturbations is aimed to explore the inter-basin modes depicted in Fig. 15c.

Exp. SH-A, SH-B and SH-C
In this set of experiments the model was restarted from the reference state and forced by additional small positive perturbations in zonal wind stress, as defined in Eq. (9). The wind stress perturbations were linearly increased from 0 at t = 0 to the specified strength at the end of 20 years. Afterward, the wind stress perturbations were kept constant and the model run for additional 100 years.
The reason of running the experiments for 120 year is as follows. Adjustment of wind-driven circulation is carried out by wave like perturbations, mostly the baroclinic Rossby waves, which move quite slowly at high latitudes, typically on the order of a few centimeters per second. Due to the large interfacial friction imposed in the model, the adjustment time scale of wind-driven circulation in our model ocean is primarily determined by the high latitude basin crossing time of first baroclinic Rossby waves. The theoretical speed for the baroclinic Rossby waves is The southern tips of the continents in the model ocean is 45°S, where the Coriolis parameter and beta (15) c = βg ′ H/f 2 in our beta-plane model is f = 1.0395 × 10 −4 /s, β = 2.1 × 10 −11 /m/s. Since g ′ = 0.01 m/s 2 , assuming H = 800 m gives a typical wave speed of 0.0168 m/s at high latitudes. For a beta-plane model ocean of 360 grids with a zonal grid size of 110 km, the corresponding time for the first baroclinic Rossby waves to travel through the southern tip of the South Pacific basin in the model is about 31 years, and the corresponding time for the waves to travel around the whole Southern Oceans in the model is around and 75 years. Thus, running the model for additional 100 years after the wind stress perturbations reaching its final amplitude is long enough for the circulation to reach a quasi-equilibrium state.
We begin with Exp. SH-A, the case of a weakened equatorial easterly. As shown in Fig. 18a, due to weakening of the equatorial easterly, the zonal mean equatorial thermocline moves upward and the warm water above the equatorial thermocline is pushed poleward. Thus, warm water volume in the equatorial band declines with time; since the total amount of warm water in the model ocean is conserved, warm water volume at middle/high latitudes increases.
The adjustment of wind-driven gyre in the upper layer implies an anomalous MOC. When wind stress perturbations apply to the equatorial band, the adjustment of winddriven circulation induces a negative MOC because warm water in the upper ocean is pushed southward. This anomalous MOC grows with time, and its peak value at year 20 is around 0.4 Sv; however, as soon as the wind perturbation is no longer increased, the MOC quickly diminishes, Fig. 18b.
The adjustment of wind-driven circulation also induces warm water redistribution in the vertical direction. In the reference state, thermocline near the low latitude western boundary is the deepest, reaching 956 m in the Pacific basin. Hence, negative thermocline perturbation at low latitude bands induces a negative heat content anomaly at this depth range, but heat content anomaly in the shallow layers is positive, Fig. 18c.
When wind stress perturbations apply to the 30°S latitudinal band, warm water volume anomaly south of 30°S is negative; warm water from this latitudinal band is pushed to lower latitudes, Fig. 18d. The movement of warm water induces a northward MOC, which reaches the peak value of 1 Sv at year 20. However, as soon the wind stress perturbations are level off, the MOC diminishes, and quickly drops to a quite low level, Fig. 18e.
In fact, for the last 20 years, the MOC becomes negative, as indicated by the white area on the right hand side of Fig. 18e. Since the amplitude of negative MOC is so small, it cannot be shown in the contour Fig. 18e; instead, we include a fine figure for this area of negative MOC in Fig. 19a. This figure indicates that for the last 30 year of the numerical experiment the MOC is reversed to a rather small negative value. This indicates that the solution enters into an oscillating mode. Oscillation modes in the world oceans have been discussed in many previous studies. For example, Cessi and Primeau (2001) explored the low-frequency modes in closed basins. Recent studies indicated that for the case with low friction, the oscillations associated with these eigen-modes may take a long time, on the order to multi decades or even centennial, to decay, e.g., Allison et al. (2011), Jones et al. (2011) and Samelson (2011. However, due to the strong dissipation imposed by the interfacial friction and lateral friction in our model, these oscillations are strongly damped.
In Exp. SH-B, the redistribution of warm water in the ocean leads to a downward shifting of warm water. As a result, warm water volume below 850 m increased, but it is reduced above this depth, Fig. 18f.
In Exp. SH-C, the positive zonal wind stress perturbations apply to the 60°S latitude band, the wind-induced Ekman pumping rate in the subtropical basins is greatly enhanced. Consequently, this induces large positive warm water volume anomaly north of 50°S, and volume anomaly  Fig. 18g.
The movement of warm water induces a northward anomalous MOC, which reaches the peak value of 1.2 Sv at 50°S in year 20. However, as soon the wind stress perturbations are leveled off, the MOC diminishes, and quickly drops to a quite low level, Fig. 18h.
Similar to Exp. SH-B, for the last 20 years in Exp. SH-C the MOC becomes negative, as indicated by the white area on the right hand side of Fig. 18h; a refine figure for this area of negative MOC is shown in Fig. 19b.
The redistribution of warm water in the ocean leads to a downward shifting of warm water. As a result, warm water below 850 m increased, but it is reduced above this depth. In particular, warm water volume decline is maximum near the 200 m level, corresponding to the strong shoaling of the front in the modeled ACC, Fig. 18g, i.
Since the solutions vary greatly with time, it is also meaningful to examine the final states of these three experiments at the end of the 120 year experiments. The meridional distribution of the volume anomaly is shown in Fig. 20a. In Exp. SH-A, the volume anomaly is negative north of 30°S, but it is positive south of 30°S. The meridional shifting of warm water implies a MOC. Although the MOC averaged over the 120 experiments is much smaller than its peak value at year 20, it is still quite sizeable. For Exp. SH-A, the MOC minimum (averaged over 120 year run) is −0.059 Sv, Fig. 20b. Due to the anomalous MOC, there is a poleward heat flux, with a mean value of nearly −1.8 TW averaged over the 120 years of the experiment. In comparison, the volume anomaly in Exp. SH-B is much larger and it is of opposite sign, the blue curve in Fig. 20a. The MOC averaged over 120 year reaches its peak of 0.19 Sv around 30°S. This MOC also carries a sizeable poleward heat flux with its peak value of 5.7 TW, blue curve in Fig. 20b.
In Exp. SH-C, there is a negative volume anomaly in the periodic channel, giving rise to a MOC (maximum values of 0.23 Sv) and poleward heat flux (maximum value of 7.4 TW), red curve in Fig. 20b.
The anomalous MOC and poleward heat flux diagnosed from these experiments are much smaller than the corresponding values inferred from long term mean observations. For example, Talley (2013) put the estimate of the global overturning cell associated with bottom water formation in the world oceans at 29 Sv, and the associated poleward heat flux on the order of 0.1-0.2 PW (1 PW = 10 15 W). Nevertheless, these transient MOC and poleward heat flux suggest that adiabatic adjustment of wind-driven circulation may contribute to a substantial portion of variability in the MOC and poleward heat flux diagnosed from observations or numerical simulations of the world oceans.
The adjustment of wind-drive circulation also induces the vertical redistribution of heat content. At the end of 120 year experiment, the heat content anomaly all appears in the form of first baroclinic modes, Fig. 20c. In Exp. SH-A, the basin mean heat content anomaly is positive above 850 m; it is negative below 860 m, and reaches the minimum of −5.1 × 10 19 J/m at depth of 910 m. In Exp. SH-B, the basin mean heat content anomaly has signs opposite to Exp. SH-A: it is negative above 800 m; it is negative below 820 m, and reaches the maximum of 13.5 × 10 19 J/m at depth of 870 m. In Exp. SH-C, the basin mean heat content anomaly pattern is similar to Exp. SH-B; however, the upper branch looks quite different: this negative branch reaches its minimum of −9.3 × 10 19 J/m at the depth of 220 m; the heat content profile has a zero crossing at the depth of 610 m, and it reaches its maximum of 17.1 × 10 19 J/m at depth of 880 m.
The basin-scale redistribution of heat content in the vertical direction implies a vertical heat flux, defined in Eq. (14). For Exp. SH-A, the vertical heat flux peak is 1.1 TW (0.005 W/m 2 ); for Exp. SH-B, it is −2.6 TW (−0.012 W/m 2 ); for SH-C, it is −4.9 TW (−0.024 W/m 2 ), which is much stronger than previous two cases, Fig. 20d.
The baroclinic mode structure shown in Fig. 20c is for the profile averaged over the whole model ocean; but the heat content profile in each sub-basin can be quite different, depending on the spatial distribution of wind perturbations.
Due to assumption of adiabatic motions, the basin mean heat content anomaly must appear in the form of baroclinic modes. For each sub-basin, however, the net heat content anomaly is non-zero due to inter-basin shifting of warm water. Thus, heat content anomaly in a sub-basin can contain a barotropic component. In Exp. SH-A, heat content anomaly in both the Atlantic and Indian basins (blue and green curves in Fig. 21a) has a single positive lobe below 750 m; heat content anomaly in the Pacific basin (black curve in Fig. 21a) is positive between 750 m and 830 m, Dashed black lines denote the mean warming rate in the upper ocean inferred from observations but it is negative below 840 m. The warm water from the deep part of the Pacific basin is pushed toward the ACC and piled up at shallower levels, red curve in Fig. 21a. Heat content anomaly in the ACC is zero above 160 m, and it is positive below, except for the depth below 910 m, where it has a rather small negative lobe (not visible in this figure). Note that above 750 m, heat content anomaly in all sub-basins is zero; hence the global heat content anomaly profile (magenta curve) overlaps with heat content anomaly profile in the ACC. Below this depth, the global heat content profile is dominated by the contribution from the Pacific basin (black curve). The heat content profiles in each sub-basin have quite different baroclinic structures in the present case. Thus, our results demonstrate that heat content anomaly profile induced by adiabatic motions of the wind-driven circulation in the Southern Oceans can have complex baroclinic structure.
In Exp. SH-B, the heat content anomaly is quite different from Exp. SH-A. The maximum depth (>950 m) of the thermocline in the reference state is located near 30°S, where wind stress perturbations can induce large positive heat content anomaly in all three sub-basin around the depth of 900 m. On the other hand, heat content anomaly is mostly negative above 760 m, Fig. 21b. Thus, the baroclinic structure of global heat content is in the form of a first baroclinic mode. If wind stress perturbations apply to the latitude band around 20°S, the global heat content anomaly appears in the form of a third baroclinic mode (figure not shown).
In Exp. SH-C wind stress perturbations apply to the 60°S latitude band where thermocline is much shallower. The transport of warm water from high latitudes to lower latitudes creates a positive heat content anomaly below 800 m. Although the basin mean heat content anomaly is in the form of a first baroclinic mode, heat content anomaly profiles in all three sub-basins have no zero-crossing, and they seems to be a combination of a barotropic mode and a second baroclinic mode, but heat content anomaly in the ACC appears in a form close to a first baroclinic mode, Fig. 21c. The transported warm water is mostly originated from the ACC at much shallower levels. Hence, the heat content anomaly above 630 m in the ACC is negative, with a peak near the 210 m level, indicating that the southern edge of the ACC loses a lot of warm water at the shallow level (red curve in Fig. 21c). The combination of heat content anomaly from these four sub-basins creates a first baroclinic mode with shape peaks at 200 m and 890 m level, magenta curve in Fig. 21c.
It is interesting to compare the amplitude of the baroclinic modes inferred from our simple model with observations. According to Abraham et al. (2013), over the period of 1970-2012 the planetary heat storage in the upper 700 m is estimated at 27 × 10 19 J/m. Accordingly, the baroclinic modes of heat content inferred from these experiments is smaller, but may be comparable with the mean warming rate inferred from observations. The corresponding time evolution of volume anomaly is each sub-basin is shown in Fig. 22. It is well known that the adjustment of wind-driven circulation in a closed basin is carried out through wave motions, in particular the first baroclinic mode of Rossby waves. As discussed above, the time scale of adjustment for the model ocean is estimated at 75 years, and the volume anomaly ratio (instantaneous volume anomaly/final volume anomaly) gradually reaches the final value of 1 unit after 80-100 year integration, low panels of Fig. 22.
Although the total amount of warm water at low/middle latitudes declines, the situation in each sub-basin is different. In Exp. SH-A, the Pacific basin loses warm water quickly, but the Atlantic basin actually gains warm water slowly. The Indian basin also loses warm water during the first 20 years, but it starts to gain warm water afterward and ending up with more warm water at the end of the 120 year run, Fig. 22a.
In this case, wind forcing anomaly along the Equator and low latitudes leads to a decline of the east-west slope of the thermocline at low latitudes. These signals reach the western boundary and form the coast-trapped Kelvin waves, which move toward the Equator. They propagate eastward along the Equator. After reaching the eastern boundary, these waves reflect and bifurcate into the poleward propagating waves. On their poleward propagation path along the eastern boundary, these waves gradually shed their energy and mass, forming the westward Rossby waves which carry the signal through the ocean interior, e.g., Huang et al. (2000).
In this experiment, the Pacific Basin sets the pace of the adjustment because of its large size. Only after the adjustment of Pacific basin is nearly complete, the final signals propagate downstream, i.e., eastward, leading to the completion of adjustment in the Atlantic basin, and then finally the Indian basin. In fact, the warm water volume anomaly in the Indian basin reaches 95 % of its final value (at year 120) only after 74.7 years. Note that warm water volume anomaly in the Atlantic basin actually overshoots the final value reached at year 120, lower part of column in Fig. 22a. The oscillations shown in model solutions discussed here are similar to the oscillatory solutions discussed in previous study (e.g., Cessi and Paparella 2001;Cessi et al. 2004). Due to the selection of parameters, however, these oscillations are strongly damped.
In Exp. SH-B, all sub-basins gain warm water in the final state, except the ACC which loses warm water. The corresponding time evolution of vertical volume anomaly is shown in Fig. 22b. The total volume of warm water in the Pacific basin actually overshoots and then turns back to the final value around year 80, lower part of Fig. 22b.
In Exp. SH-C, warm water in the upper ocean in the ACC band is pushed northward; thus, volume anomaly in the Indian, Pacific and Atlantic basin increases, but it is greatly reduced in the ACC, Fig. 22c. Since wind stress perturbations apply to rather high latitudes, the blockage of continents does not play much important role and thermocline adjustment in all basins are mostly synchronized and nearly completed within 60 years and without much delay in individual basins, as shown in the previous cases.

Exp. SH-D, SH-E, SH-F, SH-G and SH-H
In these experiments, the model was restarted from the reference state and forced by wind stress with small perturbations where �τ = 0.02 N/m 2 is the amplitude, x = 3300 km, y = 1100 km, (x 0 , y 0 ) are the longitude and latitude of the center of wind stress perturbations. The wind stress perturbations were linearly increased from 0 at t = 0 to the full scale of the specified strength at the end of 20 years. Afterward, the wind stress perturbations were kept constant and (16) τ x′ = �τ e −((x−x 0 )/�x) 2 −((y−y 0 )/�y) 2 the model run for additional 100 years. The aim of carrying out these experiments is to explore the inter-basin modes associated with warm water adjustment due to wind stress perturbations applied to individual basins.

Exp. SH-D and SH-E
In this set of experiments, positive zonal wind stress perturbations apply to the equatorial band. In Exp. SH-D wind stress perturbations applies to the Indian basin, (x 0 , y 0 ) = (60°E, 0°S); in Exp. SH-E wind stress perturbations applies to the Pacific basin, (x 0 , y 0 ) = (180°E, 0°S). Similar experiment was carried out for the case with zonal wind stress perturbations applied to the Atlantic basin; however, results from such an experiment are rather similar to Exp. SH-D; hence, the corresponding results are not included here.
These wind stress perturbations induce warm water volume decline at low latitudes (with the peak near 10°S) and warm water increase at high latitudes, Fig. 23a. The pattern of changes in the circulation is similar in both cases; however, when wind stress perturbations apply to the Pacific basin (Exp. SH-E), the changes in the circulation are much The meridional transport of warm water implies an anomalous MOC and poleward heat flux. When wind stress perturbations apply to the Pacific basin, the amplitude of the MOC (−0.017 Sv) and poleward heat flux (−0.52 TW) are much larger than the case when wind stress perturbations apply to the Indian basin (with the corresponding value of −0.0073 Sv and −0.23 TW), Fig. 23b. In addition, heat content anomaly in the vertical direction also changes in response. When wind perturbations apply to the Indian basin, the heat content anomaly appears in the form of a second baroclinic mode, black curve in Fig. 23c. However, when wind stress perturbations apply to the Pacific basin, the heat content anomaly appears in the form of a first baroclinic mode, but with much large amplitude, which is clearly due to fact that the Pacific basin is much larger than the Indian basin. As a result, the corresponding vertical heat flux in Exp. SH-E is also much larger than that in Exp. SH-D, Fig. 23d.
Changes in the vertical stratification in each sub-basin are different for these two cases. When wind stress perturbations apply to a sub-basin, thermocline there shoals, indicted by the large negative heat content anomaly in the corresponding sub-basin. In Exp. SH-D, heat content anomaly in the Indian basin is positive above 800 m; it is negative for depth of 800-940 m, Fig. 24a. On the other hand, in the Pacific and Atlantic basins it is positive over the depth of 750-950 m; in the ACC it is mostly positive over the depth of 170-900 m, with an extremely small negative lobe over the depth of 910-960 m (not visible in Fig. 24a). The contributions from all sub-basins give rise to a global heat content anomaly in the form of a second baroclinic mode, indicated by the magenta curve in Fig. 24a. On the other hand, when wind stress perturbations apply to the Pacific basin (Exp. SH-E), warm water volume anomaly is much stronger (Fig. 24b), basically doubled the amplitude in the previous case. Since wind stress perturbations apply to the Pacific basin, they induce a large negative volume anomaly in the Pacific basin over the depth of 850-960 m, but it has a small positive value over the depth of 750-840 m. In both the Indian and Atlantic basins, heat content anomaly is non-negative, and it has positive value over the depth of 750-950 m.
Heat content anomaly in the ACC has a positive lobe over the depth of 170-900 m. (There is a very weak negative heat content anomaly over the depth of 910-940 m, not visible in Fig. 24b). The contributions from these four sub-basins generate a global heat content anomaly in the form of a first baroclinic mode, magenta curve in Fig. 24b.
The time evolution of warm water volume anomaly for these two cases is shown in Fig. 25. When wind stress perturbations apply to the Indian basin (Exp. SH-D), the warm water in this basin is reduced (red curve in the upper panel of Fig. 25a) and pushed downstream to the other basins. With warm water directly from the upstream basin, the adjustment in the Pacific basin is the fastest, nearly completed within 30 years. The adjustment in the Indian basin is completed secondary, and that of the Atlantic basin is the last because it is located at the end in the downstream direction. Before the thermocline in the Atlantic basin is close to its final state, the circulation and thermocline in the Indian basin cannot be completed, and this is part of In Exp. SH-E, wind stress perturbations apply to the Pacific basin; thus, warm water in this basin is reduced (blue curve in the upper panel of Fig. 25b) and pushed downstream to other basins. With warm water directly from the upstream basin, the adjustment in the Atlantic basin is the fastest, it reaches the final value before year 30; after overshooting, it approaches the final state first. The adjustment in the Pacific basin and the ACC is completed secondary. The adjustment of Indian basin is the last because it is located at the end in the downstream direction.

Exp. SH-F, SH-G and SH-H
In this set of experiments, wind stress perturbations apply to the 30°S latitude band. In Exp. SH-F, the wind stress perturbations apply to the Indian basin, (x 0 , y 0 ) = (60°E, 30°S). Changes in the circulation are opposite to those in Exp. SH-D and Exp. SH-E. First of all, warm water is now transported from high latitudes (south of 30°S) to low latitudes, Fig. 26a.
Meridional transport of warm water in the reduced gravity model is equivalent to a MOC and poleward heat flux in the ocean. In this case, there is a northward MOC, Fig. 26b. Averaged over the 120 years of model run, the mean MOC rate is 0.024 Sv, and the mean poleward heat flux is 0.75 TW.
The meridional transport of warm water from middle/high latitudes creates baroclinic modes of heat content anomaly, Figs. 26c and 27. When wind perturbations apply to the Pacific basin, the global heat content anomaly appears in the form of a second baroclinic mode; when wind perturbations applied to either the Indian or Atlantic basins, the global heat content anomaly appears in the form of a first baroclinic mode, Fig. 26a, c. The amplitude of heat content anomaly and the corresponding vertical heat flux for the case with wind perturbation applied to the Pacific basin is much larger than for the cases with wind perturbation applied to the Indian or Atlantic basins, Fig. 26d.
When wind stress perturbations apply to the Indian basin, a negative heat content anomaly is created in this basin over the depth of 740-950 m (Fig. 27a), but its amplitude in the Indian basin is only half of that in Exp. SH-D when wind stress perturbations apply to the same longitude but along the Equator. On the other hand, the heat content anomaly in the Pacific basin is positive over the depth of 750-950 m, and its amplitude is double of that in Exp. SH-D. Heat content anomaly in the Atlantic basin is in the form of a first baroclinic mode, with a small negative value over the depth of 750-780 m, and a relatively larger positive value over the depth of 790-950 m.
Heat content anomaly in the ACC has a negative lobe at depth of 170-760 m; however, over depth of 770-890 m, it is positive; and for the depth of 910-950 m, it has an extremely small negative value (not visible for the scale used in this figure). The contributions from all sub-basins give rise to a global heat content profile in the form of a first baroclinic mode, magenta curve in Fig. 27a.
The time evolution of warm water volume in these subbasins is shown in Fig. 28a. Both the Indian basin and the ACC lose warm water, but the Pacific and Atlantic basins gain warm water. The adjustment in the ACC leads the process and completes the first. Indian basin is last one to complete.
In Exp. SH-G, wind stress perturbations apply to the Pacific basin. Warm water volume anomaly north of 30°S is positive. Because the size of the Pacific basin is more than In the ACC the major part of heat content signals is negative over depth of 170-750 m, red curve in Fig. 27b. Over depth of 760-900 m, it is positive. There is a very weak negative segment over depth of 910-940 m (its magnitude is too small to be visible in Fig. 27b). The contributions from these four basins combine into the global heat content anomaly in the form of a second baroclinic mode, magenta curve in Fig. 27b.
The time evolution of volume anomaly is shown in Fig. 28b. Similar to Exp. SH-F, the adjustment in the ACC leads that in the whole model ocean; it is nearly completed by year 35. The adjustment of the Atlantic basin is the second, and that in the Indian basin is the third. It is interesting to note that adjustment in the Pacific basin is the last one to finish. In fact, although the final volume anomaly in the Pacific is negative, it starts positive first. After the wind stress perturbations reach its full amplitude at year 20 and increase no more, the volume anomaly in the Pacific basin starts to increase. It enters the regime of positive value only after year 32.2. In comparison, volume anomaly ratio in the ACC crosses the line of 95 % level at year 30.8. It is clear that the adjustment in the Pacific basin is completed only after the completion of adjustment in all other basins. The dynamic detail of this adjustment remains unclear at this time, and it is left for further study.
In Exp. SH-H, wind stress perturbations apply to the Atlantic basin. The results from this experiment are rather similar to Exp. SH-F with wind stress perturbations applied to the Indian basin. As shown in Fig. 26, the red curves and the black curves are almost the same.
In the vertical direction, the basic features of the heat content signals are quite similar to those in Exp. SH-F. Since wind stress perturbations now apply to the Atlantic basin, heat content anomaly there is negative value below 740 m, green curve in Fig. 27c. Heat content anomaly in the Indian basin is positive and its amplitude is larger than that in the Pacific basin, because the Indian basin is directly downstream from the wind stress perturbations. Heat content anomaly in the ACC is quite similar to that in Exp. SH-F. Similar to that in Exp. F, contributions from individual basins give rise to a first baroclinic mode of global heat content anomaly, magenta curve in Fig. 27c.
The time evolution of volume anomaly in each basin is shown in Fig. 28c. The general feature is somewhat similar to Exp. SH-F. Of course, in this case, volume anomaly in the Atlantic basin and the ACC is negative. Furthermore, the amplitude of volume anomaly is now much smaller than in Exp. SH-F.

Discussion
Using a simple reduced gravity model, we carried out several sets of numerical experiments. These experiments demonstrated that adiabatic movement of warm water in the upper ocean induced by decadal wind stress perturbations can lead to many important dynamical consequences.
One of the most important phenomena is the redistribution of warm water in three-dimensional space. In particular, heaving induced by adjustment of wind-driven circulation can lead to vertical heat content anomaly in the form of the baroclinic modes. Since the structure of the modes depends on the two-dimensional shape of the thermocline and the wind stress perturbations, these modes may appear in different forms.
These baroclinic modes may be used to interpret the time evolution of heat content anomaly diagnosed from observations or computer generated climate datasets. In particular, since the Southern Oceans are connected through the ACC, changes of stratification observed in a local region may be caused by change of the local wind or change of wind stress in a remote area of the world oceans.
Note that heat content variability inferred from our model is at least one order of magnitude smaller than the mean warming rate inferred from observations. It is reasonable to expect that such a simple model may not be able to accurately produce the heat content variability; nevertheless, heat content variability induced by adiabatic motions in the ocean can have amplitude not negligible in diagnosing the general long term trend of climate warming or cooling.
Our model showed that if negative zonal wind stress perturbations are imposed along the Equator, the corresponding heat content anomaly is cooling in the upper ocean and warming in the deep layers. A combination of such a baroclinic mode of heat content anomaly with the general trend of warming over the whole water column may be a cause of the hiatus of global SST and the accelerating warming in the deep layers reported in many recent studies.
In addition, heat content anomaly profile reproduced in Exp. 2H-B indicates that the deep layers may be cooled down due to adiabatic motions. Although there is a general trend of global warming, the deep ocean cooling also occurs over certain time periods; such a phenomenon may seem odd, but our model results showed that it is dynamically quite conceivable.
For climate study, another important consequence of water mass redistribution is the transient MOC, poleward heat flux and vertical heat flux. Results obtained from our simplified model indicate that they can reach quite nonnegligible levels, and may consist of sizeable component of the transient MOC and poleward heat flux diagnosed from in situ observations and numerical simulation of climate change. In particular, our results indicated that there is an equivalent vertical heat flux induced by adiabatic motions of water masses. The contribution to vertical heat content anomaly and the associated vertical heat flux have not yet been thoroughly explored up till now; thus, our analysis may stimulate further study in this direction.
In this study we did not include the case with the upper layer outcropping; hence, heat content anomaly discussed here is limited to the subsurface layers only. It is, however, a straightforward step to carry experiments including the upper layer outcropping. In fact, the numerical model used in this study is based on the so-called positive-definite scheme, which can handle the case with outcropping and volume conservation of the moving layer, as discussed by Huang (1987).
In our simple reduced gravity model, the hierarchy of stratification and the wave motions are greatly simplified. For example, the barotropic Rossby waves are excluded, and this may distort the adjustment process, especially for short time scales. Furthermore, representing the continuous stratification in the ocean with one moving layer of constant density is certainly an over-simplification. As such, the MOC, poleward/vertical heat flux and heat content profiles inferred from this study should be viewed as the first step in quantifying the corresponding components in the ocean. Our initial testing based on a two-moving-layer model suggested that the MOC and poleward heat flux inferred from a reduced model tend to exaggerate these fluxes. To obtain more accurate pictures one should run models based on more realistic formulation by including more density layers. Whether larger amplitude oscillations in the MOC, poleward heat flux and vertical heat flux exist in more realistic model simulation and how to identify such anomalous field require further analysis of climate datasets generated from computer simulations.
In addition, it is well know that circulation in the ACC is closely linked to the eddy motions. As such, results from numerical simulations are rather sensitive to the resolution of the numerical models. There is a well-known phenomenon, the so-called eddy saturation. For example, Hallberg and Gnanadesikan (2006), Meredith and Hogg (2006) showed that in low-resolution models, the transport of the ACC increases with the amplification of the zonal wind stress; however, for the eddy-permitted models, the transport of the ACC reaches a plateau with the increase of wind stress, and would no longer be enhanced even with the further increase of zonal wind stress. Thus, results obtained from our simplified reduced gravity model are not expected to portrait the dynamics of the Southern Oceans accurately.
Nevertheless, we hope results obtained from our simple model can enlighten the fundamental physics related to the change of stratification in the world oceans induced by heaving. Further studies based on eddy-resolving models should provide much more accurate information related to such important issues.