Sea State Decadal Variability in the North Atlantic: A Review

Abstract

Long-term changes of wind-generated ocean waves have important consequences for marine engineering, coastal management, ship routing, and marine spatial planning. It is well-known that the multi-annual variability of wave parameters in the North Atlantic is tightly linked to natural fluctuations of the atmospheric circulation, such as the North Atlantic Oscillation. However, anthropogenic climate change is also expected to influence sea states over the long-term through the modification of atmospheric and ocean circulation and melting of sea ice. Due to the relatively short duration of historical sea state observations and the significant multi-decadal variability in the sea state signal, disentangling the anthropogenic signal from the natural variability is a challenging task. In this article, the literature on inter-annual to multi-decadal variability of sea states in the North Atlantic is reviewed using data from both observations and model reanalysis.

1. Introduction

It is now well-established that the sea state displays variability on inter-annual to decadal timescales. Various datasets based on observations or on numerical models have been used to study this variability. The advantages and disadvantages of each dataset will be reviewed in this article. This slow variability is due to the fluctuations of the atmospheric circulation, such as the North Atlantic Oscillation, that modifies the trajectories of extra tropical storms in the North Atlantic (NA, a glossary of acronyms is given in

Table 1. Glossary of acronyms used in this review.

WEF	Wave Energy Flux	TC	Tropical Cyclone
GCOS	Global Climate Observing System	VOS	Visual Observations of Sea state
NA	North Atlantic	ENSO	El Niño–Southern Oscillation
JFM	January February March	EOF	Empirical Orthogonal Function
SLP	Sea Level Pressure	PC	Principal Component
CCA	Canonical Correlation Analysis	EA	Eastern Atlantic
SCA	Scandinavian Pattern	EA/WR	East Atlantic/West Russia Pattern
PNA	Pacific-North American	TNH	Tropical North Hemisphere
SAM	Southern Annular Mode	WEPA	West Europe Pressure Anomaly
WW3	Wave Watch 3	NDBC	National Data Buoy Center
MEDS	Marine Environmental Data Section	CFSR	Climate Forecast System Reanalysis
SSM/I	Special Sensor Microwave/Imager	CCI	Climate Change Initiative

) and thus the associated sea state. The relationship between the surface atmospheric field and the sea state is, however, not straightforward, because sea states consist of locally (wind wave) and remotely (swell) forced waves. In particular, sea states at mid- and low latitudes in the Atlantic are strongly influenced by swell originating from higher latitudes, forced by extra tropical storms, and might also be influenced by swell originating as far as the Southern Ocean [1]. In addition, the same wind intensity might result in different wave amplitudes depending on the duration and area over which the wind is blowing. Therefore, the wind direction, the storm trajectory, and the presence of land or sea ice on the path of the storms all impact on the wave generation process. Moreover, local effects, such as shoaling, refraction, diffraction, and wave–current interactions may also modify the wave field along its propagation. Despite these difficulties, numerous studies, for example, [2,3,4] have established the possibility to statistically link the slow sea state variations to the slow atmospheric circulation variations, as we shall review in this article. Given that the sea state variations are tied to atmospheric variations, it is expected that changes in atmospheric circulation due to anthropogenic climate change will also lead to changes in the sea state climate.

The remainder of this manuscript is organized as follows. In Section 2, we begin by reviewing the different datasets used in the literature to study the inter-annual to decadal variability of sea states: visual and instrumental in situ observations, satellite altimetry, seismic data, and model reanalysis. In Section 3, we review our current understanding on the link between sea state and atmospheric variability by first reviewing the various statistical methods used to establish this link, and then by describing climate modes’ effects on the North Atlantic sea wave height variability. In Section 4 we review the current understanding of the projected change in sea states due to climate change. In Section 5, we summarise and conclude.

2. Observing Sea State Decadal Variability in the North Atlantic

In this section, we review the long-term sea state dataset that has been used to investigate decadal wave climate variability in the North Atlantic, and we address their advantages and limitations for this purpose. For a more exhaustive review on wave-observing systems, please see Ardhuin et al. [5].

2.1. Visual and Instrumental In Situ Observations

Early evidence of decadal variability in the wave climate was found in the North Atlantic and the North Sea in the late 1960s based on both visual and instrumental in situ observations of significant wave heights [6,7]. Visual observations of sea state (VOS) conditions reported every six hours from hundreds of ships cruising across the North Atlantic were used to produce synoptic wave charts, revealing large inter-annual variability and a general upward trend of the mean $H_s$ between 1970 and 1982 [8]. These results were instrumentally confirmed from the aggregated measurements of shipborne wave recorders installed on the weather ships regularly visiting Ocean Weather Stations in the North Atlantic [9]. Today, both visual and instrumental in situ observations remain extremely valuable sources of historical wave information required to investigate decadal variability in the wave climate. Visual observations from Voluntary Observing Ships have been collected since 1853 until today without changes in observational practice [10]. However, it is only since 1950 that both the wind sea and swell components have been systematically reported. One of the main limitations of visual observations is the coding precision, which is 0.5 m for wave heights (and 1 s for periods and 10 degrees for directions), well above the 0.1 m measurement uncertainty in the wave height requirement of the Global Climate Observing System (GCOS hereafter) for sea state Essential Climate Variables [11]. Other sources of uncertainties are related to incorrectly recorded and unprofessionally observed wave parameters or errors when digitizing old logbooks. After quality control procedures are applied, only 30% of the initially available data remains [10]. Another limitation of VOS data is the spatial sampling which is constrained to the major navigation routes across ocean basins. However, in the North Atlantic, given the dense marine traffic between the European and American continents, a large part of the basin is covered with sea state conditions reported approximately 10 to 50 times per month over 2 × 2° bins, and up to 350 times per month along the US and European coastlines (see VOS density map on Figure 1).

Regarding historical instrumental observations, most of the Ocean Weather Stations regularly visited by weather ships during the second half of the 20th century were discontinued during the 1980s, as unmanned weather buoys represented a much cheaper opportunity to continuously record met-oceanic conditions. Still, these stations represent the locations with the highest number of collocated instrumental and visual observations in the deep North Atlantic (NA hereafter) ocean (see Figure 1). Since the 1980s, in situ wave buoys have become the standards for recording sea state conditions, and some of the longest time-series now cover more than 40 years, with little interruption. Several authors used these data to investigate sea state decadal variability [3,12,13], although major issues remain when computing long-term trends. For instance, Gemmrich et al. [14] showed that modifications of the wave measurement hardware, as well as the analysis procedures resulted in step changes in the records of several wave buoys in the North-East Pacific, significantly impacting the magnitude and sometimes the directions of the trends. Moreover, moorings are regularly relocated after maintenance operations, sometimes several hundreds of km away from their nominal location, which may also impair the consistency of the long-term record. In the NA, we identified in the Copernicus Marine Service In Situ database (http://www.marineinsitu.eu/ (accessed on 10 October 2021)) 33 moored wave platforms with more than 25 years of data, with a minimum time coverage of 80% (red squares in Figure 1). These stations are mostly located in the Gulf of Mexico, along the US and Canadian east coast (operated by NDBC and MEDS), and along the Atlantic coast of Europe (operated by Met Office and Puerto del Estado).

2.2. Satellite Altimetry

Radar altimeters onboard satellite missions are certainly the most efficient way to record the significant wave height at a global scale. For instance, Woolf [15] used 10 years of altimeter data (Geosat, ERS-1, ERS-2, and TOPEX missions) to investigate the variability and predictability of the significant wave heights in the North Atlantic. In 2021, 30 years of uninterrupted multi-platform altimeter records will be available. This duration is a milestone for satellite altimetry as it represents the minimum duration required for computing climatological standards following the World Meteorological Organization recommendation [16]. In recent years, long-term multi-mission altimeter records have been analyzed to compute trends in mean and extreme $H_s$ over the altimetry era [17]. However, multi-mission altimeter records aggregate data from various platforms and instruments and may suffer from heterogeneities in the radar technologies and processing techniques, which have evolved over time, in the reference dataset used for calibration, and/or in the sampling characteristics of each mission (see [18] for a comparison of trends between different multi-platform altimetry products). Regarding the sampling issues, Jiang [19] combined the track information of altimeter missions with model reanalysis outputs to investigate the impact of altimeter undersampling on wind and wave height trends. This author showed that the increasing number of in-orbit altimeters (only one altimeter mission was in orbit in 1985, while more than five were simultaneously in orbit in the late 2010s) led to very large over-estimations of the trends of extreme values computed over the period 1985–2018. Another limitation of altimetry measurements is the reduced performance of radar measurements in the coastal zone, which makes the investigation of decadal variability of nearshore wave conditions more complicated, despite the important role of waves on a coastal sea level [20]. Indeed, surface heterogeneities in the altimeter footprint in the coastal zone can contaminate the radar waveforms and prevent the retrieval of geophysical parameters [21]. However, a number of developments are currently expanding the use of radar altimetry in the coastal zone: higher along-track resolution of radar sensors using the delay-Doppler technology is available for CryoSat-2, Sentinel-3, Sentinel-6 MF, the smaller footprint of SARAL AltiKa, new waveform retracking methods based on sub-waveform partitioning [22,23], and denoising techniques [24].

2.3. Seismic Data

Another source of sea state information at multi-decadal scales may be derived from the seismic stations deployed inland and on the ocean floor to monitor earthquake activity. Indeed, background noise with frequencies between 0.05–0.3 Hz—known as microseism—is present in every seismic record on Earth, and is predominantly caused by ocean waves. The most energetic microseism signal is around a period of 5 s, and its source is proportional to the product of wave energies propagating in opposing directions with twice that period. Moreover, acoustic resonance may amplify the microseism amplitudes at specific water depths [25]. As a result, there is no general simple relation between microseismic energy and wave heights, although increasing wave heights will generally yield increasing microseismc energy, either related to wave reflection off the coast, rapidly turning winds (such as in a tropical storm), or the presence of waves from two independent storms [26].

Despite the complex relationship that links seismic measurements to wave-induced microseism sources over thousands of kilometers, one of the greatest interests in seismic data arises from the long-term time-series that can be obtained, making possible the analysis of sea state decadal variability. Algué [27] used seismic data to monitor typhoons around the Philippines in the late 19th century, and Bernard [28] analyzed seismic data between 1910 and 1975 from 12 seismic stations, mostly located around the North Atlantic, and found near-decadal oscillations in microseism amplitudes, that he interpreted as a result of the solar cycles. Grevemeyer et al. [29] reconstructed 40 years of wave climate in the North Atlantic using records of wintertime microseisms from Hamburg, Germany. The positive trend they obtained between 1954 and 1998 was closely related to similar changes in the northwestern European storm climate and the North Atlantic Oscillation. More recently, Stopa et al. [30] used decadal time-series of seismic records in order to assess the long-term consistency in the CFSR reanalysis after applying a time-varying bias correction based on satellite altimeter measurements. Microseismic records therefore present great potential for investigating the wave climate over time periods that are currently not reachable with other sources of observations. Some difficulties remain, such as the correct inversion of the seismic measurements, the construction of consistent centennial time-series from different instruments, or the comparison between local wave measurements (e.g., in situ buoys) and seismic observations that result from a spatial integration of sources over thousands of kilometers. Longer microseismic periods, generally around 10 s, are generated in shallow water by the interaction of waves and bottom topography [31], with a microseismic spectrum that is proportional to the wave spectrum, in a way that is strongly dependent on the local bottom topography properties. Aster et al. [32] found that these longer-period microseisms gave a more robust indicator of the number of storms than the short-period microseism that is so greatly influenced by directional wave properties.

2.4. Model Reanalysis

Numerical wave modelling is currently the only means capable of providing sea state data at a global scale at the time and spatial resolutions required by GCOS for climate applications [5,11]. Moreover, thanks to decades of oceanic and atmospheric measurements and improved data assimilation techniques, several model reanalyses at multi-decadal scales have been produced and made available by met-ocean agencies, superseding the time duration of most observation products. For these reasons (and others), model hindcasts have been used extensively in the last few decades to investigate sea state decadal variability at global, regional, and local scales [5,33,34,35,36]. Examples of such numerical model reanalyses are the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-5 reanalysis [37] and the 20th century reanalysis [38] (ERA-20C), the NCEP Climate Forecast System Reanalysis [39] (CFSR), and the 20th century Climate Reanalysis [40] (20CR). A major issue in long-term model reanalysis is the changing (increasing) numbers of assimilated observations throughout the years, which may produce spurious trends. Stopa and Cheung [41] intercompared wind and wave reanalysis datasets from ECMWF ERA-I and CFSR, utilizing the same set of altimetry and buoy observations and error metrics to assess their consistency in time and space. These authors showed evidence of prominent discontinuities in the $H_s$ 95 percentiles in 1993–1994 for the Southern Hemisphere, which can be attributed to the introduction of SSM/I winds into the assimilation in 1994, and an increasing trend in both datasets after 2001, which might be due to assimilated data in the atmospheric model common to both. Meucci et al. [42] intercompared several 20th century reanalyses (that assimilate in situ observations of surface pressure and marine surface wind, but no satellite or in situ ocean wave height measurements) with a model-only integration, which did not assimilate any data. They found that reanalyses showed spurious trends in 10 m surface wind speed and significant wave height throughout the whole 1901–2010 period due to the assimilation of an increasing number of observations with changing quality over time. For the more recent altimetry period, Timmermans et al. [18] compared trends in significant wave height over 1992–2018 from two multi-platform altimeter products, the ERA5 reanalysis, and the ECMWF CY46R1 hindcast, a global wave model stand-alone run, forced by ERA5 hourly 10 m neutral winds, surface air density, gustiness, and sea ice cover, without altimeter wave data assimilation. Their comparisons revealed similar spatial distributions but significant differences in the amplitude of the trends between the considered model and altimeter data.

3. Linking Atmospheric and Sea State Inter-Annual to Decadal Variability

In this section, we review our current knowledge of the modes of atmospheric and oceanic natural low frequency variability in the North Atlantic, which may exert some control over the wave climate. We start by succinctly reviewing the most widely used methods to establish the link between the atmospheric surface pressure and wave climate variability in Section 3.1, and the sea state climatology derived from a multi-mission altimeter product is then described in Section 3.2, and the following subsections report the published link between sea state and different climate modes over the NA.

3.1. Methods

Several methods can be found in the literature to establish the relationship between the atmospheric and sea state fields, which will be briefly summarized below. The common feature between these different methods is to seek a common temporal variation between the atmospheric inter-annual and decadal variability and the wave parameters’ variability, not necessarily at the same geographical location. The wave field at a given location is indeed composed of a combination of locally generated wind sea and remotely generated swells. Sea state variability at a given point is thus not only influenced by the variability of the local wind, but also by the variability of winds over remote locations. Because large parts of the North Atlantic are dominated by swell rather than by wind waves [43], the relationship between the sea state variability and the atmospheric variability is not straightforward.

3.1.1. Climate Modes

The spatial structure of atmospheric variability follows recurrent patterns often referred to as climate modes, and since the effect of such climate modes can reach large distances (e.g., El Niño–Southern Oscillation (ENSO) influence on NA climate), climate modes are also sometimes referred to as teleconnection [44]. These modes influence the sea state climate variability, and numerous studies have described these links.

The main method used to establish the link between wave climate and low-frequency atmospheric variability consists of computing 2D maps of the temporal correlation between an established climate mode, such as the North Atlantic Oscillation index and a wave parameter, such as $H_s$, for example, [15,45]. For instance, the 2D map of the temporal correlation between the North Atlantic Oscillation (NAO) index and $H_s$ is:

$$C_{H_s,NAO}(x,y)=\mathrm{corr}(H_s(x,y,t),\mathrm{NAO}\left(t\right)),$$

(1)

where the time correlation between A and B $\mathrm{corr}(A,B)$ is:

$$\mathrm{corr}(A,B)=\frac{{\sum }_i(A_i-\overline{A})(B_i-\overline{B})}{{\left({\sum }_i{(A_i-\overline{A})}^2\right)}^{\frac{1}{2}}{\left({\sum }_i{(B_i-\overline{B})}^2\right)}^{\frac{1}{2}}},$$

(2)

where $A=(A_1,A_2,\dots ,A_{N_t})$ and $B=(B_1,B_2,\dots ,B_{N_t})$ with $N_t$ the number of time records in the dataset and where the overbar represents the time mean. These 2D maps are useful to understand the region of influence of each index and are often complemented with 2D maps of the regression coefficients to give the amplitude of the wave parameters’ anomaly that is associated with the studied index. For instance, the 2D maps of the regression coefficient between the normalized NAO and $H_S$ is:

$$R_{H_s,\mathrm{NAO}}(x,y)=\mathrm{std}\left(H_s(x,y,t)\right)\mathrm{corr}(H_s(x,y,t),\mathrm{NAO}\left(t\right)),$$

(3)

where corr$(A,B)$ is given by Equation (2) and where the time standard deviation of A $\mathrm{std}\left(A\right)$ is:

$$\mathrm{std}\left(A\right)={\left(\sum _i{(A_i-\overline{A})}^2\right)}^{\frac{1}{2}}.$$

(4)

To illustrate this method, panel D of Figure 2 shows the regression map (called “projection” on the figure) of $H_s$ on the normalized NAO index taken from the NOAA website (shown on the bottom left panel of Figure 2). More details about the NAO index and its computation will be given in the following Section 3.3. The map shows that the $H_s$ anomaly associated with the NAO index in the NA is made of a dipole of an opposite sign with a first center at high latitudes between Greenland and Ireland, and a second weaker center at mid-latitudes off the eastern U.S. coast. Up to $1.45$ $\mathrm{m}$ of the $H_s$ anomaly is associated with a unity change in the NAO index off the coast of Ireland, which represents a ∼50% increase of the mean value.

One should be careful when individually studying several climate indices, since these indices may not be orthogonal/independent when they are obtained from different sources, such as, for instance, from SST (as for ENSO-related indices) and an Empirical Orthogonal Function (EOF) analysis of Sea Level Pressure (SLP) [46].

3.1.2. Empirical Orthogonal Functions

EOFs are one of the most commonly used tools to analyse the variability of both $H_s$ and SLP. They link spatial patterns $EO{F}_i$ with a time-series called Principal Component ($P{C}_i$), where i is the mode number. These modes are computed as eigenmodes of the covariance matrix of the detrended $H_s$ [47]. The eigenvalue associated with a particular eigenmode is equal to the variance explained by the mode. In general, only a few modes are necessary to represent most of the variance. The EOF analysis gives $N_t$ modes. The $P{C}_i$ $i=1,2,\dots ,N_t$ form an orthonormal basis such that $(P{C}_i,P{C}_j)={\delta }_{i,j}$, where ${\delta }_{i,j}$ is 0 if $i\ne j$ and 1 if $i=j$, and where $(A,B)={\sum }_{i=1}^{Nt}{a}_i{b}_j$ is the euclidean scalar product. For instance, the decomposition of the $H_s$ field in terms of EOF is:

$$H_s(X,t)=\sum _{i=1}^{N_p}EO{F}_i\left(X\right)P{C}_i\left(t\right)+\mathrm{trend}\left(H_s\right).$$

(5)

$EO{F}_i\left(X\right)$ can be seen either as the projection of the time anomalies of $H_s$ on the $P{C}_i$ vectors using the Euclidian scalar product or as the linear regression coefficient of $H_s(X,t)$ on $P{C}_i\left(t\right)$ using a least-square estimation. The $N_p$ is the number of EOFs used to reconstruct the signal and is generally much smaller than $N_t$, since only a few modes are necessary to represent the variance in $H_s$. The percentage of variance (${\mathrm{perc}}_i$) explained by each mode i is given by the ratio of its eigenvalue ${\lambda }_i$ to the sum of all ${\lambda }_i$:

$${\mathrm{perc}}_i=\frac{{\lambda }_i}{{\sum }_i^{N_t}{\lambda }_i}.$$

(6)

Comparisons between the first EOFs of $H_s$ and SLP are given in Woolf [15], Semedo et al. [43] and shown in panels A and B of Figure 2 and Figure 3 of this article.

3.1.3. Canonical Correlation Analysis (CCA)

The CCA analysis first finds the leading principal components (PCs) of each dataset ($H_s$ and SLP, for instance) separately and retains only the first few modes to explain a given percentage of the variance (typically 90%). The linear combination of these PCs that maximizes the temporal correlation between the two variables [48] is then obtained using a singular value decomposition. The time correlation between the CCA PC of $H_s$ and SLP is optimized by construction, but this optimization is generally made at the expense of the explained variance of each mode. Figure 2 shows the leading CCA mode for SLP and $H_s$ (panels E and F, respectively). The $H_s$ and SLP patterns look similar to the leading EOF patterns, but the correlation between their PC is higher (0.98 vs. 0.85 for the EOF PCs), while their explained variance is lower (47% and 36% for, respectively, the SLP and $H_s$ CCA leading mode vs. 57% and 41% for the EOF analysis). This method is also used, for instance, in Woolf [15], Kushnir et al. [33].

3.1.4. Redundancy Analysis

Redundancy analysis [15,34] is a technique that is used to associate patterns of variation in a predictor field (SLP in our case) with patterns of the predictand field ($H_s$) through a regression model. It seeks to find pairs of predictor and predictand patterns that maximize the associated predictand variance.

3.1.5. How Is the Seasonal Cycle Removed?

The most simple method to remove the seasonal cycle is to subtract the time mean over each season or over each month within a year. In Woolf [15], the seasonal cycle in monthly averaged $H_s$ and pressure fields was removed by subtracting for each month the monthly mean over the entire period of the dataset (1985–1997). PC analysis was then applied on the pooled 41 winter months (December–March) of years December 1985 to December 1997 (no data for winters 1989/1990 and 1990/1991). Shimura et al. [49] argued that the months within the same year cannot be assumed independent, and thus considered the 3 months within a year as one unit. Another method used to account for the seasonal cycle is to fit an analytical model to represent the seasonal cycle variations. For instance, Izaguirre et al. [50] used the following model to account for the influence of NAO:

$$\nu \left(t\right)={\beta }_0+{\beta }_{1}cos\left(\omega t\right)+{\beta }_{2}sin\left(\omega t\right)+{\beta }_{\mathrm{NAO}}textNAO\left(t\right),$$

(7)

where $\omega =\frac{2\pi }{year}$ represents the annual cycle frequency and ${\beta }_0$, ${\beta }_1$, ${\beta }_2$, ${\beta }_{\mathrm{NAO}}$ coefficients obtained from a multivariate linear regression.

3.2. Sea State Climatology in the North Atlantic

The sea state climatology in the NA has been described in many studies, using observations from satellites [15,51], voluntary observing ships [52], or model reanalysis [43,53]. Here, we succinctly summarize the main characteristics of this variability using the multi-mission altimeter product (L4 monthly gridded data, version 1) developed within the ESA Sea State Climate Change Initiative (CCI) project, which covers the period 1993–2018 (see Dodet et al. [54] for more information).

The time mean of $H_s$, computed for each calendar month over the period 1993–2018 using altimeter data, is shown in Figure 4. The largest values, around 6 $\mathrm{m}$, are found in a region centered at high latitudes in the middle of the basin, and for winter months, from December to March. During the summer months, the largest values of $H_s$ are still located in the same region south of Greenland, but have much smaller values with maxima around 3 $\mathrm{m}$. This climatology of monthly mean $H_s$ is in agreement with the one obtained by Woolf [15] with satellite data in a shorter period of time (1991–2000).

Figure 5 shows the inter-annual variability of $H_s$, represented as the standard deviation computed on $H_s$ for each month over the 1993–2019 period. The largest inter-annual variability is located in a region centered off the coast of Ireland during the winter months (mainly November, December, January, and February) with values up to $1.4$ $\mathrm{m}$ in February. The inter-annual variability is much weaker at low latitudes (<20° N) for all months. We can also note an increased $H_s$ variability along the southern U.S. Atlantic coast during the months of August and September, which may be related to tropical cyclone activity.

3.3. North-Atlantic Oscillation

A major source of climate variability in the northern hemisphere is the North-Atlantic Oscillation (NAO), which is a north-south dipole characterized by sea level pressure anomalies of opposite signs between mid and high latitudes with a southern center near the Azores and northern center near Iceland [55]. The pattern explains approximately 18% of the winter variability of the sea level pressure field in the Northern Hemisphere [56,57] and has no preferred time-scale of variability [58]. A positive NAO is associated with a more zonal flow regime, stronger surface westerlies at mid-latitudes over the North Atlantic, and with anomalous southerly flow over the eastern U.S. and anomalous northerly flow across the Canadian Arctic and the Mediterranean [57,59,60,61]. The NAO is linked with the variability of the North Atlantic storm track [61]. During positive NAO phases, low-pressure systems follow a trajectory that goes from the North Atlantic toward Northern Europe, and during a negative phase, the pressure gradient between the Icelandic low and Azores high centers of action is weaker, displaced to the south, and directed toward the Iberian peninsula [61].

There is no unique way to define the NAO index. Two main methods to compute this index are usually found in the literature. The first method derives the index from the pressure gradient between a northern and southern location in the Atlantic, usually Stykkisholmur/Reykjavik in Iceland and Lisbon in Portugal (see, for instance, Hurrell et al. [55], where historical meteorological records are available). The second method computes the index as the PC of the leading EOF of SLP [56]. Other authors have used compiled databases of instrumental station pressure, temperature, and precipitation measurements and proxy data, such as tree rings, in order to reconstruct the multi-centennial NAO index (see, for example, Luterbacher et al. [62]). A disadvantage of station-based indices over the EOF approach is that, since the station’s positions are fixed in time, their signal contains noise due to transient meteorological phenomena not related to NAO, while the PC time-series gives a better representation of NAO spatial pattern. However, an advantage of station-based indices is that they can be computed on much longer time-series. Figure 6 shows station-based and EOF-based indices starting in 1865 and 1950, respectively. Although there are some differences between the two indices, there is good agreement between the long-term variations. A large positive trend between 1960 and 1990 (Figure 6) has prompted many authors to investigate its origin, for example, [63,64]. Comparing the natural variability arising in control runs (i.e., without increased greenhouse gas) of climate models with an observed NAO index, they suggested that it is unlikely that this trend is only due to internally generated climate signals and that greenhouse gas forcing has contributed to the observed trend. Thus, the NAO index might contain part of the anthropogenic signal due to increased greenhouse gas, and should not be interpreted as a pure mode of internal variability.

Because the NAO explains a large part of the SLP variability over the North Atlantic region—around 40% of the total variance of winter SLP in the NA [57]—it has long been identified as a major source of sea state variability. Using measurements from wave recorders installed off the southwestern tip of England (Seven Stones lightship) and from the ocean weather station Lima (57° N, 20° W) and pressure gradient data (calculated between the Azores and Iceland) from the UK Meteorological Office, Ref. [66] was among the first to establish a relationship between $H_s$ and SLP. Kushnir et al. [33] performed a Canonical Correlation Analysis (CCA) of monthly mean SLP fields from ECMWF model reanalysis and simultaneous $H_s$ hindcast fields for the period 1980–1989, and found that the SLP mode is a perturbation resembling the NAO pattern. The corresponding mode associated with $H_s$ is negative in the northeast Atlantic, and positive along 30° N.

Recognising that hindcasts of the wave climate could lead to significant systematic errors that could potentially impact long-term trends [67], Woolf [15] performed a similar analysis on $H_s$ data obtained from carefully calibrated satellite missions (ERS-1, ERS-2, and TOPEX) for the period 1991–1999. Similarly to Kushnir et al. [33], their results show that NAO contributes to a large part of the sea state variability in the northeastern part of the Atlantic, suggesting that this is a very robust feature of the North Atlantic sea state variability.

To illustrate the impact of the NAO on the sea state variability in the Atlantic, Figure 2 shows the EOF of SLP and $H_s$ (panels A and B), the projection of SLP and $H_s$ on the NAO index (panels C and D), and the CCA leading mode between SLP and $H_s$ (panels E and F). $H_s$ and SLP were obtained from ERA-5 reanalysis and computed as a January-February-March (JFM) average over the 1980–2018 period, and the limits of the studied region are 0° N–90° N and 100° W–20° E. 90% of the total variance of SLP, and $H_s$ was retained to compute the CCA. The three methods shown here depict similar variability, which has previously been described by Kushnir et al. [33] and later by Woolf [15]. The variability of the SLP consists of two anomalies of opposite signs. When the NAO is in its positive phase, the SLP mode consists of a negative anomaly over the northeastern Atlantic, east of Greenland and north of 55° N, and of a positive anomaly south of 55°, centered close to the Azores. The corresponding $H_s$ mode is positive in the northeastern Atlantic, with a center off the coast of Ireland and negative at mid-latitudes, centered in the western part of the Atlantic. The global patterns between the three methods are very similar.

The time correlation between the NAO index and SLP and $H_s$ PC from EOF and CCA analysis is shown in the bottom-right panel of Figure 2. The correlation of the CCA SLP and $H_s$ PCs is, by construction, very high (0.98), higher than the correlation between the PC of SLP and the PC of $H_s$ (0.85) obtained from the EOF analysis. The NOAA NAO index is close to both the SLP and $H_s$ PC obtained from CCA analysis (correlation 0.95 and 0.92 for SLP and $H_s$, respectively) and also close to the SLP and $H_s$ PC obtained from EOF analysis. The percentage of variance explained by each mode is shown in

Table 2. Percentage of variance explained by each method (EOF, NAO, or EA index, CCA) for the leading mode (mode 1) and the second mode (mode 2) of SLP and $H_s$.

	Mode 1			Mode 2
	EOF	NAO	CCA	EOF	EA	CCA
SLP	57%	47%	47%	12%	6%	10%
$H_s$	41%	33%	36%	24%	13%	12%

. The maximum explained variance is, by definition, given by the EOF analysis (57% of the SLP variance and 41% of the $H_s$ variance), the CCA analysis optimizes the correlation between SLP and $H_s$ at the expense of the variance explained for each field: 47% for the SLP and 36% for the $H_s$. The percentage of variance explained by the NAO index is similar to that of the CCA (47% and 33% for the SLP and $H_s$, respectively).

NAO can also be used as a predictor of extreme wave values, as shown by Izaguirre et al. [50] and Kumar et al. [68] using satellite data and SLP in the period 1992–2006.

3.3.1. Relationship between NAO and Wind Sea and Swell

Using voluntary observing ship data (VOS) [69], Gulev and Grigorieva [52] computed the first EOF of the total significant wave height $H_s$, the wind sea height $H_s^w$, and the swell height $H_s^s$. They found noticeable differences between the spatial patterns of the three fields. The maximum explained variance for the wind sea is in the central subpolar North Atlantic, while the swell pattern has two maxima: one in the North European basin, and the other in the Greenland-Iceland-Norwegian sea. Moreover, their results suggest that the PCs of $H_s^s$ and $H_s^w$ are not highly correlated, and the correlation between the NAO index is higher with $H_s$ than with any of these two fields. However, this is at odds with the studies of Semedo et al. [43] and Martínez-Asensio et al. [36] who found, both using hindcast models, that the NAO index is well-correlated with $H_s^w$ but not with $H_s^s$.

Using CCA analysis, Gulev and Grigorieva [52] further suggested that wind waves are more correlated with the wind speed, while swell is associated with the number of deep cyclones. An increased number of cyclones in the northern European basin leads to increased swell in the northeast Atlantic.

Semedo et al. [70] used NORA10, a high-resolution regional downscaling of the ERA-40 reanalysis produced by the Norwegian Meteorological Institute to study the sea state climate in the North Sea, Norwegian Sea, Greenland Sea, and Barents Sea. They found that the correlation of DJF $H_s$ and the NAO index over 1958–2001 was high in the region, except in the Barents Sea and along the east coast of the UK. The correlation of $H_s^s$ with the NAO index is generally higher than with $H_s^w$, except in the North Sea where the sheltering effect of the British Isles is strong.

3.3.2. Relationship between NAO and Other Sea State Parameters

Bromirski and Cayan [4] showed by using a hindcast model (Wave Watch 3: WW3) over 60 years (1948–2008) that the wave power $P_w=a{H}_s^2{T}_P$, with a as a constant and $T_P$ as the peak wave period, is also mainly controlled by the NAO in the North Atlantic. Satellite altimeter data can also be used to show that the wave period $T_P$ has a significant positive correlation with the NAO index [71] off the coast of Ireland, and a negative correlation at mid-latitudes around the Azores.

Using a North Atlantic wave model (WW3) forced by wind fields from the NCEP/NCAR Reanalysis project [72] over 1953–2009, Dodet et al. [45] studied the correlation of the NAO index with the mean wave direction, MWD, and $T_P$ along with $H_s$ in the northeast Atlantic. They showed that MWD correlation with NAO is weak everywhere except south of Portugal and off the coast of Morocco, where negative correlations up to 0.7 were found, and that the correlation of NAO with $T_P$ gave significant values off the coast of Ireland. A rotation of 20° clockwise of MWD per unit of the winter NAO index was found around southern England [73]. Using WW3, Charles et al. [74] found that the positive NAO phase was associated with an increase in winter wave height, mean period, and with a northerly shift in the wave peak direction.

3.4. East Atlantic Pattern

The East Atlantic Pattern is the second mode of variability over the North Atlantic region, and was identified using an EOF analysis of SLP data [56,75]. The location of its main center of action is to the south of Ireland and west of the United Kingdom [76,77], and can be seen as a modulation of the locations and intensities of the Icelandic Low and Azores High associated with NAO. Comas-Bru and McDermott [78] showed that the combination of the Eastern Atlantic (EA) index and NAO is important in correctly describing winter climate variability in Europe. Additionally, the variability of the North Atlantic jet stream is well-described by the combination of NAO and EA [79].

Using a combination of EOF and Redundancy Analysis of SLP observations taken from the Southampton Oceanography Center of climatology [80] and $H_s$ obtained from altimeter measurements, Woolf [15] found that the EA was associated with an anomaly pattern of wave height located west of the British Isles and continental Europe. They found that this pattern was close to the pattern of the second EOF of wave height that explains 19% of the total inter-annual variance of wave height in the NA (for the winter period of 1985 to 1997). Using ERA-40 reanalysis data, Shimura et al. [49] showed that EA correlates positively with JFM averages of $H_s$ in the eastern North Atlantic and negatively in the Mediterranean Sea. Similar results were obtained by Martínez-Asensio et al. [36] and Castelle et al. [81] using hindcast model data. Martínez-Asensio et al. [36] found that the impacts of EA on MWP has a similar pattern to that obtained for $H_s$ and are up to 1.1 s per positive unit of EA. EA has also been shown to induce counterclockwise changes in MWD of 20° per positive unit EA along the European coasts and clockwise changes of 60° in the northwestern North Atlantic. Extreme wave heights are positively correlated with EA in the northwest Atlantic and are negatively correlated with Mediterranean extreme wave height [50]. In the bay of Biscay, Charles et al. [74] found correlations between all seasons of $H_s$ and EA, and a positive EA was associated to an increase of $H_s$ and the mean wave period.

Figure 3 shows the second EOF modes of JFM SLP and JFM $H_s$ computed from ERA-5 data over the period 1980–2018. The SLP and $H_s$ second EOF modes explain 12% and 24% of the SLP and $H_s$ total variance, respectively, and the correlation between the two PCs is weak (0.52). The SLP second EOF pattern is a tripole with two of the same sign anomalies at low and high latitudes (the low latitude anomaly is weak) and an opposite sign anomaly in the middle, at mid-latitudes. The $H_s$ second EOF pattern takes the form of a dipole with a first center of action in the mid-NA at 45° N and a second in the Norwegian sea.

The projection of SLP anomalies on the EA index (panel C Figure 3) gives a dipole with a weak anomaly at low and mid-latitudes and a stronger anomaly of an opposite sign south of Iceland. The projection of $H_s$ on the EA index gives a zonal anomaly along the latitude 45° N. Patterns of SLP and $H_s$ associated with the EA index are both quite different from the EOF pattern, and this is also apparent from the time correlation between the EOF PC of SLP (0.17) and EA and EOF PC of $H_s$ and EA (0.51). Because the NOAA EA index is obtained from an EOF analysis of the SLP field on the Northern Hemisphere region, we could have expected it to be similar to the second EOF PC calculated on SLP. Several reasons might explain this difference: (1) In this study, the SLP EOF were computed on a smaller region (i.e., the NA) than the entire Northern Hemisphere used by the NOAA index; (2) the definition of the NOAA EA index uses a longer time period (1950 to 2010 compared to the 1980 to 2018 in ERA-5); and (3) the NOAA index makes use of an additional processing technique (varimax) that we do not use in our EOF computation. It is interesting to note that the variance explained by the second EOF of $H_s$ (24%, see Table 2) is larger than that obtained with the EA index (13%).

Because the CCA analysis optimizes the correlation (0.97), the SLP and $H_s$ patterns obtained from this CCA analysis are quite different from the EOF patterns of SLP and $H_s$ (see Figure 3). There is no correlation between the SLP CCA PC and the EA index (0.05) nor between the $H_s$ CCA PC and the EA index (0.04). The $H_s$ CCA PC represents 12% of the JFM $H_s$ variance, which is similar to the variance explained by the EA mode, but less than the variance explained by the second $H_s$ EOF. In Appendix A, it is shown that it is possible to construct a new index with a high correlation with the second PC of $H_s$. This new index is based on the second and third PCs of winter SLP. The fact that the second EOF of $H_s$ has a better correlation with a combination of the second and third PC of SLP than with the second PC of SLP alone outlines the complex relationship between winter SLP and winter $H_s$.

To summarise this subsection, if the EA index clearly has an impact on the NA wave climate as already demonstrated by a number of studies, for example, by [15,49], this mode of climate variability does not seem to represent all the variance contained in the second EOF of $H_s$. The representation is greatly improved when using a linear combination of the second and third EOFs of SLP over the NA instead.

3.5. Other Modes of Variability

In this subsection, we review other climate modes that are described in the literature as having an influence on sea state variability in the NA. The time correlation between the JFM average of the different indices discussed in this review and the JFM average of ERA-5 $H_s$ for the period 1980–2018 is shown in Figure 7.

3.5.1. Scandinavian Pattern (SCA)

The Scandinavian pattern (SCA) was first defined by Barnston and Livezey [56] through a Rotated Principal Component Analysis analysis of monthly mean 700 mb height anomalies over the extratropical latitudes in the northern hemisphere. The pattern has the primary centre of action around the Scandinavian Peninsula, and two others of opposite signs over the northeastern Atlantic and central Siberia.

Shimura et al. [49] found a positive correlation between winter $H_s$ and SCA off the coast of western Africa and negative correlations between Iceland and Ireland. Izaguirre et al. [50] found a positive correlation between SCA and extreme wave height in the eastern part of the Atlantic (mostly around the Canary Islands).

The time correlation between ERA-5 JFM $H_s$ and the NOAA SCAN index is shown in Figure 7. Significant positive correlations were found in the Mediterranean sea and along the coast of Morocco. Negative correlations were found north of Ireland and the U.K.

3.5.2. East Atlantic/West Russia pattern (EA/WR)

EA/WR is characterized by two main large-scale anomalies of SLP located over the Caspian Sea and western Europe [56]. Depending on the region and time period used for the EOF analysis, the corresponding index is obtained as the third or fourth EOF of SLP. Shimura et al. [49] showed that winter $H_s$ has a weak but significant (at the 95% level) positive correlation with EA/WR in most of the mid-NA region, and a negative correlation was found in the Mediterranean sea and off the U.S. east coast. Martínez-Asensio et al. [36] also computed the correlation between EA/WR and NA winter wave parameters and found a significant negative correlation at low latitudes between 12° and 24°, at odds with the results of Shimura et al. [49]. This might be partly due to the two different seasonal periods used in these studies, that is, DJF means in Shimura et al. [49] and DJFM means in Martínez-Asensio et al. [36]. Izaguirre et al. [50] found a correlation between EA/WR and extreme wave height in the Mediterranean basin, particularly between the Balearic Islands and Corsica and Sardinia Islands. We found a weak significant negative correlation between the JFM EA/WR index and ERA-5 JFM $H_s$ in the Bay of Biscay and Western Mediterranean sea, and a positive correlation off the northern coast of Norway (Figure 7).

3.5.3. Pacific-North American (PNA)

The PNA pattern has alternating centers of anomalous pressure that arc northeastward through the North Pacific Ocean, through Canada, and then curve southeastward into North America [56,75]. The associated PNA index is obtained by applying a Rotated PC Analysis to SLP in the northern hemisphere. The correlation between the PNA index and winter $H_s$ is negative in the eastern part of the mid- and high-latitude NA and at low latitudes [49]. Bromirski and Cayan [4] found a positive correlation between PNA and wave power in the western Atlantic, off the U.S. coast. These patterns of correlation were also found between ERA-5 JFM $H_s$ and the NOAA JFM PNA index (Figure 7).

3.5.4. Tropical North Hemisphere (TNH)

The TNH has a center off the Pacific northwest coast of the US at 40°–50° N and 125°–140° W, the same sign center near Cuba and an oppositely signed center near the U.S. Great Lakes [56]. The correlation between the TNH index and winter $H_s$ is positive at low latitudes and at high latitudes. At mid-latitudes, the correlation is negative along the U.S. East coast and positive off the coast of Morocco [49] and Figure 7.

3.5.5. El Niño-Southern Oscillation (ENSO)

ENSO is the dominant mode of inter-annual climate variability at a global scale [82]. It is characterized by large-scale sea surface temperature (SST) anomalies in the eastern equatorial Pacific Ocean with time-scales ranging between 2 and 7 years [83]. SST anomalies linked with ENSO usually appear in the eastern Pacific, and terminate in the central Pacific. ENSO is composed of two active phases (warm El Niño and cold La Niña) determined by SST anomalies in the Equatorial Pacific Region.

Shimura et al. [49] computed the correlation between the NINO-3.4 index (defined by SST anomalies in the 5° S–5° N and 170° W–120° W region) and winter $H_s$, and found a positive correlation off the U.S. east coast, and negative correlations east of the Caribbean Islands (see also Figure 7). Bromirski and Cayan [4] found similar correlations between the NINO-3.4 index and winter wave power $P_w$, but no significant correlation between NINO3.4 and summer $P_w$. Stopa et al. [53] used the Southern Oscillation Index (SOI) to represent ENSO, this index being obtained from pressure measurements at Tahiti and Darwin (see references in [53]). Their studied wave parameters consist of the wave intensity, defined as the monthly average of the highest 10% wave height, and their occurrence was defined as the monthly records above the 90th percentile of the entire time-series. They found that El Niño was significantly associated with wave intensity mostly in the western half of the NA. Wave occurrence was associated with El Niño off the U.S. east coast and with La Niña at low latitudes mostly to the east of the Caribbean Islands, which gave a pattern similar to that obtained by Shimura et al. [49] and Bromirski and Cayan [4] with winter $H_s$ and NINO-3.4.

3.5.6. Southern Annular Mode (SAM)

The Southern Annular Mode (sometimes called Antarctic Oscillation) is an important mode of variability of the high and mid-latitudes of the Southern Hemisphere [84]. SAM can be defined as the leading PC of the 850-hPa geopotential height anomalies south of 20° S [85]. It is a zonally-symmetric mode of variability with latitudinal bands of pressure anomaly with opposite signs between the mid and high latitudes. A positive (negative) SAM is associated with negative (positive) sea level pressure anomalies at high latitudes and positive (negative) SLP anomalies at low latitudes. During positive SAM, westerlies are intensified and weakened south and north of 40° S, respectively.

Although this mode principally influences the wave climate in the Southern Hemisphere, Marshall et al. [86] has suggested a link between the SAM and NA wave climate. Using WW3 forced by a 31 year reanalysis of wind given by the Climate Forecast System Reanalysis for the period 1979–2010, they showed that SAM is significantly linked with winter $H_s$ in the North Atlantic, off the coast of continental Europe, during neutral ENSO years. When El Nino and La Nina years are included, the relation between North Atlantic $H_s$ and SAM is much smaller, suggesting that ENSO plays a key role in offsetting the link between the SAM and the Northern Hemisphere winter climate. We found significant positive correlations in the western part of the Atlantic basin at low latitudes between JFM SAM (obtained from NOAA) and ERA-5 JFM $H_s$ (Figure 7).

3.5.7. West Europe Pressure Anomaly (WEPA)

Castelle et al. [81] defined a new station-based climate index called the “West Europe Pressure Anomaly” (WEPA) computed to optimize the explanation of winter $H_s$ variability along the Atlantic coast of western Europe from surface pressure measurements. The index is based on SLP gradients measured between the stations Valentia (Ireland) and Santa Cruz de Tenerife (Canary Islands). A positive phase of WEPA is associated with an intensified latitudinal SLP gradient in the NE Atlantic. Contrary to the other indices presented in this section, which are mainly constructed to describe atmospheric variability and then correlated with sea state climate, WEPA is constructed specifically to optimize the representation of wave parameters along the European coast. Figure 7 shows that significant positive correlations between JFM $H_s$ and JFM WEPA were indeed found in the bay of Biscay and off the western coast of Europe.

To summarize this section, a large number of climate indices have been shown to be correlated with winter $H_s$ in several regions of the NA. To better understand what part of JFM $H_s$ variance is explained by the sum of all climate modes considered here, we apply a multivariate linear regression to the sum of the modes and a linear trend to each point of $H_s$:

$$H_s(x,y,t)=\sum _{i=1}^9{\alpha }_i(x,y){\mathrm{index}}_i\left(t\right)+\alpha (x,y)t+\beta (x,y)+ϵ(x,y),$$

(8)

where ${\alpha }_i$ is the coefficient given by the multivariate linear regression for ${\mathrm{index}}_i$, which is one of the nine indices studied here (i.e., NAO, EA, SCA, EA/WR, PNA, TNH, ENSO, SAM and WEPA). The $\alpha$ is the coefficient of the linear trend, $\beta$ is the intercept, and $ϵ$ is the error. The coefficient of determination $R^2(x,y)$, which is the ratio of the variance explained by the model over the total variance of $H_s$, is shown in Figure 8. High values of $R^2$, which indicate that the sum of the modes explain most of the $H_s$ variance, are found off the coast of western Europe and along the coast of Morocco, in a latitudinal band around 35° N and in the western half of the basin at low latitudes. Low values of $R^2$, associated with low levels of explained $H_s$ variance, are found in the eastern half of the basin at low latitudes, in the Gulf of Mexico, in a latitudinal band north of 20° N, east of the UK, along the northeastern coast of North America and in the south-west of Greenland.

3.6. Decadal Change in Wave Heights

A number of studies have investigated the long-term trends in various wave parameters using observations and hindcast models. Because it is not always clear whether these trends are associated with anthropogenic forcing or with the modes of atmospheric variability described above, we confine trend description to this section and dedicate the next section to the projected wave climate change under the greenhouse gas scenario obtained using simulations from climate models.

Carter and Draper [9] and Bacon and Carter [12] used instrumental records from a ship-borne wave recorder at Ocean Weather Station L and Seven Stone Light Vessel between 1962 and 1986 in the NA and found a positive trend of about 1% per year in significant wave height at these locations.

Using voluntary observing ship data, Gulev and Hasse [87] found an increase of $H_s$ of the order of 1 $\mathrm{cm} {\mathrm{year}}^{-1}$ to 3 $\mathrm{cm} {\mathrm{year}}^{-1}$ between 1964 and 1993 over the NA except for the western and central subtropics. They also showed that wind wave and swell had different trends: swell trend is similar to that of $H_s$, while wind wave has positive trends only in the central mid-latitude Atlantic. They attributed the positive trends in swell to the intensification of high-frequency synoptic processes. Similarly, Bouws et al. [88] found an increase in wave height of approximately $2.3$ $\mathrm{cm} {\mathrm{year}}^{-1}$ in the NA, using a collection of more than 20,000 hand-drawn wave charts covering the period 1960–1988. Studying 50 years of ship observations in the NA, Gulev [89] found slightly lower mean wave height trends of the order of $1.4$ $\mathrm{cm} {\mathrm{year}}^{-1}$.

Numerous studies have used wave hindcast models to compute wave height trends. Günther et al. [90] used a 40-year (1955–1994) reconstruction of the wave climate in the NA with the WAM wave model and found a large increase of 5–10 $\mathrm{cm} {\mathrm{year}}^{-1}$ for most parts of the Northeast Atlantic. Using a 109-year (1900–2008) numerical wind-wave hindcast, Bertin et al. [91] also found a positive trend in $H_s$, but much smaller, reaching 1 $\mathrm{cm} {\mathrm{year}}^{-1}$ north of 50° N in the NA. They attributed this trend to a positive trend in wind speed over the same region. In the Gulf of Mexico, Appendini et al. [92] used a 30-year hindcast to show that extreme wave heights increased at a rate of 7–8 $\mathrm{cm} {\mathrm{year}}^{-1}$ in October/September over the period of study because of an increased cyclone intensity in the last decade. However, trends computed from reanalysis should be interpreted with caution; indeed, Meucci et al. [42] showed that spurious trends in $H_s$ are present in the entire period of 1901–2010 due to an increasing number of observations and change in the data quality.

Using satellite data from various missions covering the period of 1985–2008, Young et al. [93] and Young and Ribal [17] found positive trends in wind speed over the North Atlantic, and stronger positive trends for the 90th percentile, indicating that extreme values of wind speed are increasing faster than the mean. Using an updated version of their multi-satellite sea-state database covering the period 1985–2018, Ref. [17] found weak positive trends for the mean $H_s$, of the order of $0.3$ $\mathrm{cm} {\mathrm{year}}^{-1}$ and larger positive trends for the 90th percentile of $H_s$ in the mid to high latitudes of the NA (∼$0.8$ $\mathrm{cm} {\mathrm{year}}^{-1}$). No significant trends of $H_s$ (mean or 90th percentile) were observed in the equatorial region.

Because the trends computed from multi-mission satellite data might be quite sensitive to the method used to calibrate the various satellite missions, Timmermans et al. [18] compared the JFM linear trends for the period 1992–2017 for four different datasets: Ribal and Young [94] database, the ESA Sea State Climate Change Initiative database, and two ERA5-based reanalysis and hindcast surface wave products from ECMWF. They showed that significant differences in trend magnitude and pattern in the NA existed, particularly at low latitudes, where Ribal and Young [94] database gives negative trends, while the trend is positive in the three other databases. Nonetheless, some patterns seem to appear in all databases, such as the positive trend of approximately 2 $\mathrm{cm} {\mathrm{year}}^{-1}$ off the northeastern U.S. coast and the negative trend (2–4 $\mathrm{cm} {\mathrm{year}}^{-1}$) in the northeastern Atlantic at high latitudes.

Figure 9 shows the $H_s$ JFM trends computed in the NA over the period 1993–2018 from the ESA Sea State CCI v1 database using a simple linear regression. Negative values, up to $-4$ $\mathrm{cm} {\mathrm{year}}^{-1}$, were found in the Norwegian Sea and are statistically significant at the 5% level. Strong positive trends of $H_s$, up to 4 $\mathrm{cm} {\mathrm{year}}^{-1}$, were found in the western part of the Mediterranean sea. Some significant positive trends can also be found along the eastern U.S. coast and at low latitudes.

4. Projected Climate Change Impact on Sea State Variability

In this section, we review our current understanding on the projected modification of the wave climate due to anthropogenic climate change. Because coupled atmosphere-ocean general circulation models (GCM) do not generally include wind-wave-dependent parameterizations, projection of the wave climate is performed using atmospheric forcing derived from climate models [95,96]. Two different approaches to obtain the wave field from the projected atmospheric forcing are generally found in the literature. The first consists of determining a statistical link between SLP anomalies and the SLP gradient to the wave field using, for instance, data from a present-day reanalysis (such as ERA-Interim) [97], for example, Wang et al. [2,98]. The second approach is dynamical, for example, Hemer et al. [99], and makes use of wind wave models to resolve the wave field. The advantage of the statistical approach is a lower computational cost, whereas the dynamical approach is expected to better resolve the characteristics of the simulated wave field in the context of a changing climate.

Morim et al. [100] identified three dominant sources of uncertainty for these wave climate projections: the emission scenarios, the GCM employed (and thus the quality of the simulated surface wind or SLP used to obtain the wave field), and the choice of wind wave model, that is, the parameterization of the model when a dynamical approach is used, or the details of the statistical method otherwise. The confidence of the wave climate projection crucially depends on the forcing field, that is, the atmospheric circulation; however, in general, low agreement across climate models in the projected atmospheric circulation is found [101], particularly in the NA. This low agreement can be attributed to the strong natural variability as compared to the anthropogenic forced signal and to the sensitivity to the model error of the circulation response to climate forcing [101]. It is also interesting to note that the IPCC AR6 report ([102], chapter 9) considered that there is medium level of confidence in projections of changes in the mean wave climate and low confidence in the projected changes in extreme wave conditions.

4.1. Changes in Extra-Tropical Atmospheric Circulation

The most prominent feature linked with global warming in NA extra-tropical atmospheric circulation is the poleward shift of the storm tracks and the decrease of the total number of cyclonic storms, for example, Yin et al. [103,104,105,106].

The Hemer et al. [95] used results from CMIP3 multi-model ensemble wave-climate projections from the Coordinated Ocean Wave Climate Project [107] to assess changes in the wind-wave climate due to future climate change. A decrease of approximately 5% of $H_s$ in the NA for all seasons for the period 2070–2100 relative to the present climate (1979–2009) was found, as well as a decrease of the mean wave period. They explained this decrease by a projected positive trend in the North Atlantic Oscillation that is associated with a decrease in $H_s$ in the central Atlantic. Similar results were found by Wang et al. [98] and Camus et al. [108] using statistical approaches in CMIP5 and by Morim et al. [100] with an analysis of 10 independent state-of-the-art studies using atmospheric forcing fields also obtained from CMIP5.

Meucci et al. [109] investigated the evolution of extreme wave conditions by 2100 under two IPCC greenhouse gas emission scenarios (RCP4.5 and RCP8.5). In the NA, they found a decrease of the 100-year return period significant wave height ($H_s^{100}$) by 2100 in low to mid-latitudes (−5 to −15%) and an increase at high latitudes (however, this increase is not statistically significant). The magnitude of this increase and decrease pattern is stronger in the case of high-emission scenarios (RCP8.5) compared with mid-emission scenarios (RCP4.5). They further found that the decrease of $H_s^{100}$ in the low to mid-latitudes of the NA was linked with a decrease in the frequency of extreme storms. Focusing on the extreme wave energy flux (WEF) along coastlines, Mentaschi et al. [110] showed that WEF was expected to decrease almost everywhere in the NA except in the Baltic Sea, where they found a significant increase of WEF.

4.2. Changes in Tropical Cyclone Activity

Tropical cyclones (TC) can generate extreme wave conditions, which may result in catastrophic coastal damage, due to the combined effect of destructive winds, storm surge, and heavy rainfalls [111]. Because of the strong coupling between sea surface temperature and TC intensification, ocean warming induced by anthropogenic forcing is expected to impact TC dynamics [112]. Although there is still little confidence that currently observed changes in TC activity can represent detectable or attributable anthropogenic changes [113], a WMO expert team, which assessed model projections of the TC activity response to anthropogenic warming in climate models, stated with medium to high confidence that 2 °C anthropogenic global warming is projected to increase the global average TC intensity, and that very intense (category 4–5) levels will increase with a median projected change of +13% [114]. Timmermans et al. [115] forced a global spectral wave model with wind forcings from the Community Atmosphere Model at horizontal resolutions of 1.0° and 0.25° resolution, and found that differences in extreme wave height manifested most strongly in tropical cyclone (TC) regions, emphasizing the need for high-resolution forcing in those areas.

Regarding TC changes in the NA, Belmadani et al. [116] investigated future projections of NA tropical cyclone-related wave climates using a new configuration of the ARPEGE-Climate global atmospheric model on a high-resolution grid able to reproduce the distribution of tropical cyclone winds, including Category 5 hurricanes. They used historical (1984–2013, five members) and future (2051–2080, five members) simulations with the IPCC RCP8.5 scenario to drive the MFWAM (Météo-France Wave Action Model) spectral wave model over the Atlantic basin during the hurricane season. They found that future projections exhibited a large-scale decrease in average wave heights throughout the hurricane season. This is driven by a weaker and poleward-displaced anticyclone, yet concurrent increases in extreme TC-related wave heights within a large region extending from the African coasts to the North American continent were associated with stronger TC activity around Cape Verde and related changes in the wind sea. However, they also found that the largest expected changes in extreme wave heights were found near the coast of the northeastern United States, possibly related to slight poleward displacement of TC activity, and highest during the season peak phase from mid-August to mid-September.

4.3. Wave–Ice Interactions

Sea ice can be found in several marginal seas of the subpolar North Atlantic, such as the Labrador Sea, Greenland Sea, and Denmark Straight, and is known to impact wave generation and propagation [117,118]. In addition, waves propagating into sea ice may cause ice break-up and modify the ice structure and its mechanical properties [119]. Because of the declining Arctic sea ice extent during the last decades, the influence of Atlantic swells on the Arctic wave climate has increased, and the connection between these two basins has strengthened [120]. Therefore, understanding the feedback mechanisms between waves and sea ice is key for addressing the impact of climate change on the subpolar NA sea state variability. Waves have been observed to propagate hundreds of kilometers into sea ice [121], yet they suffer strong attenuation rates that vary with both wave (amplitude and frequency) and ice characteristics (thickness and mechanical properties) [118]. Stopa et al. [122] combined altimeter data and wave hindcasts to investigate past changes in the Arctic wave climate. They found that the reduction in extent of ice over the 1992–2014 time period enhanced sea states with higher wave heights, longer wavelengths, and more persistent swells. Moreover, the sea ice minimum occurred later in fall when the wind speeds increased, creating more favorable conditions for wave development. Projections of the wave climate in the Arctic under the RCP 8.5 scenario indicated that the maximum significant wave height and associated wave periods will significantly increase by the end of the century as the ice-free season lengthens and Arctic waters become more exposed to storms in autumn [120]. These projected changes in wave conditions will lead to a general increase of wave-driven erosion along the Arctic coastlines, as already observed, for instance, along the Canadian coast of the Beaufort Sea [123].

5. Conclusions

In this article, we have first reviewed the main datasets used to study the sea state on inter-annual to decadal time-scales. The VOS dataset [10] provides unique information on the historical sea state (since 1853) but suffers from low coding precision and is spatially constrained to the major navigation routes. Decadal time-series of sea state data can also be retrieved from instrumental in situ observations, but periodic changes of the wave instruments and relocations have impacted the long-term consistency of these data. Radar Altimetry has now almost 30 years of continuous data with a global scale coverage, and is thus particularly suited for the study of the long-term sea state variability. However, some limitations arise because of the low revisit time period of altimeter missions and the need to aggregate multi-mission records with efficient inter-calibration strategies to reduce the impact of data heterogeneities. Because of the availability of long-term and historical measurements, seismic data could be a valuable source of sea state observations but have rarely been exploited until now because of the complex relationship between seismic data and sea states. Apart from these observations, numerical wave modelling is widely used to study the sea state’s long-term variability. Taking advantage of decades of atmospheric and oceanic measurements, reanalysis products such as ERA-5 can provide accurate sea state information at a global scale over several decades. However, due to the assimilation of an increasing number of observations with changing quality over time, spurious trends may appear in long-term evolution derived from model reanalysis. Future progress in our understanding of sea state variability will be fostered by the exploitation of decades of spectral swell observations provided by either SAR missions or the French Chinese Ocean Satellite mission (CFOSat) dedicated to the observations of ocean waves and wind [124]. Indeed, spectral wave parameters are key components of the sea state description, yet they are still poorly constrained in model reanalysis due to a lack of observations, particularly in the southern hemisphere [125]. The analysis of historical seismic time-series also represents a promising avenue for exploring sea state conditions in the late 19th and early 20th centuries that have not been observed otherwise. Moreover, the recent explosion of in situ measurement through low-cost buoys [126] also provides potential future opportunities for tracking change, though with potential unknown issues, such as homogeneity between technologies or public/private data-sharing agreements.

Maps of standard deviation of wave height computed on each month of the altimeter dataset for the period 1993–2018 in the NA (Figure 5) show that inter-annual variability is larger during winter months and for latitudes north of 40° N. Therefore, most of the literature focuses on the winter months, and we adopted the same strategy. We then reviewed the current understanding of the relationship between the inter-annual to decadal sea state variability and that of the atmosphere. The main method described in the literature to establish the link between the two fields is to seek time correlations between significant wave height and the appropriate climate index, chosen to adequately represent the major fluctuations of atmospheric circulation. This method has proved to be useful to link the knowledge associated with these climate modes to the wave height variability. However, a question that remains unanswered is, what part of the inter-annual to decadal $H_s$ variability can be explained by the atmospheric variability on the same time-scale? As shown in Figure 8, some regions, particularly off the coast of western Europe, are well-explained by the main climate modes, NAO and EA. Others, such as the Mediterranean sea, the Gulf of Mexico, or low latitudes regions remain poorly explained by these modes. Whether it is because the relevant modes of SLP or atmospheric variability for these regions are missed, or because the inter-annual variability of $H_s$ in these regions can be accounted for only when the atmospheric variability of shorter time-scales is used, is not known. Since the relationship between $H_s$ and wind intensity is not linear and further complicated by the effect of swell propagation, there is indeed no reason to believe that the $H_s$ inter-annual and longer variability is entirely explained by the atmospheric variability with the same time resolution.

It is now well-established, with both observations and models, that the NAO is the main driver of the NA inter-annual winter sea state variability. A positive NAO is indeed associated with a positive $H_s$ anomaly in the eastern half of the basin at latitudes higher than 45° N and with a negative $H_s$ anomaly in the western part of the basin and centred at 30° N. The three different statistical methods presented in Figure 5 (namely SLP PC and $H_s$ PC correlation, SLP and $H_s$ projection on the NAO index and CCA analysis) give similar results, and therefore give confidence in the robustness of the link between NAO and sea state variability. The EA, which is the second mode of variability of winter SLP in the Northern hemisphere, is associated with a positive anomaly of $H_s$ centred at 45° N. Unlike what is found for the NAO, the different statistical methods give somewhat different results and suggest that neither the EA index, nor the second PC of SLP over the NA, entirely capture the second mode of winter $H_s$ variability. It is however possible to construct a new index, based on the second and third PCs of winter SLP over the NA, that accurately describe the second mode of variability of $H_s$ (Figure A1). This emphasises the fact that the link between both SLP and $H_s$ long-term variability is not straightforward. We have, then, reviewed the effects of the main climate indices that are known to affect the NA $H_s$ winter climate: SCA, EA/WR, PNA, TNH, ENSO, SAM, and WEPA. Note that this last mode (WEPA) follows a different approach because it is specifically computed to maximise the explained variance of $H_s$ off the western coast of Europe.

Most of the climate modes studied in this review represent the SLP variability. The link between $H_s$ and the wind intensity modes on these time-scales is almost everywhere and stronger than between $H_s$ and the climate modes (not shown). Statistical models that links $H_s$ with the atmospheric field generally use both the SLP and the SLP gradients which are proportional to the geostrophic wind intensity, for example, in Wang et al. [98]. One possible way to improve our understanding of the atmospheric and $H_s$ co-variability at inter-annual time-scale is thus to use a combination of these atmospheric fields rather than the SLP alone. It is also interesting to note that when the SLP EOFs are computed over the North Atlantic rather than over the whole Northern Hemisphere as for the NOAA climate modes, the $H_s$ variance explained by these modes increases (see Table 2). However, the advantage of the use of standard climate indices rather than other statistical methods, such as the ones presented in Section 2, is that each of these climate indices and associated pattern have generally been well-studied in terms of dynamic, climate impacts, and variability. However, they are not explicitly optimised to represent the maximum $H_s$ variance, and thus might miss some $H_s$ variability. Indeed, Figure 8 shows that the JFM $H_s$ is not completely captured by the climate modes in substantial regions of the NA. The regions where the variability is well-captured by the climate modes are located off the coast of western Europe, along the NA storm track. Whether the relatively weak explained $H_s$ variance in some regions can be attributed to the fact that the appropriate climate modes are missed or to the absence of link between slow variability of SLP and $H_s$ in these regions is not known and would necessitate further studies.

In the last section, we reviewed our current understanding of the sea state projection under future climate change in the NA. Three main sources of sea state climate change are identified for the NA. The first deals with changes in extra-tropical atmospheric circulation, such as the expected storm track poleward shift. In response, $H_s$ is expected to decrease at the end of the century compared to present day climate. However, it was also shown by [100] that these projections are sensitive to the accuracy of the simulated atmospheric circulation at the end of the century for which low agreement between different climate models is generally found [101]. The second identified source of sea state climate change is the increase in tropical cyclone activity by the end of the century that lead to increased extreme wave heights in the NA [116]. The third identified source is the wave–ice interactions. Because wave are attenuated by the presence of ice, the projected decrease of sea ice extent in Arctic waters will mechanically increase the wave height in this region. The increasing extent of ice-free waters is also expected to lead to an increase of wave height because of a larger fetch [127].