# Improving the thermocline calculation over the global ocean

Emmanuel Romero, Leonardo Tenorio-Fernandez, Esther Portela, Jorge Montes-Aréchiga, Laura Sánchez-Velasco

► To cite this version:

Emmanuel Romero, Leonardo Tenorio-Fernandez, Esther Portela, Jorge Montes-Aréchiga, Laura Sánchez-Velasco. Improving the thermocline calculation over the global ocean. Ocean Science, 2023, 19 (3), pp.887-901. 10.5194/os-19-887-2023. hal-05312354

HAL Id: hal-05312354

https://hal.univ-brest.fr/hal-05312354v1

Submitted on 13 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Improving the thermocline calculation over the global ocean

Emmanuel Romero<fnref n="1" />, Leonardo Tenorio-Fernandez<fnref n="1" />,<fnref n="2" />, Esther Portela<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, Jorge Montes-Aréchiga<fnref n="6" />, and Laura Sánchez-Velasco<fnref n="1" />

**Correspondence:** Leonardo Tenorio-Fernandez (ltenoriof@ipn.mx)

Received: 12 August 2022 – Discussion started: 29 August 2022

Revised: 28 April 2023 – Accepted: 28 April 2023 – Published: 22 June 2023

**Abstract.** According to the typical thermal structure of the ocean, the water column can be divided into three layers: the mixed layer, the thermocline and the deep layer. In this study, we provide a new methodology, based on a function adjustment to the temperature profile, to locate the minimum and maximum depths of the strongest thermocline. We first validated our methodology by comparing the mixed layer depth obtained with the method proposed here with three other methods from previous studies. Since we found a very good agreement between the four methods we used the function adjustment to compute the monthly climatologies of the maximum thermocline depth and the thermocline thickness and strength in the global ocean. We also provide an assessment of the regions of the ocean where our adjustment is valid, i.e., where the thermal structure of the ocean follows the three-layer structure. However, there are ocean regions where the water column cannot be separated into three layers due to the dynamic processes that alter it. This assessment highlights the limitations of the existing methods to accurately determine the mixed layer depth and the thermocline depth in oceanic regions that are particularly turbulent such as the Southern Ocean and the northern North Atlantic, among others. The method proposed here has shown to be robust and easy to apply.

## 1 Introduction

In most of the ocean, a typical vertical temperature profile shows maximum temperature at the surface, due to solar radiation, and can be divided into three main layers according to the thermal structure of the ocean: (i) the mixed layer, where the turbulence generated by atmospheric processes homogenizes the temperature and distributes heat throughout the layer; (ii) the thermocline, the layer with the strongest stratification, which separates the upper mixed layer from the deep layer of the ocean; and (iii) the deep layer, where the temperature is practically invariant over time and relatively constant from the lower thermocline to the seafloor. This three-layer structure is similar to the fundamental vertical density structure of the world ocean (Sallée et al., 2021), where the central layer is the pycnocline.

In most of the ocean, the temperature exerts the main control on the density of the water column. Exceptions to this are mainly found in polar regions, where temperature is very low, and the seawater density is mostly determined by salinity (de Boyer Montégut et al., 2004; Yamaguchi and Suga, 2019; Sallée et al., 2021), and in the so-called barrier layer (BL) regions (Lorbacher et al., 2006). The latter are regions where the mixed layer depth (MLD) is determined by a halocline. In these regions, the MLD based on temperature (the isothermal layer) is deeper than MLD based on density profiles (the

<footnote n="1">¹Departamento de Oceanología, Instituto Politécnico Nacional-Centro Interdisciplinario de Ciencias Marinas (IPN-CICIMAR), Av. IPN s/n, La Paz, B.C.S., 23096, Mexico</footnote>

<footnote n="2">²CONAHCyT, Consejo Nacional de Humanidades Ciencia y Tecnología, Av. Insurgentes Sur 1582, Col. Crédito Constructor, Alcaldía Benito Juárez, C.P. 03940, Mexico City, Mexico</footnote>

<footnote n="3">³School of Environmental Sciences, University of East Anglia, Norwich, UK</footnote>

<footnote n="4">⁴Univ. Brest, Laboratoire d'Océanographie Physique et Spatiale, CNRS, IRD, Ifremer, Plouzané, France</footnote>

<footnote n="5">⁵Institute for Marine and Antarctic Studies, University of Tasmania, Hobart 7001, Australia</footnote>

<footnote n="6">⁶Departamento de Física, Universidad de Guadalajara, Gral. Marcelino García Barragán 1421, Olímpica, 44430 Guadalajara, Jal, Mexico</footnote>
isopycnal layer). In the opposite case, when the isothermal layer is shallower than the MLD, the vertical compensation between salinity and temperature causes compensated layers (CLs) located below the MLD (de Boyer Montégut et al., 2004).

The thermocline is the ocean layer where the temperature changes sharply with depth as compared to the upper and lower layers (Fiedler, 2010). Consequently, the thermocline depth is often defined as the depth of the maximum vertical temperature gradient. The characteristics and depth of this layer vary spatially. At low latitudes, due to relatively higher temperatures in the upper water column, the stratification is high, and permanent thermoclines are relatively strong and thin. In contrast, at high latitudes, where there is generally little difference between the surface and deep layer temperature, the thermocline is generally weaker and deeper (Webb, 2021). The strengthening of the upper thermocline at mid latitudes during summer, when net heat flux at the surface is positive, and wind mixing is low, is known as the seasonal thermocline (Sprintall and Cronin, 2001). Due to cooling, wind-driven mixing and a well-stratified thermocline, the mixed layer is deeper in winter (Sprintall and Cronin, 2001). In tropical and polar regions the seasonal changes are weak.

Other classifications of the thermocline have been proposed from a machine learning approach. For instance Jiang et al. (2017) classify thermoclines on the basis of their form as positive, inverse or mixed thermoclines as well as multi-thermoclines. The forms that originate from this classification could be related to the temperature inversions that occur at the base of the BL and in the polar regions (de Boyer Montégut et al., 2004; Dong et al., 2008) and by the double-diffusive staircase stratification events (Timmermans et al., 2008; Toole et al., 2011).

The MLD (which is also the top of the thermocline) as well as the maximum thermocline depth (MTD) and thermocline strength plays a key role in determining the vertical distribution of many physical and ecological parameters (Fiedler, 2010). The thermocline is a physical gradient that plays a key role in climate variability and ocean-atmosphere interactions (Chu and Fan, 2019). The thermocline strength affects buoyancy, heat budgets, circulation and exchange of properties. Its depth is associated with the habitat and abundance of zooplankton organisms (Southward and Barrett, 1983; Ruvalcaba-Aroche et al., 2022) and is also an ecological boundary for pelagic organisms. Strong temperature changes can set habitat distributions, and the thermocline often corresponds to gradients in nutrients (nutricline), oxygen (oxycline) or other limiting factors. The thermocline thickness also affects the intensity of the primary production. Particularly at the poles, where the thermocline is weaker, enhanced mixing distributes nutrients throughout the water column. In contrast, in equatorial and tropical regions, the strong thermocline prevents nutrient-rich water from the deep layer from reaching the surface (Webb, 2021). Observations of

tracer concentrations and model simulations suggest a connection between the equatorial thermocline and mid-latitude ventilation regions (Harper, 2000).

The strength of the ocean stratification has strong implications for the ventilation of the interior ocean and the injection of traces as oxygen and carbon (Sallée et al., 2012; Portela et al., 2020). Therefore, the knowledge and monitoring of the depth and strength of the thermocline are particularly relevant in the context of ocean warming and its effects on the pelagic ecosystem. Recent global studies have found an overall increase in the global ocean stratification (Li et al., 2020) as well as important regional variability in the global pycnocline trend over the past decades (Sallée et al., 2021). For instance, the thickness of the equatorial and tropical thermocline is enhanced under ocean warming because the surface layer warms more and faster than the lower layers (Yang and Wang, 2009).

Previous regional studies have identified a shallowing and strengthening thermocline in the western Pacific Ocean (Vecchi and Soden, 2007) and in the equatorial Pacific (Zelle et al., 2004). Additionally, modeling studies have suggested that important changes in the Pacific Ocean, such as rising sea levels and temperatures, affect the structure of the thermocline from the subtropics to the tropics (Landerer et al., 2007; Overland and Wang, 2007). These changes in thermocline in the western tropical Pacific are considered to influence the properties of El Niño–Southern Oscillation (ENSO; Luo et al., 2009), which has strong climatic and socio-economical consequences at the basin scale because the organisms change their distribution and abundance.

Different methodologies have been proposed to locate the maximum depth of the mixed layer (e.g., de Boyer Montégut et al., 2004, 2007; Lorbacher et al., 2006; Holte and Talley, 2009) and the strength and trends of the ocean stratification (Yamaguchi and Suga, 2019; Li et al., 2020; Sallée et al., 2021). However, little effort has been dedicated to identifying and mapping the MTD on a global scale. This depth can be delimited empirically by locating the rapid temperature change in the profile; other studies are based on, for instance, calculating the thermocline gradient using the exponential leap-forward gradient method (Chu and Fan, 2017) or the maximum curvature point method (Jiang et al., 2016) or using a matrix to calculate the temperature gradient strength of each point and filtering those points that meet the thermocline standard ( > 0.2 °C m⁻¹) (Jiang et al., 2017). Additionally, Fiedler (2010) has compared different methods to estimate the mixed layer depth, thermocline depth and thermocline strength. In his study, the method that gave the best results was the variable representative isotherm (VRI). This method locates the thermocline from the base of the mixed layer to the depth at which temperature has dropped halfway toward the deep-water temperature at 400 m (Fiedler, 2010). Despite few studies applying the above methods to particular regions of the ocean, to the best of our knowledge, there are no studies addressing the MTD on a global scale. The
methodologies mentioned above to locate the MLD use data from the profilers of the Argo program. Argo is an international program that measures the ocean water column using a fleet of autonomous profilers, which move along ocean currents and measure the water column by making profiles from a depth of 2 km to the surface (Argo, 2022a).

This paper proposes a simple and efficient methodology to locate the minimum and maximum depth of the thermocline and its thickness, making an adjustment of the sigmoid function to the temperature profiles. Locating these depths helps to conduct research on thermocline-related ocean warming, and through the proposed methodology, it will be possible to conduct local and global studies on changes in ocean thermal structure through time and space.

In this study, we first describe the proposed method to calculate the MLD and MTD; then we compare the results with other methods found in the literature; and finally we calculate the thickness and strength of the thermocline to obtain the climatologies of the mixed layer depth, the maximum thermocline depth, the thermocline thickness and the thermocline strength index.

## 2 Data collection

For all the diagnostics carried out in this study we used the Argo dataset. We downloaded the snapshot of January 2022 (Argo, 2022b), and we used the profiles already evaluated by the delayed-mode quality control (DMQC) from January 1998 to December 2021 (more than 2 million) that have been classified as “good” or “probably good” data.

We selected pressure, temperature and salinity profiles from the core Argo floats, which typically sample down to 2000 m. We then transformed the in situ temperature and practical salinity into conservative temperature (Θ) and absolute salinity (Sₐ) using the definition of the Thermodynamic Equation of Seawater – 2010 (IOC et al., 2010).

## 3 Methodology

Previous studies (e.g., de Boyer Montégut et al., 2004; Holte and Talley, 2009) have proposed different methodologies to calculate the MLD on a global scale. Despite the existence of methods to calculate MTD (such as those compared in the revision study of Fiedler, 2010), these have been evaluated with a limited amount of data and in relatively small tropical and subtropical areas, therefore excluding profiles of high latitudes. Here we propose to use a new method based on the sigmoid function adjustment in the temperature profile to localize the MTD. Our method takes advantage of the characteristics that this function shares with the typical temperature profiles in most of the ocean: a straight line that represents the homogeneity of the MLD, a diagonal that represents the rapid increase or decrease in temperature with depth (changing the sign of the function) in the thermocline, and a straight line that represents the little variability in the temperature of the deep ocean (Fig. 1).

To locate the MTD, we computed the vertical maximum of the contribution of temperature to the Brunt–Vaisala frequency squared (i.e., maximum of N²ₜ) to locate the most stratified point from the temperature profile. We assume that this point is within the thermocline, as the most stratified point of the water column given by N² is inside the pycnocline (IOC et al., 2010). N²ₜ is given by Eq. (1), where g is the gravitational acceleration, ρ is the density, α is the coefficient of thermal expansion, and ΔΘ is the difference between conservative temperatures of vertically adjacent seawater parcels separated in pressure by ΔP (IOC et al., 2010).

$$N_T^2 = g^2 \rho \frac{-\alpha^{\Theta} \Delta \Theta}{\Delta P} \qquad (1)$$

Schematically, most of the temperature profiles at all latitudes have a shape similar to the sigmoid function (S shape); for this study we used the logistic function shown in Eq. (2), where a is the steepness of the curve, and b is the value of the midpoint of the function, also known as the inflection point.

$$f(x) = \frac{1}{1 + e^{-a(x-b)}} \qquad (2)$$

To perform the function adjustment, we first locate the greatest absolute value of N²ₜ, and we take the temperature profile from the surface to its depth multiplied by 2. In this way, we reduce the data from the deep layer but make sure not to exclude the isothermal layer or the thermocline. The sigmoid function presents central symmetry with respect to its inflection point; from this point, in both directions, the sigmoid presents a diagonal line, a curve and a straight line. Given these characteristics, by fitting the sigmoid function, we seek to fully represent the mixed layer with a straight line, locate the inflection point in the center of the thermocline and consequently represent the thermocline with the diagonal line.

First, we evaluate the direction of the vertical temperature change. To do this, we compare the temperature value near the surface against the deeper one; if the value closest to the surface is greater, the profile decreases with depth; otherwise it increases. If the temperature decreases with depth, the sigmoid function is inverted by multiplying it by -1; then we normalize the temperature data between 0 and 1.

Next, nonlinear least squares is used to fit the function to obtain the optimal values of the parameters a and b. Once these parameters are obtained, it is possible to approximate the temperature values at any depth above the sigmoid. Despite the central symmetry that the sigmoid function presents, the nonlinear fit of least squares allows the fit to place one straight line shorter than the other one (without losing its shape), thus losing the symmetry and placing the inflection point in the center of the thermocline, regardless of whether or not it coincides with the greatest value of N²ₜ. We assess
the goodness of the fit with the coefficient of determination (R²); this coefficient provides information on how well the adjusted function approximates the real data, 1 being the best adjustment.

Once the sigmoid has been fitted to the temperature profile, we can determine the MLD and MTD by scrolling through the function. The temperature at a depth of 10 m resulting from the adjustment of the function is taken as a reference and is denormalized, that is, it is transformed again to be represented as a function of depth. The MLD is then determined as the depth where the potential temperature is 0.2 °C higher (or lower) than the reference temperature at 10 m (de Boyer Montégut et al., 2004). To locate the MTD, we used the same procedure but going upwards in the function; in this case we take the reference temperature where the deep layer should be located, and we look for the difference of 0.2 °C by decreasing the depth through the function. Because the method is based on a single nonlinear function adjustment, we can have a precision of even centimeters. The procedure explained above can be seen in Fig. 1 and can be used through the script developed by Romero et al. (2022). To visualize the profiles of Fig. 1 up to 2000 m depth, see Fig. S1 of the Supplement.

This methodology was applied to each of the DMQC Argo profiles, and consequently we provide the monthly average of R² in a 2° × 2° grid as a proxy to know the regions of the ocean where the proposed methodology is reliable.

We have validated the method by comparing our results for the MLD with other existing methods. To do so, the MLD of each profile of the dataset was calculated in four different ways: (i) with the proposed method, (ii) following the methodology of Holte and Talley (2009) and with the methodology of de Boyer Montégut et al. (2004), (iii) using the density threshold, and (iv) using the temperature threshold (hereinafter, we refer to the former three methods as HT09, B04D and B04T respectively). HT09 performs an evaluation of several criteria (calculated from temperature, salinity and density separately) to determine the MLD for each profile, while B04D uses a threshold of 0.03 kg m⁻³ compared to the reference value at 10 m depth in the density profile and B04T an absolute difference of 0.2 °C compared to the reference value at the same depth but in the temperature profile, to locate the MLD. To compare the four methods and compute the MLD, the ocean was divided into regions following the reference of the Working Group I contribution to the Sixth Assessment Report (AR6-WGI) (Iturbide et al., 2020) of the Intergovernmental Panel on Climate Change (IPCC), and the regional monthly average of the MLD was calculated. Regions with fewer than 10 averaged values were not taken into account. Finally, the Spearman correlation was calculated between the results of the four methodologies.

To carry out our computations we averaged all profiles available for each climatological month in 2° × 2° cells. The choice of the 2° × 2° cells responds to a compromise between keeping reasonable resolution and enough data in each cell

for each climatological month. With these data, we then obtained climatologies of the MLD produced by each methodology described above including the one proposed here. Once the calculation of the MLD was validated, the monthly climatologies of the MTD, the thickness and the strength of the thermocline were obtained. The thermocline strength was calculated using the thermocline strength index (TSI), defined as ΔΘΔh⁻¹ (°C m⁻¹), where ΔΘ and Δh are the differences in temperature and depth between the MLD and MTD (Yu et al., 2010).

## 4 Results

The methodology developed here was applied to each of the DMQC Argo profiles marked as “good” or “probably good” data. As a preliminary assessment of the adequacy of the sigmoid function fit, we performed a first visual scan of random temperature profiles at different ocean latitudes. Figure 1 shows, illustratively, some examples of different temperature-profile-adjustment situations with different characteristics and geographical locations, where the MTD and the MLD computed with the method proposed here are indicated.

The temperature profiles shown in Fig. 1a–e were taken from high latitudes in the Southern Hemisphere to high latitudes in the Northern Hemisphere, while the profiles in Fig. 1f–j are located in regions where thick BLs are found (see de Boyer Montégut et al., 2007). In the profiles from Fig. 1a–e the temperature drops in the thermocline as the depth increases, while the profiles in Fig. 1f and h–j show temperature increase with depth. In both cases our methodology seems to accurately determine the MTD. Profiles in Fig. 1f and g show the greatest variability in N_T², but the quality of our adjustment differs between them. In panel f, despite the high variability in the deep layer, the methodology correctly determines the MLD and MTD. However, in panel g, high variability occurs from the end of the isothermal layer, and our methodology cannot perform the adjustment of the function correctly. In the same way, the B04T, D_σ and VRI methods failed to correctly locate the MLD and the thermocline (as shown in Fig. S1 in the Supplement). To illustrate the precision of our method and to identify regions where it should be applied with caution due to the variability in the temperature profiles, we provide a map of the monthly average of R² (Fig. 2).

In general terms, the adjustment of the sigmoid function is very good (with R² ≥ 0.9) at low latitudes and mid-latitudes. However, the cells with red and gray colors should be taken with caution. These present R² < 0.3 and < 0.7 respectively, which indicates that the adjustment of the sigmoid function was poor or not optimal. The worst adjustments correspond to the core of the Antarctic Circumpolar Current in the Southern Ocean, the North Pacific and the western North Atlantic. These are regions where the stratification of the water col-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Location of the MTD (dashed red line) and the MLD (dashed black line) in temperature profiles (blue dots): (a) 52.95° S, 90.05° W, on 23 January 2003; (b) 25.13° S, 93.47° W, on 12 January 2013; (c) 1.90° S, 126.07° W, on 25 August 2013; (d) 20.02° N, 41.14° W, on 15 December 2015; (e) 49.00° N, 174.69° W, on 13 December 2017; (f) 60.00° S, 116.86° W, on 12 August 2015; (g) 55.42° S, 162.63° W, on 12 August 2020; (h) 63.23° N, 54.20° W, on 8 February 2010; (i) 56.07° N, 174.91° W, on 20 February 2014; (j) 61.84° N, 54.27° W, on 1 February 2016. Goodness of fit is shown at the top of each profile with $R^2$.</figcaption></figure>

umn is dominated by salinity; there are temperature inversions and/or present strong currents and associated turbulent dynamics. In general terms, in the regions where the adjustment was worse, it was not as good in winter months.

The results of the preliminary evaluation showed that both the visual examination (not shown) and the values of $R^2$ indicate that our methodology correctly locates the MTD and the MLD at different latitudes. After this first step, we carried out the validation against other methodologies.

## 4.1 Atlas of the mixed layer depth

The monthly climatology of the MLD computed with the proposed method (Fig. 3) reproduces the spatial patterns and the seasonal variability in the mixed layer well, as shown in previous studies (e.g., de Boyer Montégut et al., 2004; Holte et al., 2017). It captures the regions with the deepest (northern North Atlantic and Southern Ocean) and the shallowest values (tropical and subtropical areas of both hemispheres) and their magnitudes.

The MLD shows strong seasonality as well as hemispheric asymmetry, mainly in the subtropical and subpolar regions. In the Northern Hemisphere, in summer months, the mixed layer is generally shallower than 50 m, while in late winter,
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Monthly average of $R^2$.</figcaption></figure>

it reaches climatological mean values over 1000 m in some regions of the North Atlantic basin such as the Labrador Sea and the Nordic Seas. In the Southern hemisphere, the MLD is generally deeper than in the Northern Hemisphere, and it is dominated by the signal of the Antarctic Circumpolar Current. The mixed layer in this region varies between 75–100 m depth in summer and around 500 m depth in winter, mainly in the Indian and Pacific basins. At tropical and subtropical latitudes, the MLD is generally very shallow, varying below 15 m in summer up to 150 m depth in winter.

In general terms, our climatology agrees with those of Holte et al. (2017), de Boyer Montégut et al. (2004), HT09, B04D and B04T (Figs. S2–S5), but there are some differences. The comparison between our method and HT09 (Fig. S2) shows no net underestimation or overestimation, and the relative difference between the two methods is less than |25%| over most of the ocean. The HT09 method gives slightly shallower mixed layers than our method at subpolar latitudes and deeper mixed layers in tropical and subtropical regions. On the other hand, B04D and B04T generally overestimate the MLD with respect to our method (Figs. S2 and S3). B04D has the greatest differences (above 50%) com-

pared to our method, especially in the Northern Hemisphere from July to September, while B04T presents its greatest differences at tropical and subtropical latitudes throughout the year. To make a more quantitative comparison, the MLDs computed with the four different methods were averaged within the reference regions of AR6-WGI, as these were designed for regional synthesis (Fig. 4).

The red-delimited regions are fully continental regions and were not used for our analysis. The regions with fewer than 10 averaged values were also excluded. The averages of the MLD computed with each method in each region were plotted for a representative month of each season (Fig. 5).

The MLD shows good agreement between the four methods in most regions. In general terms, the method proposed here is in better agreement with HT09 and B04T, with B04D exhibiting the greatest differences. In February and May, the B04D method seems to overestimate the MLD in the regions of northeastern North America and the Arctic Ocean (regions 03 and 48; Fig. 5a and d). These are polar regions containing semi-enclosed seas in the Northern Hemisphere. It is likely that these regions exhibit particular dynamics that complicate the detection of the MLD with a global-threshold-based
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Climatology of the MLD estimated from individual profiles.</figcaption></figure>

method. One possible explanation is that coastal regions are
generally worse sampled but also that the computation of
the MLD can be complicated by a number of coastal pro-
cesses such as river discharges or shallow bathymetry, among
others. Moreover, the traditional delta density criterion of
0.03 kg m⁻³ has been suggested to underestimate the MLD
in polar regions, as demonstrated in the study by Peralta-
Ferriz and Woodgate (2015), where it was found that a bet-
ter criterion for these regions is 0.1 kg m⁻³. Interestingly, the
agreement between the four methods is also good in the re-
gions where our adjustment was not considered to be good
(R² < 0.7; Fig. 2). In other complicated regions such as the
one around Greenland (region 01) there are some differences
between the four methods, but the one proposed here agrees
with B04T in February and HT09 in August and gives an in-
termediate MLD value with respect to two other methods (as
in May and November). The Spearman correlation between
the results was calculated (Fig. S5) and showed high corre-
lation between the results of all the methodologies, with the
results of HT09 and B04T being the most correlated (0.98),
followed by the proposed method and B04T (0.95). Finally,

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Reference regions of AR6-WGI. The regions used for comparison are those delimited in black.</figcaption></figure>

all the correlations between B04D and the rest of the method-
ologies showed values close to 0.90.

All this suggests that although our method is not perfect
in highly dynamical regions, it gives results that compare
well with other broadly used methods to detect the isother-
mal layer and the MLD, even in salinity-dominated regions.
Moreover, this highlights the possible deficiency of all existing methods in detecting the MLD in these regions.

## 4.2 Climatology of the maximum thermocline depth, thickness and strength

Once the proposed methodology is validated with the calculation of the MLD, we computed the monthly climatology of the MTD (Fig. 6).

As expected, the shape of the MTD follows that of the MLD, but the hemispheric asymmetry is not so evident. The subtropical and subpolar regions of the North Atlantic as well as the Southern Ocean exhibit the deepest thermoclines of the ocean. Similarly to the MLD, the deepest thermoclines are found in late winter: March–April in the Northern Hemisphere and September–October in the Southern Hemisphere. In the Northern Hemisphere, during the summer the MTD is generally no deeper than 100 m, while in winter, it reaches depths greater than 1000 m in the same regions where the MLD reaches its maximum values (Fig. 3), extending to the Gulf Stream in winter and spring months. In the Southern Hemisphere, the mean climatological MTD for the summer and winter months is similar to those of the opposite hemisphere. MLDs are relatively deep in the Antarctic Circumpolar Current in the Southern Ocean. During winter and spring, the MTD reaches values deeper than 500 m in its core, reaching more than 1000 m in localized areas. The deeper thermoclines in the Southern Ocean, in winter, coincide with the seasonality of the zonal band where the inertial horizontal kinetic energy in the mixed layer is larger (Flexas et al., 2019). This energy is injected by the relatively strong winds during this season, which is also related to the relatively deep MLD (as shown in Fig. 3).

During summer, the MTD rarely exceeds 50 m depth in the Northern Hemisphere and 75 m in the Southern Hemisphere. In these same regions, the climatologies of the thermocline thickness (Fig. 7) mainly presented values below 50 m.

The thermocline thickness follows a similar pattern to the MLD and the MTD. As for the MLD and MTD (Figs. 3 and 6 respectively), the climatology of the thermocline thickness (Fig. 7) presents a marked seasonality at subtropical and subpolar latitudes. As expected, the thickest thermoclines ($> 500$ m) are found in late winter and in the regions of lower stratification: March–April in the North Atlantic basin and September–October in the Southern Ocean, particularly in the region dominated by the Antarctic Circumpolar Current. Away from these regions, the seasonal variability in the thermocline thickness is low in tropical regions, where it varies between 150 and 250 m depending on the region. The thinnest thermoclines are observed in summer at subtropical/subpolar latitudes of both hemispheres (July–August in the Northern Hemisphere and January–February in the Southern Hemisphere).

Finally, the climatology of the thermocline strength (Fig. 8) maintains the seasonality at subtropical and subpolar

latitudes with a $|TSI| < 0.1$. On the other hand, the tropical eastern Pacific and tropical eastern Atlantic have TSI $> 0.1$ throughout the year, as do the North Pacific and North Atlantic, but from June–November. The Black Sea shows the highest values of TSI ($> 1$), while the lowest values ($< -0.1$) are scattered in subpolar regions and in the discharge region of the Ganges and Brahmaputra rivers (December–February).

## 5 Discussion

In this study we propose a new method to locate the strongest thermocline that lies just below the MLD of the water column. This method is based on the adjustment of a sigmoid function (in this case the logistic function) that relies on the principle that the thermal structure of the ocean consists of three main layers: the mixed layer, the thermocline and the deep layer of the ocean. Although not all temperature profiles have the S shape mentioned above throughout the ocean, the proposed method is based on the absolute maximum point of $N_T^2$ and the nonlinear least squares to fit the sigmoid function and place the diagonal line of the function in the thermocline. Since the most stratified point of the temperature profile is used to place the sigmoid, the method locates the strongest thermocline, which in most cases will coincide with the seasonal one.

The proposed method, due to the shape of a typical temperature profile in the ocean, also allows us to determine the MLD, and we were therefore able to validate it. The climatology of the MLD generated in our study (Fig. 3) is in good agreement with those provided by de Boyer Montégut et al. (2004) and Holte et al. (2017). The four methods compared in this study (based on temperature and density thresholds) reproduce the magnitude, the spatial variability and the seasonal cycle of the MLD throughout the global ocean in a similar and consistent way for most of the ocean regions. This is because over most of the ocean the stratification is dominated by the temperature (de Boyer Montégut et al., 2004). Moreover, in the few regions where all three methods disagreed, the one proposed here was mostly in line with the results of B04T and HT09, while B04D exhibited extreme MLD values. The latter is explained because the proposed method and B04T use the temperature profile and the same threshold to calculate the MLD. Moreover, HT09 uses a method that combines different thresholds, including a temperature-based one. The coincidences between the proposed method and HT09 can also be partly explained by the fact that both methodologies use function adjustments. The study of HT09 is based on the homogeneity of the MLD to perform a linear function adjustment in addition to performing another linear function adjustment in the thermocline, while the proposed methodology uses this same MLD feature to determine how deep the sigmoid function will be adjusted. The regions with the greatest differences between methodologies were found
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of methodologies to locate the MLD in (a) February, (b) May, (c) August and (d) November.</figcaption></figure>

in a small region of the Arctic Ocean and the Labrador Sea (Fig. 5a and d). There, the method of B04D strongly differed from the others; we attribute this to the fact that staircase stratification has been reported in these areas (Timmermans et al., 2008; Toole et al., 2011), defined as BL and CL regions (de Boyer Montégut et al., 2004). Indeed, our preliminary results have shown inverse thermoclines, multi-thermoclines and mixed thermoclines (Jiang et al., 2017) in these regions (not shown). In this sense, we provide some examples of how our method could accurately locate the BL in some cases (Fig. S1 in the Supplement). In these cases, the MLD calculated by the 0.2 °C threshold (B04T) coincides with that calculated with our methodology. Conversely, in other cases, where the isothermal layer was highly variable (Fig. S1g and i), our method was unable to locate the thermocline. The BL was also located in these profiles using the $D_σ$ and $D_{T-02}$ criteria defined in de Boyer Montégut et al. (2007), showing that when the thickness of the BL is of the order of hundreds of meters, the proposed method locates the inverse thermocline produced by the inversion of temperature that is found below the isothermal layer. In contrast, in thinner BLs, our method locates the thermocline below the lower limit of the BL. Sprintall and Tomczak (1992) define the BL as the distance that separates the MLD from the MTD; however, when this BL is thick, the proposed methodology locates the inverse thermoclines that are within this barrier and not those that could be below $D_{T-02}$.

All previous studies compared here have used Argo data in their climatologies. However, de Boyer Montégut et al. (2004), being an older study, had less data, while Holte et al. (2017) used real-time quality control (RTQC) data in a 1° × 1° grid (in contrast to the 2° × 2° grid used by us and de Boyer Montégut et al., 2004). This smaller-sized grid is not optimal to be used with the amount of DMQC data available in this snapshot (Argo, 2022b), as about a third of the cells would contain fewer than three values for monthly averages. Using RTQC data would increase the amount of data available for our computations. However, Argo recommends using only DMQC for scientific research, since the RTQC tests are automated and may contain bad data, as explained in the manuals (Argo Data Management Team, 2022) and even using the best quality control flag, as shown in Romero et al. (2021). Therefore, using these data could cause the erroneous computation of the MLD and MTD.

Using the density to estimate the MLD usually gives good results, since it depends on temperature and salinity. However, the density can show vertical compensation below the well-mixed layer (de Boyer Montégut et al., 2004), causing deeper MLDs calculated from density thresholds. Although the B04D and HT09 methods use density profiles for the calculation of the MLD, this does not give good results with the methodology proposed in this paper. The adjustment of the sigmoid function is based on the typical shape of the temperature profile. In order to calculate the pycnocline with a
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Climatology of the MTD estimated from individual profiles.</figcaption></figure>

similar methodology to the one proposed here, it would be
necessary to find another function that better fits the den-
sity profile since, despite also being able to be represented
in three layers, the typical density profiles present an inclina-
tion along the entire profile, which makes it difficult to fit a
conventional sigmoid function. In the polar and some subpo-
lar regions of both hemispheres, where salinity is the major
contributor to the density gradient, it dominates the stratifi-
cation. In these cases the thermocline and pycnocline may
differ significantly, and this is why MLD calculations based
on temperature profiles differ considerably from those based
on density in these regions.

The calculations of the climatologies of the MTD, the
thermocline thickness and the thermocline strength were not
compared with the calculation of any other method since no
method of calculating these parameters was found that works
on a global scale. Helber et al. (2012) mention that the transi-
tion layer thickness (TLT) used in their study may encompass
the entire thickness of the thermocline, and in fact their TLT
climatology presents some coincidences with our climatol-
ogy of the thermocline thickness (Fig. 7). The most notori-
ous coincidences are in the regions of the northeast Pacific

Ocean, the northern North Atlantic Ocean and the Antarc-
tic Circumpolar Current in the Southern Ocean, where the
thermocline thickness and the TLT show values greater than
350 m in the winter months, in addition to the marked season-
ality that is present in both tropical and subtropical regions.
The climatology of the thermocline strength (Fig. 8) was cal-
culated through the TSI; this index indicates the steepness
of the thermocline (Duka et al., 2021). The further TSI is
from 0, the less steep the slope and therefore the stronger the
thermocline. The strongest thermoclines found in the Black
Sea are associated with thin thermoclines (15–20 m) between
warm surface waters and cold intermediate waters (20–8 °C)
(Akpinar et al., 2017), which produce small slopes. Nega-
tive values of TSI are caused by inverse thermoclines. These
were found mainly in subpolar regions and in the Ganges and
Brahmaputra river discharge; these regions present TSI close
to 0, which means steep slopes and therefore weak thermo-
clines. The formation of intermediate-strength thermoclines
(i.e., 0.2 < TSI < 0.8) in the North Pacific and North At-
lantic coincides mainly with the months (July–September)
when there are no BLs in these regions (de Boyer Montégut
et al., 2007). In contrast, from January–March, when the BLs
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Climatology of the thermocline thickness estimated from individual profiles.</figcaption></figure>

are thicker (de Boyer Montégut et al., 2007), weak thermoclines and regions with inverse thermoclines are shown.

To compare the location of the MTD, the VRI method was applied to the profiles shown in Fig. S1 and only gave good results with those located at tropical latitudes. As shown in Fig. S1, far from the tropics, the calculation of the MLD and MTD with the VRI method does not give good results, in addition to not considering the inverse thermoclines (Fig. S1f-j). For these reasons, the performance in fitting the sigmoid function was used to validate the method.

In different scientific areas, $R^2$ is used as a goodness-of-fit measure for sigmoid functions (e.g., Cao et al., 2019; Bhatgal et al., 2014; Ritz and Spiess, 2008; Liu and Saint, 2002; Van der Graaf and Schoemaker, 1999). Through this measure, our method showed generally good performance in the adjustment of the sigmoid function to the temperature profiles, with the exception of a few regions: the North Pacific Ocean, the northern North Atlantic Ocean, the Arctic Ocean and the core of the Antarctic Circumpolar Current in the Southern Ocean. While the interpretation of the MLD and the MTD in these regions has to be made with caution, it is noteworthy that the values of $R^2$ are not a direct indicator

of the precision of the method to calculate the MLD and the MTD. Rather, this index shows the goodness of the sigmoid function fit to the temperature profile. Precisely, in the most problematic regions mentioned above, $R^2$ is lower than 0.7 in some months (Fig. 2). However, the three methods used for the MLD calculation give very similar values (Fig. 5) even in these regions. This suggests that it is not a particular shortcoming of our adjustment. Previous studies such as those by Peralta-Ferriz and Woodgate (2015) and Pellicherro et al. (2017) have shown the variability in the calculation of the MLD in these regions depending on the methods and thresholds used. This has evidenced that these are complex, highly dynamical regions (i.e., turbulent regions with important eddy activity), where the estimation of the MLD in a reliable way is a complicated task. In this sense, it has been suggested that in the Southern Ocean the MLD calculation is less accurate than in regions at lower latitudes, where the water column is strictly temperature-stratified (Dong et al., 2008). The low values of $R^2$ shown in Fig. 2 are due to the abrupt changes in temperature in the profiles measured in these regions. These abrupt changes might be related to well-known processes taking place in certain regions of the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Climatology of the thermocline strength estimated from individual profiles.</figcaption></figure>

ocean such as (i) double-diffusive staircase stratification in the Arctic (Timmermans et al., 2008; Toole et al., 2011) and (ii) the temperature inversions due to the influence of salinity that present different vertical structures in the Southern Ocean (Dong et al., 2008).

The results of our adjustment also evidence the regions of the ocean where the water column exhibits a typical vertical thermal structure in three layers and the regions where, due to their dynamics, the structure of the water column cannot be divided into these three layers. The efficiency provided by the proposed method for the calculation of the MLD and the MTD allows local to global studies to be performed. For example, in the context of ocean warming, the differences in these layers could be compared over different timescales to analyze the changes in the water column and detect areas of the ocean where the thermocline has changed its depth, thickness or strength over time and could therefore be considered to be a parameter of the potential effects on the pelagic ecosystem and socio-economic repercussions.

## 6 Conclusions

In this study, we present a methodology to locate the minimum and maximum depths of the strongest thermocline, its thickness, and its strength by adjusting the sigmoid function to the temperature profiles in the global ocean. This methodology can be applied in those areas of the ocean where the water column can be divided into three layers according to its thermal structure. Our methodology gave good results in its validation against three other broadly used methodologies in the global ocean. The MLD computed with the four methods showed a high correlation, even in regions where the coefficient of determination suggested a poor adjustment. This suggests that it is not a particular shortcoming of our method, but rather a general difficulty in determining the limits of the three typical oceanic thermal layers in highly turbulent regions.
*Code availability.* The methodology presented in this study was developed in Python 3.7 and is licensed under a Creative Commons Attribution 4.0 International License. The source code is available at https://doi.org/10.5281/zenodo.6985561 (Romero et al., 2022). The latest package version is v1.1.0.

*Data availability.* These data were collected and made freely available by the international Argo program and the national programs that contribute to it (https://argo.ucsd.edu, last access: January 2022; https://www.ocean-ops.org, last access: January 2022). The Argo program is part of the Global Ocean Observing System. The data used are from the snapshot of January 2022 (https://doi.org/10.17882/42182; Argo, 2022b).

*Supplement.* The supplement related to this article is available online at: https://doi.org/10.5194/os-19-887-2023-supplement.

*Author contributions.* ER developed the methodology described in this work and carried out the analyses. All co-authors contributed to the conceptualization and design of the study, the interpretation of the results, and the preparation of the article.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

*Acknowledgements.* We are grateful to CONACYT for granting scholarship no. 780669 to Emmanuel Romero. We appreciate that these data were collected and made freely available by the international Argo program. We also thank the Centro Interdisciplinario de Ciencias Marinas (CICIMAR) for their institutional support and acknowledge the critical comments from the reviewers.

*Financial support.* This research has been supported by CONACYT (scholarship no. 780669).

*Review statement.* This paper was edited by A.J. George Nurser and reviewed by two anonymous referees.

## References

Akpinar, A., Fach, B. A., and Oguz, T.: Observing the subsurface thermal signature of the Black Sea cold intermediate layer with Argo profiling floats, Deep Sea Res. Pt. I, 124, 140–152, https://doi.org/10.1016/j.dsr.2017.04.002, 2017.

Argo: Argo, https://argo.ucsd.edu/ (last access: January 2022), 2022a.

Argo: Argo float data and metadata from Global Data Assembly Centre (Argo GDAC) [data set], https://doi.org/10.17882/42182, 2022b.

Argo Data Management Team: Argo user's manual, https://doi.org/10.13155/29825, 2022.

Bhogal, A. A., Siero, J. C., Fisher, J. A., Froeling, M., Luijten, P., Philippens, M., and Hoogduin, H.: Investigating the nonlinearity of the BOLD cerebrovascular reactivity response to targeted hypo/hypercapnia at 7T, NeuroImage, 98, 296–305, https://doi.org/10.1016/j.neuroimage.2014.05.006, 2014.

Cao, L., Shi, P.-J., Li, L., and Chen, G.: A New Flexible Sigmoidal Growth Model, Symmetry, 11, 204, https://doi.org/10.3390/sym11020204, 2019.

Chu, P. C. and Fan, C.: Exponential leap-forward gradient scheme for determining the isothermal layer depth from profile data, J. Oceanogr., 73, 503–526, https://doi.org/10.1007/s10872-017-0418-0, 2017.

Chu, P. C. and Fan, C.: Global ocean synoptic thermocline gradient, isothermal-layer depth, and other upper ocean parameters, Sci. Data, 6, 119, https://doi.org/10.1038/s41597-019-0125-3, 2019.

de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., and Iudicone, D.: Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology, J. Geophys. Res.-Oceans, 109, C12003, https://doi.org/10.1029/2004JC002378, 2004.

de Boyer Montégut, C., Mignot, J., Lazar, A., and Cravatte, S.: Control of salinity on the mixed layer depth in the world ocean: 1. General description, J. Geophys. Res.-Oceans, 112, C06011, https://doi.org/10.1029/2006JC003953, 2007.

Dong, S., Sprintall, J., Gille, S. T., and Talley, L.: Southern Ocean mixed-layer depth from Argo float profiles, J. Geophys. Res.-Oceans, 113, C06013, https://doi.org/10.1029/2006JC004051, 2008.

Duka, M. A., Shintani, T., and Yokoyama, K.: Thermal stratification responses of a monomictic reservoir under different seasons and operation schemes, Sci. Total Environ., 767, 144423, https://doi.org/10.1016/j.scitotenv.2020.144423, 2021.

Fiedler, P. C.: Comparison of objective descriptions of the thermocline, Limnol. Oceanogr.-Meth., 8, 313–325, https://doi.org/10.4319/lom.2010.8.313, 2010.

Flexas, M. M., Thompson, A. F., Torres, H. S., Klein, P., Farrar, J. T., Zhang, H., and Menemenlis, D.: Global Estimates of the Energy Transfer From the Wind to the Ocean, With Emphasis on Near-Inertial Oscillations, J. Geophys. Res.-Oceans, 124, 5723–5746, https://doi.org/10.1029/2018JC014453, 2019.

Harper, S.: Thermocline ventilation and pathways of tropical-subtropical water mass exchange, Tellus A, 52, 330–345, https://doi.org/10.3402/tellusa.v52i3.12269, 2000.
Helber, R. W., Kara, A. B., Richman, J. G., Carnes, M. R., Barron, C. N., Hurlburt, H. E., and Boyer, T.: Temperature versus salinity gradients below the ocean mixed layer, J. Geophys. Res.-Oceans, 117, C05006, https://doi.org/10.1029/2011JC007382, 2012.

Holte, J. and Talley, L.: A New Algorithm for Finding Mixed Layer Depths with Applications to Argo Data and Subantarctic Mode Water Formation, J. Atmos. Ocean. Technol., 26, 1920–1939, https://doi.org/10.1175/2009JTECHO543.1, 2009.

Holte, J., Talley, L. D., Gilson, J., and Roemmich, D.: An Argo mixed layer climatology and database, Geophys. Res. Lett., 44, 5618–5626, https://doi.org/10.1002/2017GL073426, 2017.

IOC, SCOR, and IAPSO: The international thermodynamic equation of seawater – 2010: Calculation and use of thermodynamic properties, p. 196, 2010.

Iturbide, M., Gutiérrez, J. M., Alves, L. M., Bedia, J., Cerezo-Mota, R., Cimadevilla, E., Cofiño, A. S., Di Luca, A., Faria, S. H., Gorodetskaya, I. V., Hauser, M., Herrera, S., Hennessy, K., Hewitt, H. T., Jones, R. G., Krakovska, S., Manzanas, R., Martínez-Castro, D., Narisma, G. T., Nurhati, I. S., Pinto, I., Seneviratne, S. I., van den Hurk, B., and Vera, C. S.: An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets, Earth Syst. Sci. Data, 12, 2959–2970, https://doi.org/10.5194/essd-12-2959-2020, 2020.

Jiang, B., Wu, X., and Ding, J.: Comparison of the calculation methods of the thermocline depth of the South China Sea, Mar. Sci. Bull, 35, 64–73, 2016.

Jiang, Y., Gou, Y., Zhang, T., Wang, K., and Hu, C.: A Machine Learning Approach to Argo Data Analysis in a Thermocline, Sensors, 17, 2225, https://doi.org/10.3390/s17102225, 2017.

Landerer, F. W., Jungclaus, J. H., and Marotzke, J.: Regional Dynamic and Steric Sea Level Change in Response to the IPCC-A1B Scenario, J. Phys. Oceanogr., 37, 296–312, https://doi.org/10.1175/JPO3013.1, 2007.

Li, G., Cheng, L., Zhu, J., Trenberth, K. E., Mann, M. E., and Abraham, J. P.: Increasing ocean stratification over the past half-century, Nat. Clim. Change, 10, 1116–1123, https://doi.org/10.1038/s41558-020-00918-2, 2020.

Liu, W. and Saint, D. A.: Validation of a quantitative method for real time PCR kinetics, Biochem. Biophys. Res. Commun., 294, 347–353, https://doi.org/10.1016/S0006-291X(02)00478-3, 2002.

Lorbacher, K., Dommenget, D., Niiler, P. P., and Köhl, A.: Ocean mixed layer depth: A subsurface proxy of ocean-atmosphere variability, J. Geophys. Res.-Oceans, 111, C07010, https://doi.org/10.1029/2003JC002157, 2006.

Luo, Y., Rothstein, L. M., and Zhang, R.-H.: Response of Pacific subtropical-tropical thermocline water pathways and transports to global warming, Geophys. Res. Lett., 36, C00A02, https://doi.org/10.1029/2008GL036705, 2009.

Overland, J. E. and Wang, M.: Future climate of the north Pacific Ocean, EOS T. Am. Geophys. Un., 88, 178–182, https://doi.org/10.1029/2007EO160003, 2007.

Pellicherro, V., Sallée, J.-B., Schmidtko, S., Roquet, F., and Char-rassin, J.-B.: The ocean mixed layer under Southern Ocean sea-ice: Seasonal cycle and forcing, J. Geophys. Res.-Oceans, 122, 1608–1633, https://doi.org/10.1002/2016JC011970, 2017.

Peralta-Ferriz, C. and Woodgate, R. A.: Seasonal and interannual variability of pan-Arctic surface mixed layer properties from 1979 to 2012 from hydrographic data,

and the dominance of stratification for multiyear mixed layer depth shoaling, Prog. Oceanogr., 134, 19–53, https://doi.org/10.1016/j.pocean.2014.12.005, 2015.

Portela, E., Kolodziejczyk, N., Vic, C., and Thierry, V.: Physical Mechanisms Driving Oxygen Subduction in the Global Ocean, Geophys. Res. Lett., 47, e2020GL089040, https://doi.org/10.1029/2020GL089040, 2020.

Ritz, C. and Spiess, A.-N.: qpcR: an R package for sigmoidal model selection in quantitative real-time polymerase chain reaction analysis, Bioinformatics, 24, 1549–1551, https://doi.org/10.1093/bioinformatics/btn227, 2008.

Romero, E., Tenorio-Fernandez, L., Castro, I., and Castro, M.: Filtering method based on cluster analysis to avoid salinity drifts and recover Argo data in less time, Ocean Sci., 17, 1273–1284, https://doi.org/10.5194/os-17-1273-2021, 2021.

Romero, E., Tenorio-Fernandez, L., Portela, E., Montes-Aréchiga, J., and Sánchez-Velasco, L.: romeroqe/mld-mtd: MLD & MTD, [code], https://doi.org/10.5281/zenodo.6985561, 2022.

Ruvalcaba-Aroche, E. D., Sánchez-Velasco, L., Beier, E., Barton, E. D., Godínez, V. M., Gómez-Gutiérrez, J., and Martínez-Rincón, R. O.: Ommastrephid squid paralarvae potential nursery habitat in the tropical-subtropical convergence off Mexico, Prog. Oceanogr., 202, 102762, https://doi.org/10.1016/j.pocean.2022.102762, 2022.

Sallée, J.-B., Matear, R. J., Rintoul, S. R., and Lenton, A.: Localized subduction of anthropogenic carbon dioxide in the Southern Hemisphere oceans, Nat. Geosci., 5, 579–584, https://doi.org/10.1038/ngeo1523, 2012.

Sallée, J.-B., Pellicherro, V., Akhoudas, C., Pauthenet, E., Vignes, L., Schmidtko, S., Garabato, A. N., Sutherland, P., and Kuusela, M.: Summertime increases in upper-ocean stratification and mixed-layer depth, Nature, 591, 592–598, https://doi.org/10.1038/s41586-021-03303-x, 2021.

Southward, A. J. and Barrett, R. L.: Observations on the vertical distribution of zooplankton, including post-larval teleosts, off Plymouth in the presence of a thermocline and a chlorophyll-dense layer, J. Plankton Res., 5, 599–618, https://doi.org/10.1093/plankt/5.4.599, 1983.

Sprintall, J. and Cronin, M. F.: Upper Ocean Vertical Structure, https://doi.org/10.1006/rwos.2001.0149, 2001.

Sprintall, J. and Tomczak, M.: Evidence of the barrier layer in the surface layer of the tropics, J. Geophys. Res.-Oceans, 97, 7305–7316, https://doi.org/10.1029/92JC00407, 1992.

Timmermans, M.-L., Toole, J., Krishfield, R., and Winsor, P.: Ice-Tethered Profiler observations of the double-diffusive staircase in the Canada Basin thermocline, J. Geophys. Res.-Oceans, 113, L04601, https://doi.org/10.1029/2008JC004829, 2008.

Toole, J. M., Krishfield, R. A., Timmermans, M.-L., and Proshutinsky, A.: The ice-tethered profiler: Argo of the Arctic, Oceanography, 24, 126–135, http://www.jstor.org/stable/24861307, 2011.

Van der Graaf, P. and Schoemaker, R.: Analysis of asymmetry of agonist concentration-effect curves, Journal of Pharmacological and Toxicological Methods, 41, 107–115, https://doi.org/10.1016/S1056-8719(99)00026-X, 1999.

Vecchi, G. and Soden, B.: Global Warming and the Weakening of the Tropical Circulation, J. Climate, 20, 4316–4340, https://doi.org/10.1175/JCLI4258.1, 2007.
Webb, P.: Introduction to Oceanography, Roger Williams University, http://rwu.pressbooks.pub/webboceanography (last access: January 2022), 2021.

Yamaguchi, R. and Suga, T.: Trend and Variability in Global Upper-Ocean Stratification Since the 1960s, J. Geophys. Res.-Oceans, 124, 8933–8948, https://doi.org/10.1029/2019JC015439, 2019.

Yang, H. and Wang, F.: Revisiting the Thermocline Depth in the Equatorial Pacific, J. Climate, 22, 3856–3863, https://doi.org/10.1175/2009JCLI2836.1, 2009.

Yu, H., Tsuno, H., Hidaka, T., and Jiao, C.: Chemical and thermal stratification in lakes, Limnology, 11, 251–257, https://doi.org/10.1007/s10201-010-0310-8, 2010.

Zelle, H., Appeldoorn, G., Burgers, G., and van Oldenborgh, G. J.: The Relationship between Sea Surface Temperature and Thermocline Depth in the Eastern Equatorial Pacific, J. Phys. Ocean., 34, 643–655, https://doi.org/10.1175/2523.1, 2004.