World Ocean Thermocline Weakening and Isothermal Layer Warming

: This paper identiﬁes world thermocline weakening and provides an improved estimate of upper ocean warming through replacement of the upper layer with the ﬁxed depth range by the isothermal layer, because the upper ocean isothermal layer (as a whole) exchanges heat with the atmosphere and the deep layer. Thermocline gradient, heat ﬂux across the air–ocean interface, and horizontal heat advection determine the heat stored in the isothermal layer. Among the three processes, the e ﬀ ect of the thermocline gradient clearly shows up when we use the isothermal layer heat content, but it is otherwise when we use the heat content with the ﬁxed depth ranges such as 0–300 m, 0–400 m, 0–700 m, 0–750 m, and 0–2000 m. A strong thermocline gradient exhibits the downward heat transfer from the isothermal layer (non-polar regions), makes the isothermal layer thin, and causes less heat to be stored in it. On the other hand, a weak thermocline gradient makes the isothermal layer thick, and causes more heat to be stored in it. In addition, the uncertainty in estimating upper ocean heat content and warming trends using uncertain ﬁxed depth ranges (0–300 m, 0–400 m, 0–700 m, 0–750 m, or 0–2000 m) will be eliminated by using the isothermal layer. The isothermal layer heat content with the monthly climatology removed (i.e., relative isothermal layer heat content) is calculated for an individual observed temperature proﬁle from three open datasets. The calculated 1,111,647 pairs of (thermocline gradient, relative isothermal layer heat content) worldwide show long-term decreasing of the thermocline gradient and increasing of isothermal layer heat content in the global as well as regional oceans. The global ocean thermocline weakening rate is ( − 2.11 ± 0.31) × 10 − 3 ( ◦ C m − 1 yr − 1 ) and isothermal layer warming rate is (0.142 ± 0.014) (W m − 2 ).


Introduction
The thermocline with a strong vertical gradient is the transition layer between the vertically quasi-uniform layer from the surface, the isothermal layer (ITL), and the deep-water layer [1,2]. Intense turbulent mixing near the ocean surface causes formation of ITL. The thermocline resists the turbulent mixing from the ITL and limits the ITL deepening and heat exchange between the ITL and deeper layer due to its strong vertical gradient. The mean thermocline gradient (G) directly affects the ITL warming or cooling [3][4][5][6].
The oceanic climate change is represented by temporal variation of the upper layer ocean heat content (OHC), which is the integrated heat stored in the layer from the sea surface (z = 0) down to a below-surface depth. Here, z is the vertical coordinate. The below-surface depth can be either a fixed depth or the ITL depth (h). With h as the below-surface depth, the upper layer OHC is called the ITL heat content (H ITL ). Both G and H ITL are important in ocean prediction [6]. Table 1. Several reported trends (warming rate) of fixed depth ocean heat content anomaly.
Here, the Earth coordinate system is used with r = (x, y) for the longitudinal and latitudinal and z in the vertical (upward positive). The second one is the 0.25 • × 0.25 • gridded climatological monthly mean temperature fieldsT c (t c , r c , z c ) of the NOAA World Ocean Atlas 2018 (WOA18) (https://www.nodc.noaa.gov/OC5/woa18/) [20] with t c representing the month, r c the horizontal grid point, and z c the WOA standard depth. The third one is the global ocean thermocline gradient, ITL depth, and other upper ocean parameters calculated from WOD CTD and XBT temperature profiles from 1 January1970 to 31 December 2017 (NCEI Accession 0173210) (https://data.nodc.noaa.gov/cgibin/iso?id=gov.noaa.nodc:0173210) [21], which was established using the optimal schemes [22][23][24]. The third dataset contains 1,202,061 synoptic G(t, r) and h(t, r) for the world oceans [25]. Each data pair has a corresponding observed temperature profile T(t, r, z) in the first dataset.
The decadal variations of statistical parameters for the global (G, h) data are depicted in

Describing Mixed Layer Dynamics
The connection between G (unit: • C m −1 ) and the ITL heat content H ITL (unit: J m −2 ) can be found using the conservation of heat in the ITL [1,[3][4][5]26] H ITL = ρ 0 c p hT I where v (unit: m s −1 )is the bulk ITL horizontal velocity vector; h is the isothermal layer depth (unit: m); T I is the ITL temperature (unit: • C); Q (unit: W m −2 ) is the net heat flux at the ocean surface (upward positive); w e (unit: m s −1 ) is the entrainment velocity; Λ is the step function taking 0 if w e < 0, and 1 otherwise due to the second law of the thermodynamics; ∆T is the temperature difference between ITL and thermocline; T th (t, r) is the thermocline temperature (unit: • C). It is reasonable to assume that the temperature difference is proportional to the thermocline gradient, ∆T = γG, with the proportionality γ (unit: m). Using conservation of mass in the ITL gives the H ITL equation whereT I is the temperature below the ITL; the last term in the right-hand side of (2) occurs only the mixed layer dynamics and thermodynamics are included. For the entrainment regime (w e > 0, Λ = 1), the thermocline gradient G is related negatively to ∂H ITL /∂t. Strong (weak) thermocline gradient causes the decrease (increase) in H ITL . Such an entrainment effect does not show up with the upper ocean heat content with fixed depth ranges commonly used in oceanography such as 0-300 m, 0-400 m, 0-700 m, 0-750 m, and 0-2000 m.

Methodological Procedure
For each observational CTD/XBT temperature profile T(t, r, z) from the first dataset WOD18 [22], a pair of [G(t, r), h(t, r)] is identified from the third dataset with the same time t and horizontal location r. We use the 3D linear interpolation to calculate a climatological profileT c (r, z) at the observational location (r, z) from four neighboring grid points of the temperature fieldT c (t c , r c , z c ) from the second dataset WOA18 [23] for the same month.
The ITL depth (h) was determined from an individual temperature profile using the optimal fitting and exponential leap-forward gradient scheme ( Figure 1). It contains four steps: (1) estimating the ITL gradient (near-zero), (2) identifying the thermocline gradient G, (3) computing the vertical gradient at each depth (non-dimensionalized by G), and (4) determining h with a given threshold (or user input) to separate the near-zero gradient layer (i.e., the ITL) and the non-zero gradient layer (i.e., the thermocline). Figure 1 illustrates the procedures of this method. Interested readers are referred to [22][23][24] for detailed information.
For an observational temperature profile T(t, r, z) the corresponding climatological monthly mean profileT(r, z) is obtained. With h(t, r), we calculate the difference T(t, r, z) −T(r, z) first and then the ITL heat content (H ITL ) with the mean monthly variation removed (i.e., relative ITL heat content), where c p = 3985 J kg −1 • C −1 is the specific heat for sea water; ρ 0 = 1025 kg m −3 is the upper ocean characteristic density. Let be the in-situ and corresponding climatological monthly ITL temperature. With (3) and (4), we have

Time Series of (HITL, G) for Global and Regional Oceans
To further explore the temporal variability for global and regional oceans, we represent 1,111,647 (HITL, G) pairs into the format of [HITL(m, τi, rj), G(m, τi, rj)] with m = 1, 2, …, 12 the time sequence in months; and τi = 1970, 1971, …, 2017 the time sequence in years with i = 1, 2, ..., N (N = 48), then we divide the world oceans excluding the Arctic Ocean into 12 regions (see the first column in Table 4), and average the synoptic data pairs (HITL, G) within each region and the year to obtain 13 time-series of [<HITL(τi)>, <G(τi)>], including global ocean, with the time increment of a year. Differences among the 13 time-series shows the spatial variability of the yearly variation.
The trends of the time series are determined using linear regression (taking <G(τi)> for illustration),

Time Series of (HITL, G) for Global and Regional Oceans
To further explore the temporal variability for global and regional oceans, we  Table 4), and average the synoptic data pairs (HITL, G) within each region and the year to obtain 13 time-series of [<HITL(τi)>, <G(τi)>], including global ocean, with the time increment of a year. Differences among the 13 time-series shows the spatial variability of the yearly variation.
The trends of the time series are determined using linear regression (taking <G(τi)> for illustration),

Time Series of (H ITL , G) for Global and Regional Oceans
To further explore the temporal variability for global and regional oceans, we represent 1,111,647 The trends of the time series are determined using linear regression (taking <G(τ i )> for illustration), Appl. Sci. 2020, 10, 8185 6 of 13 The trend b is calculated from the data <G(τ i )>, where The error in determining the trend b is estimated by the t-test. For a given significance level α, the confidence interval for the real trendb is given by where t α/2 is a value of the t distribution with (N − 2) degrees of freedom.

Decadal Variability of Global (H ITL , G)
In analyzing global spatial and temporal variability, we divided the 1,111,647 synoptic data pairs of [H ITL (t, r), G(t, r)] into six temporal periods: 1971-1980, 1981-1990, 1991-2000, 2001-2010, 2011-2017. For each time-period, we average the (H ITL , G) data in 4 • (longitude) × 2 • (latitude) cell if the number of data pairs is larger than 5, and identify the cell with no data if the number of data pairs is less than 5. Figure 3 shows decadal variations of the global distributions of (H ITL , G). The number of data points are on the right panels. The white grid cells denote no data or a number of data pairs less than

Trends of (H ITL , G) for Global and Regional Oceans
The 13 time series pairs [<H ITL (τ i )>, <G(τ i )>] show the overall temporal (yearly) variability for the global and individual regional ocean. The significant level to estimate the trend of the time series is set as α = 0.05. The trends of G are all negative (i.e., ∂G/∂τ < 0) in the global ocean [(−2.11 ± 0.31) × 10 −3 • C/(m × yr)] (Figure 4) as well as in all the 12 regions (Figures 5 and 6).
The robust weakening trend of the thermocline gradient (∂G/∂τ < 0) and the warming trend of the isothermal layer (∂HITL/∂τ > 0) are found in global and regional oceans except the Southern Ocean, where slightly decreasing in the relative heat content (−0.0272 W m −2 ). The data are very sparse in the Southern Ocean with a total number of observations of 7942 ( Figure 3, Figure 5f, and Table 4) and make the trends in the Southern Ocean not robust. The highest warming rate (0.160 ± 0.028 W m −2 ) of HITL in the Indian Ocean (Table 4)   The robust weakening trend of the thermocline gradient (∂G/∂τ < 0) and the warming trend of the isothermal layer (∂H ITL /∂τ > 0) are found in global and regional oceans except the Southern Ocean, where slightly decreasing in the relative heat content (−0.0272 W m −2 ). The data are very sparse in the Southern Ocean with a total number of observations of 7942 (Figures 3 and 5f, and Table 4) and make the trends in the Southern Ocean not robust. The highest warming rate (0.160 ± 0.028 W m −2 ) of H ITL in the Indian Ocean (Table 4) is coherent with the recent results on the fastest rate of warming in the tropical Indian Ocean among tropical oceans [27].

Comparison in the OHC Calculation between the Isothermal Layer and Fixed-Depth Layer
The dynamical characteristics are different between the isothermal layer and fixed-depth layer (Figure 7). In the isothermal layer, the vertical velocity does not cause the temperature advection since the vertical temperature gradient is near zero. The turbulent heat flux at the base of the mixed layer (z = −h) is downward or zero [3][4][5]. However, such simple dynamic features do not exist in the fixed-depth layer (0-100 m, 0-300 m, or 0-700 m), where the vertical temperature advection is not zero. The heat flux at the base of the fixed-depth layer can be either downward or upward. This causes the OHC computed for an individual water column represented by the corresponding temperature profile having the simple mixed layer dynamics (no vertical advection and downward or zero heat flux at z = −h) with the isothermal layer (i.e., HITL) and more complicated dynamics

Comparison in the OHC Calculation between the Isothermal Layer and Fixed-Depth Layer
The dynamical characteristics are different between the isothermal layer and fixed-depth layer (Figure 7). In the isothermal layer, the vertical velocity does not cause the temperature advection since the vertical temperature gradient is near zero. The turbulent heat flux at the base of the mixed layer (z = −h) is downward or zero [3][4][5]. However, such simple dynamic features do not exist in the fixed-depth layer (0-100 m, 0-300 m, or 0-700 m), where the vertical temperature advection is not zero. The heat flux at the base of the fixed-depth layer can be either downward or upward. This causes the OHC computed for an individual water column represented by the corresponding temperature profile having the simple mixed layer dynamics (no vertical advection and downward or zero heat flux at z = −h) with the isothermal layer (i.e., H ITL ) and more complicated dynamics (nonzero vertical advection and either downward or upward heat flux) with the fixed-depth layer (H 0-D , D = 100 m, 300 m, 700 m, 2000 m . . . ). This makes mean value or trend of H ITL more representative statistically than H 0-D . The basic dynamic characteristics of the fixed-depth layer also cause diversification in the warm trends of the OHC anomaly (see Table 1).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 13 (nonzero vertical advection and either downward or upward heat flux) with the fixed-depth layer (H0-D, D = 100 m, 300 m, 700 m, 2000 m …). This makes mean value or trend of HITL more representative statistically than H0-D. The basic dynamic characteristics of the fixed-depth layer also cause diversification in the warm trends of the OHC anomaly (see Table. 1).

Figure 7.
Comparison between ocean heat content (OHC) for the isothermal layer and fixed-depth layer. Table 3 shows that the global isothermal layer depth has the decadal mean and standard deviation from (38.  [11], and 0.26 ± 0.02 W m −2 during 1975-2009 [14]. The warming rate increases with depth for the fixed-depth layer [11]. The isothermal layer is more effective for the warming than the fixed-depth layer since HITL is comparable to H0-300 m but with the isothermal layer depth much smaller than 300 m.

Conclusions
World ocean thermocline weakening with the rate of (−2.11 ± 0.31) × 10 −3 (°C m −1 yr −1 ), along with ITL warming with the rate of (0.142 ± 0.014) W m −2 , were identified from the three open datasets NOAA/NCEI WOD18, WOA18, and the global ocean thermocline gradient, ITL depth, and other upper ocean parameters calculated from WOD CTD and XBT temperature profiles from 1 January 1970 to 31 December 2017 (NCEI Accession 0173210) (https://data.nodc.noaa.gov/cgibin/iso?id=gov.noaa.nodc:0173210) using the new definition of the upper ocean heat content (i.e., ITL heat content HITL). Such two trends occur in the global as well as regional oceans except the Southern Ocean, where slightly decreasing in the HITL was identified. Weakening thermocline reduces the resistance to the ITL deepening that causes a thick ITL with more heat stored in the ITL (ITL warming). Thermocline weakening may play roles in global climate change in addition to the greenhouse effect from the atmosphere and in phytoplankton growth because strong thermocline inhibits the nutrient pumping. Besides, the global sea surface temperature is directly related to the ITL heat content through the ocean mixed layer dynamics. It is more reasonable to use ITL heat content than the heat content with fix-depth range in climate change studies.  Table 3 shows that the global isothermal layer depth has the decadal mean and standard deviation from (38.  [11], and 0.26 ± 0.02 W m −2 during 1975-2009 [14]. The warming rate increases with depth for the fixed-depth layer [11]. The isothermal layer is more effective for the warming than the fixed-depth layer since H ITL is comparable to H 0-300 m but with the isothermal layer depth much smaller than 300 m.

Conclusions
World ocean thermocline weakening with the rate of (−2.11 ± 0.31) × 10 −3 ( • C m −1 yr −1 ), along with ITL warming with the rate of (0.142 ± 0.014) W m −2 , were identified from the three open datasets NOAA/NCEI WOD18, WOA18, and the global ocean thermocline gradient, ITL depth, and other upper ocean parameters calculated from WOD CTD and XBT temperature profiles from 1 January 1970 to 31 December 2017 (NCEI Accession 0173210) (https://data.nodc.noaa.gov/cgi-bin/iso?id=gov.noaa.nodc: 0173210) using the new definition of the upper ocean heat content (i.e., ITL heat content H ITL ). Such two trends occur in the global as well as regional oceans except the Southern Ocean, where slightly decreasing in the H ITL was identified. Weakening thermocline reduces the resistance to the ITL deepening that causes a thick ITL with more heat stored in the ITL (ITL warming). Thermocline weakening may play roles in global climate change in addition to the greenhouse effect from the atmosphere and in phytoplankton growth because strong thermocline inhibits the nutrient pumping. Besides, the global sea surface temperature is directly related to the ITL heat content through the ocean mixed layer dynamics. It is more reasonable to use ITL heat content than the heat content with fix-depth range in climate change studies.