Factors Contributing to the Long-Term Sea Level Trends in the Iberian Peninsula and the Balearic and Canary Islands

1. Introduction

Sea level rise is one of the most important threats caused by climate change on coastal areas. Despite the uneven distribution of tide-gauges with long time series, global mean sea level has been reconstructed during the twentieth and twenty first centuries by means of Reduced Space Optimal Interpolation [1,2,3,4,5] and Optimal Statistical Interpolation [6]. These works and the recent use of altimetry data have revealed that global mean sea level increased at a rate of 1.7 mm/yr [1.3 to 2.2 mm/yr] from 1901 to 2018, and this rate has accelerated to 3.7 mm/yr [3.2 to 4.2 mm/yr] for the period 2006-2018 [7,8]. The main causes of sea level rise were the ocean warming and the melting of glaciers, and of Greenland and Antarctica Ice Sheets, [9,10,11]. However, changes in the global mean sea level could not be representative of local or regional ones [8,12], and the assessment of such changes is of paramount importance for establishing adaptation strategies [13].

Changes in local mean sea-level can be decomposed into a mass component (manometric sea-level change) and a density component (steric sea-level change) [14,15]. The latter can be further decomposed into thermosteric and halosteric sea-level changes. Although the halosteric contribution is considered to be negligible on a global scale [14,16], it could have a great importance on a local or regional one [17], and in concentration basins such as the Mediterranean Sea [15,18,19]. The mass component of sea level can be inferred subtracting the steric contribution to the observed sea level, or recently, since the launch of the GRACE mission in 2002, from gravimetry measurements. Nevertheless, the gravimetry data have a coarse resolution of 300 km [20] and therefore cannot be used for local or even regional studies. For instance, in the case of the Mediterranean Sea, [19] found that these data could be used to infer the mass component for the whole Mediterranean Sea, but were not suitable for smaller spatial scales.

The main goal of this work is to analyse the long-term trends of relative sea level change (RSLC) on a local scale in the north, west and south Atlantic coasts of the Iberian Peninsula, the Canary Islands, and the Spanish Mediterranean coasts, including the Balearic Islands. Long time series from tide-gauges are used for the analysis of the period 1948-2019, and both, tide-gauge data, and altimetry data (in this case relative to the reference ellipsoid; geocentric sea-level change, [14]) are compared for the more recent period 1993-2019. Previous works have corrected the observed sea level for the atmospheric forcing by means of barotropic 2D models forced by realistic pressure and wind fields [21,22,23,24]. The thermosteric and halosteric contributions to RSLC have been calculated from vertical profiles of temperature and salinity compiled in different data bases or from reanalysis projects [15,25,26]. In this work we follow a different approach, already used in [27]. Monthly time series of atmospheric pressure, zonal and meridional components of the wind, and thermosteric and halosteric contributions were calculated from different databases. Then a multiple linear regression was used to determine which factors contributed to the sea level variability and to estimate their contributions to the observed sea level trends. Section 2 describes the data used and the methodology applied. Section 3 presents the results, and finally a discussion and summary of the main results are presented in section 4.

2. Materials and Methods

2.1. Sea Level Data

Monthly time series of Revised Local Reference (RLR) sea level were downloaded from the Permanent Service for Mean Sea Level (PSMSL; https://www.psmsl.org/data/obtaining/ [28]). We downloaded the 84 time series corresponding to the tide gauges distributed along the coasts of the Iberian Peninsula, the Balearic Islands, and the Canary Islands. These time series were updated until 2018 or 2019 in most of the cases. The main objective of this work was to study long-term changes in sea level and to determine those factors that contribute to these changes. Therefore, only time series starting in 1990 or earlier were selected for their analysis. The only exceptions to this rule were the Las Palmas, and Palma time series which started in 1991 and 1997 respectively. These shorter series were included in order to have some information about the sea level variability in the Balearic and Canary Islands where long time series are very scarce. Black dots in Figure 1 show the position of the 17 tide gauges selected for their analysis.

The PSMSL carries out a quality control of data and labels those data that are suspicious. Such data were discarded in the present work. Isolated spikes with changes larger than 250 mm were also rejected. The selected time series presented some gaps that were filled by interpolation with natural cubic splines, when they were shorter than three months [1]. Larger gaps were filled by means of linear regression on nearby tide gauges ([27] and supplementary material). Open circles in Figure 1 show the location of tide gauges whose time series were not long enough to be considered for their analysis, but could be used for filling the gaps in the long ones. In some cases, the open circles coincide with the black dots. This is due to the existence of two tide gauges in the same location operated by different institutions. In these cases, the longest time series was taken, while the shortest one was used for completing or even extending the long one.

There are four exceptions to the data processing method described above. In the case of the Cádiz time series, we used the recovered record by [29], after a very thorough work of data archaeology in which some problems of levelling were detected and corrected. Alicante and Santander time series were those reconstructed by [30]. Finally, in the case of Tenerife, we used the time series reconstructed by [31]. The latter extends until 2012 and was completed until 2019 using the data from the PSMSL. In the case of Alicante, there are two different tide gauges with available long time series. One is located in the inner harbour, and the other one is in the outer one. Finally, we analysed the series corresponding to the outer harbour because it had a larger extension, and used that from the inner harbour for filling some gaps. A detailed description of the regression analysis used to fill the gaps and statistical information is presented in supplementary material (including Figures S1, S2, and Table S1). These time series have different starting dates, and extend to 2018 or 2019.

Sea level data from tide gauges are affected by the effect of Glacial Isostatic Adjustment [32,33]. The GIA contribution to the RSL at each location was obtained from the files available at the PSMSL [33,34,35]. They are presented in the first column of Table 1.

2.2. Altimetry Data

Relative Sea Level data from tide gauges are affected by Glacial Isostasy Adjustment (GIA; [32,33]) and should be corrected for this effect. Besides GIA, other land movements can also impact RSLC (for instance, tectonic subsidence/uplift or groundwater extraction etc. [9]) and should be corrected by means of Global Navigational Positioning System (GNPS; [36,37]). However, no GPS measurements were available for most of the tide gauges used in this work. Altimetry data are also influenced by GIA, but not by other land movements. Therefore, comparison of sea level trends estimated from tide gauges and altimetry data could provide some information on the existence of such movements.

Altimetry gridded sea level anomalies were downloaded from the Copernicus Marine Service (https://resources.marine.copernicus.eu/products). We used the product Global Ocean gridded L4 Sea Surface Heights and derived variables (SEALEVEL_GLO_PHY_L4_MY_008_047). This dataset includes monthly sea level data referred to a reference ellipsoid and has a spatial resolution of 0.25º x 0.25º in latitude and longitude. We selected the grid points closest to each of the 17 tide gauges considered in this work and obtained the corresponding monthly time series. These time series extend from 1993 to 2021. Unlike tide gauge observations, altimetry data were corrected for the effect of the atmospheric forcing by means of a barotropic 2D model forced by atmospheric pressure and wind fields [24].

2.3. Temperature and Salinity Data

Various sources of temperature and salinity data have been used to compute the density component of sea level. First, gridded temperature and salinity data were obtained from the Met Office Hadley Centre observations datasets (version EN.4.2.1; [38]). This data set offers monthly data on a 1º x 1º grid with 42 vertical levels. Monthly temperature and salinity profiles were obtained for the closest grid points to each of the 17 tide gauges selected for their analysis. Monthly time series of steric (Δη_ST), thermosteric (Δη_T), and halosteric (Δη_H) sea-level change, and mass of salt contributions to RSLC were estimated for each of these 17 grid points following [15,19,27]:]:

(1.1)

(2.2)

(3.3)

(4.4)

Where ρ_S is sea surface density, η is the sea surface elevation above an arbitrary reference level, α and β are the thermal expansion and haline contraction coefficients, and ρ, S and T are the density, salinity and temperature along the water column. We finally obtained 17 monthly time series of steric, thermosteric, halosteric and mass of salt components of the sea level variability corresponding to the 17 tide gauges selected for their analysis. These time series extended from 1940 to 2019.

For comparison, a similar gridded data set has been used, available through the NCAR/UCAR Research Data Archive [39,40].

Both data bases described above are interpolated onto regular grids and thus share similar characteristics and limitations, especially in coastal waters. Therefore, also individual continuous profiles have been used when they were available and consistent over time. In particular, in the case of the Mediterranean region, the thermosteric and halosteric contributions were also calculated using the data from the seasonal monitoring program RADMED [41]. The RADMED monitoring program consists of four campaigns per year, including CTD profiles in all the oceanographic stations. Such stations extend along the Spanish Mediterranean shelf, slope and deep waters from Málaga to Barcelona, including the Balearic Islands. (see [41] for the details of this monitoring program).

2.4. Atmospheric Variables

Atmospheric forcing is one of the causes of sea level variability. This effect can be separated from others by 2D numerical models forced by realistic atmospheric pressure and wind fields [22,23,24,42]. In this work, we analysed this forcing statistically. For this purpose, we downloaded monthly time series of atmospheric pressure and the U and V components of the wind from the reanalysis of the National Centre for Environmental Prediction/ National Centre for Atmospheric Research (NCEP/NCAR, [43]). U and V stand for the zonal and meridional components of the wind respectively, being U positive when directed eastwards and V when directed northwards. This data set has a spatial resolution of 2.5º x 2.5º. As in previous cases, the grid points closest to the position of the tide gauges were selected. Finally 17 monthly time series were obtained for each variable (P, U, V) and each tide gauge. These time series extended from 1948 to 2019.

2.5. Data Processing

The long-term changes in RSLC were analyzed for the 17 locations marked as black dots in Figure 1. For each location the following monthly time series were compiled: RSL from tide gauges, altimetry sea level from altimeters, steric, thermosteric and halosteric sea-level change (η_ST, η_T, and η_H) calculated from temperature and salinity profiles, atmospheric pressure (P), and the zonal (U) and meridional (V) components of the wind from reanalysis data. All these time series were de-seasoned removing the mean seasonal cycle. To do so, for each variable, the data were grouped for each of the 12 months of the year, and mean values were calculated. These 12 values represent the climatological or average seasonal cycle for each variable. Then, the seasonal cycles were subtracted to the monthly times series and finally time series of anomalies or residuals were obtained with the same length than the original ones.

2.6. Statistical Model

The statistical approach followed in this work is the same used in [27]. To summarize it here, we considered that the de-seasoned sea level anomalies (η) depended on the variability of the meteorological forcing (P, U, V) and the thermosteric and halosteric sea level change (η_T, η_H). Hereafter, the sea level anomalies will be named as the dependent variable, and P, U, V, η_T, η_H will be the predictors or independent variables, to which the seasonal cycle has also been subtracted. Sea level anomalies also contain a part which is independent of the predictors (O). Hence, we can model the sea level variability as:

(2)

O represents all the other factors that are not accounted for in the other terms in (2). This includes the local mass addition component of sea level (manometric sea level change: freshwater and salt), land movements (including GIA), and errors or noise.

Predictors will vary on different time scales. We assumed that such variability could be decomposed into a linear trend (representing the average change over time scales longer than the length of the time series), and shorter time scales that include monthly, inter-annual, and decadal variability:

(3.1)

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

Where t is time, and z_P, z_U, z_V, z_T and z_H are the monthly time series of predictors de-seasoned and de-trended. The part of η not accounted for by the atmospheric forcing and the steric contribution can also be affected by a linear trend. In fact, the GIA contribution can be simply modeled as a linear trend for many applications. Other contributions such as the addition of fresh water and salt will be made of a linear trend and a term (z₀) expressing the monthly and inter-annual variability which we cannot calculate directly.

Substituting expressions (3.1) to (3.6) into (2), the sea level variability can be expressed as:

(4.1)

Being:

(4.2)

The first step was to estimate the linear trends for the five de-seasoned predictors to obtain de-seasoned and de-trended time series. Then, coefficients a, b, b₁, b₂, b₃, b₄, b₅ were estimated by stepwise forward linear regression of the de-seasoned sea level on time and on the five de-seasoned and de-trended predictors (for tide gauges and for altimetry). Notice that the stepwise regression does not include necessarily all the predictors and only those predictors that contribute significantly to the variance of the independent variable will be finally selected (a detailed description of this method can be seen in [44], for instance).

Once the coefficients of the selected predictors have been determined, the contribution of each predictor to the linear trend of η can be estimated as the product of that coefficient times the slope of the linear trend in equations (3) (b₁b_P, b₂b_U, b₃b_V, b₄b_T, b₅b_H). After correcting b for the GIA [33], b_O should represent the mass addition contribution (unless there are other vertical land movements). It can be estimated as b minus the sum of the previous contributions. Confidence intervals in the 95 % confidence level were calculated for each coefficient in expressions (3) and (4). Using these confidence intervals, the uncertainty for each contribution to the sea level trend was calculated using the formula for the error propagation.

It is worth noting that the contribution of the mass of salt to sea level was not included as a potential predictor. This will be further discussed below.

3. Results

3.1. Sea Level Linear Trends

Black lines in Figure 2 show the time series of de-seasoned sea level residuals from tide gauges and from altimetry measurements. This figure and its continuation correspond to the northern and western Atlantic coasts of the Iberian Peninsula and the Canary Islands. Figure 3 shows similar results for the Gulf of Cádiz (southern Atlantic coast of the Iberian Peninsula) and the Spanish Mediterranean coast. The left columns in Figure 2 and Figure 3 are the results from tide gauge data and the period 1948-2019. The central columns correspond to the tide gauge data for the period 1993-2019, and the right columns present the results for the altimetry data and the period 1993-2019. Each row corresponds to a different location. The linear trends estimated for each time series, and their confidence intervals in the 95% confidence level, are inserted in these figures. These linear trends are corrected for the effect of GIA and they are also presented in the second column of Table 1 and Table 2. The GIA contribution to the RSLC for each location is presented in the first column of table 1. Trends and confidence intervals in all the tables are presented with two significant figures when they are ≤ 25, and they are rounded to one single figure when they are >25. Table 1 shows those significant results obtained from tide gauge data corresponding to the period 1948-2019. Table 2 shows the results from tide gauge and altimetry data for the period 1993-2019.

Considering the period 1948-2019, with the only exception of Gibraltar, all tide gauges showed positive trends ranging from 0.59± 0.16 (Arrecife) to 2.29 ± 0.22 mm/yr (A Coruña). For the period 1993-2019, most tide gauges showed, as expected, positive trends that ranged between 1.3 ± 0.8 mm/yr (Leixoes) and 4.7 ± 0.7 mm/yr (Tarifa and Gibraltar). The linear trends estimated from altimetry data had a lower dispersion, ranging between 2.4 ± 0.4 mm/yr at Algeciras and Ceuta, and 4.1 ± 0.4 at Málaga (Note that altimetry trends at Palma are lower than altimetry trends in other locations, but they were estimated for the period 1997-2019). Initially, the sea level trend from the Vigo tide gauge presented a negative value. The comparison with the nearby tide gauges, and with altimetry results, makes this result suspicious. Apart from the time series used in this work, there is another tide gauge in Vigo, named as Vigo II in the PSMSL. Initially this series had not been used for the reconstruction of the Vigo time series as it was not necessary. Their analysis for the period 1998-2019 showed positive trends as in all the other cases. For this reason data for Vigo after 1998 were discarded and reconstructed by regression on the redundant tide gauge (Vigo II). Once corrected the suspicious period, the Vigo trend yielded a positive value for the periods 1948-2019 and 1993-2019 (see Table 1 and Table 2). Therefore, these will be the results discussed hereafter. The results from Gibraltar are also suspicious or simply erroneous and will not be considered for further discussion.

3.2. Linear trends of the atmospheric forcing, and the thermosteric and halosteric sea-level change.

Columns 3 to 7 in Table 1 and Table 2 show the linear trends for atmospheric pressure (#3), zonal or U(#4) and meridional or V(#5) components of the wind, and the thermosteric (#6) and halosteric (#7) sea-level change. These trends are presented only in those cases that are statistically significant in the 95 % confidence level. Table 1 corresponds to the period 1948-2019, and Table 2 is for the period 1993-2019.

For the period 1948-2019, atmospheric pressure increased with a trend of 0.02 ± 0.01 mbar/yr. The U and V components of the wind did not experience significant trends in most of the cases. The thermosteric contribution increased in all the locations. Nevertheless, the halosteric contribution was negative (with the only exceptions of Santander and Arrecife), significantly decreasing the final value of the steric contribution to the sea level long-term variability (see Table S2 in supplementary material for the trends of the steric component). Notice that the steric contribution is not strictly the sum of the thermosteric and halosteric ones. Nevertheless, the comparison of the steric component estimated by means of equation (1.1) and the sum of (1.2) and (1.3) showed that this is a valid approximation and the nonlinear terms in the equation of state of sea water were negligible for this purpose.

3.3. Linear Model

The de-seasoned sea level time series from tide gauge and altimetry data were regressed on time, the atmospheric variables, and the thermosteric and halosteric sea-level change, according to equation 4. The regression was performed for both periods 1948-2019, and 1993-2019. The model selected different variables in each case, depending on the location, the period of time, and the type of sea level data considered. Table 3 shows the coefficients for the selected predictors for the regression corresponding to the period 1948-2019. In this case the results were estimated only from tide gauge data. The multiple correlation coefficient (R) is also included. The square of this coefficient expresses the percentage of variance of the dependent variable explained by the linear model. Table 4 presents the coefficients for the linear model and the multiple correlation coefficients for the regression corresponding to the period 1993-2019. In this case two different sets of coefficients were obtained, one for the tide gauge data, and another one for the altimetry data.

The atmospheric pressure was always selected as a predictor for the sea level variability measured by means of tide gauges. These coefficients were always negative and close to -10 mm/mbar for the two periods of time analyzed. The U and V components of the wind also had a significant influence on the sea level variability in most of the tide gauges. The magnitude and sign of these coefficients depended on the location. The along coast component of the wind produces changes in the RSL according to the expected upwelling or downwelling processes. Therefore, positive values of the zonal component of the wind (westerlies) induce a sea level decrease in the southern coast of the Iberian Peninsula, which is evidenced by negative coefficients associated to the predictor U at Cádiz, Tarifa, Algeciras and Málaga (see values associated to tide gauges in Table 3 and Table 4). On the contrary, this coefficient is positive for the Santander tide gauge. In this case the positive U component induces downwelling and sea level rise, according to the orientation of the coast. Along the west coast of the Iberian Peninsula, it is the meridional component of the wind that produces the upwelling processes. In these cases the coefficients of the predictor V were positive, as positive values of the wind (southerlies) induce downwelling and sea level rise.

Despite the fact that altimetry data were corrected for the atmospheric forcing, there is a statistical contribution of the atmospheric pressure and the U and V components of the wind to the sea level variability inferred from altimetry data. Nevertheless, the coefficients were lower than in the case of the tide gauge data, especially in the case of the atmospheric pressure.

The thermosteric and halosteric contributions to sea level were selected by the linear model in some cases, not showing any clear pattern. It is noteworthy that these coefficients were always lower than 1 in those cases when these predictors entered the model (see discussion section). It could be argued that the thermosteric and halosteric contributions are negatively correlated. In the case that these two variables were highly correlated (co-linearity), the opposite effects of these two predictors on sea level could prevent these variables from being selected in the linear model. For this reason the halosteric component was regressed on the thermosteric one. The percentage of the variance explained ranged between 11 % and 47 %, with a mean value for the 17 tide gauges of 32 %. Furthermore, the Variance Inflation Factor (VIF) was calculated for each predictor to check for a possible co-linearity. This statistic was never over 2, and therefore the co-linearity of the predictors was acceptable. Furthermore, the linear model was also applied using the meteorological forcing (P, U, V) and the steric component of sea level (without decomposing it into its thermosptheric and halosteric components).The steric component was never chosen when both the thermosteric and halosteric components had not been selected in the linear model of equation (4) (see Table S3). These calculations were also repeated for all the geographical areas using the thermosteric and halosteric time series obtained from the NCAR/UCAR Research Data Archive. In the case of the Mediterranean Sea, the regression was also repeated using the temperature and salinity profiles from the monitoring program RADMED. The use of these different data sets did not improve the correlation of the thermosteric and halosteric contributions with the observed RSL.

The selected linear models explain a large fraction of the sea level variance for the two periods analyzed and both for tide gauge and altimetry data (see R values in Table 3 and Table 4). Considering the contribution of each predictor to the sea level variability (coefficients in the linear model, Table 3 and Table 4) and the linear trends of such predictors (Table 1 and Table 2) we calculated the contribution of each predictor and the mass addition contribution to the observed sea level trends (Table 5 and Table 6). The only clear contribution of the atmospheric forcing to the long-term changes of sea level was observed during the period 1948-2019 when a negative contribution (atmospheric pressure increase) was observed at all the tide gauges analyzed (Table 5). The U and V components of the wind did not contribute significantly to the sea level trends, with the only exceptions of Santander and Málaga. The thermosteric contribution for this period was positive in some tide gauges, but in those cases it was partially compensated by a negative halosteric contribution. Notice the negative thermosteric contribution at Leixoes (Table 5). During the period 1993-2019 (Table 6) no contribution of the atmospheric pressure was observed and there was a negative contribution of the V component of the wind at some locations. The positive contribution of the thermosteric component of sea level was observed at several locations, whereas the halosteric contribution was positive or negative depending on the location and the data set used.

3.4. Mass of Salt and Freshwater Contributions

The mass contribution was estimated subtracting the thermosteric, halosteric and atmospheric contributions to the observed sea level trends corrected for GIA (column #2 in Table 5 and Table 6). It can be considered that the sea level trend caused by the addition of freshwater can be estimated subtracting the mass of salt contribution, estimated by means of equation (1.4), to the mass contribution. The mass addition contribution experienced a clear increment during the period 1993-2019, compared to the period 1948-2019. Nevertheless, these figures should be taken very cautiously because of the difficulties for estimating the steric contribution. This will be discussed in detail in the next section.

4. Discussion and Conclusions

The atmospheric correction by means of the linear regression model yields the expected results for the pressure and the alongshore component of the wind, in those cases where the coast has a clear orientation. Assuming that the sea reaches the equilibrium state after a change of pressureΔP, the response of sea level should be:

(5)

The mean surface density in the area of study ranges between 1025.9 kg/m³ at the Canary Islands, and 1027.1 kg/m³ at the northernmost tide gauges of the Mediterranean Sea. According to these values, the relation between the sea level variations and those of the atmospheric pressure should be between -9.94 and -9.95 mm/mbar, that can be rounded to -10 mm/mbar and which is the well-known inverse barometer effect (see for instance [14]). Taking into account the uncertainty of the coefficients calculated in the linear regression, most of such coefficients are not different from the theoretical value (column #3 in Table 3 and Table 4).

It is more difficult to predict quantitatively the response of sea level to the wind variability. Nevertheless, our results show a qualitative agreement with the expected upwelling and downwelling processes along the northern and southern coasts of the Iberian Peninsula, as well as on its Atlantic coast. Westerly (positive) and northerly (negative) winds induce upwelling (decrease of sea level) on the southern and Atlantic coasts respectively. Therefore, the coefficients of the linear regression were negative in the first case and positive in the second one (see columns #4 and #5 in Table 4 and Table 4). In the northern coast of the Iberian Peninsula, westerly winds were responsible for downwelling (increase of sea level), yielding a positive value for the coefficient of the U component of the wind at Santander. The numerical values for these coefficients showed that the sea level change for m/s of variation of the wind is comparable to the effect of sea level change for each mbar of pressure.

The atmospheric contribution to the sea-level change is frequently estimated by means of 2D barotropic circulation models [21,22,23,24]. The final validation of these results is usually done calculating the variance reduction of the tide-gauge time series, or calculating the correlation coefficients between observed sea level and that predicted by the 2D barotropic models. [23] found a correlation ranging between 0.6 and 0.7 depending on the tide gauge considered, and [42] found correlations as high as 0.8. If the linear model was applied, considering only the atmospheric forcing, the multiple correlation coefficients ranged between 0.49 and 0.78, with the only exception of Arrecife, where it was 0.37. In any case, the statistical model used in this work should be only considered as a complementary approach to the numerical modeling, as the response of sea-level to the atmospheric forcing could be far from linear in some cases, depending on the topography and geometry of the coast. The explanation of the coefficients of the thermosteric and halosteric contributions to sea level presents more difficulties. The sea level variability can be decomposed into a mass and a steric contribution. The latter can be divided into thermosteric and halosteric ones. This decomposition can be expressed by means of equation (6) (see [15] for a detailed discussion of this expression and the interpretation of the different terms of it):

(6)

δm is the mass in a water column of area δA. The first term in the right hand side of equation (6) is the change in sea level produced by the change of mass per unit of area. The second and third terms are the thermosteric and halosteric contributions as defined in (1.2) and (1.3). Let us consider that we know the monthly change of sea level at any tide gauge or position where altimetry data are available. If we also know the change that the temperature and salinity have experienced along the water column for the same location and month, we could calculate their contributions to the observed sea level variability. In other words, according to equation (6), the coefficients that relate the thermosteric and halosteric terms and the observed sea level should be 1. If we wanted to calculate the mass component of the sea level variability we should simply subtract both contributions from the sea level change. However, the coefficients for these two predictors were always below 1 in the linear regression. Furthermore, in some cases the forward stepwise regression model did not select some of these predictors for explaining the variance of sea level. This simply indicates that, in those cases, the time series of thermosteric and halosteric contributions were not significantly correlated with the observed sea level. As explained in the results section, this cannot be attributed to the possible co-linearity of the thermosteric and halosteric sea-level changes. Nor it can be attributed to the data base used, as similar results were obtained using the NCAR/UCAR Research Data Archive and the RADMED project data. Therefore there are two possible explanations for this result.

First, the temperature and salinity data used for determining the steric contribution for each tide gauge or altimetry grid point, corresponded to a large area of 1º x 1º. This large area does not necessarily represent the local conditions of the tide gauges. If the steric component of sea level was spatially homogeneous, it would not be a problem to use this large geographical area of 1º x 1º in our calculations. On the contrary, we cannot be sure that the same kind of sea-atmosphere interaction occurs in coastal and open-sea waters, and the 1:1 relationship would not hold necessarily. We should be careful interpreting this result. It does not mean that the equation (6) is not right, or that equations (1.2) and (1.3) do not represent the steric contributions to the sea level variability. It would simply mean that the available time series of temperature and salinity profiles, represent the real conditions at the open sea, but they are not representative of the local changes occurring at the costal tide gauges.

This explanation has another problem. The altimetry sea level also corresponds to an open-sea area similar to that representative of the temperature and salinity profiles used for the calculation of the steric contribution. Neither in this case the η_T and η_H predictors were significantly correlated with the altimetry sea-level in all the locations, and in those cases when these predictors were selected by the linear model, the coefficients were also lower than 1, as in the case of tide-gauge data.

A second explanation is that the available temperature and salinity data do not yield reliable estimations of the steric contributions, neither for the open-sea areas, nor for the coastal ones. This problem has already been evidenced by [19] for the case of the Mediterranean Sea. These authors pointed out that the different available data bases did not allow obtaining consistent estimations of the steric component of the sea level and of the mass of salt contribution. As a consequence of this, consistent estimations of the contribution of the mass of freshwater could not be obtained. This problem arises from the scarcity of temperature and salinity data along the water column which makes very difficult to calculate the monthly, inter-annual and long term variability of the heat and salt content in the upper layer of the sea, where there is a very large natural variability [45,46]. [47] sub-sampled the results from a numerical model at the same times and locations where real temperature and salinity profiles were available. These authors interpolated the sub-sampled data onto a regular 3D grid on a monthly basis. Linear trends estimated from these interpolated data were not able to capture the real long-term variability of the simulated data, evidencing the limitations of the present gridded climatologies.

Taking into account all the problems related to the scarcity of data, we considered that the monthly thermosteric and halosteric contributions calculated in the present work could not be the real ones, but still could have a certain correlation with them. Once again the coefficients relating the observed sea level to the thermosteric and halosteric contributions would not be 1 and we could not subtract directly such contributions to calculate the addition of mass.

Either the first explanation is true, or the second one, the degree of correlation between the observed sea level and the available steric contributions was determined by the time series themselves, and the stepwise correlation analysis. In those cases in which there was a significant correlation, the coefficients that related the observed sea level to η_T and η_H and the linear trends calculated for these predictors were used to estimate the thermosteric and halosteric contributions to the linear trends of sea level. In such cases, the salinity profiles were also used to estimate the contribution of the mass of salt. Notice that as expression (1.3) is not directly used to estimate the contribution of the halosteric component, it would not be consistent to use equation (1.4) for the calculation of the mass of salt. If (1.4) is divided by (1.3), and neglecting the variability of β along the water column compared to the changes in ΔS, the ratio of both contributions would be

(7)

being ρ₀ and β₀ some average or reference values. Therefore, the mass of salt contribution was estimated as the halosteric contribution multiplied by the factor (7). It could also be argued that the term (1.4) could have been included as another predictor in the linear model and the mass of salt could have been estimated directly by the linear regression. In such case the two predictors given by equations (1.3) and (1.4) would be proportional and the resulting system of equations would be ill-conditioned.

In those cases in which none of η_T and η_H were significantly correlated with the sea level, those predictors were not used to estimate the mass contribution of sea level. We have to emphasize that this does not mean that there is no steric contribution to the observed sea level, nor contribution of the mass of salt. It simply means that the available time series of steric contribution do not resemble the real ones and therefore cannot be used for this purpose. Figure 4A (period 1948-2019) and 4C (period 1993-2019) show the mass contribution of sea level, once corrected for the atmospheric forcing. This contribution includes the addition of freshwater and salt. Figure 4B and 4D show the contribution of the mass of salt (red triangles) and the contribution of the mass of freshwater. In those figures corresponding to the 1993-2019 period, both the tide gauge (black lines) and the altimetry (blue lines) analyses have been included.

Considering the long period 1948-2019, only tide gauge data were available. Gibraltar sea level trend was negative (see column #2 in Table 1) and so was the contribution of mass addition and freshwater (Figure 4A and B). This is not a reliable result and indicates that some sort of leveling or any other problem affected to this tide gauge during this period. For this reason, this tide gauge will be excluded in the following discussion. The sea level, averaged for the whole area, (corrected for the effect of GIA) increased at a rate of 1.58 ± 0.19 mm/yr. This result is coincident with the trend observed for Global Mean Sea Level (GMSL) from 1901 to 2018, which is 1.7 mm/yr [7,8]. During this period, the atmospheric pressure had a positive trend in all the tide gauges analysed which induced a decrease of the sea level. However, the behavior of the thermosteric and halosteric contributions and that of the addition of mass of salt had different values for the different geographical areas analyzed. Averaging the results for the Atlantic coast of the Iberian Peninsula (northern and western coasts), the RSLC was 2.05 ± 0.21 mm/yr. It was observed a warming and salting of the water column which produced positive and negative trends for the thermosteric and halosteric sea-level changes. Nevertheless, there was no significant correlation between these terms and the observed RSLC. Hence, the thermosteric, halosteric and mass of salt terms did not contribute significantly to the RSLC and the addition of mass and that of freshwater were the same and larger than the observed RSLC (2.3 ± 0.4 mm/yr).

This is an unrealistic result. The water column has warmed during the period 1948-2019 (see column #6 in Table 1). Although the halosteric term had a negative trend that could partially counterbalance the thermosteric one, it is also associated to a positive contribution of the term of mass of salt. Therefore, the warming of the oceans, that have absorbed the 90 % of the heat stored by the Earth Climate System [48], should have contributed to the RSLC as observed on a global scale [9]. The lack of correlation between our η_T, η_H time series and the RSL do not allow us to use them for the calculation of their contribution to the sea level trends. It could be argued that the available thermosteric and halosteric time series are not able to reproduce the monthly and inter-annual variability of sea level, but they still could capture its long-term variability. In that case, we could simply subtract the thermosteric and halosteric trends to the sea level trend to obtain that of the mass component. This would be equivalent to accept that the coefficients relating η_T and η_H to sea level are lower than 1 for the monthly and inter-annual time scales, but they are equal to 1 for the long time scales. Similarly, the mass of salt could be directly subtracted to the mass component to assess the change of mass of freshwater. When these calculations were carried out, the results were unrealistic, with negative freshwater contributions in some cases.

In the case of the Canary Islands, for this period of time (1948-2019) the RSL increased at a rate of 1.09 ± 0.14 mm/yr, with a mass contribution of 1.3 ± 0.3 mm/yr. Once again the analysis of the steric and salt contributions did not yield significant results. Finally, in the case of the southern Atlantic Iberian coast and that of the Mediterranean Sea, the sea level trend was 1.35 ± 0.18 mm/yr, which is similar to that obtained for the Western Mediterranean by [6] (1.2 ± 0.14 mm/yr). The thermosteric and the mass of salt contributions were positive and the halosteric one was negative yielding contributions of 1.6 ± 0.4, 0.15 ± 0.05 and 1.5 ± 0.4 mm/yr for the mass, mass of salt and mass of freshwater contributions respectively. These results are more in agreement with other works that estimate that the addition of fresh water on a global scale is higher than 1 mm/yr [49] or 1.31 mm/yr [9]. Nevertheless, it should be taken into account that these latter estimations correspond to the 1990s decade and the first decade of the present century, whilst our results are calculated for a longer period of time.

The discussion above suggests that the available temperature and salinity data are not suitable for the analysis of the monthly, inter-annual, and long-term variability of the steric and mass of salt components of sea level. This would support those results in [19,47]. [46] also showed that the available temperature and salinity data collected in the Mediterranean Sea during the second half of the twentieth century and the first decades of the twenty first one, could not capture the long-term trends of these variables at the upper layer of the sea, because of the scarcity of data combined with the large natural variability of this layer. Despite the data scarcity problem, it should be considered that these tide gauges are not provided with GNSS receivers. As no altimetry data are available for comparison during this period, other problems such as vertical land movements, or changes in the location of the instruments, that have not been properly documented, could not be discarded.

Both tide gauges and altimetry data are available for the period 1993-2019. The comparison of these results shows that trends estimated from altimetry data are quite homogeneous (see Table 5 and Table 6 and blue lines in Figure 4). When averaging for the whole area of study, the RSL increased at a rate of 2.8 ± 0.3 mm/yr with mass, salt and freshwater contributions of 2.8 ± 0.5, 0.10 ± 0.09, and 2.7 ± 0.6 mm/yr. The results obtained from the analysis of tide gauge data are much more variable reflecting possible local effects or even errors in the instrumentation. However, these effects are canceled when the results are averaged, and the trends estimated for RSL, mass, salt, and freshwater contributions are similar to those calculated from altimetry data: 2.5 ± 0.7, 2.7 ± 1.1, 0.25 ± 0.25 and 2.5 ± 1.4 mm/yr respectively.

During this recent period there are no contributions from the atmospheric forcing. Although these variables were significantly correlated with the observed sea level (Table 4), the atmospheric pressure did not show any long-term trend, whereas the U and V components of the wind increased in some cases and decreased in others. The agreement between the altimetry and tide gauge results confirm the acceleration of the sea level trends during recent decades, but once again the observed sea level was not significantly correlated to the thermosteric and halosteric components in most of the cases. Consequently, the linear model is not able to estimate the steric contribution, nor it can estimate the mass, salt, and freshwater ones. Nevertheless, it is worth considering the particular case of L'Estartit. It is well known that the WMED has undergone a warming process along the twentieth century which has accelerated during the beginning of the twenty first one. Therefore, the thermosteric component should contribute not only to the monthly and inter-annual variability of the sea level variability, but also to its long-term trend (see column #6 in table 6), as observed on a global scale [7,8,9]. On the other hand, this sea has also suffered an intense salinity increase [18,46,50] and there should be a negative contribution of the halosteric component and a positive one of the mass of salt (column #7 in table 6 and Figure 4D). In this case all the predictors, both those representing the atmospheric forcing (P, U, V) and the thermosteric and halosteric components, contributed to the observed sea level variability (see Table 4 and Table 6). Furthermore, the linear model has a multiple correlation coefficient of 0.88, which means that the model is able to explain the 77 % of the sea level variance. Hence, we can be confident that the different contributions to the long-term trends of the sea level are accurately estimated. These results indicate that the sea level increased at L'Estartit at a rate of 2.7 ± 0.8 mm/yr since 1993. The atmospheric forcing is significantly correlated with the RSL and explains part of the monthly and inter-annual variability, but as these variables did not experience any long-term trend during the period 1993-2019, they did not contribute to the sea level linear trend. The thermosteric and halosteric contributions were positive and negative, as expected, and the resulting mass addition had a positive trend of 3.3 ± 1.9 mm/yr. The mass of salt contribution was 1.9 ± 1.0 mm/yr and the freshwater contribution was 1.4 ± 2.9 mm/yr for the tide-gauge data, and 1.3 ± 2.0 for the altimetry ones. The large uncertainty in the latter component arises from the formula for the expansion of errors. As the trend for the mass of freshwater is derived from the subtraction of the mass of salt from the mass component, its uncertainty is the sum of both errors. However, the obtained value is close to those recently reported, based on the melting of glaciers and Greenland and Antarctic ice sheets [9,49] who estimated this contribution as 1.31 mm/yr. Furthermore, in the case of L'Estartit, the estimations of the sea level linear trends, and its different components, are almost the same when calculated from tide gauge or altimetry data. Our hypothesis is that the good behavior of the statistical model for L'Estartit arises from the quality of the data. First, the tide gauge started operating in recent decades (1990) with no location changes or any other known problem. Beside this, its location is very close to an area very intensively sampled by different monitoring programs [51], and has received considerable attention because of its proximity to the area of formation of Western Mediterranean Deep Water.

In summary, averaging for all the tide gauges, it can be stated that the sea level increased at a rate of 1.58 ± 0.19 mm/yr from 1948 to 2019. The large dispersion of these results, based on the analysis of tide-gauge data, makes it difficult to estimate differences between the three large geographical areas analyzed: Atlantic Iberian Peninsula, Canary Islands, and Southern Peninsula and Spanish Mediterranean Sea. The trends of sea level accelerated during the period 1993-2019. The results for the three areas and both for tide-gauge and altimetry data are very similar (indistinguishable within the uncertainty level) and range from 2.3 ± 0.8 to 3.0 ± 0.3 mm/yr. The available data do not allow us to estimate the thermosteric and halosteric contributions in a reliable way, and therefore the mass and salt components cannot be estimated. According to previous works, this seems to be the result of the data scarcity. On the contrary, the results from L'Estartit, which is located in a well sampled area, allowed us to estimate the positive and negative contributions of the thermosteric and halosteric components of sea level, and the mass of salt. The mass of freshwater at this location increased at a rate of 1.4 ± 2.9 mm/yr for tide gauge data, and 1.3 ± 2.0 mm/yr for the altimetry data, in agreement with recent observations derived from glacier ice melting. These results evidence once more the importance of the monitoring systems of the oceans for estimating the different contributions to the present sea level rise.