Background & Summary

Ocean has upmost layer with near-zero vertical gradient such as isothermal layer (ITL) for temperature (T), or mixed layer (ML) for density (ρ) and underneath layer with strong vertical gradient such as thermocline (for T), or pycnocline (for ρ)1. Temperature and salinity (S) are observed in oceanography. Usually, the Thermodynamic Equation of Seawater-2010 (https://www.teos-10.org/) is used to compute ρ from (T, S) data. Thermocline or pycnocline with strong vertical gradient is the transition layer between the vertically quasi-uniform layer from the surface (ITL or ML), and the deep-water layer. The vertically quasi-uniform ITL and ML are caused by intense turbulent mixing near the ocean surface. Such mixing is driven by shear due to surface wind stress and by convection due to heat loss from ocean to atmosphere. The thermocline (pycnocline) resists the turbulent mixing from the ITL (or ML) and limits the ITL (or ML) deepening and heat (or water mass) exchange between the ITL (or ML) and deeper layer due to its strong vertical gradient. The thermocline gradient GT (or pycnocline gradient GD) directly affects such exchange2,3,4.

The ITL and ML provide dynamic-thermodynamic links and mediates the exchange of momentum, heat, and moisture between the atmosphere and the oceans; and hence plays a key role to affect weather and climate. Such exchange depends on an important parameter, i.e., the ocean isothermal layer depth (ILD) hT [or mixed layer depth (MLD) hD)], which determines the heat content (or freshwater content) and mechanical inertia of the layer. Temporal variability of hT (or hD) is caused by many processes occurring in the isothermal (mixed) layer such as surface forcing, lateral advection, internal waves, etc., ranging from diurnal, seasonal, to interannual variability2,3,4. Spatial variability of hT (or hD) is evident from less than 20 m in summer to more than 500 m in winter in subpolar latitudes5. Therefore, determination of the four parameters (hT, hD, GT, GD) from (T, ρ) profiles becomes important. Note that the two depths (hT, hD) are not necessary the same with the occurrence of barrier layer if (hT > hD)6,7,8,9,10,11,12,13,14,15 (Fig. 1a), and compensated layer if (hT < hD)15.

Fig. 1
figure 1

Schematic illustration of (T, ρ) profiles for occurrence of (a) barrier layer when ILD > MLD, and (b) compensated layer when ILD < MLD.

Full size image

Let a T-profile (or ρ-profile) starting from the ocean surface down to depth zk be represented by T(zk) [or ρ(zk)]. k = 1, 2, …, K, with z1 the surface and zK the bottom of the profile. The corresponding vertical gradient is represented by GT(zk) [or GD(zk)]. The vertical gradient for upper ocean ITL-thermocline (or ML-pycnocline) is depicted by

$$\left\{\begin{array}{l}{G}_{T}({z}_{k})\approx 0\;{\rm{for}}\;{z}_{k} > -\,{h}_{T}\;{\rm{(in}}\;{\rm{ITL),}}\quad \quad \quad {G}_{D}({z}_{k})\approx 0\;{\rm{for}}\;{z}_{k} > -\,{h}_{D}\;{\rm{(in}}\;{\rm{ML)}}\\ {G}_{T}({z}_{k})\;{\rm{or}}\;{G}_{D}({z}_{k})\;{\rm{evident}}\;{\rm{for}}\;{z}_{k} < -\,{h}_{T}\;{\rm{(in}}\;{\rm{thermocline)}}\;{\rm{or}}\;{z}_{k} < -\,{h}_{D}\;{\rm{(in}}\;{\rm{pycnocline)}}\end{array}\right.$$
(1)

Two types of methodology, single near-zero gradient and double gradients, are available to determine hT and hD from vertical T and ρ profiles. The single near-zero gradient method requires either the deviation of T (or ρ) from its value near the surface (i.e., reference level) to be smaller than a certain fixed value, such as 0.8 °C9, 0.2 °C15,16 to 0.1 °C17 for T, (0.03 kg/m3)15, (0.125 kg/m3)18, to (0.05 kg/m3)19 for ρ, or the near-zero vertical gradient to be smaller than a certain fixed value, such as (0.025 °C/m)6,10, (0.02 °C/m)20, to (0.015 °C/m)21 for GT(zk). To eliminate or reduce such uncertainty, the split-and-merge (SM)22, maximum curvature23,24,25, optimal linear fitting26, and maximum angle27 methods have been developed. The double gradient method is based on the transition of a near-zero gradient in the ITL (or ML) to an evident gradient in the thermocline (or pycnocline). The double gradient method with the exponential leap-forward gradient (ELG)28,29 was used in this study.

After (hT, hD) are determined from individual (T, ρ) profiles, the barrier (or compensated) layer depth, ocean heat content (OHC) for ITL (HITL)30, freshwater content (FWC) for ML (FML) were calculated. Note that HITL (or FML) is the vertical integration of temperature (salinity) profile from the surface (z = 0) down to the base of the ITL (z = −hT) [or ML (z = −hD)] rather than to fixed depths such as (H700, F700) for the upper 700 m. Obviously, HITL and FML are new and different from traditionally defined OHC (or FWC) in the oceanographic community with fixed depth intervals from the ocean surface such as 0–150 m (in the Indian Ocean)31, 0–300 m32, 0–400 m33, 0–700 m34, 0–750 m35, 0–2000 m36, deep layer OHC such as below 2000 m37, and the full layer OHC38. Interested readers are referred to two excellent review papers39,40.

From vertical gradient data below ITL (or ML) GT(zk), zk < −hT [GD(zk), zk < −hD], other thermohaline parameters can be identified from each T-profile (ρ-profile), such as maximum temperature (density) gradient, thermocline (pycnocline) depth, temperature (density) at thermocline (pycnocline) depth. The quality of the (hT, hD) data are estimated by the identification index29 (i-index, see Error Estimation Section). Altogether, this dataset, containing 17 thermohaline parameters with i-index, is located at the NCEI website (Data Citation 1) for public use.

Global climatological (annual mean and monthly mean) data of ocean thermohaline parameters can be established through two approaches: (1) analysing climatological (T, S) profiles such as earlier version of WOA18 to obtain climatological (hT, hD) data5,9,13,41, and (2) analysing observational profiles to get synoptic thermohaline parameters29,30 and then using optimal interpolation42, Kalmen filter43, or optimal spectral decomposition44 to produce gridded climatological thermohaline parameters. At present, we cannot estimate how big the difference is if using these two approaches. We can estimate only after the two approaches have been used. In this study, we take the first approach to derive climatological thermohaline parameter data from the NOAA/NCEI World Ocean Atlas 2018 (WOA18)45 annual and monthly mean temperature and salinity (S) profiles with regular 102 vertical levels (Table 1).

Table 1 Standard vertical depths of WOA18 data.
Full size table

Furthermore, on the base of the single gradient method (i.e., near-zero gradient in the mixed layer), serval global climatological datasets of ocean (hT, hD) were produced from the earlier version of WOA18 such as the NOAA/NCEI Mixed Layer Depth Data5, the Naval Research Laboratory Mixed Layer Depth (NMLD) Climatologies41, and comprehensive Mixed Layer Data for the Indian Ocean13. All these datasets don’t include any parameters below the ITL/ML. On the base of the double-gradient method, this dataset45 contains more thermohaline parameters from WOA18 such as thermocline gradient (GT), pycnocline gradient (GD), ITL heat content (HITL), mixed layer fresh-water content (FML), maximum thermocline gradient (GTmax), thermocline depth (hth), temperature at thermocline depth (Tth), maximum density gradient GDmax), pycnocline depth (hpyc), density at pycnocline depth (ρpyc), and (hT – hD) (barrier layer if positive or compensated layer if negative). Among them, ITL heat content (HITL), and mixed layer fresh-water content (FML) are different from the commonly used fixed-depth heat and freshwater contents.

Methods

These methods can be found in our related work28,29.

Main part of pycnocline (or Thermocline)

The ρ-profile [ρ(zk)] is taken for illustration. Let the depths corresponding to ρmin and ρmax be z1 and zK. The vertical density difference, Δρ = ρmax − ρmin, represents the total variability. Theoretically, the variability is 0 in ML and large in pycnocline beneath the ML. It is reasonable to identify the main part of the pycnocline between the two depths: z(0.1) and z(0.7), with the difference to ρmin as 0.1Δρ and 0.7Δρ (Fig. 2), respectively. Here, ρ(z(0.1)) = ρmin + 0.1Δρ, and ρ(z(0.7)) = ρmin + 0.7Δρ.

Fig. 2
figure 2

Illustration for determination of z(0.1) and z(0.7) for ρ-profile.

Full size image

Pycnocline and thermocline gradients

The data between z(0.1) and z(0.7) is rearranged into [ρi, i = 0, 1, 2, …, I] with [ρ0 = ρ(z(0.1)), ρI = ρ(z(0.7))]. Vertical gradients are calculated between ρi (i = 1, 2, …, I) and ρ0,

$${G}_{Di}=-\frac{\rho ({z}_{(0.1)})-{\rho }_{i}}{{z}_{(0.1)}-{z}_{i}},\quad {\rho }_{i}=\rho {{\rm{(z}}}_{i}{\rm{),}}\;i=1,2,\ldots ,I$$
(2)

Their median

$${G}_{D}={\rm{Median}}\left\{{G}_{D1},{G}_{D2},\ldots ,{G}_{DI}\right\}$$
(3)

is used to represent the characteristic gradient for the pycnocline, and it is simply called the pycnocline gradient (GD). Similarly, the same procedure id used to obtain the thermocline gradient (GT). The annual mean (GT, GD) maps show gradients are stronger in tropical regions (20°S–20°N) than middle and high latitudes. In the low latitudes, the gradients are stronger in the eastern than western Pacific and Atlantic. The January and July mean (GT, GD) maps show stronger seasonal variability in the northern hemisphere than in the southern hemisphere (Fig. 3).

Fig. 3
figure 3

Annual, January, and July mean (GT, GD) maps with the left panels for GT and the right panels for GD. The areas in open oceans with the white colour indicate low quality of the identification. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts).

Full size image

ELG for determining MLD (ILD)

Go back to the original profile shown in Fig. 2. Let the number of the data points between z1 and z(0.7) be Ng, and let N = [log2(Ng)] with the bracket indicating the integer part of the real number inside. N is much smaller than Ng. Starting from z1, the (N + 1) exponential leap-forward gradients (ELGs) are calculated at depth zk [between z1 and z(0.7)] (Fig. 4)

$${D}_{n}\rho ({z}_{k})=\frac{\rho \left({z}_{k}\right)-\rho \left({z}_{k+{2}^{n}}\right)}{{z}_{k}-{z}_{k+{2}^{n}}},\;n=0,1,2,\ldots ,N$$
(4)

where Dn is the difference operator: \({D}_{1}\rho ({z}_{k})=[\rho ({z}_{k})-\rho ({z}_{k+2})]/({z}_{k}-{z}_{k+2}),\) \({D}_{2}\rho ({z}_{k})\,=\,\)\([\rho ({z}_{k})-\rho ({z}_{k+4})]/({z}_{k}-{z}_{k+4}),\) etc. . The averaged value among (N + 1) gradients [D0ρ(zk), D1ρ(zk), …, DNρ(zk)] is computed by

$${\widetilde{G}}_{D}\left({z}_{k}\right)=-\frac{\mathop{\sum }\limits_{n=0}^{N}{D}_{n}T({z}_{k})}{N+1}$$
(5)

which represents the gradient effectively at the depth zk with capability to filter out noises in the gradient calculation28,29.

Fig. 4
figure 4

Illustration of the exponentially leap-forward gradient (ELG) method.

Full size image

Since \({\widetilde{G}}_{D}\left({z}_{k}\right)\approx 0\) if zk in the mixed layer; \({\widetilde{G}}_{D}({z}_{k})={G}_{pyc}\) if zk in the pycnocline, it is reasonable to use (an order of smaller gradient in mixed layer than in pycnocline),

$${\widetilde{G}}_{D}\left({z}_{k}\right)/{G}_{pyc} < 0.1,\quad {z}_{k}\;{\rm{in}}\;{\rm{mixed}}\;{\rm{layer}}$$
(6)

to identify hD (similarly hT). The annual mean (hT, hD) maps (upper panels in Fig. 5) show deeper ILD and MLD in the western than eastern Pacific and Atlantic in low latitudes (30°S–30°N). The January and July mean (hT, hD) maps (middle and lower panels in Fig. 5) show stronger seasonal variability with deeper (hT, hD) in January (July) in the northern (southern) hemisphere, deeper (hT, hD) in the Gulf Stream and Kuroshio areas in January, deeper (hT, hD) in east of Southern Pacific in July, and deeper hT in the Circumpolar current west of the Drake Passage in July.

Fig. 5
figure 5

Annual, January, and July mean (hT, hD) maps with the left panels for hT and the right panels for hD. The areas in open oceans with the white colour indicate low quality of the identification. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts). The red bull-eye in the southeast Pacific in the right-lower panel indicates the area of hD = 150 m.

Full size image

Barrier (Compensated) layer

Barrier layer occurs if hT > hD, with the barrier layer depth (BLD) of (hT − hD). Compensated layer occurs if hT < hD, with the compensated layer depth (CLD) of (hD − hT)15. We generate the difference (hT − hD) data to identify the occurrence of barrier or compensated layer, with annual, January, July mean (hT − hD) maps. The occurrence of barrier layer is much often than the occurrence of compensated layer with the ratio of 23562/170 (~139) from the annual data. Evident barrier layer occurs in the extra-tropical (20°–30°N, 10°–28°S) eastern Pacific Ocean, tropical (0°–8°S) western Pacific west of New Guinea, and low latitudinal Brazilian coast in the annual mean; Kuroshio and Gulf Stream extension regions in January; and south (30°–40°S) Southern Atlantic Ocean in July (Fig. 6).

Fig. 6
figure 6

Annual, January, July mean (hT − hD) maps. The areas in open oceans with the white colour indicate low quality of the identification. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts).

Full size image

ITL temperature and ML density

The ITL temperature is the vertically averaged temperature from the sea surface (z = 0) down to ILD (z = − hT),

$${T}_{ITL}=\frac{1}{{h}_{T}}\underset{-{h}_{T}}{\overset{0}{\int }}Tdz$$
(7)

The ML density is the vertically averaged ρ from the sea surface (z = 0) down to MLD (z = −hD)

$${\rho }_{ML}=\frac{1}{{h}_{D}}\underset{-{h}_{D}}{\overset{0}{\int }}\rho dz$$
(8)

Figure 7 shows the annual, January, July mean (TITL, ρML). Note that (TITL, ρML) are not the same as the sea surface temperature and density. The variables (TITL, ρML) follow the ocean mixed layer dynamics. For example, TITL satisfies following equation due to the ITL layer heat balance46

$${h}_{T}\frac{\partial {T}_{ITL}}{\partial t}+{h}_{T}{\bf{V}}\,\bullet \,\nabla {T}_{ITL}+\Lambda {w}_{e}\Delta T+\nabla \,\bullet \,\underset{-{h}_{T}}{\overset{0}{\int }}\widehat{{\bf{V}}}\widehat{T}dz=\frac{{Q}_{0}-{Q}_{-{h}_{T}}}{\rho {c}_{p}},\quad \Delta T\equiv {T}_{ITL}-{T}_{th}$$
(9)

where V is the vertically averaged horizontal velocity from the surface to −hT; \(\widehat{{\bf{V}}}\) and \(\widehat{T}\) are deviation from the vertical average; Q0 is the net surface heat flux adjusted for the penetration of light below the ITL; \({{\rm{Q}}}_{-{{\rm{h}}}_{{\rm{T}}}}\) is the vertical turbulent diffusion at the base of ITL; we is the entrainment velocity; Λ is the Heaviside unit function taking 0 if we < 0, and 1 otherwise due to the second law of the thermodynamics; ΔT is the temperature difference between ITL and thermocline; Tth is the thermocline temperature.

Fig. 7
figure 7

Annual, January, July mean (TITL, ρML) with the left panels for TITL and the right panels for ρML. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts).

Full size image

ITL heat content and ML freshwater content

The ITL heat content is the integrated heat stored in the layer from the sea surface (z = 0) down to ILD (z = − hT),

$${H}_{ITL}={c}_{p}\underset{-{h}_{T}}{\overset{0}{\int }}\rho (z)T(z)dz$$
(10)

where cp = 3,985 J kg−1 °C−1, is the specific heat for sea water. The annual, January, and July mean HITL show that HITL is higher in low latitudes than in middle and high latitudes and is higher in the western than eastern Pacific and Atlantic Oceans within the low latitudes (left panels in Fig. 8). The ITL heat content HITL represents the warming/cooling of the ITL only. However, the fixed-depth OHC represents warming/cooling of ITL, combined ITL-part of thermocline, or combined ITL-thermocline-part of deep layer depending on the selection of depth with various warming trends. For example, the warming trend is estimated as 0.64 ± 0.11 W m−2 for OHC (0–300 m) in 1993–200832, and 0.20 W m−2 for OHC (0–700 m) in 1955–201236. The global ocean ITL warming rate from 1970 to 2017 is identified as (0.14 W m−2)47.

Fig. 8
figure 8

Annual, January, July mean (HITL, FML) maps with the left panels for HITL and the right panels for FML. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts).

Full size image

The ML freshwater content (FML) is the integrated heat stored in the layer from the sea surface (z = 0) down to MLD (z = − hD)48,

$${F}_{ML}=\underset{-{h}_{D}}{\overset{0}{\int }}\left[1-\frac{S(z)}{{S}_{ref}}\right]dz,\;\;{S}_{ref}={\rm{35}}\;{\rm{psu}}$$
(11)

The annual mean FML is mostly positive in the North Pacific Ocean, the Southern Ocean, the Indian Ocean except the Arabian Sea, and mostly negative in the Atlantic Oceans. Four evident negative FML areas are in the Arabian Sea, central and eastern South Pacific Ocean between equator to 25°S, the North Atlantic Ocean from the Caribbean Sea to the west coast of Spain, and the South Atlantic Ocean between the equator to 30°S (upper right panel in Fig. 8).

Seasonal variability is evident with stronger positive FML in the North Pacific Ocean, weaker negative FML in the South Pacific Ocean, stronger negative FML in the North Atlantic Ocean, weaker negative FML in the South Atlantic Ocean in January than in July (central and lower right panels in Fig. 8). Strongest positive FML (>5 m) occurs in the northern (35° N–60° N) North Pacific Ocean in January. Strongest negative FML (<−5 m) occurs in the central to eastern subtropical North Atlantic Ocean in January, and in the eastern subtropical South Pacific Ocean and the subtropical South Atlantic Ocean in July. Note that the ML freshwater content (FML) dataset is also new since it is different from the fixed-depth freshwater content.

Thermocline and pycnocline

Maximum gradients can be obtained from vertical gradients of (T, ρ) between z(0.1) and z(0.7) calculated by Eq. (2),

$${G}_{T\max }={\rm{Max}}\left({G}_{T1},{G}_{T2},\ldots ,{G}_{TI}\right),\;{G}_{D\max }={\rm{Max}}\left({G}_{D1},{G}_{D2},\ldots ,{G}_{DI}\right)$$
(12)

which are used to represent the depths of thermocline (hth) and pycnocline (hpyc)

$${G}_{Ti}\left({z}_{i}=-\,{h}_{th}\right)={G}_{T\max },\;{G}_{Di}\left({z}_{i}=-\,{h}_{pyc}\right)={G}_{D\max }$$
(13)

The temperature at the thermocline depth is defined as the thermocline temperature

$${T}_{th}=T\left(z=-\,{h}_{th}\right)$$
(14)

The density at the pycnocline depth is defined as the pycnocline density

$${\rho }_{pyc}=\rho \left(z=-\,{h}_{pyc}\right)$$
(15)

Annual mean GTmax and GDmax have evident latitudinal variability with large values in the tropical regions and small values outside the tropical region (upper two panels in Fig. 9). In the tropical Pacific, Atlantic, and Indian Oceans, the large values are in the central and eastern parts. Seasonal variability of (GTmax, GDmax) is strong in the Northern Hemisphere and weak in the Southern Hemisphere. In the Northern Hemisphere, the difference is small between annual and January mean (GTmax, GDmax) but large between annual and July mean (GTmax, GDmax). In July, very strong (GTmax, GDmax) occur in the northeastern Asian and American coastal regions, and strong (GTmax, GDmax) appear in the Kuroshio extension and Gulf Stream regions (middle and lower panels in Fig. 9).

Fig. 9
figure 9

Maps of annual, January, July mean (GTmax, GDmax).

Full size image

Annual mean hth and hpyc have evident spatial variability with large values in low latitudes (30°S–30°N) with shallowing of (hth, hpyc) in the eastern Pacific and Atlantic Oceans (upper two panels in Fig. 10). Seasonal variability of (hth, hpyc) is strong in middle and high latitudes with deep (hth, hpyc) in (30°N–60°N) in January and deep (hth, hpyc) in (30°S–60°S) in July (middle and lower panels in Fig. 10).

Fig. 10
figure 10

Maps of annual, January, July mean (hth, hpyc) with the left panels for hth and the right panels for hpyc. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts). The “slice-like” feature may be caused by the vertical resolution of 25 m from 100 m to 500 m depth in the WOA18.

Full size image

Annual mean Tth and ρpyc have evident spatial variability with warm Tth and light ρpyc in low latitudes (30°S–30°N) and cold Tth and dense ρpyc in middle and high latitudes (30°S–60°S, 30°N–60°N) (upper two panels in Fig. 11). Seasonal variability of (Tth, ρpyc) is generally weak except in subtropical oceans where warmer Tth (light ρpyc) in the Sothern Hemisphere in January and in the Northern Hemisphere in July. Furthermore, cold Tth and dense ρpyc appears in North Atlantic Ocean near Greenland in January (middle and lower panels in Fig. 11).

Fig. 11
figure 11

Maps of annual, January, July mean (Tth, ρpyc) with the left panels for Tth and the right panels for ρpyc. The land and the areas of low quality of the identification are represented by white with the land enclosed by black curves (i.e., coasts). The “slice-like” feature may be caused by the vertical resolution of 25 m from 100 m to 500 m depth in the WOA18.

Full size image

Global statistics of the thermohaline parameters

The identified thermohaline parameters are on the grid points. We use the standard Matlab codes to calculate with area-weighted mean, standard deviation, skewness, and kurtosis for each parameter. Table 2 shows the statistical characteristics of thermohaline parameters derived from the annual mean (T, ρ) profiles. These values can be treated as the overall climatological values of global thermohaline parameters, such as 51.4 m for isothermal layer depth, 38.7 m for mixed layer depth, 14.3 m for barrier layer depth, 99.4 m the thermocline depth, 73.9 m for the pycnocline depth, 0.0638 °C/m for thermocline gradient, and 0.0212 kg/m4 for pycnocline gradient.

Table 2 Statistical characteristics of thermohaline parameters derived from the annual mean (T, ρ) profiles.
Full size table

Data Records

This global ocean climatology of thermohaline parameter dataset is publicly available at the NOAA/NCEI data repository as a NetCDF file, which includes data citation, dataset identifiers, metadata, and ordering instructions. The dataset is located at the NOAA/NCEI website (https://doi.org/10.25921/j3v2-jy50).

Technical Validation

From the identified parameters (hT, GT) and (hD, GD), fitted temperature and density profiles [\(\widehat{T}({z}_{k}),\widehat{\rho }({z}_{k})\)] can be constructed using zero gradient in the ITL and ML, and linear gradient (GT, GD) in the thermocline and pycnocline,

$$\widehat{T}({z}_{k})=\left\{\begin{array}{l}{T}_{ITL}\quad \quad {\rm{for}}\;{z}_{k}\ge -\,{h}_{T}\;{\rm{(in}}\;{\rm{ITL)}}\\ {T}_{ITL}+({z}_{k}+{h}_{T}){G}_{T}\quad {\rm{for}}\;{z}_{k} < -\,{h}_{T}\;{\rm{(in}}\;{\rm{thermocline)}}\end{array}\right.$$
(16a)
$$\widehat{\rho }({z}_{k})=\left\{\begin{array}{l}{\rho }_{ML}\quad \quad {\rm{for}}\;{z}_{k}\ge -\,{h}_{D}\;{\rm{(in}}\;{\rm{ML)}}\\ {\rho }_{ML}-({z}_{k}+{h}_{D}){G}_{D}\quad {\rm{for}}\;{z}_{k} < -\,{h}_{T}\;{\rm{(in}}\;{\rm{pycnocline)}}\end{array}\right.$$
(16b)

The methodology is evaluated through the comparison between the fitted profiles [\(\widehat{T}({z}_{k}),\widehat{\rho }({z}_{k})\)] and the WOA18 profiles [\(T({z}_{k}),\rho ({z}_{k})\)].

We take T-profile for illustration. The fitted \(\widehat{T}({z}_{k})\) profile is represented by two lines with the first one (near zero gradient in the ITL) from the top to the ITL base (circles in Fig. 12) and the second one (non-zero gradient in the thermocline) from the ITL base to the bottom of the profile. The quality in determination of (hT, GT) is identified by the sum of the square error (SSE) between the fitted \(\widehat{T}({z}_{k})\) and WOA18 T(zk) at the WOA18 vertical levels (see Table 1),

$$SSE=\sum _{k}{\left[\widehat{T}({z}_{k})-T({z}_{k})\right]}^{2}$$
(17)
Fig. 12
figure 12

The i-index to represent the quality of the ITL depth determination. (a) ITL not existence or ITL in existence, but not identified (hT = 0). (b) Identified ITL depth shallower than the real ITL depth. (c) Perfectly identified ITL depth. (d) Identified ITL depth deeper than the real ITL depth. (e) Identified ITL depth at z(0.7), i.e., thermocline not existence. Here, the ITL depth is marked by a circle. The thick line in each figure is a single linear function fitted to the temperature profile data.

Full size image

Such double-gradient fitting has errors in the ITL represented by the sum of the square error in ITL (SSEITL) and in the thermocline represented by the sum of the square error in thermocline (SSETH). The whole T-profile data fitted to a single linear function (thick lines in Fig. 12) represents the maximum error since it disregards the existence of ITL and thermocline. Such a maximum error is represented by the total sum of the square error (SSET). The identification index (called i-index) is defined by29

$${I}_{ITL}=\sqrt{1-\left(SS{E}_{ITL}+SS{E}_{TH}\right)/SS{E}_{T}}.$$
(18)

If an ITL exists but is not identified (hT = 0) (Fig. 12a), SSEITL = 0, SSETH = SSET; which gives IITL = 0. If the identified hT is shorter than the real one (Fig. 12b), SSEITL = 0, SSETH > 0, SSETH < SSET; which leads to 0 < IITL < 1. If the identified hT is the same as the real one (Fig. 12c), SSEITL = 0, SSETH = 0, SSET > 0; which makes IITL = 1. If the identified hT is deeper than the real one (Fig. 12d), SSEITL > 0, SSETH = 0, SSETH < SSET; which leads to 0 < IITL < 1. If the identified hT reaches the bottom of the thermocline (Fig. 12e), SSETH = 0, SSEITL = SSET; which gives IITL = 0.

Thus, the value of IITL represents the quality of the determination of hT from T-profile: (a) IITL = 1 for no error (Fig. 12c), (b) IITL = 0 for 100% error with no ITL identified but actual existence of ITL (Fig. 12a) or identified hT at the bottom of the thermocline (Fig. 12e,c) 0 < IITL < 1 for identified hT shallower than the actual hT (Fig. 12b) and for identified hT deeper than the actual hT (Fig. 12d).

The annual, January, July mean maps and histograms of the i-index for ILD show high quality of the method. The global average i-index is 0.849 for the annual mean, 0.865 for the January mean, and 0.882 for the July mean (left panels in Fig. 13). The annual, January, July mean maps and histograms of the i-index for MLD show quality of the method. The global average i-index is 0.841 for the annual mean, 0.858 for the January mean, and 0.865 for the July mean (right panels in Fig. 13).

Fig. 13
figure 13

Maps and histograms of the i-index to determine (hT, GT) from the annual, January, July mean temperature profiles (left panels) and (hD, GD) from the annual, January, July mean density profiles (right panels). Here, top panels are for the annual mean; middle panels are for the January mean; and lower panels are for the July mean.

Full size image