Observed winter salinity fields in the surface layer of the Arctic Ocean and statistical approaches to predicting large-scale anomalies and patterns

Abstract

Significant salinity anomalies have been observed in the Arctic Ocean surface layer during the last decade. Our study is based on an extensive gridded dataset of winter salinity in the upper 50 m layer of the Arctic Ocean for the periods 1950–1993 and 2007–2012, obtained from ~20 000 profiles. We investigate the interannual variability of the salinity fields, identify predominant patterns of anomalous behavior and leading modes of variability, and develop a statistical model for the prediction of surface-layer salinity. The statistical model is based on linear regression equations linking the principal components of surface-layer salinity obtained through empirical orthogonal function decomposition with environmental factors, such as atmospheric circulation, river runoff, ice processes and water exchange with neighboring oceans. Using this model, we obtain prognostic fields of the surface-layer salinity for the winter period 2013–2014. The prognostic fields generated by the model show tendencies of surface-layer salinification, which were also observed in previous years. Although the used data are proprietary and have gaps, they provide the most spatiotemporally detailed observational resource for studying multidecadal variations in basin-wide Arctic salinity. Thus, there is community value in the identification, dissemination and modeling of the principal modes of variability in this salinity record.

INTRODUCTION

The Arctic Ocean is very sensitive to changing environmental conditions. Its surface layer is a key component of the Arctic climate system, which serves as the dynamic and thermodynamic link between the atmosphere and the underlying waters (Carmack, 2000). Thermohaline characteristics of the surface layer are markedly influenced by atmospheric and sea-ice processes, and wind and buoyancy forcing on this important layer ultimately impact the entire upper ocean (Cronin and Sprintall, 2009). The rejection of salt during sea-ice formation strongly impacts upper ocean salinity, so that the stability and development of the ice cover are closely associated with the thermohaline properties of the upper ocean, such as the depth of the mixed layer and halocline. In this context, the Arctic Ocean surface layer is a critical indicator of climate change (Toole and others, 2010).

Here, salinity is chosen as the main characteristic of thermohaline structure variations of the Arctic Ocean surface layer because, at high latitudes, it mainly determines the density structure (Weyl, 1968; Morison and Smith, 1981; Walin, 1985). The thermohaline structure of the Arctic Ocean surface layer has undergone significant changes in recent years (Macdonald and others, 2005). Of particular interest is the great salinification of the surface layer of the Eurasian and Makarov Basins in the early 1990s – a phenomenon unprecedented in the record back to 1950 (Fig. 1). One hypothesis for this is that the increase of Arctic atmospheric cyclone activity in the 1990s led to a large change in the salinity in the Eurasian Basin through changes in river inflow, and increased brine formation due to changes in Arctic sea-ice formation (Dickson, 1999; Polyakov and others, 2008). The other reason for salinification is the influence of Atlantic waters (AW), which by 2007 became warmer by ~0.24°C than they were in the 1990s. Observations show that increases in Arctic Ocean salinity have accompanied this warming as it was associated with significant shoaling of the upper AW boundary and weakening of the upper ocean stratification in the Eurasian Basin as well. That led to facilitated exchange between AW and the upper layer (Polyakov and others, 2010). However, recent observations also show that the upper ocean of the Eurasian Basin was appreciably fresher in 2010 than it was in 2007 and 2008 (Timmermans and others, 2011).

Fig. 1. Temporal changes in salinity averaged over the depth range 5–50 m. Dashed curves show salinities from PIOMAS data. Grids with spatial resolution 200 × 200 km were obtained as the result of interpolation and reconstruction (see section ‘Field reconstruction and interpolation’) of bottled and CTD data.

In addition, there have been observations of surface-layer freshening in the Canada Basin. Jackson and others (2012) emphasized that processes related to warming and freshening the surface layer in this region had transformed the water mass structure of the upper 100 m. With these changes, energy absorbed during summer can enter the deepening winter mixed layer and melt sea ice (Timmermans, 2015).

The problem of variability of Arctic Ocean salinity is challenging from a theoretical perspective. For example, Lique and others (2009) performed an analysis of the variability of Arctic fresh water content informed by a global ocean/sea-ice model. The authors uncovered important spatial contrasts in the influence of velocity and salinity fluctuations on ocean fresh water transport variability. They conclude that variations of salinity (controlling part of the Fram fresh water export) arise from the sea-ice formation and melting north of Greenland. Jahn and others (2012) compared simulations from ten global ocean–sea-ice models of Arctic fresh water, and concluded that improved simulations of salinity variability are required to advance understanding of liquid fresh water export.

Improving the representation of the salinity distribution is crucial. However, inclusion, representation and parametrization of a number of processes is required (Steele and others, 2001; Komuro, 2014). For example, in many global ocean–sea-ice models, salt rejection in the first level of the ocean model during ice formation leads to an unrealistic salt distribution in the mixed layer and below (Nguyen and others, 2009).

The transfer of fresh water and sea ice from the Arctic Ocean to the North Atlantic are significant components of global ocean circulation (Haak and others, 2003; Gelderloos and others, 2012). Thus, the investigation of the variability of the surface layer (the upper 50 m) can make a significant contribution to understanding ocean-climate feedback. In particular, abrupt changes in surface-layer salinity may lead to critical transitions in the patterns of global ocean circulation, such as convection shut downs and climate disruptions (Hall and Stouffer, 2001; Gelderloos and others, 2012). A robust conceptual statistical model may help to describe the features of anomalies in salinification of the Arctic Ocean, which are key players in the formation of surface-layer salinity patterns. In this case, investigation of the structure of patterns and quality of anomalies leads to a better understanding of possible critical transitions in the patterns of global ocean circulation. Variations of Arctic Ocean surface-layer salinity have complex spatial and temporal structures, which are affected by many external factors. Our aim is to distinguish the most significant factors that led to recent changes in surface-layer salinity patterns.

Here, we present a statistical model for Arctic Ocean salinity fields based on multiple linear regression analysis, which builds on ideas presented in prior studies (e.g., Timokhov and others, 2012). This statistical model of variability of Arctic Ocean winter salinity in the 5–50 m layer is used as a method of reconstruction of observed winter salinity fields presented in Pokrovsky and Timokhov (2002). The model is based on an empirical orthogonal function (EOF) decomposition of the salinity data (e.g., Hannachi and others, 2007), and a multiple correlation analysis of the time series associated with the first three leading modes, or principal components (PCs); see Appendix Box 1 for a schematic diagram of the model. The contribution of atmospheric factors and hydrological processes in the spatial distribution of surface-layer salinity was interpreted by determining the structure of the multiple correlation equations. The variability patterns and relationships identified through the statistical analysis and modeling inform a conceptual model for principal drivers of Arctic salinity.

METHODS

Data set

This study is based on a collection of more than 9800 instantaneous temperature and salinity profiles, with data available at the standard levels (5, 10, 25, 50, 75, 100, 150, 200, 250, 300, 400, 500, 750, 1000 and so on every 500 m), collected between 1950 and 1993. The data were obtained from the Russian Arctic and Antarctic Research Institute (AARI) database, which was also used in the creation of the Joint US–Russian Atlas of the Arctic Ocean for winter (Timokhov and Tanis, 1997). The AARI database was first introduced by Lebedev and others (2008). This is complemented by the data made available over the period 2007–2012 from the expeditions of the International Polar Year (IPY) and afterward, which consist of Conductivity Temperature Depth (CTD) and eXpendable Conductivity Temperature Depth (XCTD) data, as well as data from the Ice-Tethered Profiler (ITP)-buoys (more than 14 600 stations in total). The average vertical resolution of all these profiles was 1 m. In areas where observations were missing, temperature and salinity data were reconstructed in a regular grid for the period 1950–1989 as detailed in the next subsection. The time period of 1994–2006 was not included in the analysis as there are too few data to construct reliable gridding fields. Thus, the working database is represented by grids with 200 km horizontal spacing, covering the deep part of the Arctic Ocean (with depths of more than 50 m). According to Treshnikov (1959), Rudels and others (1996, 2004) and Korhonen and others (2013), the average thickness of the Arctic Ocean mixed layer for the winter season is ~50 m. A description of the data sources for other physical parameters, used as predictors for the statistical model, can be found in Table 1.

Table 1. Predictors used for the approximation of PCs

The database used in this study belongs to the Oceanography Department of the Arctic and Antarctic Research Institute and it is not freely available. To mitigate the related issues, we provide additional data description. Table A1 in the Appendix contains a list of the expeditions and number of stations that were used for reconstruction and gridding of salinity fields. Figure A1 shows the overall observation density and the year associated with each observation. The data exhibit a spatiotemporal non-uniformity that is undesirable but expected given the logistical challenges associated with recovering long-term observations of Arctic salinity. While this dataset has gaps and is proprietary, we believe it provides one of the most spatiotemporally detailed observational resources for studying multidecadal variations in basin-wide Arctic salinity. This manuscript is motivated by the potential to advance understanding through identification, dissemination and modeling of the principal modes of variability in these long-term salinity observations.

Gridded fields of surface winter salinity were compared with fields from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS; Zhang and Rothrock, 2003; Lindsay and Zhang, 2006) for their overlapping period 1978–1993 (dashed curves, Fig. 1). PIOMAS is a coupled ice-ocean model which uses data assimilation methods for ice concentration and ice velocity. Forced by atmospheric observations, its output is available for 1978 to near present and is widely used as a reference for Arctic variables with limited long-term observations including salinity and sea-ice thickness. Maps of long-term means for both datasets are similar (Figs 2a and A2a), with a correlation coefficient of 0.88. Nevertheless, PIOMAS data generally provide higher salinity values for the Amerasian Basin (Canada Basin together with Makarov Basin) for the overlapping period (Fig. 1). The associated variance maps are also significantly correlated with each other (correlation coefficient R = 0.36; statistical significance level p = 0.05) but exhibit some salient differences (Fig. A3). In particular, a high variance zone along the Lomonosov Ridge is prominent in the AARI dataset, but is absent in PIOMAS data. PIOMAS instead features several centers of high variance along the shelf.

Fig. 2. The average salinity field (a) and first three modes of the average salinity field decomposition for the layer 5–50 m: (b), (c), (d) – first, second and third modes, respectively, for the periods 1950–1993 and 2007–2012.

To test for artifacts from the data gaps and the interpolation procedure (reviewed in the next subsection), we make several comparisons across methods and to independent datasets in the subsections to follow. For example, we also performed the EOF analysis with and without the additional 2007–2012 period, and report only modest change to the resulting modes of variability (section ‘Decomposition of surface-layer salinity fields by EOF’). In section ‘Decomposition of surface-layer salinity fields by EOF’, we also compare EOFs from the AARI data to those from PIOMAS.

Field reconstruction and interpolation

To provide temporal and spatial continuity, we have unified existing datasets using reconstruction and gridding. The technique of computing gridded fields for the period from 1950 to 1993 was described by Lebedev and others (2008), and is summarized here.

These techniques are based on the method of ocean field reconstruction, proposed by Pokrovsky and Timokhov (2002). This method, which was used to obtain gridded salinity fields, is given by

(1)$$\eqalign{z_i = & z_i^{\left( r \right)} + e_i,\; \; {\langle} z_iz_j{\rangle} = \sigma _{x_ix_j},\; \; {\langle}z_ie_i{\rangle} = 0, \cr {\langle}e_i{\rangle} = & 0,\; \; {\langle}e_ie_j{\rangle} = \delta _{ij}\sigma _{\rm e}^2 = \sigma _{{\rm e}_{ij}}^2.} $$

Here z(t, x) is the measured value of an oceanographic variable (e.g., temperature or salinity), and is a random function of time t and spatial coordinates x; i and j are the nodes of the irregular data grid; the notation <…> denotes the ensemble average of a value. We can write the observed value of z(t, x) as the sum of a true value z ^(r)(t, x) of the oceanographic variable and an observational error e(t, x). In addition, we introduce $\sigma _{x_ix_j}$ as a std dev. of spatial coordinates and $\sigma _{{\rm e}_{ij}}^2 $ as a std dev. of errors. We also propose that z _i^(r) has spatial correlations to the oceanographic parameters; that a systematic error is not identified; a std dev. of error (σ _e²) does exist; and $\delta _{ij} = \left\{ {\matrix{ {0,\; \; {\rm if}\; i \ne j} \cr {1,\; \; {\rm if}\; i = j} \cr}} \right.$ is the Kronecker δ.

Biorthogonal decomposition of the oceanographic variable can help to identify the connections between spatial and temporal distributions within the data:

(2)$$z\left( {t_j,\; x_i} \right) = \mathop \sum \limits_k c_k^j f_k\left( {x_i} \right) + e(t_j,x_i),$$

where f _k(x _i) is the spatial EOF, and c _k^j is the calculated coefficient or so-called k ^th PC.

As the next step, we approximate the EOF through linear combinations (with coefficients b _kl) of convenient analytical functions P _l(x _i) (e.g., polynomials, splines, etc.):

(3)$$f_k\left( {x_i} \right) = \mathop \sum \limits_l b_{kl}P_l\left( {x_i} \right).$$

Thus, the modified biorthogonal decomposition can be written as

(4)$$z\left( {t_j,\; x_i} \right) = \mathop \sum \limits_k d_l^j P_l\left( {x_i} \right) + e(t_j,\; x_i),$$

where $d_{l}^{j} = \mathop \sum {b}_{kl}\;c_k^j $.

The main goal of this spectral analysis method is to estimate the coefficients of spectral decomposition $\{ c_k^j \} $ and {b _kl}, in order to identify dominant modes of behavior. In this case, we rewrite formula (2) in the following matrix form:

(5)$$Z = F \, \cdot \, C + e,$$

where · denotes matrix multiplication.

The matrix F is formed by the values of the EOF, the matrix Z is composed of the totality of the measurement data at the points of the observation net x _i, and the matrix e is filled out by observational error values.

The system of linear equations (5) with respect to the unknown coefficients $c_k^j$ can be solved on the basis of the a priori statistical information (1) with the use of the standard estimation of the least squares method. A formula for the estimation of the matrix of the unknown coefficients C was obtained in (Pokrovsky and Timokhov, 2002), and is written

(6)$$\hat C = (F^{\rm T} \, \cdot \, K_{\rm e}^{ - 1} \, \cdot \, F + K_{\rm c}^{ - 1} )^{ - 1}F^{\rm T} \, \cdot \, K_{\rm e}^{ - 1} \, \cdot \, Z^{\prime},$$

where X ⁻¹ and X^T denote the inverse and transpose of a matrix X, respectively. The covariance matrix of errors of the expansion coefficients K _c is a diagonal matrix composed of eigenvalues of the covariance matrix K _z. Here, the matrices K _z and K _e are covariance matrices of the std dev. of spatial coordinates and the std dev. of errors, respectively.

In order to obtain an estimate of the unknown variables at the nodes of the regular grid $\tilde x_i$, it is necessary to interpolate the EOF into the corresponding nodes of the grid and obtain a new matrix $\tilde F$ of the EOF, a new matrix $\tilde C$ of decomposition coefficients, and new matrix $\tilde e$ of observational errors. Using the matrix $\tilde F$ obtained in this way and the estimates of the coefficients $\hat C$ from formula (6), from the matrix relationship

(7)$$\tilde Z = \tilde F \, \cdot \, \tilde C + \tilde e,$$

we obtain an estimate of the unknown parameters at the nodes of a regular grid. Simultaneously with the salinity fields in the nodes of a regular grid, we can also calculate the covariance matrices of the errors of estimations obtained from the following equations:

$$K_{\hat{\rm C}} = \left( {I + \left( {K_{\rm c} \, \cdot \, \left( {{\tilde F}^{\rm T} \, \cdot \, \tilde K_{\rm e}^{ - 1} \, \cdot \, \tilde F} \right)} \right)} \right)^{ - 1} \, \cdot \, K_{\rm c},$$

(8)$$K_{\hat{\rm Z}} = \tilde F \, \cdot \, K_{\rm z} \, \cdot \, \tilde F^{\rm T},$$

where I is the identity matrix and $\tilde K_{\rm e}^{ - 1} $ is the covariance matrix of the observation error expanded over the regular grid x _i.

This approach is a combination of singular value decomposition and statistical regularization. These coefficients (modes) can be linked to the real physical processes that influence salinity as presented in section ‘Statistical approach’. Preparation of the average salinity field data for 2007–2012 consisted of several stages as detailed in the Appendix. First, we checked the data for random errors. Then, we used linear interpolation and assimilated the real plane with the field data through the virtual plane of data. Next, we constructed an interpolation via a grid of nodes (separately for each plane). The gaps in the data for uncovered sites were filled with climatic values from the Joint US–Russian Atlas of the Arctic Ocean (Timokhov and Tanis, 1997). Using an older climatology to fill in recent (2007–2012) salinity fields is appropriate here as the areas where the most significant changes are observed, i.e., Canadian and Eurasian Basins, have sufficient data coverage to lessen concerns about potential inaccuracies (Fig. A1).

Statistical approach

Here we describe the approaches to data analysis which were used for physical interpretation of our statistical model. Polyakov and others (2010), Rabe and others (2011) and Morison and others (2012) have emphasized that the thermohaline structure of the surface layer has undergone significant change over the last decade. However, there is still room for discussion on what physical processes led to these changes or what the future trends may be.

The analysis of the variability of the surface layer (including salinity fields) of the Arctic Ocean may be based on a decomposition using EOFs. This approach is useful in our case because decomposition by EOF analysis gives modes (spatial patterns) and PC time series, which allow us to divide the variability into spatial and temporal components. Each mode describes a certain fraction of the total variance of the initial data, and the EOFs are conventionally ordered so that the first EOF explains the most variance and subsequent EOFs explain progressively less variance (Hannachi and others, 2007). The first three modes describe more than half the variance of the analyzed fields as further detailed below, which allows significant compressing of the information contained in the original data (Hannachi and others, 2007; Borzelli and Ligi, 1998). The EOF decomposition was carried out for the average salinity fields for the layer 5–50 m, yielding PCs for the periods of 1950–1993 and 2007–2012.

Multiple linear regression was used to model the PC time series to identify predictors that determined variability of the salinity fields. The regression equations can give projections of future changes because the predictors lead the salinity field by various temporal lags. The statistical model is characterized by a system of linear regression equations constructed for the first three PCs. The candidate predictors were as follows: atmospheric circulation indices (Arctic oscillation (AO) and Arctic dipole (AD) indices; see Table 1) calculated for winter (October–March to cover the period of active ice formation) and summer (July–September to cover the period of active ice melting), river runoff, the area of the ice-free surface in Arctic seas in September and water exchange with the Pacific and Atlantic Oceans. For the latter two water exchanges, we used the Pacific decadal oscillation (PDO) and Atlantic multidecadal oscillation (AMO) indices as respective proxies because of their influence on the temperature and salinity of water, which is entering through the Bering Strait and the Faeroe-Shetland Strait to the Arctic Ocean (Zhang and others, 2010; Dima and Lohmann, 2007). Atmospheric indices were averaged by different time periods within indicated winter and summer seasons. The optimal periods of averaging for a particular index were chosen on the basis of maximal correlations with PCs.

RESULTS

Decomposition of surface-layer salinity fields by EOF

As a result of EOF decomposition of the salinity fields for the 5–50 m layer, we obtained two sets of modes and PCs – one for the period of 1950–1993 (series 1) and one for the same period adding the years 2007–2012 (series 2). North's rule of thumb states the following: if the sampling error in an eigenvalue is comparable to the distance to a neighboring eigenvalue, then the sampling error of the EOF will be comparable to the size of the neighboring EOF (North and others, 1982). Based on this rule, the first three modes were accepted for further analysis as physically significant. The first three modes obtained by the decomposition of series 1 describe more than 55% of the total dispersion of the initial fields. The first three modes of series 2 describe more than 61% of the total dispersion. Nevertheless, the first modes for both decompositions have very similar shapes. The only differences are a more distinct dipole structure between the Canada Basin and Eurasian Basin, and positive EOF loading (instead of negative in the first mode of series 2) over much of the Nordic seas that appears in the first mode of series 2. As the Nordic seas region is the pathway of AW in the Arctic Basin (Karcher and others, 2007), we assume that the change of sign of EOF values is associated with increased temperature and salinity of AW inflow and subsequent salinification of the Eurasian Basin (Polyakov and others, 2010; Beszczynska-Möller and others, 2012). Thus, the modes obtained by decomposition in series 1 cannot take into account the essential features of the distribution of surface-layer salinity fields associated with the salinification of the Eurasian Basin. Therefore, for further analysis we used the PCs and modes obtained upon decomposition in series 2 (Fig. 2).

As a point of comparison for these EOFs, Appendix Fig. A2 presents a similar EOF analysis using model output from the PIOMAS. The spatial pattern of PIOMAS mean salinity (Fig. A2a) reasonably resembles the pattern shown for our data in Figure 2a. EOFs of PIOMAS salinity after detrending (Figs A2b–d) repeat the main features of corresponding results in Figure 2b–d, despite incomplete overlap over the analysis periods. In particular, for both datasets, the leading EOF for salinity features a prominent dipole between the Canada Basin and Eurasian Basin (Figs 2b and A2b). However, the leading EOF of PIOMAS salinity has a more patchy structure. In particular, there are negative centers of action situated along the Siberian shelf break, along with freshened areas (Fig. A2a) that are not as clearly pronounced in the AARI salinity data (Fig. 2a). The second EOF of salinity features a negative center of action elongated along the Lomonosov Ridge surrounded by positive loading strongest along Siberian shelf (Figs 2c and A2c). Less agreement is seen for the third EOF (Figs 2d and A2d), which is perhaps unsurprising in moving toward modes accounting for less variance.

The linear regression equation for the PCs

We present here a statistical model of interannual variability of Arctic Ocean surface-layer salinity. This research builds on already established approaches used by Pokrovsky and Timokhov (2002), specifically their reconstruction of salinity fields applying modified EOF methods.

We suggest some additions to improve on the ideas presented in previous research by Timokhov and others (2012). For example, the analysis presented is based on a dataset updated for the period 2007–2012, which is important for understanding the physical processes during the dramatic recent changes in Arctic sea ice. The area under consideration was extended and includes the Nordic seas and part of the Siberian shelf with depths of more than 50 m. This expanded domain allowed consideration of the areas strongly influenced by river runoff and ice processes, and makes visible their effects on the structure of the EOFs. Also, in contrast to our previously published research (Timokhov and others, 2012), we do not use the previous values of the PCs (history) as predictors for linear regression, which simplifies the physical interpretation of the equations obtained.

A set of external factors having the most correlation with salinity values (the predictors as described in section ‘Statistical approach’) have been defined based on the results of correlation analysis. As a result of linear regression analysis, we obtained empirical equations for the first three PCs (see Table A2 in Appendix). The structure of these equations can be explained through the sets of factors that simulate the effects of both atmospheric and hydrological processes.

Thus, the predictors used can be divided into two groups. The first group includes atmospheric circulation indices and reflects the influence of atmospheric processes. The second group corresponds to hydrological processes: river runoff into Arctic seas, inflows through the Bering Strait and the Faeroe-Shetland Strait, which were characterized by the PDO and AMO indices, and the areas of open water in the Arctic seas in September. Predictors were included in the equations with different time shifts (lags). The value of the time shift was 1–10 years and was chosen to maximize correlations of predictors with the PCs as noted above.

The contribution of each group to the explained variance in PC₁ through PC₃ can be calculated based on the magnitude and sign of the regression coefficients of corresponding predictors included in that particular group. In this case, hydrological processes have a dominant contribution to the explained variance of all PCs (from 86 to 61%). Atmospheric factors (i.e., AO and AD) contribute from 14 to 39%.

DISCUSSION AND SUMMARY

The first mode of the surface-layer salinity decomposition (EOF₁) displays an out-of-phase relationship between salinity anomalies in the Canada and Eurasian Basins (which includes the Nansen and Amundsen Basins) and Makarov Basin (Fig. 2b).

In the late 1980s, as a consequence of surface air temperature rising, the atmospheric circulation regime in the Arctic began to change (Steele and Boyd, 1998; Proshutinsky and others, 2009; Morison and others, 2012). Degradation of the Arctic anticyclone, shifting of the pressure pattern in the Beaufort Sea counterclockwise from the 1979–1992 pattern (Morison and others, 2000) and strengthening of the dipole pressure pattern (Overland and others, 2010) were observed. According to Wang and others (2009), large values of the AD indices (higher than 0.6 std dev.) could be a primary reason for the historical record lows of sea-ice extent in the summers of 1995, 1999, 2002, 2005 and 2007. In addition, in the late 1980s, the inflow of warm and highly saline AW into the Arctic Basin increased (Morison and others, 2000; Polyakov and others, 2010). Observed shoaling of the AW upper boundary, together with a decrease of static stability in the halocline layer, led to an increase in upper layer temperature and salinity in the Eurasian Basin (Polyakov and others, 2010). At the beginning of this century, the heat flux through the Bering Strait to the Chukchi Sea increased (Woodgate and others, 2010). Comparatively warm and fresh (salinity range 31< S < 32) summer Pacific waters, due to their low density, were able to inject heat close to the ocean surface (Stigebrandt, 1984) and enhance ice melting in the Canada Basin (Shimada and others, 2006) which led to decreasing surface-layer salinity in this region.

These observations allow us to suggest that salinity differences between the Canada Basin and Eurasian Basin, which became more pronounced in recent years (Fig. 1), are the consequence of these processes. Our suggestion is supported by the regression equation for PC₁ (Table A2) from which we see that PC₁ is a function of AMO, PDO, open water area in the East Siberian and Chukchi seas (that can be considered as an indirect indicator of fresh water inflow from these seas to the Arctic) and summer AO index. The time lags must be related with the time that it takes for Atlantic and Pacific waters to reach the Arctic Basins and become involved in associated circulation.

According to Karcher and others (2002), travel time for the propagation of anomalies in AW is 5–10 years from the Nordic Seas to the Eurasian Basin. This is reflected by AMO time lags. Bourgain and Gascard (2012) revealed that the warm signal from the Bering Strait propagated in the interior of the Canada Basin over 4–5 years. A proxy of this process is the PDO. The travel times for the Siberian river water from the river mouths to the shelf edge are estimated to be 2–5 years (Schlosser and others, 1994; Karcher and others, 2002). These time periods are in good agreement with time lags of the statistical model predictors, which may vary from 3 to 10 years (Table A2).

EOF₂ exhibits opposite polarity of salinity anomalies in the central Arctic Ocean and near-slope areas (Fig. 2c). Spectral analysis of the associated PC₂ revealed a 9-year cycle (periodogram is not shown here), which we associate with shifts between cyclonic and anticyclonic circulation regimes (Proshutinsky and Johnson, 1997; Rigor and others, 2002). The regression equation for PC₂ demonstrates its dependence on the summer AO and AD indices. Thus, during an anticyclonic regime of atmospheric circulation, fresh surface waters tend to flow to the center of the Arctic basin and negative salinity anomalies form there. At the same time, along the slopes there is upwelling of AW and positive salinity anomalies occur (Proshutinsky and Johnson, 1997). During a cyclonic circulation regime, reversed momentum forcing should likewise produce positive salinity anomalies in the central Arctic and negative salinity anomalies along the slopes.

The contribution of each predictor in the variability of a particular PC was evaluated as:

(9)$$I = \displaystyle{{\sigma _{i} \cdot \alpha} \over {\mathop \sum \nolimits (\sigma _i \cdot \alpha )}} \cdot 100\%, $$

where σ _i is the std dev. of the predictor and α is the regression coefficient of the predictor. According to this formula, contributions of the predictors for PC₂ (Table A2) were calculated. PDO and river runoff from the Laptev, East Siberian and Chukchi Seas make the largest contribution to the variability of PC₂ (33.8 and 26.9%, respectively) with slightly weaker effects from AO and AD (20.5 and 18.9%).

EOF₃ is represented by a field with multicore structure (Fig. 2d). The positive centers of action spread from the Beaufort Sea over the North Pole to the Kara Sea and are surrounded by negative centers of action. This kind of distribution is associated with an AD (Wu and others, 2006). The winter AD index accounts for ~22% of the variability of PC₃. The most distinct negative cores are located along the shelf of the Laptev and Chukchi Seas and also near Greenland. In our statistical model, associated variations are accounted for by river runoff from the Laptev, East Siberian and Chukchi Seas, the PDO index, and the AMO index, which account for 20, 25.8 and 32.4% of the variability of PC₃, respectively (Table A2).

All predictors included in the regression equations (with particular time lags and averaging periods) are statistically independent, i.e., they are not significantly correlated with each other, except for AMO (−10) and PDO (−3). These indices have a slight positive correlation (R = 0.33), but this is not a concern because they have different regions of influence and are associated with different proxies.

Time series of PC₁₋₃ show a mixture of interannual and quasidecadal oscillations (Fig. 3). Based on the configurations of the EOFs (Fig. 2), the regression equations (Table A2) and results of the spectral analysis of the PCs, we may assume that large-scale surface-layer salinity anomalies (with periods longer than 20 years) are the result of water exchange effects. The shorter period (8–9 years) variations appear to be determined by atmospheric circulation processes. Also, interannual variations occur due to interannual variability of both atmospheric and hydrological processes.

Fig. 3. The actual (black line) principal components and calculated principal components (red dashed line) with the help of the equations of linear regression. Correlation coefficients between the calculated time series of PCs and actual PCs are: r(PC₁) = 0.88; r(PC₂) = 0.73; r(PC₃) = 0.55.

The derived equations in the Appendix (Table A2) describe the first three PCs for the period 1950–2012. Calculated with these equations, the modeled PCs agree well with the values of the PCs directly derived from the decomposition of salinity fields via EOF analysis (Fig. 3).

Theoretically, the salinity fields for 1994–2006 can be reconstructed using this statistical model. We noted above that this period had gaps in observational data. The salinity fields of 1995, 2000 and 2005 were chosen for reconstruction to demonstrate the capabilities of the statistical model. The results were compared with PIOMAS model data. Although reconstructed fields have high significant correlations with the PIOMAS fields (correlation coefficients are 0.84, 0.88 and 0.81), in the Amerasian Basin they show lower salinity values than the PIOMAS data. Differences in this region may reach 2 psu. In the Eurasian Basin, specially over the Lomonosov Ridge, the reconstructed salinity values are higher than the PIOMAS data with differences of up to 1.5 psu (for 2005) (Fig. 4).

Fig. 4. Maps of differences of salinity fields reconstructed with the statistical model and those from PIOMAS data.

We also applied our statistical model to the reconstruction of salinity fields for 2013–2014, extending beyond the data record in order to develop a retrospective forecast (sometimes referred to as a hindcast). As a result, we obtained salinity fields that correspond to the trends observed in recent years. This preserved the freshening in the Canada Basin as well as salinification of the surface layer over the Lomonosov Ridge (compared with average surface-layer salinity for 1961–1990) (Fig. 5). According to our modeled values for 2013–2014, freshened water from the Beaufort Gyre should have moved westward along the Canadian Continental Slope in 2013. There are also negative salinity anomalies observed in the Arctic seas along the Siberian shelf. These processes were able to freshen the Eurasian Basin slightly so that in 2014 positive salinity anomalies over the Lomonosov Ridge were lower than in 2013. To demonstrate the quality of the forecast, we compared the salinity field for 2013 with the corresponding gridded field of observational data and PIOMAS data (Fig. 5). As we were not able to find enough data in winter 2014 to produce a reliable gridded field, comparison for this year was not conducted. In both cases, the reconstructed values are lower in the Canadian Basin and higher in the Eurasian Basin. However, our results are closer to the observational data than to the PIOMAS data, as differences in the first case are not larger than 0.7 psu (Fig. 5d), and in the second case, they are as large as 2.5 psu (Fig. 5f). In any case, PIOMAS data always show higher salinity values in the Amerasian Basin and lower values in the Eurasian Basin than prognostic or reconstructed fields. Some discrepancies from PIOMAS are expected for at least two reasons. First, PIOMAS salinity fields are simulated values not directly guided by salinity data assimilation. Second, the statistical model we use here is based on a subset of the EOFs and thus does not represent all the observed salinity variance as further discussed below.

Fig. 5. Reconstructed salinity fields for the layer 5–50 m in 2013 (a) and 2014 (b); actual salinity field for the layer 5–50 m in 2013 (c) (from AARI data), difference between actual salinity field and reconstructed one for 2013 (d); PIOMAS salinity field for 2013 (e) and difference between PIOMAS salinity field and reconstructed one for 2013 (f).

This method of salinity reconstruction may suffer from inaccuracies due to the higher frequency variability of the calculated PCs. The model may not reliably generate PCs for short-term time series and does not reflect properly the effect of short-term processes such as seasonal variations of ice melting/formation or river runoff, although the trend in variability of all three PCs is reproduced correctly. Therefore, the model can be used for tracking long-term processes of the structural transformation of salinity fields.

Validation of the model was carried out by calculating an error of reconstruction for the surface-layer salinity fields. The difference between the actual and reconstructed salinity fields (ε) was determined as a percentage by the following formula:

(10)$$\varepsilon = \left( {\sigma \left( {S_{\rm f} - S_{\rm c}} \right)/\sigma \left( {S_{\rm f}} \right)} \right) \cdot 100\%, $$

where σ is the std dev., S _f is the actual salinity and S _c is the calculated salinity.

Twelve surface-layer salinity fields from the time series under consideration (fields for the years 1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 2007, 2009 and 2011) were reconstructed, using modeled values of the PCs. The years were chosen at approximately equal intervals in order to reflect the different stages of the salinity field evolution through all decades. The average error of reconstruction for the chosen fields was 18.4%. As the first three EOF modes describe 61% of the variability of the initial fields, the error obtained is less than the variance not covered by the first three EOFs. Thus, we have a system of regression equations (statistical model) that may skillfully reproduce long-term salinity anomalies. The rest of the surface-layer salinity variance captured by higher order EOFs (~39%) is likely explained by short-term and probably local processes such as ice formation and cascading in polynya regions (Ivanov and Watanabe, 2013), deep convection or mixing with the AW upper boundary (Ivanov and others, 2012).

Thus, we have identified various patterns in Arctic Ocean surface-layer salinity fields using a reliable statistical model. In addition, we have found anomalies in the salinity fields which have occurred in the past, and conclude that more than 60% of surface-layer salinity variability is related to long-term processes ranging from water exchange with neighboring oceans (with periods longer than 20 years) to interannual variations of atmospheric circulation, ice pack characteristics (such as ice extent and volume export) and river runoff. The other nearly 40% is due to short-term and local processes, such as changes in polynya activity or local mixing processes. Our findings again raise questions about non-linearities in global ocean circulation, particularly in the Arctic Ocean, which is strongly connected with Earth's climate system. In the future, information obtained about these anomalies may be helpful in determining whether Arctic Ocean salinity, and related oceanographic phenomena, have reached a critical threshold.