## Abstract

The ocean’s circulation is analyzed in Absolute Salinity *S*_{A} and Conservative Temperature Θ coordinates. It is separated into 1) an advective component related to geographical displacements in the direction normal to *S*_{A} and Θ isosurfaces and 2) into a local component, related to local changes in *S*_{A}–Θ values, without a geographical displacement. In this decomposition, the sum of the advective and local components of the circulation is equivalent to the material derivative of *S*_{A} and Θ. The sum is directly related to sources and sinks of salt and heat. The advective component is represented by the advective thermohaline streamfunction . After removing a trend, the local component can be represented by the local thermohaline streamfunction . Here, can be diagnosed using a monthly averaged time series of *S*_{A} and Θ from an observational dataset. In addition, and are determined from a coupled climate model. The diathermohaline streamfunction is the sum of and and represents the nondivergent diathermohaline circulation in *S*_{A}–Θ coordinates. The diathermohaline trend, resulting from the trend in the local changes of *S*_{A} and Θ, quantifies the redistribution of the ocean’s volume in *S*_{A}–Θ coordinates over time. It is argued that the diathermohaline streamfunction provides a powerful tool for the analysis of and comparison among ocean models and observation-based gridded climatologies.

## 1. Introduction

The global ocean circulation is of great importance for the earth’s climate system through its heat transport and the cycling and storage of carbon. Understanding the underlying physical processes that drive the global ocean circulation is essential for our understanding of the earth’s climate system and our ability to model it accurately. However, observing and modeling the global ocean circulation remains a challenge.

A component of the ocean circulation is buoyantly or density-driven due to thermohaline forcing and is often referred to as thermohaline circulation. Here, thermohaline forcing refers to boundary freshwater and heat fluxes and salt and heat fluxes by diffusive mixing. In this paper, we develop a description of the thermohaline circulation using its unambiguous relationship to the thermohaline forcing. The use of thermodynamic coordinates for understanding the relationship between the circulation and thermodynamic forcing was first made clear by Walin (1982). He showed that the area-integrated surface heat flux between two outcropping isotherms could be used to calculate the diathermal advection in a steady-state ocean.

Speer and Tziperman (1992) applied Walin’s framework to isopycnals instead of isotherms, which together with many subsequent studies (Marshall 1997; Marshall et al. 1999; Nurser et al. 1999; Sloyan and Rintoul 2000, 2001; Iudicone et al. 2008; Badin and Williams 2010; Nikurashin and Ferrari 2013) allowed an understanding of ocean circulation driven by thermohaline forcing. Using this framework, Nurser and Marsh (1998) showed that the total diapycnal transport can be expressed as a sum of a streamfunction difference and nonsteady component. The diapycnal transport can be fully expressed as a streamfunction only when the nonsteady component is small.

Streamfunctions have been widely used as a diagnostic to study ocean circulation. The advantage of scalar streamfunctions lies in their ability to represent the complex three-dimensional and time-varying ocean circulation in two dimensions, allowing new insight into the ocean circulation. A classic example is perhaps the meridional overturning streamfunction Ψ_{λz}, representing the zonally averaged ocean circulation in latitude *λ* and depth *z* coordinates. Because ocean currents tend to follow isopycnal surfaces rather than surfaces of constant depth, Döös and Webb (1994) calculated a streamfunction in latitude and potential density *σ* coordinates Ψ_{λσ}. They identified circulation cells related to the Atlantic overturning, the subtropical gyres, and the Antarctic Bottom Water circulation and showed that the Deacon cell, which is a strong circulation cell visible in Ψ_{λz}, almost entirely disappeared in *λ*–*σ* coordinates. From observations, Marsh (2000) used Walin’s framework to derive a streamfunction in latitude and density *ρ* coordinates Ψ_{λρ}. A more thermodynamic streamfunction was presented by Nycander et al. (2007), who constructed a streamfunction in density and depth coordinates Ψ_{ρz}. This study showed that certain ocean circulation cells are either mechanically forced, buoyantly forced, or forced by a combination of both. This clearly illustrates that streamfunctions provide a different representation of the ocean circulation for each choice of coordinates.

Simultaneously, Zika et al. (2012) and Döös et al. (2012, hereafter referred to as ZD12) constructed the thermohaline streamfunction, a streamfunction in salinity *S* and temperature *T* coordinates Ψ_{ST}. The Ψ_{ST} is essentially a two-dimensional extension of Walin’s framework, applied in *S*–*T* coordinates. However, ZD12 calculated Ψ_{ST}, using a model’s three-dimensional velocity field instead of using the thermohaline forcing. Therefore, their Ψ_{ST} represented only advection normal to *S*–*T* isosurfaces and not necessarily net water mass transformation rates. ZD12 identified similar circulation cells and obtained estimates for circulation time scales and diapycnal freshwater and heat transports. A recent study by Zika et al. (2013) shows the influence of the magnitude of the Southern Ocean winds with circulation in different coordinates, including *S*–*T* coordinates.

Although ZD12 were the first to calculate Ψ_{ST}, Speer (1993) used salt and heat fluxes, in Walin’s framework, and projected the water mass transformation rates in *S*–*T* coordinates. This provided a distinction between transformation of different water masses and their locations. They emphasized that they were unable to distinguish if this transformation was associated with either local changes in density without actual motion of water or due to advective fluxes. Kjellsson et al. (2014) analyzed the atmospheric circulation in dry and moist isentropic coordinates to construct the hydrothermal streamfunction. They considered both the effect of Eulerian advection and local changes of properties and found that the latter was negligible for these coordinates.

Here we derive the complete ocean circulation in *S*–*T* coordinates (section 2), namely, the diathermohaline circulation. We show that the diathermohaline circulation is composed of an advective component and a component due to local changes in *S*–*T* values, without fluid motion in geographical space. The advective component is due to ocean circulation normal to isohaline and isothermal surfaces and can be expressed as an advective thermohaline streamfunction (section 3). The local component is due to local changes of the *S* and *T* values, expressed as geographical displacements of isohalines and isotherms. This results in changes in the ocean volume distribution in *S*–*T* coordinates and can be expressed as a local thermohaline streamfunction, after removing a trend (section 4). The sum of the advective and local thermohaline streamfunctions is the diathermohaline streamfunction and represents the nondivergent diathermohaline circulation in *S*–*T* coordinates and is the novel aspect of this paper (section 5). The analysis of a numerical model and observational climatology in *S*–*T* coordinates shows that the local thermohaline streamfunction is a significant component of the ocean circulation in *S*–*T* coordinates and cannot be ignored (section 6). We suggest that the diathermohaline streamfunction provides a tool for the analysis of, and comparison among, ocean models and observation-based gridded climatologies (sections 7 and 8).

## 2. Circulation in *S*_{A}–Θ coordinates

*S*

_{A}–

Here, we construct a formalism to calculate circulation in Conservative Temperature Θ and Absolute Salinity *S*_{A} coordinates. Conservative Temperature is proportional to potential enthalpy (by the constant heat capacity factor ), which represents the heat content per unit mass of seawater (McDougall 2003; Graham and McDougall 2013). Absolute Salinity is measured on the Reference Composition Salinity Scale (Millero et al. 2008) and represents the mass fraction of dissolved material in seawater (g kg^{−1}) (IOC et al. 2010; McDougall et al. 2012). The framework presented here could in principle be applied to any two coordinates defined by two conservative tracers.

### a. Thermohaline forcing

The *S*_{A}–Θ values of a fluid parcel define their position in *S*_{A}–Θ coordinates. For the fluid parcel to change its position in *S*_{A}–Θ coordinates, it must change its *S*_{A} and/or Θ values. The Conservative Temperature represents the ocean’s heat content. Processes by which the heat content of a fluid parcel can be changed are 1) heat fluxes at the ocean's boundary and 2) diffusive mixing (Graham and McDougall 2013). There are three processes by which one can change the *S*_{A} of a fluid parcel: 1) surface freshwater fluxes, 2) diffusive mixing, and 3) biogeochemical processes. The latter is assumed to be negligible (IOC et al. 2010). Hence, all changes of a fluid parcel’s *S*_{A} and Θ and motion in *S*_{A}–Θ coordinates are driven by boundary freshwater and heat fluxes, and salt and heat fluxes by diffusive mixing, that is, thermohaline forcing. Using continuity, it can be shown that for volume *V*, this results in the following conservation equation for *S*_{A}

and Θ

Here the velocity **u** = (*u*, *υ*, *w*), and and (*f*_{Θ} and *d*_{Θ}) are the net convergence of salt (heat) due to boundary fluxes and diffusive mixing processes, respectively.

Surface radiative and sensible heat fluxes cause movement in the Θ direction (Fig. 1). Precipitation results in movement in the negative *S*_{A} direction and evaporation results in movement in the positive *S*_{A} direction and simultaneously movement in the negative Θ direction due to the related latent heat loss (Fig. 1). Mixing acts to reduce both the stratification and the *S*_{A}–Θ range of the fluid parcels. Hence, the initial range enclosed by an isolated fluid parcel in *S*_{A}–Θ coordinates will never be increased by mixing processes. In the absence of external forcing, loss of salt, heat, or mass, the initial range will eventually be reduced to a point in *S*_{A}–Θ coordinates. Mixing depends on the magnitude and orientation of geographical gradients of *S*_{A} and Θ (Fig. 1). In geographical coordinates, mixing occurs as downgradient tracer diffusion due to epineutral mesoscale eddies or as small-scale vertical turbulent diffusion (Gent et al. 1995; Griffies 2004). However, the latter is in fact small-scale isotropic (rather than vertical) turbulent diffusion (McDougall et al. 2014).

### b. The diathermohaline velocity

A fluid parcel that moves in *S*_{A}–Θ coordinates must cross isohalines and/or isotherms, leading to a diahaline or diathermal velocity and related volume transport. For example, consider a volume *V*, with uniform Θ. If the Θ increases, this volume will be displaced, leading to a diathermal volume transport, in the positive Θ direction in *S*_{A}–Θ coordinates. Equivalently, if the *S*_{A} of *V* increases, this volume will be displaced, leading to a diahaline volume transport in the positive *S*_{A} direction in *S*_{A}–Θ coordinates. We reserve the use of the prefix “dia” for the net displacement through a surface (for a discussion on a diasurface velocity, refer to Griffies (2004, 138–141).

To calculate the diahaline or diathermal volume transport we first define as the diahaline velocity with which a fluid parcel moves through an isohaline in the positive *S*_{A} direction (in geographical coordinates, in the **∇***S*_{A} direction). Equivalently, we define as the diathermal velocity with which a fluid parcel moves through an isotherm in the positive Θ direction (in geographical coordinates, in the **∇**Θ direction). With this definition we can construct the diathermohaline velocity vector:

where and are perpendicular in *S*_{A}–Θ coordinates, but may have any angle relative to each other in Cartesian coordinates (Fig. 2). We define and relate the diathermohaline velocity with boundary fluxes and mixing by dividing (1) by |**∇***S*_{A}| and (2) by |**∇**Θ| to obtain the diahaline velocity

and the diathermal velocity

Here is the component of the fluid’s velocity with which a fluid parcel is geographically advected normal to the isohalines, and is the component of the fluid’s velocity, with which a fluid parcel is geographically advected normal to the isotherms. We define as the velocity of the displacement of an isohaline in geographical coordinates, normal to the isohalines, due to local changes of *S*_{A}. Finally, is the velocity of the displacement of an isotherm in geographical coordinates, normal to the isotherms, due to local changes of Θ. This derivation results in the diathermohaline velocity vector:

where and (Fig. 3). To analyze the diathermohaline circulation (6), it is convenient to express , using the diathermohaline streamfunction . However, only a velocity that is nondivergent in two coordinates can be represented by a streamfunction. To obtain the diathermohaline streamfunction we will first represent by the advective thermohaline streamfunction as previously obtained by ZD12. We will show that , after subtracting a trend, can be represented with the local thermohaline streamfunction . The sum of both the advective and local thermohaline streamfunction is the diathermohaline streamfunction.

## 3. The advective thermohaline streamfunction

Based on an averaging technique developed by Nurser and Lee (2004), Ferrari and Ferreira (2011) defined a meridional streamfunction for an arbitrary tracer *C*_{1}, where represents the advective meridional transport of *C*_{1} for an instantaneous velocity field *υ*. Zika et al. (2012) generalized this technique to compute a mean streamfunction for any two tracers *C*_{1} and *C*_{2}. This represents advection of fluid parcels in *C*_{1}–*C*_{2} coordinates. Applied to an instantaneous velocity field for a Boussinesq ocean such that **∇** · **u** = 0 and for tracers *S*_{A} and Θ, this results in

Here is the area integral over all Θ′, smaller than or equal to Θ on a surface of constant *S*_{A}. Equivalently, is the area integral over all , smaller than or equal to *S*_{A} on a surface of constant Θ, and is a time average. The advective thermohaline streamfunction is given by and represents the time-averaged thermohaline component of a nondivergent geographical ocean circulation, in *S*_{A}–Θ coordinates, as previously obtained by ZD12.

## 4. The local thermohaline streamfunction

We separate into a cycle and a trend term. We can then construct using the detrended component of .

For a conceptual description, imagine a control volume containing a motionless fluid with fixed volume *V*_{control} and uniform *S*_{A}–Θ values that (uniformly) change over time. Changing the *S*_{A}–Θ properties of the fluid in *V*_{control} over time leads to a volume transport in *S*_{A}–Θ coordinates. We now consider different changes of the *S*_{A}–Θ values of *V*_{control} to analyze and understand the resulting motion in *S*_{A}–Θ coordinates.

### a. The diathermohaline trend

First, consider heating the control volume from time *t*_{1} to *t*_{2}. When the fluid with volume *V*_{control} is heated, it may cross several isotherms and the related volume transport through these isotherms, in positive Θ direction, will be *V*_{control}/Δ*t*, where Δ*t* is a constant time step (Fig. 4a, cycle 1). Equivalently, a change in *S*_{A} causes the properties of *V*_{control} to cross isohalines, resulting in a similar volume transport in the positive *S*_{A} direction (Fig. 4a, cycle 2). Although the fluid in the control volume has not moved in geographical space, its *S*_{A}–Θ values have changed, resulting in motion in *S*_{A}–Θ coordinates. The resulting volume transport can be calculated by analyzing the total changes of the *S*_{A}–Θ values of the fluid, rather than analyzing the displacements of the isohalines and isotherms (see the appendix).

A net displacement of fluid in *S*_{A}–Θ coordinates can be understood as a diathermohaline trend (Fig. 4b, cycle 4). However, if this displacement is continued over multiple time steps, we can obtain both a diathermohaline trend and a time-averaged steady-state and nondivergent component, which we will refer to as the diathermohaline cycle.

### b. The diathermohaline cycle

The diathermohaline cycle in *S*_{A}–Θ coordinates takes place when changes of *S*_{A}–Θ values in the control volume returns to their initial *S*_{A}–Θ values over time, but the heating and cooling (freshening and salinification) takes place at a different *S*_{A} (Θ), such that the net transports through isohalines for a certain Θ range and isotherms for a certain *S*_{A} range do not exactly cancel out (Fig. 4b, cycle 5). This cycle results in a net volume transport of *V*_{control}/Δ*t*, across isohalines, for different Θ ranges, and net volume transports across isotherms, for different *S*_{A} ranges. Note that any cycle in which there are (simultaneous) changes of *S*_{A}–Θ values of the fluid in the control volume over time, which causes a transport across exactly the same isohalines and isotherms when moving away from and returning back to its initial state, results in no net transport over time (Fig. 4a, cycles 1, 2, and 3 from *t*_{1} to *t*_{3}).

For example, consider a volume in a motionless ocean, enclosed by a pair of isotherms and a pair of isohalines that change position with time (Fig. 5). For illustrative purposes we consider only the dynamics of three volumes that are within our grid indicated with *a*–*d*: 1) the volume *V*_{Θ}, of which Θ changes; 2) the volume , of which *S*_{A} changes; and 3) the volume , of which both *S*_{A} and Θ changes. We have the following sequence of events: (i) the (Θ + *d*Θ) isotherm changes position, warming the fluid in both *V*_{Θ} and ; (ii) the (*S*_{A} + *dS*_{A}) isohaline changes position, increasing the salinity of the fluid in both and ; (iii) the (Θ + *d*Θ) isotherm returns to its initial position, cooling the fluid in *V*_{Θ} and back to their original temperatures; and finally (iv) the (*S*_{A} + *dS*_{A}) isohaline returns to its initial position, freshening the fluid in and to their initial states. Averaging over time, we find that the diathermohaline volume transport through the isotherm (isohaline) due to the transport of the fluid in *V*_{Θ } in *S*_{A}–Θ coordinates cancels [Fig. 4a, cycle 1 (2)]. However, has undergone cyclic motion in *S*_{A}–Θ coordinates that does not cancel out over time, thus leaving net diathermohaline volume transports through isohalines and isotherms (Fig. 4b, cycle 5). Note that, for this particular example, the circulation in *S*_{A}–Θ coordinates is established without contributing to the geographical ocean circulation, only to the circulation in *S*_{A}–Θ coordinates.

### c. The local thermohaline streamfunction

We have argued that diathermohaline volume transports can be decomposed into a diathermohaline cycle and trend. The diathermohaline trend represents net changes in the ocean’s volume distribution in *S*_{A}–Θ coordinates. A trend can be the result of 1) the net change of the geographical position of the bounding isohalines and isotherms induced by a net change in local *S*_{A}–Θ values, and 2) a nonzero integral of **u** in the direction normal to the bounding isohalines and isotherms, over the surface area enclosing a volume.

The latter is only of interest for the construction of the advective thermohaline streamfunction, but is not an issue for an instantaneous velocity field in a Boussinesq ocean. Hence, we focus on removing a potential trend in , which does not allow us to represent with a streamfunction. The trend in is due to a time-averaged trend in local changes of *S*_{A}–Θ:

Note that, although the integral results in a constant, |**∇***S*_{A}| and |**∇**Θ| may change in time and space. This results in . Here is the net velocity of the movement of isohalines in geographical space, in the direction normal to themselves, due to a local trend of *S*_{A} over time period Δ*t*, and is the net velocity of isotherms in geographical space, in the direction normal to themselves, due to a local trend of Θ. We can now use (8) to obtain the nondivergent component of the diathermohaline velocity, which can be represented by the local thermohaline streamfunction :

Here, represents the thermohaline volume transport due to cyclic movement of isohalines and isotherms due to local changes in *S*_{A} and Θ. Cyclic refers to isohalines and isotherms (or *S*_{A} and Θ) moving away from their initial state and then subsequently returning back to it, within a given period of time Δ*t* (e.g., Fig. 4b, cycle 5, and Fig. 5). To calculate , we use a detrended ocean hydrography (from either a model or observations). The local changes in *S*_{A}–Θ values need not to be periodic. Depending on the period Δ*t*, they are expected to be dominated by different physical processes as for example, diurnal cycles, seasonal cycles or climate modes such as El Niño–Southern Oscillation (ENSO).

## 5. The diathermohaline streamfunction

In the preceding sections, we have shown that the diathermohaline velocity consists of two components: one related to geographical advection of fluid parcels and one due to local changes of both *S*_{A} and Θ. Following ZD12, we have expressed the advective component as an advective thermohaline streamfunction . In addition, we have expressed the component due to local changes in *S*_{A} and Θ, as the local thermohaline streamfunction . The sum of both thermohaline streamfunctions allows us to construct the diathermohaline streamfunction . This can be shown using only the component of the diathermohaline velocity that is cyclic in *S*_{A}–Θ coordinates, that is, detrended diathermohaline velocity . Making use of (4) and (5) we obtain

Here, represents the cyclic or nondivergent component of the steady-state geographical ocean circulation in *S*_{A}–Θ coordinates. Note that is a diatransport and requires water mass transformation to occur. Therefore, is solely forced by, and can be directly related to, the thermohaline forcing. If there is no thermohaline forcing, we obtain . However, as and are not diatransports, geographical advection of water parcels may lead to cycles and an opposing advective cycle that does not require water mass transformation (and therefore no thermohaline forcing), allowing for . Therefore, circulations shown by and , may also contain a signal driven by geographical advection that is not quantified. This limitation illustrates that one should be careful when interpreting the circulations shown by and as being caused by thermohaline forcing only. Note that the diathermohaline volume transport due to also requires a part of the available thermohaline forcing. This needs to be taken into account when providing an interpretation of .

## 6. Application to model- and observational-based hydrography

In the following, we calculate , , and using output from a climate model and also calculate from an observation-based climatology.

### a. Description of used model and observations

We use the final 10 yr of a 3000-yr spinup simulation of the University of Victoria Climate Model (UVIC). This model is an intermediate complexity climate model with a horizontal resolution of 1.8° latitude by 3.6° longitude, 19 vertical levels, and a 2D energy balance atmosphere [Sijp et al. (2006), specifically their case referred to as GM]. The ocean model is the Geophysical Fluid Dynamics Laboratory Modular Ocean Model, version 2.2 (MOM2), with Boussinesq and rigid-lid approximations applied (Pacanowski 1996). We use the monthly averaged potential temperature *θ*(°C) and practical salinity *S*_{P}.

The observation-based climatology we use is the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Atlas of Regional Seas (CARS) 2009. This is a global, except for the boundary at 75°S, high-resolution (0.5° × 0.5° grid spacing, 79 vertical levels) seasonal atlas of in situ temperature (i.e., *T*), practical salinity (i.e., *S*_{P}), and several other properties (Ridgway et al. 2002; Ridgway and Dunn 2003). At higher latitudes, observations are limited and biased to the summer state. In summary, CARS provides a gridded, time-continuous standard year value of *S*_{P} and *T* on each grid point. To be able to compare the results with UVIC, we use monthly averaged values.

We have used the Thermodynamic Equation of Seawater—2010 (TEOS-10) software (IOC et al. 2010; McDougall and Barker 2011) to convert both the UVIC output and the CARS data to *S*_{A} and Θ. The practical details of the calculations of the different variables from both model and observations are provided in the appendix.

### b. The thermohaline volume transports and diathermohaline trend

The diathermohaline volume transport in *S*_{A}–Θ coordinates for volumes with a certain *S*_{A} and Θ grid size can be calculated and compared for the advective and local thermohaline streamfunction and the diathermohaline trend (Fig. 6). The diathermohaline trend shows maximum values on the order of 1 Sverdrup (Sv; 1 Sv ≡ 10^{6} m^{3} s^{−1}), which is small compared to those of the advective and local streamfunctions, as expected from an equilibrated climate model. Note that, using the Student's *t* test, we found that none of these trends are significant within the 95% level. Hence, for UVIC represents UVIC’s diathermohaline circulation and can be directly related to UVIC’s thermohaline forcing. As CARS is constructed as a standard year, no diathermohaline trend can be calculated.

Climate trends (e.g., long-term warming of the ocean) lead to a permanent change in the positions of isohalines and isotherms, thereby shifting the ocean’s volume toward a different state. Diathermohaline trend terms that do not stem from a climate change complicate the interpretation of the diathermohaline trend. For example, the time period over which we average to analyze the diathermohaline trend is chosen to include an integer amount of closed cycles of the dominant cyclic motions in *S*_{A}–Θ coordinates. Climate variability modes that do not fit in the selected time period will cause a net shift of the ocean’s volume in *S*_{A}–Θ coordinates. Rounding errors due to time and spatial averaging of model or observational data, numerical diffusion, and other possible sources of error also lead to a diathermohaline trend in *S*_{A}–Θ coordinates.

### c. The advective thermohaline streamfunction

The is only calculated for UVIC as CARS does not provide **u**. As we calculate using the same model output as Zika et al. (2012), we only provide a short summary of our results as they are similar to ZD12.

The three dominant circulation cells shown by (Fig. 7) can, in conjunction with the knowledge of ocean basin’s *S*_{A}–Θ properties, be related to a geographical location and associated thermohaline processes. We distinguish (i) a high-latitude cell, related to processes in the Arctic region and the Southern Ocean such as Antarctic Bottom Water formation; (ii) a tropical cell related to the equatorial shallow wind-driven circulation and surface freshwater and heat fluxes; and (iii) the global cell that spans a large range in *S*_{A}–Θ coordinates that can be interpreted as a globally interconnected ocean circulation.

### d. The local thermohaline streamfunction

We have calculated for UVIC (Fig. 7) and CARS (Fig. 8). The majority of the circulation of is explained by near-surface dynamics, as the ocean abyss is not expected to have significant cyclic changes in *S*_{A}–Θ values on subdecadal time scales. For an interpretation of , we will use the ocean’s surface *S*_{A}–Θ values to relate the circulation cells of shown in *S*_{A}–Θ coordinates with geographical locations (Fig. 9).

An anticlockwise-rotating cell is found in an area in *S*_{A}–Θ coordinates associated with equatorial surface waters (Fig. 9). It shows strong freshening due to precipitation, associated with the seasonal variability of the surface freshwater and heat fluxes related to the position and strength of the intertropical convergence zone (ITCZ), and is named the precipitation cell.

A clockwise-rotating cell is mainly related to evaporation over high *S*_{A}-valued water masses (Fig. 9), which are most likely the formation regions of South Pacific and Indian Ocean Subtropical Mode Water (Hanawa and Talley 2001) and high-salinity evaporation areas such as the Mediterranean, Persian Gulf, Red Sea, and tropical Atlantic, and is named the evaporation cell.

There is a clear separation in *S*_{A}–Θ coordinates between the evaporation *E* and precipitation *P* cells and the Southern Hemisphere (SH) and North Pacific (NP) radiative cells, which are dominated by radiative heating and cooling (Fig. 9). The anticlockwise-rotating SH radiative cell stretches over a wide temperature range in *S*_{A}–Θ coordinates. This is due to summer radiative heating and winter radiative cooling in the Southern Ocean and adjacent ocean basins at subtropical latitudes (Fig. 9). The curved shape of the SH radiative cell can be explained by a change in the *E*–*P* fluxes from being dominated by freshening (precipitation dominating evaporation and possible cryospheric processes) at cold temperatures, to being dominated by salinification (evaporation dominating precipitation) at slightly warmer temperatures. As a result the streamfunction changes direction on the *S*_{A} axis, leading to the observed shape. The latter also explains the anticlockwise-rotating NP radiative cell, due to a combination of the cryospheric processes and seasonal radiative heat fluxes. The NP radiative cell is separated from the SH radiative cell by the fresh North Pacific surface waters.

The clockwise-rotating cryosphere cell represents a circulation at high latitudes. We remind the reader that CARS is based on limited high-latitude observations, biased toward summer and extending only to 75°S. This might reduce the magnitude and spread of the cryosphere cell in *S*_{A}–Θ coordinates and complicate the interpretation. However, we suggest the cryosphere cell represents cryospheric processes as we find that warming (cooling) is associated with freshening (salinification), with the minimum influence of seasonal radiative heat fluxes (Fig. 9).

Comparing from CARS (Fig. 8) with that from UVIC (Fig. 7) shows that the model has a reduced spread in *S*_{A}–Θ coordinates. An incorrect model representation of the extremes in the thermohaline coordinates may be caused by 1) excessive near-surface mixing (too much damping), 2) underestimation of the surface freshwater and heat fluxes (not enough forcing), 3) heating and cooling (freshening and salinification) taking place at the same salinity (temperature) because of incorrect modeling of the salt (heat) fluxes resulting in cancellation (Fig. 4a, processes 1, 2, and 3), and 4) an incorrect representation of advection in the model, influencing the unquantified, advective-driven component of .

The reduced spread of the UVIC precipitation cell suggests an incorrect simulation of the ITCZ freshwater and heat fluxes. The evaporation cell also shows a reduced magnitude and spread at higher salinities suggesting a problem with the formation and properties of the South Pacific and Indian Ocean Subtropical Mode Water and the tropical Atlantic surface waters.

Areas covered by the evaporation and precipitation cells in *S*_{A}–Θ coordinates are expected to diverge in time on the *S*_{A} axis, if relatively salty water increases its salinity and relatively freshwater reduces its salinity (Durack et al. 2012). Such dynamics can be analyzed from climatological trends and are related to changes in the surface freshwater and heat fluxes of these particular areas.

In UVIC, both SH and NP radiative cells have merged in *S*_{A}–Θ coordinates, suggesting that the surface freshwater pool in the North Pacific is not correctly modeled and coincides with the SH salinities. Differences in the SH surface freshwater fluxes between CARS and UVIC suggest that the transition between the evaporation- to precipitation-dominated regime for decreasing temperatures for the SH radiative cell is not captured by the model. Finally, the UVIC cryosphere cell has an increased spread and magnitude compared to CARS, which has a bias toward the summer and lacks Arctic data.

## 7. Discussion

We have presented the derivation of the diathermohaline streamfunction , which is composed of an advective thermohaline streamfunction and a local thermohaline streamfunction . The latter is calculated from both a model and observations. The circulation shown by (Fig. 7) enables physically unconnected volumes in the geographical ocean (e.g., in different ocean basins), bounded by similar isohalines and isotherms, to be represented by one single grid in *S*_{A}–Θ coordinates, reducing the three-dimensional spatial and temporal ocean circulation to a two-dimensional time-averaged circulation in *S*_{A}–Θ coordinates, that is, driven by thermohaline forcing. The framework presented here could in principle be applied to any two coordinates defined by two conservative tracers.

To drive the observed circulation of , freshwater and heat fluxes are required. Consequently, part of the available freshwater and heat fluxes will not be available to drive the circulation observed in the . A quantitative insight in the relative importance of and is obtained by calculating the diapycnal salt [*F*_{Salt}(*ρ*_{r})] and heat [*F*_{Heat}(*ρ*_{r})] transports for reference potential density *ρ*_{r}. Following Zika et al. (2012), this is given by

where we have applied an integration along lines of constant reference potential density *ρ*_{r}. The salt flux can be recalculated to an equivalent freshwater flux [*F*_{FW}(*ρ*_{r}) = *F*_{Salt}(*ρ*_{r})/*S*_{ref}], when dividing it by a reference salinity (*S*_{ref} = 35). The resulting diapycnal freshwater and heat transports, for both and , emphasize that requires a significant amount of the total available freshwater and heat fluxes and cannot be ignored when analyzing the ocean’s diathermohaline circulation (Fig. 10). Hence, only the combination of both and provide a diagnostic for modelers to analyze and understand the magnitude, structures, and types of diathermohaline circulation in their models. Finally, in contrast to most other streamfunctions, can also be calculated from observation-based products.

We have used surface *S*_{A}–Θ properties to provide understanding of the observed circulation of . Hieronymus et al. (2013, manuscript submitted to *J. Phys. Oceanogr.*) represented the different components of thermohaline forcing (boundary salt and heat fluxes and the epineutral and dianeutral salt and heat diffusion components) as water mass transformation vectors in *S*_{A}–Θ coordinates. Their study provides a tool to understand in more detail how the diathermohaline circulation provided by is related to the different components of the thermohaline forcing.

Ballarotta et al. (2013) showed that temporal and spatial averaging applied to construct the overturning streamfunction significantly influences the observed circulation and requiring care with the interpretation. Zika et al. (2012) showed that the parameterized eddies have significant influences on . The study of the sensitivity of , , and to spatial and temporal averaging of the model’s *S*_{A}, Θ, and **u** output is beyond the scope of this research. Hence, one can only interpret our results to be valid on seasonal to decadal time scales.

## 8. Conclusions

We have introduced the diathermohaline circulation, driven by boundary freshwater and heat fluxes and mixing processes (thermohaline forcing). The time-averaged nondivergent component of the diathermohaline circulation in *S*_{A}–Θ coordinates is quantified by the diathermohaline streamfunction . Here, is composed of an advective and local component. The advective component is related to the geographical ocean circulation normal to isohaline and isothermal surfaces, as recently introduced in the literature (ZD12) and here quantified by the advective thermohaline streamfunction . The local component is quantified by the local thermohaline streamfunction . Here represents geographical displacements of isohalines and isotherms that, over time, return to their initial position, as a result of detrended local changes in *S*_{A}–Θ values. The term is expressed as the sum of and and provides a model diagnostic that is related directly to thermohaline forcing. The time-averaged unsteady component of the diathermohaline circulation is referred to as the diathermohaline trend and represents permanent, time-averaged water mass changes, related to thermohaline forcing.

Currently, can only be calculated using models as it requires global knowledge of the three-dimensional velocity field (i.e., **u**). As presented can be calculated using model output and an ocean hydrography, allowing for a comparison between model and observed . In future work we will present a technique using an inverse method that allows us to calculate from observational air–sea freshwater and heat fluxes and an ocean’s hydrography, in combination with a representation of ocean-mixing processes [Groeskamp et al. (2014, manuscript submitted to J. Phys. Oceanogr.)].

## Acknowledgments

SG was supported by the joint CSIRO–University of Tasmania program in quantitative marine science (QMS) and the CSIRO Wealth from Oceans flagship and through the Office of the Chief Executive (OCE) Science Team Postgraduate Scholarship Program. BMS was supported by the Australian Climate Change Science Program, jointly funded by the department of Climate Change and Energy Efficiency and CSIRO. JDZ is supported by the U.K. National Environment Research Council.

We thank two anonymous reviewers for their helpful comments that have improved the paper. We thank Willem Sijp for providing the UVIC output and associated discussion, Ken Ridgway for providing useful discussion about the CARS climatology, and we acknowledge numerous discussions with Paul Barker. We thank Kristoffer Döös, Jonas Nycander, Johan Nillson, and Leo Maas for valuable discussions on this work. Finally we thank Walter Munk for inspiration and encouragements to continue this study in an “unfashionable” area of oceanographic research.

### APPENDIX

#### Calculating from Models or Observations

##### a. Discretization processes

We calculate as a sum of and . The streamfunctions are defined on discrete *S*_{A}–Θ intervals. The *S*_{A}–Θ values are at grid points (*x*_{i}, *y*_{j}, *z*_{k}, and *t*_{n}). Here we have *i* = 1–*I*, *j* = 1–*J*, *k* = 1–*K* and *n* = 1–*N*, where *I*, *J*, and *K* indicate the number of grid points in the east, north, and downward direction, respectively, and *N* is the number of time steps. Around each value (*x*_{i}, *y*_{j}, *z*_{k}, and *t*_{n}) one can define six interfaces between adjacent tracer locations, which together enclose a volume Δ*V*_{i,j,k}. We omit the *n* index as Δ*V*_{i,j,k} does not change in time, as the data are on a fixed geographical grid. At each grid point (*x*_{i}, *y*_{j}, *z*_{k}, and *t*_{n}), we have a and Θ_{i,j,k,n} value, located at the center of Δ*V*_{i,j,k}.

In addition to the geographical discretization there is also a discretization in *S*_{A}–Θ coordinates. Consider volume Δ*V*, bounded by a pair of isotherms that are separated by 2*d*Θ and a pair of isohalines that are separated by 2*dS*_{A}. The volume’s Θ ranges between Θ ± *d*Θ and *S*_{A} ranges between *S*_{A} ± *dS*_{A}. While Δ*V* may have any shape in Cartesian coordinates, it covers a square grid in *S*_{A}–Θ coordinates (Fig. 2).

A position in *S*_{A}–Θ coordinates is given by . Here *m* = 1–*M* and *p* = 1–*P*, and *M* and *P* indicate the number of grid points in the *S*_{A} and Θ direction, respectively. We find the isotherms to be at Θ^{p} ± *d*Θ and the isohalines to be at . We argue that *d*Θ and *dS*_{A} must be small enough to distinguish between different water masses in the ocean’s interior and large enough to contain a significant amount of the ocean volume in a grid. For the circulation to be described by an equal resolution in both *S*_{A}–Θ coordinates directions, we adjust the ratio between *d*Θ and *dS*_{A} to be about *dS*_{A} = 0.05 g kg^{−1} and *d*Θ = 0.75°C.

##### b. Calculating from ocean hydrography products

Here we provide the practical method to calculate , which is used to calculate . Let us define as the time-averaged (denoted by the overbar) thermohaline volume transport vector due to the movement of isohalines and isotherms, where is the time-averaged volume transport in the positive *S*_{A} direction through the area of the surface of constant *S*_{A}, which has a Θ range between Θ ± *d*Θ, due to the movement of isohalines induced by local changes of *S*_{A} and Θ and is given by

Here integrates over the area that has a constant salinity of *S*_{A} and Θ range of Θ ± *d*Θ. Equivalently, is the time-averaged volume transport in the positive Θ direction through the area of the surface of constant Θ, which has a *S*_{A} range between *S*_{A} ± *dS*_{A}, due to the movement of isotherms, and is given by

Let it be clear that the size of depends on the grid sizes *dS*_{A} and *d*Θ and is not a thermohaline streamfunction difference, because may include a trend. To calculate from the hydrography, we use the fact that when isohalines and isotherms move in space and time, they change the distribution of the ocean’s volume in *S*_{A}–Θ coordinates. Therefore, to calculate , we will address the change of the ocean’s volume in *S*_{A}–Θ coordinates in time.

In summary, one tracks the path of Δ*V*_{i,j,k} in time, in *S*_{A}–Θ coordinates, and assigns the related volume transport to the crossed isotherms (isohalines) in between the associated *S*_{A} (Θ) range. Repeating this for each time step, for all the defined volumes in the ocean and taking the time average of all contributions to obtain . One needs only the time series of volumes *S*_{A}–Θ values, such that an ocean hydrography is sufficient to calculate , applying the following steps:

Consider Δ

*V*_{i,j,k}and the associated and Θ_{i,j,k,n}, at each grid point (*x*_{i},*y*_{j},*z*_{k}, and*t*_{n}).Then consider Δ

*V*_{i,j,k}and the associated and Θ_{i,j,k,n+1}, at each grid point (*x*_{i},*y*_{j},*z*_{k}, and*t*_{n+1}). Hence, it is the same volume at the same position, but analyzing the tracer values at the next time step, that is,*t*=*t*_{n+1}rather than*t*=*t*_{n}*.*If and/or Θ

_{i,j,k,n}≠ Θ_{i,j,k,n+1}, and as a result of this change, Δ*V*_{i,j,k}has moved to a different*S*_{A}–Θ grid, then there will be motion of Δ*V*_{i,j,k}in*S*_{A}–Θ coordinates in time. The associated volume transport can be calculated as with Δ*t*= (*t*_{n+1}−*t*_{n}). If the change in*S*_{A}and/or Θ is not sufficiently large to move the*S*_{A}–Θ value to the next grid, then this change does not contribute to the volume flux across isotherms or isohalines. With very small values for*dS*_{A}and*d*Θ (i.e., a high resolution), this becomes less important.The volume transport , associated with the displacement of Δ

*V*_{i,j,k}in*S*_{A}–Θ coordinates in time, has to be assigned to the correct*S*_{A}interval of the isotherms that Δ*V*_{i,j,k}has crossed and to the correct Θ interval of the isohalines that Δ*V*_{i,j,k}has crossed. We find these intervals by applying the shortest route method. This method uses a straight line in*S*_{A}–Θ coordinates between the grid that Δ*V*_{i,j,k}occupies at*t*=*t*_{n}and the grid that Δ*V*_{i,j,k}occupies at*t*=*t*_{n+1}(Fig. A1). The volume transport is then assigned to the Θ interval (*S*_{A}interval) on the isohalines (isotherms) that are crossed by this straight line, which is the shortest route in*S*_{A}–Θ coordinates. In case of only warming this is a vertical line, crossing only isotherms. In case of salinification this is a horizontal line, crossing isohalines only. When there is a change in both*S*_{A}and Θ, there will be a volume transport across both isohalines and isotherms. It should be clear that the motion of one Δ*V*_{i,j,k}in*S*_{A}–Θ coordinates at a succession of times can cause a transport assigned to multiple intervals on isohalines and isotherms.

Note that the above has not removed the trend from . To calculate we require

The calculation of is basically similar to that of . By replacing with in (A1) and (A2) and then repeating steps 1–5 as for , while we consider only the motion of Δ*V*_{i,j,k} in *S*_{A}–Θ coordinates between the initial state at *t* = *t*_{1} and the final state at *t* = *t*_{N}, this provides the net transport over the considered time, that is, the diathermohaline trend.

##### c. Calculating from an ocean model

To calculate one requires the three-dimensional velocity from a model. Following Zika et al. (2012), we recognize that tracers in UVIC are advected by the sum of an Eulerian mean velocity **u**^{EM} and quasi-Stokes velocity **u**^{GM}. Hence, we use

We can use this to calculate as follows:

Consider is the volume transport through a grid interface

*A*_{i}in the positive*x*direction [**e**_{1}= (1, 0, 0)] and the associated and Θ^{i,j,k,n}at (*x*_{i},*y*_{j},*z*_{k}, and*t*_{n}).Consider and Θ

_{i+1,j,k,n}at (*x*_{i+1},*y*_{j},*z*_{k}, and*t*_{n}).If and/or Θ

_{i,j,k,n}≠ Θ_{i+1,j,k,n}, and the resulting combination of*S*_{A}–Θ values occupies a different*S*_{A}–Θ grid, then the volume transported through*A*_{i}, that is, , is normal to isotherms and/or isohalines. The volume transport will then be assigned to the correct*S*_{A}interval of the crossed isotherms and to the correct Θ interval of the crossed isohalines, using the shortest route method.

##### d. Calculating and

As and are no longer continuous, neither will and be continuous. A volume transport is defined as the total transport through a Θ range (*S*_{A} range) on an isohaline (isotherm). This is exactly between the “corners” of an *S*_{A}–Θ grid. The coordinates of the corners of the grid, with at the center, are given by for the top-right corner, for the bottom-right corner, for the top-left corner, and for the bottom-left corner.

For a grid with at its center, we define diahaline and diathermal transports in positive Θ and *S*_{A} directions. Thus, the diahaline transport takes place at and in between Θ^{p} ± *d*Θ. Thus, the diahaline volume transport is given by . Equivalently, the diathermal volume transport in the positive Θ direction, for the grid with at its center, at Θ^{p} + *d*Θ in between salinity range , is given by (Fig. A2).

As is calculated using a summation of that only adds up at discrete locations, which are defined to be exactly those corner points. Hence, we can only define at the corner points. Calculating using our discretization, leads to

This summation reduces to an integration again when Δ*S*_{A} and ΔΘ are taken as infinitesimally small. A similar exercise obtains :

## REFERENCES

*Fundamentals of Ocean Climate Models.*Princeton University Press, 518 pp.

*Ocean Circulation and Climate: Observing and Modelling the Global Ocean,*International Geophysics Series, Vol. 77, Academic Press, 373–386.

**71,**916–928, doi:.

*J. Phys. Oceanogr.,*doi:, in press.

*International WOCE Newsletters,*No. 31, WOCE International Project Office, Southampton, United Kingdom, 36–39.