Deep-reaching acceleration of global mean ocean circulation over the past two decades.

Ocean circulation redistributes Earth's energy and water masses and influences global climate. Under historical greenhouse warming, regional ocean currents show diverse tendencies, but whether there is an emerging trend of the global mean ocean circulation system is not yet clear. Here, we show a statistically significant increasing trend in the globally integrated oceanic kinetic energy since the early 1990s, indicating a substantial acceleration of global mean ocean circulation. The increasing trend in kinetic energy is particularly prominent in the global tropical oceans, reaching depths of thousands of meters. The deep-reaching acceleration of the ocean circulation is mainly induced by a planetary intensification of surface winds since the early 1990s. Although possibly influenced by wind changes associated with the onset of a negative Pacific decadal oscillation since the late 1990s, the recent acceleration is far larger than that associated with natural variability, suggesting that it is principally part of a long-term trend.

INTRODUCTION

Earth has experienced rapid warming for decades, driven by increased emissions of greenhouse gases (). As the major heat reservoir for the climate system, the ocean has chronicled this long-term warming trend from the surface to great depths (–). However, these trends are conflated with multidecadal variability (–), for which ocean circulation plays an important role. Therefore, an accurate assessment of changes in ocean circulation is key to understanding global climate change (–).

Regional ocean circulation has responded differently to historical greenhouse forcing in recent decades. Subtropical western boundary currents, the main oceanic conduits for heat transport from the equator to the poles, have warmed () and intensified since the 1900s due to enhanced wind shear between low and high latitudes (). However, the regional boundary current systems have shown diverse tendencies. For example, there is little evidence that the Kuroshio has markedly intensified over the past two decades (, ). Although broadened, the Agulhas Current has not intensified because of increased eddy activity since the early 1990s (). On the contrary, the Pacific Ocean shallow overturning cells have accelerated in response to intensified Pacific trade winds since the early 2000s, contributing to a recent climate warming hiatus () and leading to an increased leakage of heat and fresh water into the Indian Ocean via the Indonesian Throughflow (, , ). Nonetheless, all these trends in regional current systems do not necessarily reflect trends in the global ocean current system.

Because of complex internal dynamics, regional trends usually include strong signals from natural and/or internal variability, more so than global mean variables (). Therefore, assessing the trend of global mean ocean circulation, which averages out short-term regional variations, offers a simpler perspective to understanding global climate change.

However, because of a lack of systematic direct ocean current observations (), the issue of whether there is a trend emerging over in the global mean ocean circulation remains elusive. Here, guided by in situ observations, reanalysis products, and model experiments, we show that a substantial acceleration of global mean ocean circulation has occurred over the past two decades.

RESULTS

Deep-reaching acceleration of the global mean ocean circulation

We refer to “ocean circulation” as the horizontal movement of seawater in terms of a monthly average. For ocean circulation, the usable energy is mainly in the form of kinetic energy (KE), available gravitational potential energy (APE), available internal energy, and a residual term related to dissipation of KE due to molecular and small-scale diffusion and subgrid processes (). Among these terms, KE describes the energy that an object has because of its motion, while APE in the ocean is the potential energy and is dependent on an object’s mass and vertical distance from a certain reference level (). Thus, KE directly describes the movement of seawater and is a suitable index for measuring the intensity of the ocean circulation. In contrast, the available internal energy makes only a small contribution and negative energy to the total potential energy (). We calculate the total global integral of KE [total KE (TKE)] using absolute geostrophic velocity data derived from Argo profiling observations (AGVA) (), five reanalysis/assimilation products based on observations (Fig. 1), and 12 numerical simulations (fig. S1). See Materials and Methods and table S1 for details on the datasets and methods.

Fig. 1

TKE anomalies and linear trends in various datasets.

(A to F) Thick (thin) color lines denote monthly TKE′ with (without) 61-month running mean in Estimating of Circulation and Climate of the Ocean, phase 2 (ECCO2), European Centre for Medium-Range Weather Forecasts (ECMWF) ORA-S3, Global Ocean Data Assimilation System (GODAS), AGVA, Geophysical Fluid Dynamics Laboratory Ensemble Coupled Data Assimilation (GFDL ECDA), and ECMWF ORA-S4 datasets, respectively. Black lines are the linear trends of the 61-month low pass–filtered series. (G) The ensemble-mean series of the TKE′ from the six reanalysis datasets. (H) Linear trends of global ocean mean KE anomaly before 1990 [Period A: ECMWF ORA-S3 and ORA-S4 1959–1990, GFDL ECDA 1960–1990, and National Oceanic and Atmospheric Administration (NOAA) GODAS 1980–1990] and after 1990 (Period B: ECMWF ORA-S3 1991–2011, ECCO2 1991–2013, GFDL ECDA 1991–2012, NOAA GODAS 1991–2011, AGVA 2005–2010, and ECMWF ORA-S4 1991–2013). SEs are shown in blue bar, and trends that are statistically significant at the 99% confidence level are denoted in blue circles, while nonsignificant trends are signposted with red “×”. SE of TKE trend in individual product is defined as the SD of the detrended time series, while that of ensemble mean is defined as the SD of trends from six reanalysis.

TKE anomaly (TKE′) is calculated by subtracting the time mean from the monthly series of TKE and then low pass–filtered with a cutoff period of 61 months. Despite the discrepancy in absolute magnitude of TKE′ among these data products, a common pattern emerges: a clear shift in the multidecadal trend of TKE′ (Fig. 1). Before 1990, TKE′ varied slowly but increased sharply after about 1990 (Fig. 1, A to G). The TKE′ trends from the reanalysis datasets range from about 0.3 × 1017 to 1.8 × 1017 J decade−1 in part due to their differences in terms of spatial domain, vertical depth, resolution, and time span. Despite these differences, all the datasets indicate a common increasing trend. The ensemble mean of the linear trends in TKE′ is approximately (−0.01 ± 0.16) × 1017 J decade−1 before 1990, which is not statistically significant. The ensemble mean trend over the post-1990 period is approximately (1.37 ± 1.10) × 1017 J decade−1 and is significantly different from zero at the 99% confidence level. This rate of change is equivalent to an increase of (15 ± 12)% decade−1 over the climatological mean of 9.05 × 1017 J. The difference in trends before and after 1990 may be partly due to the fact that some reanalysis products [e.g., Estimating of Circulation and Climate of the Ocean, phase 2 (ECCO2) and AGVA] with larger trends were only used for the post-1990 period. A modified Mann-Kendall test () is used to estimate the confidence level of the trends in the time series of TKE′ that accounts for the autocorrelation of TKE′ within each of these products. Thus, the global circulation system as a whole has been accelerating substantially since the last decade of the 20th century.

To examine the spatial pattern of the acceleration of the global ocean circulation, we calculated the linear trend of the upper 2000-m mean KE′ for the period of 1991–2011 using the ensemble mean of the European Centre for Medium-Range Weather Forecasts (ECMWF) Ocean Analysis/Reanalysis System 4 (ORA-S4), ORA-S3, ECCO2, Global Ocean Data Assimilation System (GODAS), and Geophysical Fluid Dynamics Laboratory Ensemble Coupled Data Assimilation (GFDL ECDA) (Fig. 2A). Most of the ocean circulation systems are shown to have accelerated (Fig. 2A). About 76% of the upper 2000 m in the global ocean shows an increasing trend in KE′, with 28% showing significant increases at the 95% confidence level.

Fig. 2

Horizontal and vertical distribution of long-term KE change.

(A) Linear trend of oceanic KE averaged over the upper 2000-m layer (shaded color, unit in 103 J m−2 decade−1) during 1991–2011 from the ensemble mean of the ECMWF ORA-S4, ORA-S3, ECCO2, GODAS, and GFDL ECDA. Area where statistical significance is above the 99% confidence level is highlighted by black dots. (B) Vertical distribution of linear trend (blue line) of ensemble and global mean KE′ and the linear trend in percentage (embedded red line, relative to the climatological KE at each depth). Black dots indicate the 99% confidence level, and shaded area denotes the error in the trend. (C) Time-depth plot of ensemble mean KE′ (color shading and black contour lines) integrated over the global ocean. KE′ is calculated by subtracting the time means over 1991–2011 from the monthly series and low-pass filtering with a cutoff period of 25 months. (D) As in (C), but for the AGVA during 2005–2010.

The increasing trend in TKE′ is mainly attributed to an increase in KE′ in the tropics (Fig. 2A). Figure S2 shows the vertically integrated KE′ (5 to 1868 m) and its trend within the tropical zone (averaged over 20°S to 20°N), as well as meridional sections in the Pacific Ocean (180°E), Atlantic Ocean (30°W), and Indian Ocean (90°E). In all three ocean basins, there is an acceleration of the mean circulation at most latitudes (fig. S2, B to D). The acceleration is particularly prominent in the tropical oceans, especially the tropical Pacific Ocean, where a much stronger increase in KE′ is seen than higher latitudes. This indicates a substantial enhancement of the tropical current system (Fig. 2 and fig. S2). The Southern Ocean also shows positive trends but is marked by an alternating positive-negative zonal pattern (fig. S2, B to D), indicative of the shifting, or broadening and narrowing, of the Antarctic Circumpolar Current in the Southern Ocean (Fig. 2A). Separately, these results have been supported by previous observational studies in these regional ocean basins (, –).

The acceleration of global circulation extends into the deep oceans on a planetary scale. Figure 2 (C and D) shows the temporal evolution of the vertical structure of global ocean horizontal mean KE′ (derived from the ensemble mean of the reanalysis products separately for the AGVA data product that is based on in situ observations), while Fig. 2B illustrates the vertical distribution of ensemble mean KE trend. A notable feature is that the positive KE′ has increased in the upper 2000 m at a rate of the order of 1013 J m−1 decade−1, and over time, the increased KE′ extends deeper into the water column at a rate in the range of 1011 to 1012 J m−1 decade−1, which is significant at the 99% confidence level. A clear extension of positive KE′ into the abyssal oceans deeper than 1000 m was observed after the mid-1990s, reaching as deep as 3000 m in the Southern Ocean (fig. S3, D and H). Horizontal mean KE′ profiles in individual ocean basins and latitudes show similar features (fig. S3). The increasing trend relative to the mean KE is of the same order from the sea surface to the bottom with a rate of about 9% decade−1 of the climatological mean at each depth (Fig. 2B). The horizontal distribution of long-term KE trends in individual products is depicted in fig. S4, and time-depth plots of global mean KE′ in individual products are presented in fig. S5. A comparison between individual products suggests that the trends and evolution of TKE from individual products are similar to the ensemble mean.

Variance of KE′ in the Atlantic Ocean, Indian Ocean, and Southern Ocean is dominated by a statistically significant trend, but the Pacific Ocean shows interdecadal variability (fig. S3). Variability in the global mean KE′ over the upper 200-m layer at interdecadal (fig. S3) and interannual time scales (Fig. 2C) rides on the long-term increasing trend. Result from the AGVA, which is available for a short period between 2005 and 2010, also points to a deep-reaching increasing trend of the global mean KE′ (Fig. 2D). Negative phases of multidecadal variability occurred during the 1960s and 1980s–1990s, while positive phases occurred during the 1970s and after the late 1990s in the Pacific Ocean (fig. S3). The global ocean KE′ shares some of the interdecadal fluctuation of the Pacific Ocean KE′ (fig. S3, A and B), suggesting that the Pacific Ocean might be a main source of the interdecadal oscillation in the global KE. However, the increasing trend during the past two decades is much greater than the natural variability and commences from the early 1990s, which does not match the late 1990 or early 2000 onset of the Pacific decadal oscillation (PDO) (Fig. 1G).

A change in the observing system occurred in the early 2000s when broadscale upper ocean profile observations, mainly obtained using the expendable bathythermographs (XBTs), were superseded by the Argo era. A question arises as to whether the increasing density and frequency of Argo profiles influence the TKE trend. Assimilation of Argo salinity is found to cause a decrease in halosteric height after the introduction of Argo (). However, the time series of TKE′ for the global upper 700-m ocean (hence excluding the new observations below 700 m provided by Argo) shows a similar variation to the multidecadal trends when compared to the full-depth integral (fig. S6). In particular, the multidecadal increasing trend in TKE′ begins in the early 1990s, rather than the 2000s when the Argo observation system began. The robustness of the TKE trend in the reanalysis products is further discussed in Supplementary Text.

Another obvious change in the observing system is the introduction of precision altimeter data that began to be assimilated into many of the reanalysis products in 1992. However, the altimeter data have continued to be assimilated into the reanalysis products in the past two decades over which the trend in TKE has continued. Therefore, the change of observational platforms may not be responsible for the TKE increase.

Long-term trends of TKE are also investigated with simulations from the Coupled Model Intercomparison Project phase 5 (CMIP5) historical experiments. In the period before 1990, 8 of the available 12 CMIP5 models show a decreasing trend of TKE′. After 1990, 8 of the 12 models, as well as the ensemble mean, suggest an increasing trend of TKE′, but this trend is not statistically significant (fig. S1). We further discuss the differences between the reanalysis products and CMIP5 experiments in Discussion.

Energy sources of the acceleration of the global mean ocean circulation

Wind forcing is the major source of the mechanical energy of the global ocean circulation (, ). Analysis using six relatively independent products (see datasets and processing in Materials and Methods) confirms the role of wind forcing in the change of the energetics. All the products show a remarkable increase in the global mean sea surface wind speed over the past two decades. The linear trend of the global mean sea surface wind speed is 0.12 ± 0.09 m s−1 decade−1 since 1990, i.e., 1.9% decade−1 relative to the climatological mean of about 6.2 m s−1 (Fig. 3A).

Fig. 3

Global mean sea surface wind speed and wind work.

(A) Global mean wind speed (m s−1) at 10 m from various wind products. Each product is presented with a 13-month smoothed monthly time series superimposed with a 25-month smoothed time series. The thick blue line is the ensemble mean of the six wind products. The thick black line denotes the linear trend of the ensemble mean wind speed over the overlapping period of the six wind speed products (1985–2010) at the 99% confidence level. (B) Distribution of linear trend in wind work (10−3 W m−2 decade−1) during 1991–2011 from the ECMWF ORA-S4 assimilation. Area with confidence level higher than the 99% level is highlighted by black dots. (C) Comparison between globally integrated wind work (red) and ORA-S4 TKE′ (blue). All the time series are low pass–filtered (13-month running mean). Wind work is calculated using wind stress and 5-m current products from ECMWF ORA-S4 (see Materials and Methods for details).

To determine the contribution of wind forcing to acceleration, we calculate the wind work and its linear trend over the global ocean (Fig. 3B). The global mean wind work shows a statistically significant trend of about 2.04 × 10−4 W m−2 decade−1 during 1959–1989 (i.e., about 7.2 × 1010 W decade−1 for a total ocean surface area of 3.5 × 1014 m2) and 4.56 × 1014 W m−2 decade−1 (i.e., about 1.60 × 1011 W decade−1 globally) during 1990–2011. That is, the wind work contributes to an increase in energy input by 1020 J in two decades, which is much larger than the TKE′ increase (Fig. 1). The detrended global mean wind work leads detrended TKE′ by about 2 months with a correlation coefficient of 0.64 that is significant at the 99% confidence level (Fig. 3C).

Over most of the ocean basins, the spatial pattern of increased wind work is generally consistent with that of KE′ (recall Figs. 2A and 3B). The increasing wind work in the global tropical oceans, especially the tropical Pacific Ocean, is much stronger compared to extratropical regions. Wind work in the Southern Ocean overall shows a positive trend, but this is most likely due to a similar meridional shift as seen in the trend in KE′ shown in Fig. 2A. Differences exist in spatial pattern between the wind work trend and the TKE trend, but this is reasonable because the surface wind energy input is redistributed by oceanic processes such as horizontal advection, oceanic waves, and mesoscale eddies. Overall, these results indicate that the increased energy input by wind forcing of O(1020 J decade−1) is a major source of the TKE increase and has contributed in accelerating the global ocean circulation system over the past two decades (Figs. 1 and 3C). The wind-TKE′ relationship is also examined by comparing the multiproduct ensemble mean wind speed and TKE′ (fig. S7). Detrended global mean wind speed leads detrended TKE′ by about 5 months with a statistically significant correlation coefficient of 0.55 at the 99% confidence level. The time lag between the global mean TKE′ and wind speed is a result of the superposition of rapid barotropic and slow baroclinic responses of ocean circulation to surface wind forcing ().

Enhanced energy input due to stronger wind work not only directly contributes to the TKE′ increase but also contributes to TKE′ trends through C(APE, KE), i.e., the conversion from APE to KE (fig. S7). We estimated the global means of C(APE, KE) and absolute vertical velocity with the ECMWF ORA-S4 and compared them with the wind work. With the enhancement of wind work since the 1990s, the absolute vertical velocity is increased, the C(APE, KE) process is amplified, and the ocean currents obtain more energy from APE (fig. S7). Subsequently, this additional APE has also acted to increase the intensity of ocean currents since the 1990s. The increase in C(APE, KE) is of O(1020 J decade−1), which is the same as that of the wind energy input. The dissipation of KE due to molecular and small-scale diffusion of momentum is difficult to determine using available data; however, it is expected to also play an important role in the energy balance during the acceleration of global ocean circulation ().

PDO influence versus the long-term global warming trend

Previous studies have suggested that global warming weakens the ocean circulation system, especially in the global tropics. This is largely because the Walker circulation, especially the branch over the Pacific Ocean, would weaken in response to a faster warming in the eastern equatorial Pacific than in surrounding regions (). However, here, we find from multisource datasets that, at least for the past two decades since the 1990s, the global mean ocean circulation has accelerated because of an enhanced wind work input, particularly the tropical ocean circulation.

Both the TKE′ and Pacific KE′ show clear interdecadal fluctuations over the past few decades. Recent multidecadal enhancement of global wind forcing and ocean circulation was likely contributed by the phase change of interdecadal variability in the Pacific Ocean (, ). A recent study based on multiple independent observational records, including satellite measurements, indicates that there is a common strengthening trend of the Pacific Walker circulation over recent decades and that internal variability plays a dominant role in the strengthening trend of the Pacific Walker circulation (). During the recent negative PDO phase, wind work and wind speed strengthen over the tropical oceans and the Southern Ocean compared to during the positive PDO phase (fig. S8), contributing to acceleration of the global tropical ocean circulation.

A comparison between the PDO index, global mean wind speed, and TKE suggests that there is no significant relationship between the PDO index and the global mean wind speed or between the PDO index and TKE (fig. S8). However, the detrended TKE and PDO index have a correlation coefficient of −0.5, which is significant at the 99% confidence level. Thus, the multidecadal change in wind field and global ocean circulation is most likely a combined result of both a longer-term trend related to increased emissions of greenhouse gases and natural decadal variability. A comparison between ensemble mean TKE′ and its detrended component suggests that the SD of the detrended ensemble mean TKE′ is only about one-third of the increasing trend since 1990 (fig. S8). Further, although CMIP5 models project a weakening of the Pacific Walker Circulation (, ), the global mean sea surface wind shows an increasing trend that is significant at the 99% confidence level in both the historical and Representative Concentration Pathway (RCP) 8.5 runs of the CMIP5 models (fig. S9). The ensemble mean of global mean wind speed from the RCP8.5 runs shows a much stronger increasing trend than the historical runs (fig. S9). Therefore, the remarkable increasing trend of TKE′ during the past two decades is likely a part of the long-term trend that is emerging.

DISCUSSION

We have found a strong acceleration in the global mean ocean circulation over the past two decades. The acceleration is deep-reaching and particularly prominent in the global tropical oceans and can be attributed to the planetary intensification of surface winds since the 1990s, as a result of superimposition of a dominant longer-term trend and a relatively small contribution from natural variability.

We compared the global mean wind speeds from CMIP5 historical runs with that from six reanalysis products of wind over the same period (fig. S9). We find that, during their common period 1985–2005, the reanalysis winds show clear and consistent increasing trends, and the ensemble mean of reanalysis winds has an increasing trend of about 0.075 m s−1 decade−1 that is significant above the 99% confidence level. By contrast, the global mean wind speed in the ensemble mean of CMIP5 historical runs shows a very weak trend, about 0.006 m s−1 decade−1, much weaker than the wind trend in reanalysis products. Therefore, trends of TKE in the CMIP5 historical runs are much weaker than the reanalysis products. The reason for the discrepancy between CMIP5 historical run and reanalysis products is not well understood.

However, the ensemble mean of global mean wind speed from the RCP8.5 displays a significant increasing trend in the future, indicating that an increase in greenhouse gas emissions will eventually lead to acceleration of global mean wind and hence global mean ocean circulation (fig. S9). Thus, although the observed trend during the past two decades may also be driven by natural variability, for example, decadal variability associated with the PDO, the increase in greenhouse gas emissions is projected to eventually generate what is observed over the past two decades.

Nonetheless, uncertainties remain primarily since long-term direct observations of ocean currents are lacking for a direct verification. Although the increasing trend of TKE during the past two decades is a common feature of the reanalysis and data products we used, differences in spatial resolutions and spatial coverages give rise to a difference in TKE and their trend magnitude. In addition, the available in situ observations are mainly confined to the upper 2000-m layer, and, as shown in our results, the data-void abyssal ocean is likely to be important. Thus, intensive observations that monitor the deep global ocean circulation are urgently needed not only for understanding past conditions but also for reducing uncertainty in future projections of the global ocean circulation.

MATERIALS AND METHODS

Reanalyses and data products

1) The AGVA absolute geostrophic velocities of the quasi-global with a horizontal resolution of 1° and 29 regular pressure levels spanning the upper 2000 db were derived from quality controlled Argo profiling data (). Geostrophic currents within the equatorial region between 5°N and 5°S were discarded.

2) ECMWF ORA-S3 () datasets including monthly temperature, zonal/meridional current, and wind stress data were used. The monthly ECMWF ORA-S3 data cover the global ocean spanning 53 years from January 1959 to December 2011, with 29 vertical layers within 5 to 5250 m and a horizontal resolution of 1.4° longitude and 0.9° latitude on average.

3) ECMWF ORA-S4 () replaced the assimilation method with the variational assimilation system NEMOVAR-based Nucleus for European Modeling of the Ocean. The monthly ORA-S4 has 42 vertical layers within 5 to 5350 m and a horizontal resolution of 1° by 1°.

4) The GODAS is based on the GFDL MOM.v3 with a resolution of 1° longitude and 1/3° latitude and vertically 40 levels between 5 and 4478 m since January 1980 (). The model domain extends from 75°S to 65°N. The GODAS is forced by National Centers for Environmental Prediction (NCEP) atmospheric reanalysis 2 and is relaxed to weekly analyses of sea surface temperature and annual climatological salinity. The GODAS assimilates temperature profiles from XBTs; from moorings of the Tropical Atmosphere Ocean, the Triangle Trans Ocean Buoy Network, and the Prediction and Research Moored Array in the Atlantic; and from Argo profiling floats. The known fall-rate error of the XBT data for the years before 1990 has been corrected in the current version (www.esrl.noaa.gov/psd/).

5) The JPL ECCO2 () from NASA’s Jet Propulsion Laboratory (http://ecco2.jpl.nasa.gov/). The ECCO2 assimilates the available satellite and in situ data based on the Massachusetts Institute of Technology general circulation model with Green’s function optimization and outputs eddy-resolving ocean current field. Here, both monthly and 3-day ECCO2 datasets were applied. The ECCO2 has a resolution of 0.25° by 0.25° by 50 vertical levels within 5 to 5906 m and covers a period of 1992–2013.

6) The ECDA v3.1 from Ocean Data Assimilation Experiment was provided by National Oceanic and Atmospheric Administration (NOAA)’s GFDL ECDA () (www.gfdl.noaa.gov/ocean-data-assimilation). The GFDL ECDA v3.1 covers 1961–2012 and is based on the fully coupled climate model with a zonal resolution of 1°, meridional resolutions increasing from 1° to 0.3° in the equatorial region, and 50 vertical levels within 5 to 5316 m.

7) Numerical simulations from the CMIP5 () forced by the historical (1950–2005) radiative forcing were also used to calculate the trend of KE. These models include the following: the Beijing Normal University Earth System Model (BNU-ESM) (monthly, 50 vertical layers within 10 to 5579 m, monthly, 1° longitude and 0.96° latitude on average), the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Mk 3.6.0 (monthly, 31 vertical layers within 5 to 4800 m, 1.88° longitude and 0.93° latitude on average, two members), the CanESM2 (monthly, 40 vertical layers within 5 to 5233 m, monthly, 1.41° longitude and 0.94° latitude on average), the Flexible Global Ocean-Atmosphere-Land System Model (FGOALS) s2 (monthly, 30 vertical layers within 5 to 5243 m, monthly, 1° longitude and 0.8615° latitude on average), the GISS E2 R/R-CC (monthly, 32 vertical layers within 6 to 4887 m, monthly, 1.25° longitude and 1° latitude), the GISS E2 H-CC (monthly, 33 vertical layers within 0 to 5500 m, monthly, 1° longitude and 1° latitude), the HadCM3 (monthly, 20 vertical layers within 5 to 5192 m, monthly, 1.25° by 1.25°), the MIROC ESM/ESM CHEM (monthly, 44 vertical layers within 2.5 to 5450 m, monthly, 1.41° longitude and 0.93° latitude on average), and the HadGEM2 AO (monthly, 40 vertical layers within 5 to 5327 m, monthly, 1° longitude and 0.84° latitude on average).

8) Wind products used in this study include the Objectively Analyzed air-sea Fluxes (OAFlux), Wave- and Anemometer-based Sea Surface Wind (WASWind), TropFlux, Twentieth Century Reanalysis (20CR), NCEP Reanalysis 1 and NCEP Reanalysis 2, and the ECMWF ORA-S3 wind stresses. Of these wind products, the TropFlux covers only the global tropical oceans (30°S to 30°N), while the others have global or quasi-global coverage. The WASWind combined ship winds and wind wave height, corrected the trend bias due to increases in anemometer height, and represents reasonable trend relative to independent sea level pressure observations (). The OAFlux uses surface meteorological fields derived from satellite remote sensing and reanalysis output from NCEP and ECMWF models (). The TropFlux were computed using the bias and amplitude-corrected ECMWF interim reanalysis (ERA-I) input data plus climatological gustiness correction and evaluated against dependent and independent mooring data. The TropFlux successfully corrected the bias induced by rain under the Intertropical Convergence Zones and is found to display the best agreement with in situ data in the global tropics compared to the NCEP, NECP2, and QuikSCAT wind products (). The four reanalysis products, the 20CR (), NCEP R1 (), NCEP R2 (), and ECMWF ORA winds, are assimilations based on historical observations using different data assimilation systems. In addition, wind speed from historical (1850–2005) and RCP8.5 (2006–2100) runs of 21 CMIP5 models was used as well (fig. S9).

Calculation of KE

The KE of a water particle at time t is defined aswhere M is the mass, and u, v, and w are zonal, meridional, and vertical velocities, respectively, and x, y, and z are zonal, meridional, and vertical coordinates. In the ocean, w is usually less than 1% of the horizontal velocity (), so the KE of a water particle is reduced to

Then, the globally integrated TKE is defined as

We used a constant mean density of seawater ρ0 assumed to be 1025 kg m−3 and integrated KE over the global ocean and from surface to the bottom of the spatial grids of each corresponding dataset, giving rise to

Anomaly of TKE (i.e., TKE′) is hence defined as

The Argo-based observational data AGVA, five assimilation products, and 12 model simulations were used to calculate KE. ECMWF ORA-S4, ORA-S3, GFEL ECDA, GODAS, and ECCO2 were used to calculate an ensemble mean TKE that covers the past two decades (Figs. 1 and 2 and fig. S8B). Figures S2, S3, and S6 used ECMWF ORA-S4, because the ECMWF ORA-S4 has the longest time coverage, which is important in examining the multidecadal variability and also because major variables (e.g., ocean current, temperature, salinity, and surface fluxes) from the ECMWF ORA-S4 were available for us.

Wind work

The work performed by sea surface wind to the global oceanic general circulation is calculated as (): ∬τw ∙ vsdxdy, where τw is monthly wind stress and vs is surface current velocity. Here, we calculated the wind work applying wind stress from the ECMWF and WASWind and 5-m-depth current data from the ORA-S4.

Sources and sinks of and energy conservation equation

The KE sources and sinks that determine KE variability can be expressed by the following energy conservation equation as ()where the terms on the right-hand side respectively represent wind work, bottom friction term (τb is bottom friction and vb is near bottom velocity), conversion from APE to KE, dissipation of KE due to molecular and small-scale processes (ED), and residual term (Eres). Among these terms, the wind work and C(APE, KE) are the major sources of KE change. Conversion from APE to KE is measured by C(APE, KE), which is defined as C(APE, KE) = − g ∭ w(ρ − ρ0)dxdydz, where ρ0 is basin-scale mean density.

An increasing trend in KE of the ocean implies increased energy generation from surface forcing (e.g., wind forcing and surface heat flux) into ocean circulation. Internal energy conversion is also a potential source to increase KE. Some of the reanalysis products we used do not necessarily conserve energy due to the so-called assimilation increment, which is a statistical correction applied during sequential data assimilation time steps to correct the model state to fit the observations. ORA-S3, ORA-S4, ECDA, and GODAS do not conserve energy, and in these models, there is an additional energy source for the production of KE.

Calculation of linear trend

The linear trend was estimated by applying a one-dimensional linear regression model in the least-squares sense. All the tendency results are statistically tested by a modified Mann-Kendall test (). The Mann-Kendall test is a commonly used nonparametric trend test. Hamed and Rao () improved the Mann-Kendall test to make it suitable for autocorrelated data, and the accuracy of the modified Mann-Kendall test in terms of empirical significance level is superior to that of the original trend test without any loss of power. According to the modified Mann-Kendall test, the variance of S for autocorrelated data can be calculated by , where n is the actual number of “observations” and represents a correction related to autocorrelation in the time series . C(i) denotes the autocorrelation between the ranks of the observations ().

Statistical analysis

All the tendency results were statistically tested using a modified Mann-Kendall test (). The confidence level of trend estimated using the modified Mann-Kendall trend test infers an integral time scale, representing the autocorrelation of the time series, to produce the effective degrees of freedom used in the trend test. These statistics of the TKE′ time series calculated from various data products are shown in table S1. The effective degree of freedom is less than the sample size due to autocorrelation in the time series. Assuming that the TKE′ satisfies the t distribution, the t value and the correlation coefficient with time at the significance level of 0.05 [t(α = 0.05), i.e., γα(α = 0.05)] are also presented in table S1. The correlation coefficients between TKE′ and time γ (correlation coefficient between time series and time) are shown in the last row of table S1. The γ for various data products is greater than the γα(α = 0.05), suggesting that the TKE′ trends are significant at the 95% confidence level. However, this statistical analysis does not account for uncertainties in the TKE′ trend that are due to the lack of global oceanic observations contributing to the TKE′ estimate itself.