Global dataset of thermohaline staircases obtained from Argo floats and Ice-Tethered Profilers

Thermohaline staircases are associated with double-diffusive mixing. They are characterized by stepped structures consisting of mixed layers of typically tens of metres thick that are separated by much thinner interfaces. Through these interfaces enhanced diapycnal salt and heat transport take place. In this study, we present a global dataset of thermohaline staircases derived from observations of Argo profiling floats and Ice-Tethered Profilers using a novel detection algorithm. To establish the presence of thermohaline staircases, the algorithm detects subsurface mixed layers and analyses the interfaces in between. Of each detected staircase, the conservative temperature, absolute salinity, depth, and height, as well as some other properties of the mixed layers and interfaces, are computed. The algorithm is applied to 487 493 quality-controlled temperature and salinity profiles to obtain a global dataset. The performance of the algorithm is verified through an analysis of independent regional observations. The algorithm and global dataset are available at https://doi.org/10.5281/zenodo.4286170.

1 Introduction

Thermohaline staircases consist of subsurface mixed layers that are separated by thin interfaces. They are associated with double-diffusive processes, which in turn result from a difference of 2 orders of magnitude between the molecular diffusivity of heat and that of salt (Stern, 1960). Whenever the vertical gradients of temperature- and salinity-induced stratification have the same sign, these differences in molecular diffusivity can enhance the vertical mixing through double-diffusive convection, leading to effective diffusivities of the order of 10⁻⁴ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (Radko, 2013, and references therein).

It is still a topic of discussion how double-diffusive convection leads to the formation of thermohaline staircases in oceanic environments (Merryfield, 2000). For example, Stern (1969) argued that small-scale mixing processes trigger the formation of internal waves. On the other hand, variations in the turbulent heat and salt fluxes (Radko, 2003) or in the counter-gradient buoyancy fluxes that sharpen density gradients (Schmitt, 1994) could also lead to the formation of thermohaline staircases. Lastly, subsurface mixed layers can also arise from thermohaline intrusions (Merryfield, 2000). Although it remains unclear how these staircases arise, these studies agree that the formation of these subsurface mixed layers are related to double-diffusive processes.

Based on the Turner angle (Tu), which compares the density component of the temperature distribution with the density component of the salinity distribution, two regimes of double diffusion can be distinguished (Ruddick, 1983). Waters with $-\mathrm{90}{}^{\circ }<\text{Tu}<-\mathrm{45}{}^{\circ }$ correspond to a stratification where both temperature and salinity increase with depth and belong to the diffusive-convective regime (DC). Those with $\mathrm{45}{}^{\circ }<\text{Tu}<\mathrm{90}{}^{\circ }$ correspond to a stratification where temperature and salinity decrease with depth and belong to the salt-finger regime (SF).

Theoretical and laboratory studies have indicated that diapycnal fluxes of heat and salt in thermohaline staircases are elevated compared to the background turbulence (e.g., Schmitt, 1981; Kelley, 1990; Radko and Smith, 2012; Garaud, 2018). These results were confirmed by a tracer release experiment in the western tropical Atlantic Ocean (Schmitt, 2005). Although these enhanced fluxes were observed, the importance of these fluxes for the global mechanical energy budget remain unknown. Moreover, the vertical heat and salt fluxes in thermohaline staircases can also affect water-mass properties. In some regions, persistent thermohaline staircases with layers stretching over a few hundred kilometres have been observed (Schmitt et al., 1987; Timmermans et al., 2008; Shibley et al., 2017), which could result in significant diapycnal fluxes between water masses. For example, the double-diffusive diapycnal fluxes in the Mediterranean Sea dominate the transport between the deep water masses (Zodiatis and Gasparini, 1996; Bryden et al., 2014; Schroeder et al., 2016), and in the Arctic Ocean and Southern Ocean, an upward heat flux has been observed through staircase interfaces (Timmermans et al., 2008; Shibley et al., 2017; Polyakov et al., 2012; Bebieva and Speer, 2019).

Modelling studies that incorporated parameterizations of double-diffusive fluxes indicated that the associated double-diffusive diapycnal fluxes can reduce the strength of the global overturning circulation (Gargett and Holloway, 1992; Merryfield et al., 1999; Oschlies et al., 2003). To be able to study this with observations, we present a global dataset of the occurrence of thermohaline staircases and their properties. The dataset is based on observations from Argo floats and Ice-Tethered Profilers. In the following sections we briefly describe the raw data used to extract the dataset (Sect. 2) and the algorithm we designed to detect staircase structures (Sect. 3). The sensitivity of this detection algorithm to the chosen input parameters is assessed in Sect. 4. The dataset is verified in Sect. 5, followed by some guidelines for the use of the dataset in Sect. 7.

2 Data preparation

The dataset contains observations of autonomous Argo floats and autonomous Ice-Tethered Profilers (ITPs). The data of all active and inactive profilers are obtained from http://www.argo.ucsd.edu (last access: 14 May 2020) and http://www.whoi.edu/itp (last access: 14 May 2020) from 13 November 2001 to 14 May 2020. Details on the profilers are described in Krishfield et al. (2008) and Toole et al. (2011) for the ITP and in Argo (2020) for the Argo floats. First a quality check is performed, where a profile is excluded from analysis if it was taken by an Argo float mentioned on the grey list. This grey list contains floats that may have problems with at least one of the sensors (https://www.nodc.noaa.gov/argo/grey_floats.htm, last access: 14 May 2020). As thermohaline staircases consist of mixed layers with depths of tens of metres, we also require that profiles have continuous data up to 500 dbar with an average resolution finer than 5 dbar. Details on the origin and vertical resolution of the profiles are depicted in Table 1 and Fig. 1, in which Fig. 1b confirms that all profiles have observations deeper than 500 dbar. Furthermore, the average vertical resolution of the profiles indicates the average resolution is well below the 5 dbar that was used as a threshold (Fig. 1c). After this quality control, 487 493 vertical temperature and salinity profiles remain. Their global distribution is shown in Fig. 2.

Table 1 Number of floats and profiles in the global dataset. Profiles taken with Argo floats are categorized by the Data Assembly Centre (DAC). Profiles taken with Ice-Tethered Profilers are categorized as ITP. The percentage between brackets indicates the relative contribution to the total number of profiles in the global dataset (487 493 profiles). More details on abbreviations of DAC can be found in Argo (2019).

Figure 1 (a) Locations of observations categorized by Data Assembly Centres (DAC) when obtained by an Argo float. Profiles obtained with Ice-Tethered Profilers are indicated with ITP. (b) Cumulative fraction of profiles that reached a given pressure in 25 dbar intervals from 0 to 2000 dbar per DAC. (c) Average number of observations in 25 dbar intervals from 0 to 2000 dbar. (d) Distribution of detected mixed-layer pressures in the salt-finger (red histogram) or diffusive-convective (blue histogram) regime. (e) Number of detected mixed-layer height in the salt-finger (red histogram) or diffusive-convective (blue histogram) regime. (f) Distribution of detected mixed-layer heights in thermohaline staircases per pressure level. Panels (b, c) were obtained following Wong et al. (2020). Black lines indicate the averages in the total global dataset. More details on abbreviations of DAC can be found in Argo (2019).

Figure 2 Observation density of the number of profiles obtained from the Argo floats and Ice-Tethered Profilers after quality control (km⁻²). Observation density is binned per degree longitude and degree latitude. Empty bins indicate that no data were available at that location.

Table 2 Input parameters applied during the data preparation and the algorithm as used in this study. The sensitivity of the output of the algorithm to the input variables is discussed in the Sect. 4.

Figure 3 Schematic of a typical temperature profile with staircases, indicating the definitions of the quantities used to detect the thermohaline staircases (green: mixed layer; orange: interface).

Figure 4 Histogram of the number of detected interfaces as a function of the Turner angle (Tu) by applying a criterion for (a) conservative temperature, (b) absolute salinity, (c) potential density, and (d) all three properties given in Eq. (4) (orange shading). Each panel shows the data remaining compared to the raw interface data (grey). Vertical shaded bands correspond to Turner angles in the diffusive-convective (blue) and salt-finger (red) regimes.

Next, the profiles of the Argo floats and ITP were linearly interpolated to a vertical resolution of 1 dbar from the surface to 2000 dbar so that their data could be analysed in a consistent manner. As a result, the small steps in, for example, Arctic staircases might be missed (see Sect. 5). From these interpolated profiles we calculate several variables. Absolute salinity (S) in grammes per kilogrammes and conservative temperature (T) in degrees Celsius are computed with the TEOS-10 software (McDougall and Barker, 2011). Note that we use conservative temperature as this is more accurate than potential temperature in computations concerning heat fluxes and heat content (Graham and McDougall, 2013). We apply a moving average of 200 dbar (Table 2) to obtain the background conservative temperature and absolute salinity profiles of the water column and to compute the thermal expansion coefficient (α, ${}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) and the haline contraction coefficient (β, kg g⁻¹). A consequence of the moving average of 200 dbar is that the upper 100 dbar and lower 100 dbar of each profile are omitted in the remainder of the analysis. The Turner angle is computed using profiles that were smoothed with a moving average of 50 dbar instead of 200 dbar, which is similar to Shibley et al. (2017), following Ruddick (1983), from

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\text{Tu}={\tan }^{-\mathrm{1}}\left(\mathit{\alpha }{\displaystyle \frac{\partial T}{\partial p}}-\mathit{\beta }{\displaystyle \frac{\partial S}{\partial p}},\mathit{\alpha }{\displaystyle \frac{\partial T}{\partial p}}+\mathit{\beta }{\displaystyle \frac{\partial S}{\partial p}}\right),\end{array}$$

where the vertical gradients are approximated with a central differences scheme.

3 Detection algorithm

After the data pre-processing, we apply a detection algorithm that exploits the vertical structure of staircase profiles (Fig. 3). The benefit of using the vertical structure, instead of using assumptions based on the Turner angle, is that we can use this angle to verify the results. The detection algorithm consists of five steps. First the algorithm detects all data points that are located in the subsurface mixed layers (ML, green dots in Fig. 3) by identifying weak vertical density gradients in conservative temperature and absolute salinity. Next, the properties of any layer lying between the mixed layers (the interfaces, IF, orange dots in Fig. 3) are assessed by applying a minimum in temperature and salinity variations. Third, the height of the interface and variations within the interface are limited. The fourth step determines the regime of double diffusion (diffusive convection or salt fingers), and the fifth step is the identification of sequences of interfaces, which eventually characterizes the thermohaline staircases. The different steps of the algorithm applied to three example profiles are shown in Figs. A1–A3. In the following subsections, each algorithm step is described in more detail.

Figure 5 Histogram of the number of detected interfaces as a function of the Turner angle (Tu) by applying a criteria for (a) height, (b) maximum height, (c) inversions, and (d) all three height limitations (yellow shading). Each panel shows the data remaining compared to the interfaces detected based on the conservative temperature and absolute salinity requirements shown in Fig. 4d (orange shading). Vertical shaded bands correspond to Turner angles in the diffusive-convective (blue) and salt-finger (red) regimes.

Figure 6 Histogram of the number of detected interfaces as a function of the Turner angle (Tu) after (a) classification of the double-diffusive regime and (b) selection of sequences of the interfaces. Each panel shows the data remaining compared to the interfaces detected based on interface height requirement shown in Fig. 5d (yellow shading). Vertical shaded bands correspond to Turner angles in the diffusive-convective (blue) and salt-finger (red) regimes.

Figure 7 Example conservative temperature profiles selected by the staircase detection algorithm. They are ordered left to right by the number of steps detected. Panel (a) shows examples of increasing steps of diffusive convection. Panel (b) shows examples of the salt-finger regime.

Figure 8 Number of detected interfaces obtained with the detection algorithm for different input parameters. Each panel shows the sensitivity of the detection algorithm to one input parameter: (a) moving average window, (b) density gradient of the mixed layer, (c) density difference within the mixed layer, and (d) the maximum height of the interface. In each panel, the grey histogram corresponds to the default parameters listed in Table 2. The coloured lines correspond to the varying parameter (see legend). Shaded regions indicate Turner angles in the diffusive-convective (blue) and salt-finger (red) regimes.

Padman and Dillon (1987)Timmermans et al. (2003)Timmermans et al. (2008)Shibley et al. (2017)

Table 3 Characteristics of thermohaline staircases in Canada Basin. The region of the global dataset is confined to 75–80^∘ N, 135–145^∘ W. The observational techniques indicate whether the data were obtained from Argo floats (Argo), Ice-Tethered Profilers (ITPs), conductivity–temperature–depth (CTD) measurements, or microstructure measurements (MSs). The dominant type of thermohaline staircases is indicated by DC (diffusive convection) and SF (salt-finger) with the percentage of occurrence between brackets. Ranges of the obtained variables of the global dataset are indicated by means of percentiles 2.5 and 97.5.

Zodiatis and Gasparini (1996)Bryden et al. (2014)Buffett et al. (2017)Durante et al. (2019)

Table 4 As Table 3, but for the Mediterranean Sea (30–43^∘ N, 0–15^∘ E).

Schmitt et al. (1987)Schmitt (2005)Fer et al. (2010)

Table 5 As Table 3, but for the western tropical North Atlantic Ocean (10–15^∘ N, 53–58^∘ W).

3.1 Mixed layers

The first step of the detection algorithm is the identification of the mixed layers. Preferably, this is done by assessing a density difference relative to a reference pressure, which is the most reliable method to detect a mixed layer (Holte et al., 2017). However, in the case of a thermohaline staircase, it is necessary to detect subsurface mixed layers, because the reference pressure is unknown beforehand. To determine this reference pressure, a threshold gradient criterion is applied first (Dong et al., 2008). In this criterion, vertical density gradients are identified as a mixed layer whenever the gradients are below a certain threshold.

We apply the gradient criterion to the vertical gradients of the potential density anomaly at a reference pressure of 1000 dbar (σ₁). We used a threshold of $\partial {\mathit{\sigma }}_{\mathrm{1}}/\partial p_{\text{max}}=\mathrm{0.0005}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{dbar}}^{-\mathrm{1}}$ (Table 2), which is similar to mixed-layer gradients used by Bryden et al. (2014). Furthermore, this threshold gradient is slightly larger than the threshold used by Timmermans et al. (2008), who used 0.005 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ (which corresponds to $\partial {\mathit{\sigma }}_{\mathrm{1}}/\partial p_{\text{max}}=\mathrm{0.00036}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{dbar}}^{-\mathrm{1}}$). The threshold gradient method is applied to both conservative temperature and absolute salinity profiles, i.e.,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left|\mathit{\alpha }{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\partial T}{\partial p}}\right|\le \mathrm{0.0005}{\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\mathrm{dbar}}^{-\mathrm{1}},\\ \text{(2)}& {\displaystyle }& {\displaystyle }\left|\mathit{\beta }{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\partial S}{\partial p}}\right|\le \mathrm{0.0005}{\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\mathrm{dbar}}^{-\mathrm{1}}.\end{array}$$

Also the vertical density gradients from the combined temperature and salinity effects must satisfy this condition:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\left|{\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{1}}}{\partial p}}\right|\le \mathrm{0.0005}{\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\mathrm{dbar}}^{-\mathrm{1}}.\end{array}$$

These three conditions ensure that the vertical conservative temperature, absolute salinity, and potential density gradients are all below the threshold value. At each pressure level where all three conditions are met, the data point is identified as a mixed layer. Next, for each continuous sequence of data points, the algorithm computes the average pressure. This is then used as a reference pressure, which is required to be able to apply the mixed-layer detection.

At every reference pressure, a maximum density range is required within the mixed layers to identify the full vertical extent of each mixed layer. To allow for small variations in conservative temperature and absolute salinity in the mixed layer, but to exclude variations in the interface, we use a threshold of $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{1},\text{ML,max}}=\mathrm{0.005}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for density variations within each mixed layer (Table 2). This density range corresponds to the density range used by Holte et al. (2017) for the detection of surface mixed layers. The applied density range allows for mixed layers with heights of the order of 10 m assuming gradients of $\partial {\mathit{\sigma }}_{\mathrm{1}}/\partial p_{\text{max}}=\mathrm{0.0005}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{dbar}}^{-\mathrm{1}}$. To ensure separation between individual mixed layers, the upper and lower data points of each mixed layer are removed. Note that this results in a minimum interface height of 2 dbar, which could result in false negatives in for example the Arctic Ocean (Sect. 5)

After applying the threshold for density range, the algorithm defines each continuous set of data points as a mixed layer and computes the average pressure (${\stackrel{\mathrm{‾}}{p}}_{\text{ML}}$), average conservative temperature (${\stackrel{\mathrm{‾}}{T}}_{\text{ML}}$), average absolute salinity (${\stackrel{\mathrm{‾}}{S}}_{\text{ML}}$), mixed-layer density ratio $\left({\stackrel{\mathrm{‾}}{R}}_{\mathit{\rho }}=\mathit{\alpha }\frac{\partial \stackrel{\mathrm{‾}}{T}}{\partial p}/\left(\mathit{\beta }\frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial p}\right)\right)$, average Turner angle (${\stackrel{\mathrm{‾}}{\text{Tu}}}_{\text{ML}}$), and height (h_ML) for each mixed layer.

3.2 Interfaces: conservative temperature and absolute salinity variations

The algorithm defines an interface as the part of the water column between two mixed layers. In addition, to ensure a stepped structure, the algorithm requires that the conservative temperature, absolute salinity, and potential density variations within each mixed layer should be smaller than the variations in the interface (Fig. 3):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }max\left(\left|\mathrm{\Delta }{T}_{\text{ML,1}}\right|,\left|\mathrm{\Delta }{T}_{\text{ML,2}}\right|\right)<\left|\mathrm{\Delta }{T}_{\text{IF}}\right|;\\ {\displaystyle }& {\displaystyle }max\left(\left|\mathrm{\Delta }{S}_{\text{ML,1}}\right|,\left|\mathrm{\Delta }{S}_{\text{ML,2}}\right|\right)<\left|\mathrm{\Delta }{S}_{\text{IF}}\right|;\\ \text{(4)}& {\displaystyle }& {\displaystyle }max\left(\left|\mathrm{\Delta }{\mathit{\sigma }}_{\text{1,ML,1}}\right|,\left|\mathrm{\Delta }{\mathit{\sigma }}_{\text{1,ML,2}}\right|\right)<\left|\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{1},\text{IF}}\right|;\end{array}$$

where the subscripts 1 and 2 correspond to the mixed layer directly above and below an interface, respectively. It appears that most data points that meet these requirements (orange histograms in Fig. 4a–c) have Turner angles in the two double-diffusive regimes. This dependence of the variations in the interfaces on the Turner angle is in line with expectations that staircase-like structures are mostly found within double-diffusive regimes. In total, 28 % of all detected interfaces meet all three requirements (Fig. 4d).

3.3 Interface: height

The next step in the staircase detection algorithm is to limit the height of the interface to ensure that the mixed layers are separated from each other by a relatively thin interface (Fig. 3). We require

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{h}_{\text{IF}}<min\left(h_{\text{ML,1}},h_{\text{ML,2}}\right);\end{array}$$

i.e. the interface height is smaller than the height of the mixed layers directly above and below the interface. In total, 27 % of the interfaces that fulfilled the conservative temperature and absolute salinity requirements meet this requirement (Fig. 5a). Note that this part of the algorithm defines the top and bottom of a sequence of a staircase in a profile. Furthermore, the tallest observed interfaces are found in the Mediterranean Sea with heights up to h_IF=27 m, where they separate mixed layers of over 100 m (Zodiatis and Gasparini, 1996; Radko, 2013). To prevent false detection of large vertical interfaces of up to hundreds of metres, we limit the interface height to h_IF,max=30 dbar (Table 2, Fig. 5b). This only affects the classification of 1 % of the interfaces (Fig. 5b).

To solely detect step-like structures that are associated with the presence of thermohaline staircases, the algorithm also removes all interfaces with conservative temperature or absolute salinity inversions. This is done by limiting the number of local minima and maxima of the conservative temperature and absolute salinity allowed in each interface to two (Fig. 5c). The combination of all three interface height requirements is met by 27 % of the interfaces detected based on the conservative temperature and absolute salinity requirements discussed in the previous section (Fig. 5d).

3.4 Interface: double-diffusive regime

After the algorithm has selected all the interfaces with a step-like structure, the double-diffusive regime of each interface is assessed (Fig. 6a). In case both conservative temperature and absolute salinity of the mixed layers above and below the interface increase with pressure, the interface is classified as the diffusive-convective regime. If the conservative temperature and absolute salinity of the mixed layers above and below the interface both decrease with pressure, the interface belongs to the salt-finger regime. The algorithm detects more interfaces in the salt-finger regime (27 %) than in the diffusive-convective regime (11 %, Fig. 6a). As expected, most interfaces with diffusive-convective characteristics have Turner angles of $-\mathrm{90}{}^{\circ }<\text{Tu}<-\mathrm{45}{}^{\circ }$ (blue histogram in Fig. 6a) and most salt-finger interfaces have Turner angles of $\mathrm{45}{}^{\circ }<\text{Tu}<\mathrm{90}{}^{\circ }$ (red histogram in Fig. 6a). This implies that these interface properties are consistent with the background stratification.

3.5 Sequences of interfaces

The final step of the detection algorithm is to only select vertical sequences of at least two interfaces in the same double-diffusive regime that are separated from each other by one mixed layer (Fig. 6b). This step removes most thermohaline intrusions, as these are characterized by alternating mixed layers in the diffusive-convective and salt-finger regimes (Bebieva and Timmermans, 2017). In this final step, the algorithm also removes salt-finger interfaces and diffusive-convective interfaces outside their favourable Turner angle (compare Fig. 6a and b).

After applying this final step of the algorithm, we obtain a global dataset consisting of 166 141 interfaces in the salt-finger regime and 119 619 interfaces in the diffusive-convective regime. The distribution of the pressure levels and height of the mixed layers at these interfaces are displayed in Fig. 1. In general, mixed layers in the diffusive-convective regime are found at lower pressure levels than mixed layers in the salt-finger regime (Fig. 1d). At the same time, the height of the mixed layers in the diffusive-convective regime are smaller, which is in line with previous observations (Fig. 1e, e.g., Radko, 2013). Recall that the algorithm required a minimal interface height of 2 dbar, which implies that, following Eq. (5), the minimal mixed-layer height is 3 dbar and the detection of interfaces is cut off below these limits. Consequently, the interfaces with smaller heights are missed by the algorithm. Figure 1e indicates that this is more problematic for interfaces in the diffusive-convective regime than for interfaces in the salt-finger regime.

Examples of thermohaline staircases, which were selected based on their high number of interfaces, are shown in Fig. 7. In line with previous results (Rudels, 2015), staircases in the diffusive-convective regime (Fig. 7a) are mainly detected on the thermocline with the conservative temperature increasing with depth. These staircases are predominantly located in the Arctic Ocean at a depth between 300–400 m, which is between the warm and saline Atlantic Water and cold and fresh surface waters (Rudels, 2015). Figure 7a also indicates that the deepest mixed layer of some thermohaline staircases is located at the temperature maximum, which suggests that this lowest layer might be the result of thermohaline intrusions (Ruddick and Kerr, 2003). There, the algorithm identified a mixed layer, because temperature and salinity stratification were weak enough (see Sect. 3.1). Furthermore, both conservative temperature and absolute salinity in this mixed layer are larger than in the mixed layer above. While both are typical for a staircase in the diffusive-convective regime, the algorithm does not detect whether this mixed layer is a temperature maximum, which could indicate that it arose from thermohaline intrusions. Note that this only concerns the deepest mixed layers of the staircases and that only the characteristics of the interfaces in between mixed layers are labelled as part of a staircase by the algorithm.

Thermohaline staircases with a high number of steps in the salt-finger regime are detected on the main thermocline where the conservative temperature decreases with depth (Fig. 7b). Compared to the staircases in the diffusive-convective regime, these staircases are located slightly deeper at 400–700 m. While the locations of these staircases vary, they are located above the cold and fresh Antarctic Intermediate Water, which is observed below 700 m (Tsuchiya, 1989; Fine, 1993; Talley, 1996).

For each thermohaline staircase, characteristics of the interfaces and mixed layers, such as their conservative temperature, absolute salinity, and height, are available in the dataset. An overview of the provided variables is given in Table A1. The detection algorithm is verified by comparing our data to independent observations in three regions in Sect. 5.

4 Robustness of the detection algorithm

The algorithm requires four input parameters: the moving average window, a threshold for the maximum density gradients of the mixed layers, the maximum density difference of the mixed layers, and the maximum height of the interface (Table 2). In this section, the sensitivity of the algorithm to each input parameter is assessed (Fig. 8).

The moving average window is used by the algorithm to compute the thermal expansion coefficient (α), the haline contraction coefficient (β), and the density ratio (R_ρ). We varied the moving average window between 50 and 350 dbar to assess the sensitivity of the outcomes of this choice (Fig. 8a). We find that the varying moving average window does not result in large variations in detected mixed layers (Fig. 8a).

In contrast to the moving-average window, the detection algorithm is sensitive to the value set for the density gradient threshold of the mixed layer (Fig. 8b), which is used to obtain a reference pressure for the sub-surface mixed layers (Sect. 3.1). Not surprisingly, we detect more (fewer) interfaces when we increase (decrease) the allowed threshold density gradient. A small value allows for only the strongest mixed layers to be detected, which are usually referred to as well-defined staircases, while a large density gradient also allows for the detection of rough staircases (e.g., Durante et al., 2019). Although the number of detected interfaces depends on the value set for this density gradient, the detected interfaces remain confined to the two double-diffusive regimes, indicating a robust outcome of the algorithm for the choice of this input parameter.

Similar to the variations in the maximum density gradient, the variation in the maximum density difference allowed within the mixed-layer results in a different number of detected interfaces (Fig. 8c). The number of detected mixed layers increases when we decrease the maximum density difference allowed within the mixed layer. This effect is mostly visible in the diffusive-convective regime, as we obtained a decrease of 54 % of detected interfaces in the diffusive-convective regime compared to a decrease of 31 % of detected interfaces in the salt-finger regime in the case in which we doubled the density difference in the mixed layer ($\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{1},\text{max}}=\mathrm{10}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). This difference between the regimes is due to relatively small interface variations in the diffusive-convective regime compared to the salt-finger regime (Radko, 2013) and can be explained as follows: when a too large density difference is applied, the relatively small density gradients in the interfaces of the diffusive-convective regime are detected as mixed layers by the algorithm. Consequently, multiple mixed layers can be identified as a single mixed layer. However, if the applied density difference is too small, this could result in the detection of multiple mixed layers per staircase step.

The last input parameter of the detection algorithm concerns the interface height (Fig. 8d). As expected from Fig. 5b, variations in this input parameter do not result in large differences in the number of detected interfaces. If we omit this input parameter by setting it to infinity, we obtain a total increase in detected interfaces of 17 %.

Overall, the detection algorithm gives robust results as it predominantly detects interfaces within the double-diffusive regime (Fig. 8). In line with expectations, the detection algorithm is most sensitive to the threshold value for the maximum density gradient in the mixed layer and the density variations within the mixed layers. The four input variables allow for optimization of the detection algorithm based on the regime and characteristics of the staircases.

5 Regional verification

The characteristics of thermohaline staircases obtained with the detection algorithm are compared to those obtained from previous observational studies for three major staircase regions: the Canada Basin in the Arctic Ocean, the Mediterranean Sea, and the C-SALT region in the tropical Atlantic Ocean. An overview is given in Tables 3–5.

In the Canada Basin (75–80^∘ N, 135–145^∘ W), the algorithm detects thermohaline staircases in the diffusive-convective regime in 90 % of the profiles (Table 3). Both the occurrence and depth range are comparable to what was reported by Timmermans et al. (2008) and Shibley et al. (2017), who analysed thermohaline staircases from several Ice-Tethered Profilers, demonstrating that our detection algorithm indeed detects thermohaline staircases at the right location. Microstructure observations suggested that the thermohaline staircases in Canada Basin have interface heights of approximately h_IF=0.15 m (Padman and Dillon, 1987; Radko, 2013). Due to the vertical resolution of the profiles and the design of the algorithm (recall that the mixed layers are separated from each other by removing the upper and lower data points of the mixed layer, Sect. 3.1), the method is not capable of detecting very thin interfaces (Fig. A1). As expected from these limitations for the detection of the interface heights, the algorithm detects conservative temperature and absolute salinity steps (ΔT_IF and ΔS_IF, respectively) in the interfaces that are in the upper ranges of earlier observations (Padman and Dillon, 1987; Timmermans et al., 2003, 2008; Shibley et al., 2017). In the Mediterranean Sea, thermohaline staircases are characterized by relatively thick mixed layers that are separated by thick interfaces of up to 27 m (Zodiatis and Gasparini, 1996). In this region (30–43^∘ N, 0–15^∘ E), the detection algorithm detected thermohaline staircases with interfaces up to 21 dbar in 6 % of the profiles, which is comparable to previous observations (Table 4). An example of the detection of a Mediterranean staircase is shown in Fig. A2. We find that the depth at which the thermohaline staircases occur is underestimated by the detection algorithm. This could be explained by the fact that most Mediterranean observations are obtained by the Coriolis DAC (Fig. 1a). From this DAC, approximately 50 % of the profiles have observations that are deeper than 1000 dbar (Fig. 1b), which means that the coverage below 1000 dbar is limited in the Mediterranean Sea. Although the Argo floats, and consequently the detection algorithm, do not cover the full extent of the staircases (Fig. 1), the conservative temperature and absolute salinity steps that are found are similar to previous observations (Table 4). Note that the conservative temperature and absolute salinity steps of the staircases increase with depth (Zodiatis and Gasparini, 1996), which explains why the conservative temperature and absolute salinity steps detected by the algorithm are slightly smaller than those observed in the deeper observations (Zodiatis and Gasparini, 1996; Durante et al., 2019).

In the C-SALT region in the western tropical North Atlantic Ocean (10–15^∘ N, 53–58^∘ W), the algorithm detected thermohaline staircases in the salt-finger regime in 60 % of the profiles (Table 5). Similar to previous studies (Schmitt et al., 1987; Schmitt, 2005; Fer et al., 2010), the algorithm detected thermohaline staircases on the main thermocline (see example in Fig. A3). Again, the interface height is slightly overestimated by the detection algorithm, but the algorithm obtained conservative temperature and absolute salinity steps comparable to previous studies.

Overall, the comparison between the outcomes of the detection algorithm with previous studies indicates that the detection algorithm performs well. The small overestimation of the interface height can be attributed to the limited vertical resolution and the limitation imposed by the detection algorithm to avoid detection of false positives. Despite this overestimation, the interfaces are detected at the correct depths with conservative temperature and absolute salinity steps within realistic ranges. Therefore, we conclude that the detection algorithm is very suitable for the automated detection of thermohaline staircases in large and quickly growing datasets like the Argo float and Ice-Tethered-Profiler data.

6 Code and data availability

Both algorithm and global dataset are available at DOI: https://doi.org/10.5281/zenodo.4286170 (van der Boog et al., 2020). The algorithm is written in Python 3 and is available under the Creative Commons Attribution 4.0 License. More details on the functions and output of the algorithm are depicted in Tables A1 and A2, respectively. The structure of the algorithm is displayed in Fig. A4. The Argo Program is part of the Global Ocean Observing System (Argo, 2020).

7 Conclusions

In this study, we presented an algorithm to automatically detect thermohaline staircases from Argo float profiles and Ice-Tethered Profiles. As these thermohaline staircases have different mixed-layer heights and temperature and salinity steps across the interfaces in different staircase regions, the design of the detection algorithm is based on the typical vertical structure and shape of the staircases (Figs. 3–5). Note that by formulating the algorithm solely on this vertical structure of the staircases, we could use the Turner angle of the detected staircases for verification. Using this Turner angle, we showed that the structures are within the two double-diffusive regimes: the salt-finger regime and the diffusive-convective regime (Fig. 6).

We optimized the input of the algorithm such that it provides a global overview and limits the number of detected false positives. As a result, the regional verification in Sect. 5 indicated that the data pre-processing and data analysis have some limitations. For example, the vertical resolution of 1 dbar in the profiles is too course to capture all staircase steps in the Arctic Ocean. In the Mediterranean, the Argo floats did not dive deep enough to capture the full depth of the staircase region. However, the fact that (i) the algorithm detects thermohaline staircases at realistic depth ranges, with (ii) conservative temperature and absolute salinity steps across the interfaces, and in (iii) the same double-diffusive regime as previous studies (Tables 3–5), indicates that the algorithm itself performs well. Therefore, when considering an individual staircase region, we recommend optimizing the input variables of the algorithm for that specific region and applying the algorithm to additional data, for example high-resolution CTD or microstructure profiles, where available.

A sensitivity analysis to different input parameters showed that the results of the detection algorithm are robust; the detected staircase interfaces are confined to the double-diffusive regimes. Furthermore, the comparison between the detected interface characteristics of thermohaline staircases in three prevailing staircase regions and previous observations suggested that the detection algorithm accurately captures both double-diffusive regimes. The algorithm detected correct magnitudes of the conservative temperature and absolute salinity steps in the interfaces, which allows for adequate estimates of the effective diffusivity in thermohaline staircases.

The global dataset resulting from the detection algorithm contains properties and characteristics of both mixed layers and interfaces. Combined with their locations, these data allow for a statistical analysis of thermohaline staircases on global scales. For example, the global occurrence of thermohaline staircases could give insight into the contribution of double diffusion to the global mechanical energy budget. Moreover, the interface characteristics can be used to validate model and laboratory results on how double-diffusive mixing impacts the regional ocean circulation.

Appendix A

Figure A1 Steps of the detection algorithm applied to a profile in the Arctic Ocean, where steps are indicated on separate (a) conservative temperature and (b) absolute salinity profiles. Each profile is shifted for clarity. Similar to Figs. 4–6, an interface is not considered by the detection algorithm when the interface characteristics did not meet the requirements of a previous step. Original profile is taken from Ice-Tethered Profiler ITP64 at 137.8^∘ W and 75.2^∘ N on 29 January 2013. The details of the data preparation and the algorithm steps are discussed in Sects. 2 and 3, respectively.

Figure A2 As Fig. A1, but for a profile in the Mediterranean Sea. Original profile is taken from Argo float 6901769 at 8.9^∘ E and 37.9^∘ N on 31 October 2017.

Figure A3 As Fig. A1, but for a profile in the western tropical North Atlantic. Original profile is taken from Argo float 4 901 478 at 53.3^∘ W and 11.6^∘ N on 9 August 2014.

Figure A4 Structure of the software. Each step in the software is shown by a box. Whenever a particular step is contained inside a function, the name of the function is mentioned above the step. Details of the preprocessing of the data and the detection algorithm are discussed in Sects. 2 and 3, respectively.

Table A1 Metadata of all variables that are saved in the dataset.

Table A2 Functions used in the software.