WP data and calculation

Wind-generated waves modulate air-sea exchanges, dissipating energy, and passing momentum on to ocean turbulence and currents59. A fraction of this energy is employed for wave motion28, which results in a transport of energy by waves, known as WP. For irregular waves, the WP can be obtained from the spectral energy density function, S(f,θ), where f represents the frequency and θ represents the direction of waves60:

$$WP = \rho g{\int\!\!\!\!\!\int} {c_g} \left( {f,h} \right) \cdot S\left( {f,\theta } \right)dfd\theta$$

where c g represents the group velocity \(c_g = \frac{1}{2}\left( {1 + \frac{{2kh}}{{{\mathrm{sinh}}(2kh)}}} \right)\frac{L}{T}\), L represents the wavelength, k is the wavenumber 2π/L, h represents the water depth, and T represents the wave period. L and T are related through the dispersion equation as follows:

$$\left( {\frac{{2\pi }}{T}} \right)^2 = gk{\mathrm{tan}}h\left( {kh} \right)$$

Sustained conditions of irregular waves (i.e., sea state) are usually described by certain parameters that characterize the spectral shape and can be calculated from it. The wave height associated with the standard deviation of the surface elevation, H m0 , can be computed from the integral of the spectral density or the order-zero moment, m 0 , as follows:

$$H_{m0} = 4.004 \cdot \sqrt {m_0}$$

For the swell sea states, H m0 is comparable to the significant wave height, H S , which represents the mean value of the largest third of the wave heights in the sea state. Other variables defining a sea state are the peak or mean periods (T 02 or T 01 , respectively).

The WP for an irregular sea state can be obtained from wave spectral parameters using the following expression60:

$${\mathrm{WP}} = \frac{{\rho g^2}}{{64\pi }} \cdot T_e \cdot \left( {H_s} \right)^2$$

where T e or T −10 is known as the energy period. This parameter can be estimated from the spectral shape and other parameters. The mean period, T 01 , and T e can be related through the following relationship:

$$T_e = \frac{{m_{ - 1}}}{{m_0}} = \alpha T_{01}$$

where α depends on the spectral shape. In this work α is taken to be 0.538, which corresponds to a mean peak-enhancement factor of 3.3 on the JONSWAP spectrum14. The assumption of a constant spectral shape introduces some uncertainty into the WP estimates, but the effect is negligible because the error when estimating α is an order of magnitude smaller than that for the effects of T 01 and H S ; furthermore, the periods have less influence than H S in the WP equation. For the same reasons, correctly modeling H S is critical for an accurate assessment of the global WP.

Information on wave heights collected using buoys and satellite altimetry do not provide continuous data over space and time, and so require numerical climate reconstructions in order to study historical climate states61,62. We use numerical wave models, combined and validated with instrumental sources, to describe the global wave climate across different time periods14.

The GWP was calculated hourly for time series of significant wave heights (H S ) and mean wave periods (T S )14. WP is expressed in terms of kilowatts per meter (kw m−1) along the wave front. Global wave data containing H S at hourly intervals, the mean period and the mean directions for the period 1948–2008 were obtained from the GOW reanalysis29. The GOW reanalysis used the WAVEWATCH III model63 with NCEP/NCAR global wind and ice cover datasets61. The simulated wave parameters are later calibrated with satellite altimetry from the period 1992–2008 through statistical and directional corrections64. Hurricanes and typhoons in the satellite data are identified applying an outlier filter based on weighted least squares at a 95% significance levels and not considered in the calibration65. The calibrated data reproduce the spatiotemporal variability of the global wave climate, as demonstrated by the comparison with wave buoys and satellite altimetry, both in terms of wave heights and WP14,29.

The global long-term wave data from GOW were cross-checked with two other data sources to corroborate the global time series of GWP in terms of magnitude and temporal variability and extend the time series beyond the year 2008. First, we calculated the GWP with wave data from an independent and more recent high-resolution global reanalysis (RaA13) that covers the period 1994–2012 and uses improved parameterizations for the model WAVEWATCH III and high-resolution wind forcing (the Climate Forecast System Reanalysis, CFSR)30. RaA13 implemented and validated the improved parameterizations of Ardhuin et al.66 and found no evidence of bias in wave heights when using unbiased winds after 1994, even for high significant wave height values. Another GWP time series up to the year 2017 is calculated from the global reanalysis GOW-CFSR, which is an updated version of GOW with a resolution of 0.5 degrees globally and uses the model parameterizations in RaA13 and the CFSR wind forcing31. Additionally, satellite altimetry data from 1992–2008 (the period coinciding with that of GOW) were used to calculate the significant wave heights via satellite altimetry, and the periods were interpolated in time and space from GOW to provide a third time series of GWP. The altimeter data were produced and distributed by AVISO (http://www.aviso.ocanobs.com/)

However, CFSR winds and waves show large errors over the Southern Ocean until the 1994 compared to altimeter data30,67. Chawla et al. explained this high bias by an increase in extreme values of the wind speed67 but Ardhuin et al. also found that even the mean values are strongly biased for latitudes south of 30ºS30,68, the most energetic region of all the basins (Fig. 1). The high bias is corrected for wave data after 1994 when the CFSR reanalysis started to assimilate wind data from the Special Sensor Microwave Imager (SSM-I)62. It is possible that a correction of the CFSR wind speed histogram may be enough to correct the biases on wave parameters, but both the AaR2013 and GOW-CFSR are uncorrected before the year 1994, which induces large overestimation of WP compared to the satellite. Therefore, GOW-CFSR and AaR2013 data are only considered after the year 1994. The GOW-CFSR climatology has been validated with satellite data in the latter period showing a good agreement for mean significant wave heights with some bias in the high latitudes of the southern hemisphere31. Additionally, the three numerical datasets (GOW-NCEP, AaR2013, and GOW-CFSR) are compared with the signal of GWP calculated from satellite in the Supplementary Figure 1. The two high-resolution datasets provide a better fit than the long-term GOW-NCEP data, but both present different biases and scatters indices in terms of GWP and are therefore both used in the analysis as independent datasets.

Although numerically-generated wave data are useful to study global wave climatology, there are inhomogeneities caused by changes in the amount of assimilated observations within the forced wind fields throughout the historical time span69,70. The two numerical wave datasets used here may present regional discrepancies because they use different wind forcing data and numerical implementations; the GOW uses the NCEP/NCAR reanalysis61, while RaA13 uses the Climate Forecast System Reanalysis62 and has a higher spatial resolution (0.5° × 0.5° vs 1.5° × 1.0° in GOW). However, the discrepancies between each wind forcing dataset are minimal in the extratropical storm belts. The climatologies of wind direction and wind speed, as described by the cross-calibration and multiplatform assimilation of ocean surface wind data71,72, show that wind stress from the NCEP/NCAR reanalysis and climate forcing from the GOW reanalysis only differ from more detailed climate reanalyses in tropical zones (e.g., ERA40 or ERA-Interim). Only in the extratropical storm belts do both wind datasets show similar biases with respect to wind measurements73,74. However, the effect of wave generation on tropical zones is negligible in comparison with the effects on high-latitude generation areas45.

The NCEP/NCAR reanalysis assimilation system61 was unchanged during the reanalysis period to eliminate climate jumps, although the reanalysis is still affected by changes in the observing data70,75. There were two changes in the observing data: before and after 1957 (eras I and II) when the upper-air network was being established and that mainly affected the southern hemisphere;62 and the introduction of the global operational use of satellite soundings in 1979. It has been therefore recommended that the periods before and after 1979 are studied independently because the climatology based on the years 1979–present day is most reliable75. Following this recommendation, our analysis separates the full historical period (1948 onwards) and the satellite era (1979 onwards). The results for the eras not included in the main manuscript can be found in the Supplementary Tables 2 to 6.

Sea surface temperature

SST is the water temperature close to the surface of the ocean, usually up to 20 m below the surface of the ocean76. SST measurements became more comprehensive and diverse after 195027. Ships and buoys have been recording SST, among many other parameters, for well over 100 years, but after 1967, satellites began to remotely measure SST, with the first global SST composite created in 197077. More widespread satellite measurements of SST since 1982 have allowed for more in-depth exploration of temporal and spatial variations and enabled deeper insights into ocean atmosphere connections.

The SST anomalies used in this analysis were obtained from two sources. We used the ERSSTv3b dataset to define the global signal of upper ocean warming. The ERSSTv3b is the most recent version of the extended reconstructed SST analysis developed by NOAA35,78. The dataset is based on the International Comprehensive Ocean-Atmosphere Data Set (ICOADS). ERSST v3b does not use satellite-based SST data but in situ measurements and improved statistical methods to reconstruct temperature fields from sparse data. It provides information on a global 2º grid from 1854 to the present. The anomalies are expressed in ºC and computed with respect to the 1971–2000 monthly climatology35. We also use the Optimal Interpolation Sea Surface Temperature (OISST) dataset, which offers higher resolution than ERSSTv3b, providing data on a 0.25º global grid after 1981 and combines observations from different platforms (satellites, ships, buoys)36. This dataset is used to validate the ERSSTv3b time series and serve as a second and independent SST dataset to correlate with GWP. The global time series are calculated from monthly SST fields and downloaded from NOAA’s National Centers for Environmental Information at: https://www.ncdc.noaa.gov/data-access/marineocean-data/extended-reconstructed-sea-surface-temperature-ersst-v3b; and https://data.nodc.noaa.gov/cgi-bin/iso?id=gov.noaa.ncdc:C00844

Calculation of the global and regional time series

The hourly time series of WP were aggregated by seasons and years to calculate the regional and globally averaged time series. The global signals for both WP and SST anomalies were obtained through spatial averaging as follows:

$$Z{\mathrm{global}} = \frac{{\mathop {\sum }

olimits_i {\mathrm{cos}}({\mathrm{latitude}}_i) \cdot Z(t)_i}}{{{\mathrm{Number}}\,{\mathrm{of}}\,{\mathrm{nodes}}}}$$

where Z(t) i represents the annual or seasonal time series of the respective variable at each location in the corresponding grid.

The time series are also calculated by ocean basins and climatic regions. The ocean basins are defined by the continental limits in the Pacific, Atlantic, and Indian Oceans, and the 40 degrees south limit. The Southern Ocean is taken to be south of the 40 degrees latitudinal limit based on WP climatology and storm activity in the southern hemisphere9,13,14,33 . The regions where we calculated the regional time series corresponded to the wind-wave generation zones and were defined for each ocean sub-basin and a 30° latitudinal limit (see Fig. 8) following Alves (2006)45. These oceanic areas have been identified by the propagation footprint of waves45 and the classification of wave climate types24,46.

Calculation of long-term trends

The long-term trends were calculated for each global time series through a linear regression. The significance of the trends were checked with the Mann–Kendall (MK) test32. The purpose of this test was to statistically assess whether or not there was a monotonic upward or downward trend over time by testing if the slope of the estimated linear regression line differed from zero. The MK test is a nonparametric (in other words distribution-free) test and does not require the residuals from the fitted regression line to be normally distributed, as a parametric linear regression does.

Nevertheless, the standard P-values obtained from the MK test are based on the assumption of independence between the observations. This makes it important to check the autocorrelation in a given series and adjust the MK test if necessary. To avoid autocorrelation in the time series, we followed the approach of Wang and Swail8. The time series in question is pre-whitened (processed to make it behave statistically like white noise), and the trend is estimated by fitting with the model:

$$Y_t = a + bt + X_t$$

where a and b represent the intercept and slope, respectively, and X t is given by:

$$X_t = cX_t + {\mathrm{\varepsilon }}_t$$

The pre-whitened time series that possesses the same trend as that of the original time series (Y t ) is8:

$$W_t = \frac{{Y_t - cY_{t - 1}}}{{1 - c}}$$

To estimate W t , we follow the iterative approach presented by Wang & Swail8 as follows:

Step 1. Initial estimate of c:

The lag-1 autocorrelation of the time series Y t is taken as the first estimate of c (i.e., c 0 ).

If c 0 is less than 0.05, the effect of serial correlation is negligible, and no iteration is necessary when applying the MK test to the original time series.

Otherwise we perform the described trend analysis on the pre-whitened time series \({{W}}_{{t}} = \frac{{{{Y}}_{{t}} - {{c}}_0{{Y}}_{{{t}} - 1}}}{{1 - {{c}}_0}}\)and obtain the first estimate of b (i.e., b 0 ).

Step 2. Calculation of b 1 and c 1 :

The estimated trend component from the original data series is removed, and c is re-estimated. The new estimate of c , c 1 , is the lag-1 autocorrelation of the time series ( Y t - b 0 * t ).

If c 1 is less than 0.05, we take c 0 and b 0 as the final estimates.

Otherwise, the trend analysis is performed on the pre-whitened times series \(W_t = \frac{{Y_t - c_1Y_{t - 1}}}{{1 - c_1}}\), and a new estimate of b is obtained (i.e., b 1 ).

Step 3. Iterations to estimate b and c:

While the differences of abs( c 1 - c 0 ) and abs( b 1 - b 0 ) are larger than 1%, steps 2 and 3 are repeated, which results in b 0 = b 1 and c 0 = c 1 .

Once the appropriate estimates of c and b are obtained, the MK test is applied to the pre-whitened time series \(W_t = \frac{{Y_t - cY_{t - 1}}}{{1 - c}}\)to conclude the trend analysis.

The implementation is accompanied with a Matlab function in the Supplementary Code (see also code availability statement). For the MK test on the pre-whitened series, we adopt the implementation approach of Burkey (2006)79.

Relating changes in SST to changes in WP

Correlations were computed between the global time series, by years and seasons, and spatially over each ocean sub-basin. The correlation was assessed through the Pearson product-moment correlation coefficient, r, which is a measure of the strength and direction of the linear relationship between two variables. The statistical significance was calculated through the Student’s t-test at the 95% confidence level.

Correlations were also calculated for the time series of the non-autocorrelated residuals. The non-autocorrelated residuals were obtained after adjusting the autoregressive moving-average models to each global and regional time series in order to avoid autocorrelation effects in the statistical analysis, and to identify the existence of non-contemporaneous relationships80. The non-autocorrelated residuals ignore the effects of trends and autoregression from the original time series and, thus represent the variability (at the temporal scale) of the original time series (yearly or seasonally). Given a time series, X t , we fitted an ARMA(m, n) model by:

$$X_t - ( \propto _1X_{t - 1} + \ldots \propto _mX_{t - m}) - (\theta _1\varepsilon _{t - 1} + \ldots \theta _n\varepsilon _{t - n}) = \varepsilon _t$$

where ∝ m is the parameter of the autoregressive part of the model, θ n is the parameter of the moving-average part of the model, and ε t is the error term. Non-autocorrelation of the residuals was statistically tested with the Ljung–Box test. Table 2 provides the parameters for the adjusted models for each dataset.

Table 2 Autoregressive models. Autoregressive moving-average models parameters used for each dataset Full size table

Using the global time series of GWP and SST anomalies, linear regressions were calculated between the time series, annually and monthly, as follows: \({\mathrm{GWP}}(t) = a + b \cdot {\mathrm{SST}}(t - \Delta )\), where Δ represents a temporal lag, which was found to be zero from a lagged-correlation analysis (Supplementary Figure 3-a).

However, previous studies that use semi-empirical relationships with temperatures to study the effect of global warming on other response variables such as sea level rise48 and hurricane activity38,81 have used rates of change, instead of the value of the temperature anomalies, as a descriptor. For this reason, we calculate the yearly changes in the GWP and SST and regress this annual variability as:

$$\frac{{\partial {\mathrm{GWP}}(t)}}{{\partial t}} = a + b\frac{{\partial {\mathrm{SST}}(t - \Delta )}}{{\partial t}}$$

Each regression model is shown in Fig. 5.

In addition, information theory was used to determine the degree of MI between the WP and SST anomaly. We follow the approach undertaken by Hoyos et al.38 to study the relationship between SST and hurricane intensity. In information theory, the MI of two variables is quantified to represent the measure of independence of the two variables37. MI quantifies the distance between the joint distribution of two variables and the product of their marginal distributions. Therefore, MI measures (in bits) the independence of two variables. MI is based on the concept of entropy, which is associated with the randomness of a variable. The entropy, H(X), of variable X for random event x that occurs with a probability of f(x) is37:

$$H\left( X \right) = \mathop {\sum }\limits_x f(x) \cdot {\mathrm{log}}_2f(x)$$

The joint entropy of two variables, X and Y, measures the entropy contained in the joint system and is defined as:

$$H\left( {X,Y} \right) = \mathop {\sum }\limits_{x,y} p(x,y) \cdot {\mathrm{log}}_2p(x,y)$$

MI is then defined as:

$${\mathrm{MI}}\left( {X,Y} \right) = \mathop {\sum }\limits_{x,y} p\left( {x,y} \right) \cdot {\mathrm{log}}_2\frac{{p\left( {x,y} \right)}}{{f\left( x \right)g\left( y \right)}} = H\left( X \right) + H\left( Y \right) - H(X,Y)$$

where f and g represent the marginal distributions of each variable. If X and Y are independent, the total entropy of the system, H(X,Y), would be equal to H(X) + H(Y), in which case the MI value would be zero, indicating that X does not contain information about Y, and vice versa. MI was calculated following Hoyos et al. and our implementation was validated against the theoretical cases expressed in the Supplementary Material therein. To quantify MI, we estimate the marginal distribution of the variables (Fig. 4a, b). The product of their marginal distributions (Fig. 4c) should replicate their joint distribution (Fig. 4d) if the two variables are independent, in other words, if GWP contains no information about SST and vice versa (i.e., MI = 0). However, the distributions in Fig. 4c, d differ, which implies that there is statistical dependence between the variables.

Nevertheless, the analysis of the original GWP and SSTs time series does not enable us to directly discern whether or not these statistical relationships arise from long-term trends or short-term modes of variability (in other words, on a decadal or shorter timescale). To remedy this, we performed the MI analysis on the isolated trend/variability time series following Hoyos et al.38. The trend is removed by subtracting the least-squares linear fits and the adjusted autoregressive moving-average models to calculate the non-autocorrelated residuals. The resulting non-autocorrelated residuals also have an MI different from zero (1.05 bits for entropies of 2.86 in the SST residuals and 2.95 in the WP residuals), indicating that WP is statistically dependent on SST variability.

Climate indices

Climate index data for Nino3 and AMO for Fig. 7 were obtained from NOAA. The AMO is a mode of natural variability that occurs in the North Atlantic Ocean, with basin-wide SST changes over a period of 60–80 years82. The Niño3 index83 registers SST anomalies in the tropical Pacific (90–150ºW, 5ºS-5ºN). Historically strong El Niño events were defined based on the Oceanic El Niño Index after 1950, which uses the same region as the Niño 3.4 index: (5°N-5°S, 170°-120°W). The Oceanic El Niño Index uses a 3-month running mean; to classify events as full-fledged El Niño or La-Niña occurrences, the anomalies must exceed + 0.5 °C or -0.5 °C for at least five consecutive months. This is the operational definition used by NOAA. The time series of each climate index (Fig. 7a, c) were correlated with each local time series of WP from GOW in the global grid to provide the patterns shown in Fig. 7b, d.

Code availability

The statistical significance test for trend calculation is described and provided in a Supplementary File and can be found at: https://osf.io/pvm6c/, with the identifier https://doi.org/10.17605/OSF.IO/PVM6C. Other code used to generate the results for this project are available from the corresponding author upon reasonable request.