Trend and Variability in Global Upper ‐ Ocean Strati ﬁ cation Since the 1960s

Many studies on future climate projection point out that with progressing of global warming, upper ‐ ocean density strati ﬁ cation will strengthen over this century, and consequently, global ‐ averaged ocean primary productivity will decrease. Observed long ‐ term changes in the strati ﬁ cation to date, however, still show large uncertainties of the change itself and its driver. Focusing on the vertical difference in the emergence of the global warming signals, we used only observational pro ﬁ les to describe the spatiotemporal characteristic of long ‐ term trend and variability in the upper ‐ ocean strati ﬁ cation (de ﬁ ned as the density difference between the surface and 200 ‐ m depth). Statistically signi ﬁ cant strengthening of the strati ﬁ cation since the 1960s was detected in ~40% of the global ocean area. The global average increase in the strati ﬁ cation corresponds to 3.3 – 6.1% of the mean strati ﬁ cation. The strengthening trends considerably change depending on the regions and show dominant contribution from the tropical region. In addition to the well ‐ documented explanation of strengthening strati ﬁ cation, namely the surface intensi ﬁ cation of global warming signal, we found that changes in subsurface temperature and salinity strati ﬁ cation associated with changes in atmospheric/ocean circulations signi ﬁ cantly contribute to the long ‐ term change in the strati ﬁ cation and setting its regional difference. In midlatitude and high ‐ latitude ocean of the Northern Hemisphere, the long ‐ term trend exhibits noteworthy seasonality, which shows faster increase trend in the summer than in the winter. From the detrended time series, interannual variabilities correlated with a particular climate mode are detected in several ocean regions, suggesting that these variabilities are mainly can provide the reference of the global change in the upper ‐ ocean strati ﬁ cation to date, the results are no more than one descriptive aspect of the long ‐ term change. Further studies are needed to assess the impact of detected strati ﬁ cation change to ocean ventilation and biogeochemical process with the use of more appropriate metrics such as the potential energy of water column (Burchard & Hofmeister, 2008; Yamaguchi et al., 2019). In order to do that, globally expanded and long ‐ term sustained vertical pro ﬁ le observation of both physical and biogeochemical para-meter is essential. With the development of the Biogeochemical Argo observation network, we hope that future investigations will advance our understanding of physical ‐ biogeochemical interactions in the upper ocean.


Introduction
Upper-ocean stratification plays an important role in the climate system and in many oceanic biogeochemical processes. The strength of near-surface density stratification controls the intensity of vertical mixing (Cronin et al., 2013;Qiu et al., 2004), which in turn affects the development of the surface mixed layer (ML) and the entrainment process at the base of the ML. The ML depth directly modulates the oceanic response to atmospheric forcing and the ocean ventilation process that involves the subduction of water masses into the ocean interior, accompanied by heat, carbon, and oxygen. The upper-ocean stratification can also have a direct impact on biogeochemistry by regulating components of the upper-ocean environment crucial for biological productivity, such as light availability for photosynthesis and nutrient supply from the subsurface ocean.
As a consequence of the global warming that has already occurred, globally averaged upper-ocean thermal stratification has been enhanced due to the surface intensification of the warming signal (Rhein et al., 2013). In addition, many studies on future climate projection using climate models in the Coupled Model Intercomparison Project (CMIP) phases 3 and 5 point out that upper-ocean stratification will strengthen in this century (e.g., Cabré et al., 2015;Capotondi et al., 2012;Fu et al., 2016;Moore et al., 2018).
While the strengthened stratification may produce better light availability for the phytoplankton community, it will also prevent vertical nutrient supply to the euphotic zone (Doney, 2006). Long-term observational studies in the North Pacific/Atlantic using repeat hydrographic cruise data and/or ocean station data have reported both decreases and increase in biological productivity, partly due to change in the upper-ocean stratification (Chiba et al., 2004;Corno et al., 2007;Wallhead et al., 2014;Watanabe et al., 2005). In the context of the future projection, CMIP5 projections suggest that the global average of oceanic primary production will decrease, but there are still large uncertainties due to the range of projected changes in the upper-ocean stratification (Fu et al., 2016). Quantifying the change in the strength of the upper-ocean stratification is necessary not only for understanding the ocean response to the radiative forcing that causes global warming but also for assessing accurately the impacts on biological processes.
Although the global change in the upper-ocean heat content is elaborately measured and carefully assessed by many previous studies, vertical differences of the warming signal are less described. Observational evidence of both increasing and decreasing long-term changes in the upper-ocean stratification to date has only been determined for few individual ocean regions (Dave & Lozier, 2013;Somavilla et al., 2017) or for the global average on the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5: Levitus et al., 2009;Rhein et al., 2013). The strengthening of the stratification that has been estimated from globally averaged temperature field may not occur homogeneously, as the results of CMIP projection studies indicate spatially nonuniform changes in the stratification (Cabré et al., 2015;Capotondi et al., 2012).
Indeed, Somavilla et al. (2017) reported that the stratification north of Hawaii decreased with a large amplitude of decadal variability from the 1990s, although the sea surface temperature (SST) increased. A recent investigation over large areas of the low-latitude and midlatitude Pacific oceans using profile data from 1997 to 2010 also reported a trend of decreasing stratification (Dave & Lozier, 2013). In addition to the spatial heterogeneity of the long-term trend in the stratification, these results also suggest that the simplest relation, namely that ocean warming is intensified near the surface and will result in increases in the local stratification, does not always hold in all ocean regions and there are some other drivers of the change in stratification.
In the present study, we focused on the vertical structure of the change in the ocean temperature and salinity, that is, the change in the density stratification. The purpose of this study is to describe the long-term change in the global upper-ocean stratification and evaluate the impact of the global warming on the stratification at this present time, using as many available raw historical observational profiles as possible. We also aim to detect the interannual/decadal variability in the upper-ocean stratification in individual ocean regions from the long-term time series. The remainder of this paper is organized as follows. The data set and processing methods are explained in section 2. In section 3, we describe long-term trends in the upper-ocean stratification, its vertical structure, its seasonality, and contributions of thermal and/or haline stratification to these trends. The relationships between interannual to decadal variability of the stratification and climate mode in individual ocean regions are investigated in section 4. We finally summarize obtained results and discuss their implications in section 5.
To quantify the strength of the upper-ocean stratification, we use here a metric defined as the potential density difference between the surface and 200-m depth (Δρ200). This simple indicator is not necessarily optimal for representing detailed characteristics of the upper-ocean stratification, which consists of mostly homogeneous density in surface ML and the density change below the ML (Somavilla et al., 2017). However, as it has been widely used in both model and observational studies, it allows us to compare our results with those of previous studies, to assess their reliability. Moreover, the simplicity of Δρ200 helps to avoid influences of the changes in observational instrument and its vertical resolution (cf. Figure 1g) in analytical procedures of this study. We have checked that the results were qualitatively unchanged if the choices of the lower depth were alternated by 150-and 250-m depth. Also, the vertical difference of the long-term changes in the density field is discussed also in section 3.2.
All T/S profiles are linearly interpolated onto the surface (assigned at 10-m depth) and 200-m depth, and then the potential density and stratification (Δρ200) are calculated from individual profiles. To further control data quality, in addition to the WOD13 QC procedure for original T/S profiles, we exclude individual Δρ200 values that depart from the monthly mean value by three standard deviations in each 1°(latitude) × 1°(longitude) grid cell. Monthly climatology is obtained by averaging the quality-controlled Δρ200 values in each 1°(latitude) × 1°(longitude) grid cell (Figure 2a).
Using a linearized equation of state, we also evaluate the contributions of thermal stratification anomaly (Δρ T 200′) and haline stratification anomaly (Δρ S 200′) to the density stratification anomaly (Δρ200′), as follows: Figure 1. (a-f) Spatial distribution of T/S profiles used in this study and (g) temporal change of the number (bars; left axis) and the mean vertical resolution (black curve; right axis) of the profiles. The coloring follows the WOD13 data set categories: Ocean Station Data (OSD, brown, including low-resolution CTD/XCTD data), Conductivity-Temperature-Depth data (CTD, red, including high-resolution XCTD data), Undulating Oceanographic Recorder data (UOR, yellow), Profiling Float data (PFL, gray), Moored Buoy data (MBR, green), and Glider data (GLD, blue). The vertical resolution is defined as the number of observation layers in each profile from the surface to 200-m depth. "Number of layers = 100" in (g) means there are 100 observations from the surface to 200-m depth in a profile, and thus, the vertical resolution of the profile is approximately 2 m.
where the subscript "1" and "2" stand for surface and 200-m depth value, respectively, and, α and β are the thermal expansion and haline contraction coefficients of seawater, respectively. Note that the discrepancy between Δρ200′ and the sum of Δρ T 200′ and Δρ S 200′ is very small in the T/S and pressure ranges of this analytical procedures, so that equation (1) provides good estimates of T/S contributions (cf. Figures 2 and 4a-4c).
In general, there is the predominant seasonal cycle in the upper ocean forced by the marked seasonality in the heat/freshwater exchange with the atmosphere driving the increase/decrease in the density stratification, as well as their control in the extent of vertical mixing, especially in midlatitude and high-latitude regions. To avoid artificial trend and variability due to seasonal sampling bias, all long-term trends in this study are computed using anomalies from the monthly climatology. We calculated density stratification anomalies (Δρ200′) by subtracting 1°× 1°gridded monthly climatology of nearest grid point from individual Δρ200 values. Then, to remove mesoscale or smaller-scale signals, these anomalies are yearly averaged over 5°(latitude) × 10°(longitude) (hereafter, called annual anomaly). Long-term trends in each grid cell are calculated from the annual anomalies using a least-squares fit. We assessed the statistical significance using a Mann-Kendall rank statistic for linear trends and using Student's t test for correlation coefficients of regionally averaged time series, with estimates of the degree of freedom based on the zero-crossing correlation timescale. To investigate the relationship between the interannual variability of the stratification and climate mode in individual ocean regions, we used the Niño 3.4 index (Trenberth, 1997; https://www.esrl.noaa.gov/psd/data/ correlation/nina34.data), the Pacific Decadal Oscillation (PDO) index (Mantua et al., 1997; http://research. jisao.washington.edu/pdo/PDO.latest.txt), and the North Atlantic Oscillation (NAO) index (Hurrell, 2003; https://climatedataguide.ucar.edu/sites/default/files/nao_station_annual.txt).

Long-Term Change in the Annual Mean Stratification
Linear trends of upper-ocean density stratification (Δρ200) estimated from the annual anomalies with different starting years of the trend calculation are shown in Figure 3. In the results with relatively long analytical periods, statistically significant strengthening of stratification is detected in many regions of the global ocean (Figures 3a and 3b). The spatial pattern and the values of the trends calculated with different starting years converge to those calculated since the 1960s as the analytical period becomes longer. Remarkably, this spatial distribution (Figure 3a) is similar to that of future projections obtained by CMIP3 models (Capotondi et al., 2012, their Figure 15) and CMIP5 models (Cabré et al., 2015, their Figure 1), which can be thought of as the oceanic response to the radiative forcing that causes global warming. These suggest that we could capture the oceanic response to the radiative forcing associated with global warming that proceeded.
As shown by previous study using shorter-term observational data (Dave & Lozier, 2013;Somavilla et al., 2017), the trends since the 1990s, 2000s, and 2010s are statistically significant to a lesser extent for both strengthening and weakening (in some places) This is likely because of the prevailing long-term variability on decadal and/or interannual timescales. Hereafter, we use the trend calculated since the 1960s to discuss the long-term change in the stratification.
In Figure 4, we show the spatial distribution of the long-term change in the upper-ocean density stratification since the 1960s and its decompositions to the contribution of T/S changes. Figure 4a (same as Fig. 3a) indicates that the long-term change in the upper-ocean stratification is highly spatially nonuniform. The strongest trends in the density stratification occur in the tropical Pacific-Indian warm pool region. Statistically significant trends are also detected in the eastern tropical Pacific, subarctic and along the eastern boundary of the North Pacific, and from the equator to the subtropical North Atlantic. The density stratification trends are relatively weak and less statistically significant south of 30°S. Although the trend is positive (strengthening) in the~80% of the global ocean area, the statistically significant trends accounts for only 40% of the ocean where the trends are able to be estimated, and 50% of them concentrate in the tropical region (20°S-20°N).
Decomposing the density stratification into thermal and haline stratification components using equation (1), we estimated the contributions of the changes in thermal and haline stratification to trends in the density stratification (Figures 4a-4c). The thermal stratification trends contribute most to the density stratification trends; that is, the spatial distribution of the density stratification trend is caused mainly by changes in the vertical thermal structure. Conversely, salinity stratification trends show contributions to the density stratification trend comparable to that of the thermal stratification trends only in regions where there are significant salinity components of the density stratification in the climatological field ( Figure 2), such as the subarctic North Pacific, the tropical Pacific-Indian warm pool region, and the Southern Ocean near Antarctica. The salinity stratification changes also contribute to a reduction in the density stratification in the South Atlantic and the subtropical North and South Pacific.
In accordance with the definition of density stratification, thermal and haline stratification components can be further decomposed into density trends due to changes in temperature and salinity at the 10 and 200-m depth (Figures 4d-4i). Relatively spatially uniform negative surface density trends caused by SST warming contribute to the strengthening of density stratification ( Figure 4e). The spatial pattern corresponds well with the SST warming since the 1960s, which shows weaker warming in the central North Pacific (around 40°N) and Southern Ocean (Armour et al., 2016;Huang et al., 2019).
Strong positive density trends due to subsurface cooling also contribute to the strengthening of density stratification in the western tropical Pacific and the equatorial-side margins of the Pacific subtropical gyre (Figure 4h). Since the subsurface cooling trends are located around the western part of the tropical Pacific, it is suggested that the equatorial part of the subsurface cooling is associated with weakening of the Walker circulation from the mid-twentieth century (Tokinaga et al., 2012). In climatological equatorial Pacific, since the easterlies pile up warm water in the western part, the equatorial thermocline located at upper depth than 200 m has east-west gradient that is deeper (shallower) in the western (eastern) part. The weakening of the easterlies in the tropical Pacific flattens the east-west gradient of the equatorial thermocline and results in subsurface cooling (warming) anomalies in the western (eastern) part of the equatorial Pacific. These subsurface anomalies, therefore, contribute to strengthen (weaken) the density stratification in the western (eastern) equatorial Pacific.
Subsurface temperature warming trends stand out especially within the subtropical gyres of each basin ( Figure 4h). The subsurface decrease in density due to temperature warming contributes to partially compensate the strengthening density stratification caused by surface warming (Figures 4b, 4e, and 4h). In general, the subsurface water within the subtropical gyre has its origin in a subduction process in which the water mass formed in the winter surface mixed layer and spreads across the gyre (Luyten et al., 1983;Suga et al., 2004). The subducted water masses enhance warming of subsurface water by carrying the strong surface warming signals occurring in western boundary current regions into the subsurface (Sugimoto et al., 2017). In the subtropical North Pacific/Atlantic, the southern Indian Ocean, and the South Pacific, the strengthening of density stratification by surface warming is mitigated by comparable subsurface warming.
To understand the relative importance, we show, in Figure 5, the relationships between the contributions of sea surface density and subsurface density changes to upper-ocean stratification changes. In 55% of the global ocean area, subsurface warming mitigates (47.7%) or reverses (7.8%) the strengthening of density stratification by surface warming (third quadrant in Figure 5, left). Conversely, in 25.7% of the global ocean area, an increase in subsurface density accelerates the strengthening of density stratification due to subsurface cooling and/or salinification, complicating the traditional attribution of increased density stratification to surface warming (second, third, and fourth half of quadrant in Figure 5, left). The enhanced strengthening of the density stratification by increase in subsurface density occurs mainly in the tropical Pacific and Indian Ocean ( Figure 5, right).
The trend field of the sea surface salinity (SSS) component shows systematic spatial distribution of its contributions to the change in surface density (Figure 4f). Negative surface density trends caused by decreased salinity occur in the subarctic North Pacific, the tropical warm pool region, and the Southern Ocean. Positive  Journal of Geophysical Research: Oceans 3.2. Regional Trend and Its Vertical Structure To discuss regional and vertical differences in the characteristic of the long-term trend, we computed regional trends for nine ocean regions. Figure 6 shows the regional trends in density stratification estimated from the time series of annual Δρ200 anomaly averaged over individual ocean regions without spatial interpolation; that is, using only grid cells where data exist for regional averaging (missing grid cells are neglected). In calculating the time series, we only used years when the spatial data coverage exceeds 20% of the total area of each ocean region. We also checked following results are qualitatively unchanged if the criteria were altered from 1% to 30%. Although there are few missing years in the regional time series (Figure 6), the number of available year (N) exceeds the 80% of the analytical periods (i.e., 46 years out of 58 years), except for the Arctic Ocean.
As expected from the results in Figure 3a, significant strengthening of the upper-ocean stratification has continued since the 1960s in all ocean regions, except in the Arctic Ocean where the regional trend estimated since 1979 is statistically insignificant. Global-averaged strengthening trend estimated by area-weighted averaging of the regional trends is 0.0191 kg/m 3 /decade. Therefore, the global upper-ocean density stratification increased by~0.11 kg/m 3 over these 58 years. This increase corresponds to~6.1% of the climatological annual mean stratification (1.82 kg/m 3 ). In each region, the largest and dominant change occurs in the tropical Pacific, at a rate of 0.0454 kg/m 3 /decade, followed by the Indian Ocean. The globally averaged trend decreases by 37% to 0.0120 kg/m 3 /decade when the tropical Pacific is excluded and by 58% to 0.0081 In Figure 7, we show the vertical structures of the change in the density and the contributions from temperature and salinity changes in individual ocean regions. The large rate of strengthening of density stratification in the Tropical Pacific and Indian Ocean are attributed to extremely small or slightly positive changes in the subsurface density, in addition to strong surface warming. Conversely, in the North/South Pacific and North Atlantic, as the surface warming signals substantially penetrate to the subsurface, the strengthening trend of the upper-ocean stratification is reduced. Therefore, the subsurface response to forcing associated with global warming is also significantly important for setting regional stratification change.
Contributions from salinity changes are not negligible. Although the sign of the trend in density change follows the trend in temperature in most of the ocean regions and depths, the magnitudes of the long-term density change cannot be explained without components of salinity change, especially in the Atlantic and Arctic Ocean. In the Southern Ocean, the long-term change in density due to salinity change is larger than that due to the warming near the surface.

Seasonality of the Long-Term Trend
The upper-ocean state drastically changes with season, and the seasonality is crucial for regulating not only its physical conditions but also biogeochemical activities. In this section, we examine the seasonality of the long-term trend in the density stratification. To take account of the seasonal variability of the upper ocean, we attempt to estimate the long-term trends from the yearly time series of seasonally averaged Δρ200 anomalies. Here the seasons are defined using the maximum density stratification at each grid points. We determined the month of annual maximum of Δρ200 at each grid points from the monthly climatology (Figure 8a). And then, we defined "Season I" on each grid points as three consecutive months centered on the month of annual maximum of Δρ200. Seasons II, III, and IV are defined as the consecutive 3-month groupings following season I.
Because the number of available data for computing regional yearly time series of each season-averaged Δρ200 anomalies is inevitably smaller than that of annual average Δρ200 anomalies, ocean regions where the trends are able to be estimated based on the same criteria as the case of annual mean time series are  Table 1, we show regional trends estimated from the yearly time series of seasonally averaged Δρ200 anomalies. Note that the seasonal trends estimated from insufficient number of available year (less than 46 years) are also shown in Table 1. Focusing on the regions where the trends can be obtained for all four seasons (i.e., the North Pacific/Atlantic and the tropical Pacific/Atlantic), the strengthening trends in the density stratification differ depending on the season, and its seasonality is different between the tropics and the other regions ( Figure 9). However, it should be noted that the seasonal differences of the trends do not exceed the range of the uncertainties in these four ocean regions, except for the North Atlantic.
In the tropics of both the Pacific and Atlantic oceans, the density stratification tends to strengthen more in the Seasons I and IV than the other seasons. The strengthening trend is largest in the Season I (i.e., warm season) and smallest in the Season III (i.e., cool season) in the midlatitude and high-latitude ocean of the Northern Hemisphere. As there is predominant seasonal cycle of surface ML depth in midlatitude and high-latitude ocean, the degree of the response of the ML temperature/salinity to the same heating/freshwater forcing associated with the global warming is expected to be different among the season. Indeed, relatively thin surface MLs in warm season trap heat/freshwater forcing in the thin layer and induce larger T/S anomaly near the surface than those in cool season ( Figure 10). Consequently, these seasonal differences in the vertical distribution of the emergence of global warming signal result in the seasonality of the long-term trend in the density stratification in midlatitude and high-latitude ocean.

Detrended Variability and Climate Mode
For regions with relatively dense data coverage and spatiotemporally unbiased observations (i.e., again for the North Pacific/Atlantic and tropical Pacific/Atlantic), we investigated the relationship between the detrended variability in the density stratification and the prevailing climate mode in individual ocean regions.
In the tropical Pacific region, the upper-ocean density stratification is highly significantly correlated with the Niño 3.4 index (Figure 11c). When the El Niño (La Niña) occurs, the SST is higher (lower) in the eastern tropical Pacific, and thus, the density stratification is intensified (weakened), while in the western tropical Pacific, the subsurface cooling (warming) anomalies due to the zonally flattened (tilted) thermocline strengthen (weaken) the density stratifications. SSS freshening anomalies during the El Niño, which have the maxima in the western part (Singh et al., 2011), also contribute to strengthen the density stratification. It is well known that the El Niño and La Niña phenomena are closely related to the weakening and strengthening of the Walker circulation on an interannual timescale, named as the Southern Oscillation (Trenberth & Hoar, 1996). Thus, this correlation is interpreted as being analogous to the explanation in section 3.1 of the increasing trends of density stratification due to the weakened Walker circulation in the tropical Pacific.
As the whole of the North Atlantic, the variability in density stratification has no simultaneous correlations with any prevailing climate mode such as the NAO or Atlantic Multidecadal Oscillation. However, the detrended time series shows a statistically significant negative lagged correlation with the 2-year leading NAO index (Figure 11b). The positive NAO phase is characterized by warmer SST south of the Gulf Stream and cooler SST north of 40°N due to stronger westerly winds (Visbeck et al., 2003). The negative correlation between the density stratification and the NAO index is qualitatively consistent with the cooler (warmer) SST and stronger (weaker) near-surface mixing induced by enhanced (weakened) westerly winds to the north of 40°N during positive (negative) NAO. Conversely, the warmer SST south of the Gulf Stream during positive NAO is inconsistent with this negative correlation. The regions of the warmer SST signal in positive NAO roughly correspond to and/or include the sites of water mass formation and indicates statistically significant trend exceeding the 99% confidence level, and the value within the parenthesis is estimated from insufficient number of available years (less than 46 years).

10.1029/2019JC015439
Journal of Geophysical Research: Oceans subduction. Therefore, the variability in density stratification induced by the SST changes could be smaller in these regions because the anomalous SST signals can be carried well below the surface by subduction and thus the subsurface water also has the same signal as the surface. Indeed, Δρ200 time series averaged over the region of warm SST in positive NAO does not show statistically significant positive correlation, while the region north of 40°N has significant simultaneous negative correlation with the NAO index (data not shown).
Note also that the correlation between variation in the density stratification and NAO reaches its maximum with a NAO lead of 2 years. This feature supports the explanation of SST-driven stratification variability because cool SST anomalies induced by positive NAO extend over the wider region of the North Atlantic with a 2-year lag (Visbeck et al., 2003). Therefore, it is suggested that the lagged negative correlation of the density stratification in the North Atlantic with the NAO is largely explained by its SST variability.
Interannual variability in the upper-ocean stratification in the tropical Atlantic is also strongly affected by SST variability induced by the NAO. Its time series of Δρ200 shows statistically significant simultaneous negative correlation with NAO index. It is consistent with the third pole of NAO SST variability, which is the cooler (warmer) SST south of 20°N during positive (negative) NAO.
Averaged time series over the whole of the North Pacific do not have a statistically significant correlation with the PDO index (data not shown). However, the SST pattern of the PDO, which shows the prevailing variance on interannual to decadal timescales over the entire basin, has a spatial pattern of two poles of the opposite sign in the midlatitude region. To examine the possibility that the opposite signals cancel each other when averaging over a large area, we divide the North Pacific into two regions corresponding to the two poles of PDO SST variation: one is the southwestern North Pacific and the other includes the northern North Pacific and the region along the North American coast (Figure 6, upper map), following Newman et al. (2016). We then find a significant positive correlation in the northern and eastern regions and significant negative correlation in the Figure 9. Regional trends estimated from seasonally averaged Δρ200 anomalies in the tropical Pacific (TP), the tropical Atlantic (TA), the North Pacific (NP), and the North Atlantic (NA). Filled symbols indicate statistically significant trend exceeding the 99% confidence level and error bars show the standard errors of the trend. Figure 10. Seasonality of the long-term trend in the upper-ocean density since the 1960s in the North Pacific and North Atlantic. Red and blue lines with error bars indicate density trends and its standard error at individual depths calculated from each season-averaged density anomalies. The warm (cool) season is defined as five consecutive months centered on the Season I (Season III). Gray bars are density trend estimated from annual mean anomalies (same as Figure 7). southwestern regions covering the subtropical gyre (Figure 11a). These relationships are qualitatively consistent with cooler SSTs associated with stronger westerly winds in the southwestern region and relatively warmer SSTs in the northern and eastern region for a positive PDO phase, and vice versa. The correlation coefficients, however, are relatively small especially in the southwestern region (R = −0.43), although it is statistically significant. This suggests that subsurface variability is also an important factor on the interannual variation of the stratification.
In this study, we could not analyze the variability in stratification and its relationship to the climate mode in oceans of the Southern Hemisphere and polar regions because the observational data were not able to resolve the year-to-year variability over long periods. Regionally averaged time series over the Indian Ocean show no correlation with both Indian Ocean Dipole Mode index (Saji et al., 1999) and Niño 3.4 index, although these correlations are estimated from time series with a few gaps (cf. Figure 6).

Summary and Discussion
We investigated the long-term trends and variability in upper-ocean stratification since the 1960s only using observational profiles. We detect that the statistically significant strengthening of the upper-ocean density stratification in the~40% are of the global ocean, defined as the potential density difference between the surface and 200-m depth (Δρ200). The half of this significant strengthening concentrates in the tropical region, and the spatial distribution of this strengthening is similar to the future projections by CMIP3 and 5 climate models (Cabré et al., 2015;Capotondi et al., 2012). In the global average, the strengthening over these 58 years corresponds to an increase of 6.1% of the climatological mean stratification, including the dominant contribution from tropical regions to the global trends.
This result appears to be larger than the estimate of~4% increase in density stratification from 1971 to 2010 provided by the IPCC AR5 (Levitus et al., 2009;Rhein et al., 2013). Although there are some subtle differences in the two analyses, such as the periods and areal coverage, the main cause of the discrepancy seems to be whether the anomalies are spatially interpolated. Generally, widely used objective interpolation methods tend to attenuate anomalies, especially where data are relatively sparse, such as for the 1960s/1970s and in the Southern Hemisphere. Spatial interpolation is not used in this study, so the present estimates may represent the upper limit of the range that includes the true value. Indeed, the same analysis with the assumption that grid cells of missing data have zero anomalies showed a smaller increase in global density stratification (~3.3%).
In addition, it should be noted that in regions where the seasonal cycle cannot be resolved, estimated annual trends tend to be based on the data taken during the Season I and thus the long-term trends shown in Figures 3a and 6 could be biased to the summertime trends. Indeed, such tendency that the annual trend is similar to the trend of Season I is shown in some of ocean regions outside the tropical region (South Atlantic and Southern Ocean in Table 1 and Figure 6). Assuming a seasonal behavior similar to the remaining ocean regions outside the tropics (North Pacific/Atlantic; cf. Figure 9), the wintertime trends would be smaller, and thus, estimated annual trend in regions where the seasonal cycle cannot be resolved in this study also could be smaller.
Regional stratification trends are spatially nonuniform, but mostly everywhere positive (~80% of the global ocean area) although half of them are not significant. In addition to the well-known contribution of warming SST to strengthening the density stratification, subsurface cooling due to weakening of the Walker circulation also contributes to the strengthening in the tropical warm pool region, while the strengthening is mitigated by enhanced subsurface warming inside of the subtropical gyre. Moreover, the haline stratification changes due to the intensification of the global water cycle and changes in ocean circulation also have significant impacts on the changes in density stratification. It is revealed, therefore, that subsurface changes in temperature and salinity are also important for the estimation of the trends itself and explaining how its regional difference set.
In the midlatitude and high-latitude Northern Hemisphere ocean and the tropics of the Pacific and Atlantic, seasonality of the long-term trend in the density stratification is observed. The strengthening trend in the midlatitude and high-latitude ocean is larger in the warm season than in the cool season, resulting from larger density anomaly near the surface due to thinner ML of the warm season. Although the strengthening trend itself is smaller in the cool season, the impact of this observed stratification change on winter mixing and thus ocean ventilation, nutrient replenishment, and biological productivity still could be significant because of its weak climatological mean stratification. Considering the climatological mean stratification in the Season III (North Pacific: 0.91 kg/m 3 ; North Atlantic: 0.42 kg/m 3 ), wintertime increases in density stratification over these 58 years correspond to 9.2 ± 1.8% and 7.6 ± 2.5% of the climatological mean stratification in the North Pacific and the North Atlantic, respectively.
In individual ocean regions, the variabilities in stratification associated with a particular climate mode are found from the detrended time series. The regionally averaged variability indicates a positive correlation with the Niño 3.4 index in the tropical Pacific, a negative lagged correlation with the NAO index in the North Atlantic, and correspondence with SST variations related to the PDO and NAO in the North Pacific and tropical Atlantic, respectively. As the climate modes are closely related to variability in large-scale atmospheric circulation, oceanic circulation, and biogeochemical processes (Mantua et al., 1997;Martinez et al., 2016), these significant correlations shown in this study suggest the possibility that the variability in density stratification is a key factor to explain the mechanism of the impact of large-scale atmospheric changes on biogeochemical processes, as suggested by Behrenfeld et al. (2006).
The use of the simplest metric representing stratification in this study may have limited our ability to conduct more quantitative analysis. As mentioned above, Δρ200 does not have information about the detail of water column structure such as the ML depth, which is an important factor in air-sea interaction and biogeochemical processes. Although this study can provide the reference of the global change in the upper-ocean stratification to date, the results are no more than one descriptive aspect of the long-term change. Further studies are needed to assess the impact of detected stratification change to ocean ventilation and biogeochemical process with the use of more appropriate metrics such as the potential energy of water column (Burchard & Hofmeister, 2008;Yamaguchi et al., 2019). In order to do that, globally expanded and long-term sustained vertical profile observation of both physical and biogeochemical parameter is essential. With the development of the Biogeochemical Argo observation network, we hope that future investigations will advance our understanding of physical-biogeochemical interactions in the upper ocean. Tohoku University and the participants of Otsuchi Symposium 2018 at the International Coastal Research Center, the University of Tokyo for useful comments and discussions. Also, we would like to express our appreciation to two anonymous reviewers for giving many constructive comments. This work was financially supported through JSPS grant 17J02314. T. S. is financially supported by JSPS through grants 16H04046 and 15H02129. This work also was partly supported by activity of the Core Research Cluster of Disaster Science in Tohoku University (a Designated National University). The processed data used to estimate longterm trends and variability in this study are publicly available at https://doi.org/ 10.5281/zenodo.3466122.