1 Introduction

Periodic, short- and medium-term, as well as quickly noticeable changes in ocean level are associated with a number of processes and phenomena occurring in the atmosphere and hydrosphere. Meteorological conditions such as pressure, wind strength and direction, precipitation and evaporation cause irregular and short fluctuations in sea and ocean levels. Because the oceans have an extremely high thermal capacity, when compared to the atmosphere, the ocean temperatures fluctuate seasonally much less than the atmospheric temperature. For the same reason, when air blows over the water, its temperature tends to reach the temperature of the water rather than vice versa. Thus, maritime climates are generally less variable than climates of regions located in the interiors of the continents. The pattern of atmospheric circulation largely determines the pattern of oceanic surface circulation, which in turn determines the location and amount of heat that is released to the atmosphere. Also, the pattern of atmospheric circulation partly determines the location of clouds, which influences the locations of heating of the ocean surface.

One of the most predominant atmospheric phenomenon on the Earth is the North Atlantic Oscillation (NAO) which has been subject of scrutiny over years (e.g. Walker 1924; Bjerknes 1964; Barnston and Livezey 1987; Hurrell 1995; Marshall and Schott 1999; Hurrell et al. 2003). The NAO phenomenon forms a dominant mode of winter climate variability in the North Atlantic region, ranging from central North America to Europe, and extending its impact to northern Asia (Walker and Bliss 1932; Van Loon and Rogers 1978; Rogers and van Loon 1979). The NAO influences cover parts of the East Atlantic and West Atlantic patterns, as originally shown by Wallace and Gutzler (1981) for the winter period on the Northern Hemisphere. The NAO phenomenon consists of balance changes between relatively stable northern and southern atmospheric pressure centers (IPCC 2007). The first one is located over Greenland and called the Icelandic Low (IL), while the second one, known as the Azores High (AH), is located in mid-latitudes of the North Atlantic, i.e. between 35 and \(40^{\circ }\)N (Wanner et al. 2001). Through the fluctuations in the strength of the north-south dipole, NAO controls the strength and direction of westerly winds and storm tracks across the North Atlantic (Feistel et al. 2008). Hence, it is associated with large-scale modulations of the normal patterns of zonal and meridional heat and moisture transport (Hurrell 1995). Walker (1923) first noticed a tendency towards a simultaneous strengthening or weakening of IL and AH, and he named these changes as NAO. The positive (negative) mode of NAO occurs when IL and AH are simultaneously strong (weak). This teleconnection pattern is generally defined as the normalized sea level pressure difference between Ponta Delgada (Azores) and Akureyri (Iceland). It appears in the first eigenvector of sea level pressure field in the Northern Hemisphere (Barnston and Livezey 1987; Rogers 1990), accounting for roughly one third of the variability in sea level pressure over the North Atlantic Ocean (Deser and Blackon 1993).

To quantify the anomalous difference between the polar low and the subtropical high during the boreal winter (DJFM—from December to March), the NAO index (NAOi) was created. The first constructed NAOi involved temperature and pressure data from Europe and North America, and was calculated using a complex and multi-staged statistical procedure (Walker and Bliss 1932). The index reveals interannual variability, but it also has a tendency to persist in one phase over several years. The monthly NAOi is calculated on a basis of decadal sea level pressure measurements performed at two selected meteorological stations. The most popular versions of NAOi are: (1) the Rogers NAO index (between Ponta Delgada (Azores) and Akureyri (Iceland) stations); (2) the Hurrell NAO index (between Lisbon (Portugal) and Stykkishólmur (Iceland) stations); and (3) the Jones NAO index (between Gibraltar and Stykkishólmur stations).

Predictions of NAO are useful for the economy and environment, and therefore several dedicated forecasting methods and services have been developed. Based on the principal components of sea surface temperature (SST) changes, Saunders and Qian (2002) elaborated an empirical model to predict winter NAOi. Eshel (2003) used the North Pacific surface pressure to produce 21-month forecasts of NAO based on various autoregressive models. The Climate Prediction Center of the National Oceanic and Atmospheric Administration (NOAA) provides daily, weekly, 10-day and 14-day forecasts of NAOi along with the mean value of NAOi as well as the correlation coefficients between the observations and all the forecasts for last 120 days. The NAOi forecasts are also based on the Global Forecasts System (GFS) and the associated ensemble mean model, produced by the National Centers for Environmental Prediction (NCEP). Yet another institution which provides long-term forecasts of winter NAOi is the European Centre for Medium-Range Weather Forecasts (ECMWF). Its Monthly Forecasting (MOFC) system is based on an ensemble of 51 coupled ocean–atmosphere integrations which produce 32-day forecast every two weeks (Vitart 2004). Recently, Wang et al. (2017) developed a skillful empirical model to predict winter NAO variability, and Dobrynin et al. (2018) utilized a dynamical ensemble model to produce seasonal NAO forecasts.

Most approaches developed to forecast NAOi use predictors like SST, sea level pressure or ocean heat content, however none of them is based on SSH or sea level anomaly (SLA). It is widely known that sea level adjusts to changes in barometric pressure due to the inverted barometer effect (Wunsch and Stammer 1997; Ponte and Gaspar 1999). Sea level response to atmospheric pressure was the topic of several research, e.g. conducted by Garrett (1983), Ducet et al. (1999) and Mathers and Woodworth (2001). The SLA variation reflects changes in atmospheric patterns (Chen et al. 2014; Iglesias et al. 2017), therefore, SLA may be used as an explanatory variable to construct predictions of NAOi.

The correlation between variability of surface ocean topography and NAO phases was discussed by Esselborn and Eden (2001). The first Empirical Orthogonal Function (EOF) component of the observed seasonal SSH anomalies was found to explain 34% of variance. Pingree (2002) studied the relationship between the NAO index and the SSH along a meridional transect at 35\(^{\circ }\)N in the North Atlantic and found vivid, but delayed (one year up) correlations. Investigations into the interannual variability of SSH conducted by Cromwell (2006) led the author to the conclusion that the first principal component of SSH, which is also considered in this paper, is anti-correlated with other commonly used index over North Atlantic–East Atlantic (EA) pattern index. Similarly to Esselborn and Eden (2001), the author applied complex EOF analysis and then used wavelet analysis to determine time-varying spectral density content of the principal components. Other elements of the climate system that may improve predictability of NAO are discussed by Smith et al. (2014).

The main objective of this paper is to use SLA as the input data to issue long-term forecast of NAOi. Herein, we present, verify and validate a new approach for predicting NAOi, which integrates stochastic models implemented within the Prognocean Plus system (Niedzielski and Miziński 2013; Świerczyńska et al. 2016) with the SSH–NAOi relationship presented by Esselborn and Eden (2001). The research hypothesis is the following: “it is possible to skilfully predict NAOi with the 3-month lead time based solely on the prognoses of sea level change”.

2 Data

Three sets of data are used in this paper: SLA, SSH and winter Hurrell NAOi. The relationship between SLA and SSH is described by the formula:

$$\begin{aligned} {SLA} = {SSH} - {MSSH}, \end{aligned}$$
(1)

where MSSH is the Mean Sea Surface Height.

Both SLA (variable “sea_surface_height_above_sea_level”) and SSH (variable “sea_surface_height_above_geoid”) datasets were obtained courtesy of the Copernicus Marine Environment Monitoring Service (CMEMS) from the repository called “Global Ocean Gridded L4 Sea Surface Heights and derived variables reprocessed (1993–ongoing)”. Both datasets were processed by DUACS (Data Unification Altimeter Combination System) which is a part of the CNES (Centre National d’Etudes Spatiales) multimission ground segment SSALTO (Segment Sol multimissions d’ALTimétrie, d’Orbitographie et de localization précise). It processes data from all available altimetric missions (Jason-3, Sentinel-3A, HY-2A, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, Topex/Poseidon (T/P), ENVISAT, GFO, ERS1/2) and provides a consistent and homogeneous catalogue of products. For the SSH data a few altimeter corrections and processing algorithms were used so that the data are devoid of: the ionospheric effect, wet and dry troposphere, ocean-, pole- and solid earth tides or the sea state bias. The multimission gridded SLA time series is computed from SSH referenced to a twenty-year mean.

Spatial resolution of the data is \(1/4^{\circ }\times 1/4^{\circ }\). After Esselborn and Eden (2001), we utilize only gridded SLA and SSH data, spatially limited to 15–65°N and 0–80°W. Due to data availability (the first altimetric satellite T/P was launched in 1992), we consider data collected between 1993 and 2016, aggregated to months (12 per year) and seasons (4 per year).

The monthly and seasonal gridded data where produced by performing simple averaging for each grid of original gridded (multidimensional) time series of SLA/SSH. We created a dataset for each grid (one dimension) that consisted of SLA/SSH values for a winter month or three winter months (season), and calculated the corresponding mean value. We conducted this process for each year.

The monthly and seasonal station-based Hurrell NAOi time series were obtained from Climate Data Guide (Hurrell et al. 2017). For the purpose of this study, we consider only winter months, i.e. December/January/February (D/J/F), and winter seasons (DJF).

3 Methods

3.1 Prognocean Plus

Prognocean Plus is a real-time system for computing sea level prediction. It is an upgraded version of the Prognocean system described in detail by Niedzielski and Miziński (2013). Prognocean Plus has been developed as a new generation domain-specific service within the PL-Grid infrastructure for Polish science (http://www.plgrid.pl/en, access date: 29/03/2018). By 2018, the system was available through the dedicated real-time online service (https://prognocean.plgrid.pl/, access date: 24/07/2018). Currently, it is not operational. Prognocean Plus produced SLA forecasts based on SLA time series obtained courtesy of CMEMS. The system produced short- (up to 14 days) and long-term (up to 12 weeks) predictions of global altimetric gridded SLA time series, updated daily (Świerczyńska et al. 2016). The spatial resolution of the SLA forecasts was \(1/4^{\circ }\times 1/4^{\circ }\). Along with the prognoses, gridded statistical maps of prediction performance were produced, showing spatial distributions of Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) of forecasts. The Prognocean Plus system was implemented in R (language and environment for statistical computing and graphics) and used the Message Passing Interface (MPI) computations with OpenMPI (https://www.open-mpi.org/, access date: 24/07/2018).

3.2 Stochastic models in Prognocean Plus

The Prognocean Plus system computed the predictions of time-variable ocean topography using data-based methods. The legitimacy of the use of a statistical approach for modelling physical phenomena is widely discussed in the literature (McPhaden et al. 2006; Clark 2008; Sarachik et al. 2010). Moore et al. (2013) constructed the semi-empirical global sea level projections, compared them with the process-based projections of SLA, and arrived at the conclusion that statistical and dynamical models have comparable predictive abilities.

The following statistical models were used in the Prognocean Plus system:

  • the deterministic polynomial-harmonic model (Polynomial-Harmonic, PH),

  • a combination of the PH model and the stochastic autoregressive model (AutoRegressive, AR), briefly called PH+AR,

  • a combination of the PH model with the stochastic threshold autoregressive model (Threshold AutoRegressive, TAR), known as PH+TAR,

  • a combination of the PH model and the multivariate autoregressive model (Multivariate AutoRegressive, MAR), briefly called PH+MAR.

The PH model consists of linear trend and five basic non-tidal harmonic oscillations (annual, semi-annual, 120-, 90- and 62-day). These deterministic components can be described as a sum of polynomial or harmonic nonlinear models, the parameters of which can be estimated using the Gauss–Newton algorithm within the weighted least-squares approach. The PH-based predictions are extrapolations of the model. The differences between the observed SLA data and the simulations from the PH model are residuals. They were subsequently modelled with the use of stochastic methods, as described in detail by Niedzielski and Miziński (2013) and Świerczyńska et al. (2016).

The physical background of the PH model with five harmonic components is based on the papers by Kosek (2001), Niedzielski and Kosek (2009) and Kosek et al. (2016) who, using Fourier transform band pass filter, identified the most energetic periodic oscillations in sea level as a function of geographic location. Apart from the annual oscillation, which is associated with the annual cycle and is found to be broadband (Kosek et al. 2016), Kosek (2001) found that semi-annual, 120-, 90- and 62-day components are more energetic in the Northern Hemisphere than in the Southern Hemisphere. In addition, the presence of these oscillations has been identified in sea level variation in the North Atlantic Ocean and some of its seas, mainly in the east-equatorial Atlantic Ocean, the North American Basin, the North Sea and the Baltic Sea (Kosek 2001). Also, the annual oscillation is found to be energetic in parts of the North Atlantic Ocean, as shown for the Caribbean Sea (Niedzielski and Kosek 2009) and the western North Atlantic (Kosek et al. 2016). The 62-day component is an alias-type oscillation driven by observing the M2 and S2 tide waves by altimetric satellites with 10-day sampling interval (Wagner et al. 1994; Kosek 2001). The polynomial component in the PH model has been assumed to be linear because in a short time scale (a few years) there has been no meaningful acceleration in sea level change observed by T/P, Jason-1 and Jason-2 satellites (Niedzielski and Kosek 2011).

The AR(p) solution is widely known in the theory of one-dimensional time series (Brockwell and Davis 2002; Cryer and Chan 2008). Here, the model \(Y(t,\varphi ,\lambda )\) is fitted to the residuals of SLA data at a given location \((\varphi ,\lambda )\) to improve prediction skills, especially for a few first prediction steps. The required model parameters are calculated automatically using the Yule-Walker method (Brockwell and Davis 1996). The choice of the order p is based on minimizing the Akaike Information Criterion (AIC) (Akaike 1971). The optimal AR-based linear prediction for r steps into the future is calculated recursively on a basis of the autoregression equation.

The TAR solution is the self-excitive modified version of AR model, which allows to switch between two main equations depending on the randomly determined process in the past. When the threshold is attained, the model is automatically adjusted and the more efficient equation is chosen. Model estimation and calculation of forecasts are carried out in accordance with the procedure described in Cryer and Chan (2008), using the Time Series Analysis (TSA) package in R environment.

The MAR approach is a vector version of the AR model. It is based on cross-correlations within a multivariate time series of SLA (Lardies 1996). Here, the MAR model is implemented in a moving rectangular window (width of 1.25 degrees of longitude and height of 0.75 degrees of latitude), centred in the grid cell where prediction is needed.

3.3 Relationship between sea level changes and NAO

The NAO phenomenon reveals interseasonal, interannual (Hurrell 1995; Feldstein 2003; Benedict et al. 2004; Woollings et al. 2015) and even decadal or multidecadal (Stephenson et al. 2000; Visbeck et al. 2001; Bell and Chelliah 2006; Olsen et al. 2012) variability. Bjerknes (1964) described the ocean-atmosphere interactions with respect to North Atlantic climate variability. He stated that heat exchange in the atmosphere and sea surface temperatures play an important role in producing NAO. Based on that, Van Loon and Rogers (1978) found that there exists a significant correlation between atmospheric circulation and sea surface temperature (SST). In addition, many authors have investigated a possible role of Atlantic sea surface temperature anomalies in driving change in the surface pressure patterns resembling NAO (Ratcliffe and Murray 1970; Peng et al. 1997; Czaja and Frankignoul 1999; Rodwell et al. 1999; Rodwell and Folland 2003). Furthermore, Lozier et al. (2008) showed that changes in ocean heat content can be connected to the decadal trend in the NAO signal. The newest research conducted by Yeager et al. (2012) offered the knowledge about feedback between the NAO index and ocean heat content.

It is known that the NAO signal is present in SSH or SLA time series from the North Atlantic region (Cabanes et al. 2006). The methodical bases for quantifying the relationship between SSH and NAOi were proposed by Esselborn and Eden (2001) who utilized the Principal Component Analysis (PCA). It is a statistical procedure that allows us to transform a two-dimensional set of dependent variables to a one-dimensional set of independent variables (principal components) using EOF. The principal components are dominant patterns. The new variables represent the same amount of information as the original variables, and their total variance remains the same. However, it is redistributed among the new variables in the most “unequal” way: the first of them explains not only the majority of variance among the new variables, but also the majority of variance that a single variable can possibly explain (Jolliffe and Cadima 2016).

To perform PCA we used the “prcomp” function (Venables and Ripley 2002; R Development Core Team 2008) from the “stats” package in R environment. The above-mentioned package and the “prcomp” function are commonly used in the environmental studies when there is a need to conduct PCA (Chandler and Scott 2011; Tapolczai et al. 2018; Will and Kitaysky 2018; Al-Karkhi and Alqaraghuli 2019). The calculation is carried out by a singular value decomposition (SVD) of the (centred and possibly scaled) data matrix, not by using eigenvalues on the covariance matrix. Unbiased variances are utilized. The use of the SVD method is preferred if numeral accuracy is required.

Esselborn and Eden (2001) used the T/P altimeter SSH data provided by Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO), spanning the time interval from 1993 to 1998. The spatial resolution of the data was \(1^{\circ }\times 1^{\circ }\), and the authors limited their spatial extent to 15–65°N and 0–80°W. The gridded SSH altimetry data were referenced to the mean seasonal SSH along the satellite tracks to calculate anomalies. The seasonal Hurrel NAOi was used, and therefore in order to unify sampling intervals, the SSH time series was averaged to four seasons per year: December to February (DJF), March to May (MAM), June to August (JJA) and September to November (SON). To compare the seasonal SSH variability with NAOi, the PCA method was used to extract the first principal component of gridded SSH data (PC1SSH) for each season. Esselborn and Eden (2001) extracted only winter Hurrel NAOi values which were subsequently compared with PC1SSH, and the two variables were found to be correlated.

Although Esselborn and Eden (2001) noticed limitations of their approach due to a small number of altimetric observations available at that time, they provided possible explanation of the relationship between SSH and NAO. They presented the reasoning that the NAO-related wind stress curl controls anomalies of ocean circulation and, as a consequence, leads to the interannual changes in heat content of the North Atlantic Ocean.

3.4 New approach for forecasting NAOi

In this paper, we explore two methods for preparing input time series to perform PCA and to extract PC1 of SLA/SSH for winter seasons or months. The approaches are graphically presented in Fig. 1. The sketch consists of two parts, where part “a” presents the use of the two methods for seasonal SLA/SSH data, while part “b” corresponds to monthly SLA/SSH data. The first method, known hereinafter as “method 1”, is based on the complete time series of SLA/SSH (all seasons and months), which become input data for PCA. After analysis, the first principal components for winter seasons and months are extracted. The second method, called “method 2”, is based on subsets of SLA/SSH time series, in which only winter seasons or winter months are selected to perform PCA.

The procedure of producing the NAOi forecasts is graphically presented in Fig. 2. At least 20 years (1993–2013) of the observed SLA data are needed to fit a linear regression model. We allow four temporal resolutions (Fig. 1), and the resulting time series are processed by “method 1” and “method 2”.

In the next step, following Esselborn and Eden (2001), PCA is performed on the observed SLA time series. In “method 1”, the calculated PC1s of winter seasons or months must be extracted. Based on the calculated PC1s (4 types) and time series of NAOi (for winter seasons or months only) we can construct a regression model describing the relationship between SLA and NAOi using the following formula:

$$\begin{aligned} NAOi=a\times PC1SLA +b, \end{aligned}$$
(2)

where NAOi is expressed by the step-adjusted linear regression model describing the Hurrell winter (seasonal or monthly) NAO index, a and b are model coefficients and PC1SLA is time series of one of the four possible sets of principal component time series (PC1SLA computed using “method 1” for seasonal and monthly data—two datasets, and PC1SLA computed using “method 2” for both temporal resolutions—additional two datasets). The above-mentioned regression model is recalculated every time, when new prognoses of NAOi are needed. The length of input SLA dataset is equal to the number of the available observations at the time of forecast determination. Hence, a new set of PC1SSH data is produced at every time step. The regression model is necessary to obtain coefficients a and b, which are used in the subsequent step of the procedure.

After fitting the regression model based on the observed SLA time series, we extend this dataset by attaching a newly-predicted SLA gridded map (one map containing SLA prognoses for the next DJF season or three maps consisting of SLA predictions for the next three winter months). Such gridded predictions are computed by the Prognocean Plus system using one of three stochastic models (PH+AR, PH+TAR or PH+MAR). The observed SLA data together with the attached newest gridded forecast, prepared using the above-mentioned methods 1 and 2, becomes input data for PCA. As previously, in the case of “method 2”, it is necessary to extract PC1s for winter seasons and months.

The last step of the procedure for computing prognoses of NAOi is the use of the regression model, based on the previously calculated coefficients and the extended time series of PC1SLA, utilizing the formula:

$$\begin{aligned} p(NAOi)=a\times extPC1SLA +b, \end{aligned}$$
(3)

where p(NAOi) is the step-adjusted linear regression model predicting Hurrell winter (seasonal or monthly) NAO index, a and b are the already estimated model coefficients obtained to express NAOi, extPC1SLA is a time series composed of one of the four possible sets of principal component data (the same procedure as described above, but with the original time series of SLA extended by attaching the nearest winter prediction). The last component of extPC1SLA is calculated for the predicted gridded SLA data, attached to the observational dataset. In addition to the regression model p(NAOi), we determine its modified version \(p'(NAOi)\) produced by applying the empirical correction. It is the p(NAOi) model multiplied by a coefficient defined as a ratio between standard deviation of the NAOi time series and standard deviation of extPC1SLA data, namely:

$$\begin{aligned} p'(NAOi)=p(NAOi)\times \frac{\sigma _{NAOi}}{\sigma _{extPC1SLA}}, \end{aligned}$$
(4)

where p(NAOi) is the step-adjusted linear regression model from the first approach to predict NAOi, \(\sigma _{NAOi}\) is standard deviation of NAOi dataset, \(\sigma _{extPC1SLA}\) is standard deviation of PC1SLA time series computed based on the observed SLA with attached forecast(s). This procedure results in the same variances of NAOi time series and the NAOi data reconstructed from the \(p'(NAOi)\) regression model.

4 Results

In this study, we carried out two analyses: (1) we performed the calculations along the lines of the paper by Esselborn and Eden (2001) on a basis of SSH and SLA datasets, which allowed us to confirm the strong relationship between NAOi and SSH/SLA for the 20-year time span and to quantify the relation in question; (2) based on the dependency we used the associated regression models along with the real-time SLA predictions computed by the Prognocean Plus system. The approach was used to predict winter NAOi for four boreal winters from 2013/2014 to 2016/2017.

The length of time series included in the analyses varies between 21 and 72, depending on the selected approach, the temporal resolution (seasons vs. months) and the prognosis itself (into which season it hits). The shortest time series is for the prediction of the 2013/2014 season computed using seasonal data. In contrast, the longest time series is utilized for the prediction of the 2016/2017 season calculated using monthly data.

4.1 PC1 of SLA/SSH and NAOi

The relationships between PC1 of SSH/SLA and winter NAOi are presented in Fig. 3. For each graph we present the value of Pearson’s correlation coefficient (r) with its statistical significance (p-value). In the literature, a general appreciation for the level of correlation may include the following values: 0–0.09 indicates no correlation, 0.1–0.3 indicates a small correlation, 0.3–0.5 indicates a medium correlation, and 0.5–1 indicates a large correlation (e.g. Cohen 1988). When interpreting small correlations it is worth referring to the paper by Kozak et al. (2012) who argue that “low value of the correlation coefficient can represent incredibly strong association, while a high value can still represent weak correlation” and recommend the analysis of the context of the association.

In their study, Esselborn and Eden (2001) performed the analysis on short seasonal SSH time series using “method 1” of preparing dataset for PCA. The results of similar calculations, but using much longer time series, are shown in Fig. 3a. Their version based on monthly data is presented in graph “b”. The corresponding analysis performed using “method 2” gives higher values of Pearson’s correlation coefficient. The strongest correlation (\(r = 0.56\)) appears between DJF Hurrel NAOi and PC1SSH for DJF season (graph “c”). This confirms that there exists the strong relationship between SSH and NAOi, which can be used in the process of forecasting NAOi. The Prognocean Plus system produces predictions of SLA, therefore we conducted the same experiment with SLA data (Fig. 3e–h). The correlation coefficient between NAOi and PC1SLA in all cases is even higher than between NAOi and PC1SSH. The highest correlation value (\(r=0.68\)) occurred between winter seasonal NAOi and PC1SLA for DJF season.

4.2 Seasonal forecasts

Using two regression models described in Sect. 3.4 which utilize PH+AR-based SLA predictions, we produced four seasonal prognoses of winter Hurrell NAOi. They are presented in rows of Fig. 4. In order to keep the calibration period long enough (at least 20 years), the first NAOi forecast is issued into the DJF season 2013/2014 (graphs “a” and “b”), while the last available prediction is for the DJF season 2016/2017 (graphs “g” and “h”). Each graph consists of black solid line reflecting values of the observed NAOi, solid grey (approach 1) and dark grey (approach 2) lines which are NAOi reconstructions from two regression models. Dashed lines and squares with the greyness intensity corresponding to the models indicate values of predicted DJF NAOi. Because of the stepwise recalibration of each model, every time when forecast is produced, these parts of graphs are different as the calibration dataset extends. The model recalibration shows that both models respond properly to changes of NAOi. The presented statistics (Table 1) show that the values of Pearson’s correlation coefficient fall around 0.5, and the values of MAE or RMSE statistics fluctuate around 1.5. The empirically-corrected regression model (Eq. 4) reveals better performance and, because of the same variance as NAOi time series, it better reflects the amplitudes of NAOi changes.

Although Fig. 4 focuses on the PH+AR model only (it has the smallest misfit according to the conducted performance analysis), the skills of the remaining two deterministic-stochastic models (PH+TAR and PH+MAR) are juxtaposed in Table 1. It is apparent from the table that choice of the SLA prediction model from those available in the Prognocean Plus system has no impact on the performance of NAOi forecasts. Higher correlations are found for “method 1” (in all variants), but at the same time prediction errors are similar or worse than in case of predicting winter months of NAOi. The PH+AR model reveals the lowest values of MAE and RMSE (the lowest values of these statistics are 1.18 and 1.55, respectively, for method 1, approach 1, season 2013/2014). However, it should be clearly noted that all values of prediction performance statistics are very similar. They are often non-recognizable, with differences occurring at the 4th or 5th decimal place.

4.3 Monthly forecasts

We carried out a similar exercise, based on monthly data. In this setup, NAOi for three individual winter months is predicted. The results are presented in Fig. 5, the structure of which corresponds to Fig. 4.

The empirically-corrected regression model (Eq. 4) resolves the amplitude of NAOi better than the uncorrected model. Forecasts issued into D/J/F months for winter seasons 2015/2016 and 2016/2017 are accurate, with the NAOi amplitude correctly resolved. The usefulness of the empirical correction is well seen in Fig. 5b, where the corrected prediction curves accurately hit the true NAOi variability (predicted values correctly fitted the fluctuations in NAOi values, with absolute error smaller than 0.7). To convince the reader why the empirical correction was needed, the authors prepared Fig. 6 where the scatter plot of the mentioned case from Fig. 5b is presented. The amplitudes of the uncorrected NAOi changes determined through the reconstruction and prediction vary between − 0.5 and 0.5 (gray rectangles), while the corrected ones (dark gray squares) are between − 4 and 4 which much better reflects the real NAOi variability (− 5.5 to 4).

Table 1 shows that for monthly predictions of NAOi, the use of “method 2” leads to the highest correlation coefficient values which exceed 0.59 (winter season 2015/2016). Based on the prediction errors, “method 2” is found to produce the most accurate forecasts of NAOi for D/J/F months in winter 2015/2016. In this exercise, 69 values of PCA1 and NAOi are used.

4.4 Discussion

It is known that correlations between NAOi and its short-term predictions are very high, i.e. \(\ge 0.8\) up to 6–7 days into the future (Johansson 2007; Vitart 2014). They decrease to approximately 0.35 when lead time increases to 15 days (Johansson 2007).

Although our validation of monthly and seasonal forecasts, due to data availability, is limited to four winters (2013/2014, 2014/2015, 2015/2016, 2016/2017), we believe that the numerical exercise serves as a feasibility study to confirm the usefulness of the new methodical approach presented in this paper. Our method uses altimetric sea level change of northern Atlantic Ocean as a predictor of NAOi, and it is fundamentally different from the work of Saunders and Qian (2002) who employed SST variability to forecast NAOi. Even though the statistical robustness of our NAOi prediction validation is much worse than the robustness of the results published by Saunders and Qian (2002)—i.e. 20 years for comparing the reconstructed models with true data and 4 winters for comparing true data with predictions (our study) versus 51 years of cross-validation and 15 winters for comparing predictions with independent forecasts (Saunders and Qian 2002)—certain findings are similar. The two studies aim to anticipate NAOi in the 3-month DJF seasons. In the SST-based NAOi prediction experiment, correlations were found to vary between 0.47 and 0.63 for the cross-validation period (51 years), while in our SLA-based NAOi prediction exercise correlations between true and reconstructed NAOi (20 years of calibration with stepwise extension by one year) were equal to 0.51–0.69. Thus, SST- and SLA-based NAOi predictions have comparable skills.

Similarly, Saunders et al. (2003) obtained comparable correlations of 0.60–0.62 in numerical hindcast experiments for intervals 1972–2001 and 1987–2001. Also, Johansson (2007) found that statistically significant correlations between monthly NAOi predictions and data vary from 0.36 to 0.63, depending on month. Lower and often insignificant correlations (max 0.42) were found for seasonal predictions (Johansson 2007). For longer lead times (year-to-year) Kim et al. (2012) obtained much lower correlations (approx. 0.2), while Fan et al. (2016) presented slightly higher values (approx. 0.38–0.53). The recently developed statistical method by Wang et al. (2017) uses the multiple linear regression to forecast the NAO variability, with SST, 70 hPa geopotential height and sea ice concentration as main predictors. The authors reported high correlations between data and NAO predictions. Also, high correlations, even exceeding 0.8 for selected cases, were obtained by Dobrynin et al. (2018) who used the dynamical ensemble approach.

The performance of our NAOi predictions is either controlled by the relationship between PC1SLA and NAOi, described by Esselborn and Eden (2001), or driven by the skills of the SLA predictions issued by the Prognocean Plus system (Świerczyńska et al. 2016). Figure 3 shows that our PC1SSH–NAOi analysis for 1993–2014 leads to findings which are similar to what Esselborn and Eden (2001) discovered. It is also apparent from Fig. 3 that correlations between PC1SLA and NAOi are higher than correlations PC1SSH–NAOi.

Our approach is purely statistical, since both the relationship borrowed from Esselborn and Eden (2001) and three models implemented in the Prognocean Plus system are empirical. There are other attempts to build such statistical prediction schemes to anticipate NAOi variability. For instance, Eshel (2003) focused on the Pacific–Atlantic teleconnection and built long-term 15-month NAOi predictions using sea level pressure, not only in northern Atlantic, but in particular in northern Pacific. The correlations between one-year NAOi predictions and observations were found to be approximately equal to 0.45 for 1925–2002 cross-validation experiment.

It is rather difficult to compare the skills of our predictions of NAOi with the corresponding forecasts based on other models run by different centres. For instance, the historical ECMWF prognoses are not archived (e-mail communication dated 01/08/2018), while the fortnight ensemble predictions computed by NOAA CPC are publically available solely as graphical plots without access to historical forecasts (www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml, access date: 09/08/2018). The comparison of prediction errors and skill measures must therefore be limited to published work in which different timespans are employed.

5 Conclusions

We developed a new method for forecasting winter Hurrell NAOi based solely on SLA predictions calculated by the Prognocean Plus system. Using the procedure proposed by Esselborn and Eden (2001) we demonstrated a significant relationship between SLA data, aggregated seasonally and monthly, and seasonal and monthly winter Hurrell NAOi (Pearson’s correlation coefficient was equal to 0.69 and 0.59, respectively). Forecasts of NAOi included winter season from 2013/2014 to 2016/2017 (four seasonal and four 3-monthly predictions).

The main finding of this paper is that it is possible to produce long-term (up to 3-months) predictions of NAOi based solely on SLAs. The study shows that using predictions of SLA from the Prognocean Plus system we are able to anticipate next values of winter NAOi with prediction errors of MAE and RMSE not exceeding 1.5 and 1.9, respectively. The key findings are the following:

  • forecast accuracy using regression models is satisfying, and the best performance was achieved using SLA forecasts from PH+AR model implemented in the Prognocean Plus system,

  • differences between performance statistics indicate that the choice of SLA forecasts model has less impact than the choice of method for preparing SLA time series for PCA,

  • approach 1 of predicting NAOi fails to resolve NAOi amplitudes,

  • approach 2 (empirically-corrected model outputs), properly resolves NAOi amplitudes,

  • the calibration (reconstruction) of approaches 1 and 2 is characterized by high agreement with observations, while the predictions based on the reconstructions reveal only medium or low agreement with observations.

Fig. 1
figure 1

Scheme of preparing seasonal (part “a”) and monthly (part “b”) SLA/SSH datasets for PCA using two methods. First method uses complete seasonal or monthly data (respectively 4 gridded map per year and 12 gridded maps per year) as an input data for PCA. After analysis, first principal components for winter seasons/months are extracted. In second method, winter seasons/months are already extracted from SLA/SSH dataset, hence shorter time series are input data for PCA. After analysis, no extraction is needed

Full size image
Fig. 2
figure 2

Scheme presenting several steps in procedure of predicting winter Hurrell NAOi based solely on SLA observed and forecasted data using linear regression model based on time series of PC1SLA

Full size image
Fig. 3
figure 3

Scheme presenting relationship between principal components of winter seasons (left column) and winter months (right column) of SSH (graphs a–d) and SLA (graphs e–h) observed data computed using method 1 and 2 (alternating rows) values of PC1s are presented as grey squares with black outline, values of NAOi are presented as black triangles. On the bottom of each graph, value of Pearson’s correlation coefficient (r) with statistical significance (p-value) are presented

Full size image
Fig. 4
figure 4

Set of graphs presenting performance of two regression models in reconstructing and predicting winter seasonal Hurrell NAOi. In left column all regression models are based on the PC1SLA obtained using method 1, results using method 2 are in the right column. Each row shows results of reconstruction and prediction of one DJF season, where first row is forecasting Hurrell NAOi for winter season 2013/2014 and the last one—2016/2017. Observed NAOi is presented as solid black line. Reconstruction of NAOi using approach 1 (first regression model without correction) is an solid grey line, dark grey solid line is a performance of approach 1 (second regression model with correction). Predicted values of NAOi for each winter season are indicated as an filled square with black outline and marked centroid of the figure, joined by dashed-line. Fill of the square correspond to the line colour assigned to the forecasting approach

Full size image
Fig. 5
figure 5

Set of graphs presenting performance of two regression models in reconstructing and predicting winter monthly Hurrell NAOi. In left column all regression models are based on the PC1SLA obtained using method 1, results using method 2 are in the right column. Each row shows results of reconstruction and prediction of D/J/F months, where first row is forecasting Hurrell NAOi for winter months 2013/2014 and the last one—2016/2017. Observed NAOi is presented as solid black line. Reconstruction of NAOi using approach 1 (first regression model without correction) is an solid grey line, dark grey solid line is a performance of approach 1 (second regression model with correction). Predicted values of NAOi for each winter month are indicated as an filled square with black outline and marked centroid of the figure, joined by dashed-line. Fill of the square and colour of the dashed-lines correspond to the solid line colour assigned to the forecasting approach

Full size image
Fig. 6
figure 6

Scatter plot of the case presented in Fig. 5b. Gray triangles represent corresponding values of reconstructed and predicted values of NAOi obtained from approach 1 (before correction) and observed NAOi, where dark grey squares indicate corresponding values of reconstructed and predicted values of NAOi obtained from approach 2 (after correction) and observed NAOi. Based on the corresponding values, solid lines representing trend lines are drawn. Dashed line indicates the ideal relationship between predicted or reconstructed values of NAOi and observed ones. Correlation coefficient \(r = 0.578\) for both approaches (scaling factor does not change r)

Full size image
Table 1 Statistics of accuracy of approaches 1 (app. 1) and 2 (app. 2) using different temporal resolution (monthly and seasonal), SLA forecasting models from Prognocean Plus (PH+AR, PH+TAR and PH+MAR) and two methods of preparing SLA time series for PCA (method 1 and 2), from winter 2013/2014 till 2016/2017
Full size table