Steric Sea Level Changes from Ocean Reanalyses at Global and Regional Scales

Abstract

Sea level has risen significantly in the recent decades and is expected to rise further based on recent climate projections. Ocean reanalyses that synthetize information from observing networks, dynamical ocean general circulation models, and atmospheric forcing data offer an attractive way to evaluate sea level trend and variability and partition the causes of such sea level changes at both global and regional scales. Here, we review recent utilization of reanalyses for steric sea level trend investigations. State-of-the-science ocean reanalysis products are then used to further infer steric sea level changes. In particular, we used an ensemble of centennial reanalyses at moderate spatial resolution (between 0.5 × 0.5 and 1 × 1 degree) and an ensemble of eddy-permitting reanalyses to quantify the trends and their uncertainty over the last century and the last two decades, respectively. All the datasets showed good performance in reproducing sea level changes. Centennial reanalyses reveal a 1900–2010 trend of steric sea level equal to 0.47 ± 0.04 mm year⁻¹, in agreement with previous studies, with unprecedented rise since the mid-1990s. During the altimetry era, the latest vintage of reanalyses is shown to outperform the previous ones in terms of skill scores against the independent satellite data. They consistently reproduce global and regional upper ocean steric expansion and the association with climate variability, such as ENSO. However, the mass contribution to the global mean sea level rise is varying with products and its representability needs to be improved, as well as the contribution of deep and abyssal waters to the steric sea level rise. Similarly, high-resolution regional reanalyses for the European seas provide valuable information on sea level trends, their patterns, and their causes.

1. Introduction

1.1. Sea Level Rise as a Proxy of Climate Change and Threat from Global Warming

Sea level rise has been long used as an indicator to monitor the Earth’s energy imbalance and the associated global warming [1,2], and simultaneously represents one of the most prominent threats that climate change poses to the large human population inhabiting coastal [3]. Mean sea level rise has been recognized as the main driver for changes in sea level extremes [4,5].

At the global scale, sea level rise is mostly related to the expansion of water masses due to density variations (steric effect) and to the increase of water masses due to ice-sheet and glacier melting, and redistribution of the water cycle between different Earth system components. On the contrary, regional sea level changes depend mostly on oceanic transports of heat, salt, and water masses among different basins [6]. At long time scales (e.g., multi-decadal to centennial), other non-negligible factors exist for sea level changes, such as perturbations in the Earth’s rotation [7] and the effect of geothermal fluxes [8].

The most recent estimates of sea level rise during the altimetry period (1993–2017) consider the trend of the global mean sea level (GMSL) to be equal to 3.35 ± 0.4 mm year⁻¹, with a statistically significant acceleration over such a period equal to 0.12 ± 0.07 mm year⁻² [9]. Through a compilation of several datasets, the contributions of ocean thermal expansion, glaciers, and Greenland and Antarctica ice melting to the global mean sea level rise over the period from 1993 to present have been quantified as 42%, 21%, 15%, and 8%, respectively [10,11]. In the last decade, the water mass contribution has increased up to 70% in 2014 [12].

Sea level rise during the 21st century is expected to be larger than during the 20th century [13,14,15,16,17] even if greenhouse gas emissions stopped now [18] The regional increases may locally amplify further [19] due to, e.g., weakening of the thermohaline circulation [20], and have implications of unprecedented flooding risks in coastal areas [21], although large uncertainties are still present in the climate projections of mass-induced sea level changes [22]. It is therefore imperative to continuously assess and monitor sea level changes at both global and local scales using all possible tools that the climate science community has developed, such as in-situ and/or remote sensing-based observational products, model-based products, and products fusing models and observations (e.g., reanalysis), while taking care of the associated uncertainty estimates [23].

1.2. Methods of Estimation of Sea Level

Several approaches exist to measure or estimate sea level inter-annual variability and its components. For instance, Church et al. [24] reviewed current methodologies to estimate the steric component of the sea level. Direct estimations of GMSL change rely on in-situ profiles (XBTs, CTDs, and Argo floats) and undergo mapping procedures (i.e., objective analyses) to provide valuable information on regular grids. Examples of these analyses were initially provided by, e.g., Levitus et al. [25] and Ishii et al. [26]. Although during the recent decade the Argo floats have reached a quasi-global coverage (with the exclusion of shallow and coastal areas and high-latitude regions), such analyses generally suffer from the lack of observational samplings during the pre-Argo era, for instance, in high latitude regions [27] where global warming may be amplified [28]. Furthermore, the non-negligible impact of the choice of the mapping method and reference climatology [29] may also reduce the accuracy of objective analyses. Errors in historical instruments also affect the resulting trend significantly [30], although bias correction procedures for historical datasets are at the moment well considered [31].

Another observation-based approach of sea level estimation is the tide gauge-based sea-level reconstruction [32,33], where tide gauge data are used to estimate large-scale patterns of variability, generally based on statistical relations (e.g., empirical orthogonal functions, EOFs) computed over well-observed training periods (e.g., altimetry). Such reconstructions have the crucial merit to offer a long multi-decadal reconstruction owing to the availability of long tide gauge data records; on the other hand, different reconstructions show significant discrepancy at regional scales [33].

Remote sensing has also provided a fundamental input to sea level investigations mostly thanks to altimetry missions providing high-quality sea level anomaly dataset with a nominal precision of 2 to 5 cm since late 1992. Recently, attention has been given to climate records of sea level from altimetry data [34], aiming to produce a homogenous series of sea level anomalies by cross-calibrating the different satellite missions launched in the last decades. Their accuracy largely depends on the accuracy of the atmospheric (e.g., dry and wet tropospheric corrections) and geophysical corrections applied to the data [35]. For instance, Forget and Ponte [36] decomposed the altimetry error budget in several sources corresponding to their geophysical corrections, finding that the monthly anomaly error budget is dominated by satellite calibration and instrumental errors, and to a lesser extent by representation errors due to mesoscale activity at the global scale. At the regional scale, errors in the tidal correction in specific coastal areas and inverse barometer error in the Antarctic area, due to the lack of a rich synoptic station network, become important.

Gravimetry missions, such as the Gravity Recovery and Climate Experiment (GRACE) launched in 2002, infer the evolution of the mass component of sea level. Although gravimetry-derived equivalent sea level data requires specific filters to destripe the original gravimetric data records, they provide estimates of the mass-induced sea level variations with good accuracy except in coastal regions. Combining in-situ ocean profiles with altimetry and gravimetry data has recently allowed the sea level budget to be closed at both global and regional scales [37,38,39], and the difference between altimetry and gravimetry has allowed validation of the ocean model-based steric sea level estimates [40,41].

An alternative way to study sea level variability is through ocean general circulation models (OGCMs) that integrate in-time primitive equations of the ocean and are forced at the sea surface by atmospheric datasets, generally taken from atmosphere general circulation models (AGCMS), which are possibly complemented by assimilation of atmospheric observations. Ocean models are proven effective in reproducing regional modes of sea level variability [42,43], and helpful for understanding the causes of sea level rise, but often lack of realism [44], especially at high latitudes [45]. Coupled global circulation models (CGCMs) are other primary tools to provide centennial trend estimates [46] and sea level projections for future decades [47].

Ocean reanalyses (or syntheses) fall in between observation-based products and ocean models, as they combine information from multiple observing networks with ocean model dynamics and atmospheric forcing through multivariate data assimilation methods. Their solution is thus dynamically consistent and accounts for the information contained in the assimilated observations. Reanalyses have shown increasing maturity and reliability over the last decade [48] and are routinely used for ocean state monitoring [2] and climate investigations [49]. Ocean state estimates can be considered a subcategory of ocean reanalyses, which provide a fully consistent four-dimensional picture of the ocean [50] by optimizing the model solution over the entire reanalysis period (long assimilation time windows). Most reanalysis methods rely instead on sequential data assimilation, which can in turn introduce some discontinuity unless they are explicitly considered as budget terms [51]. Indeed, the application of analysis increments into ocean models is in general non-conservative, and may lead to some spurious high-frequency variability (i.e., at the frequency of the analysis increment application, generally daily to weekly). However, when focusing on low-frequency changes (annual or longer scales), as done in most of the present work, the discontinuities arising by data assimilation become negligible. Recent comparisons of ocean reanalyses [52] show considerable consistency between different products, even in under-sampled high-latitude regions [53], motivating their use in climate applications.

1.3. Overview of Reanalyses Performances in Capturing Sea Level Variability

Early attempts of using ocean reanalyses to quantify the sea level trend suggested that sea level rise has undergone acceleration in the past decades compared to the centennial trend. This was shown about 15 years ago by Carton et al. [54]. Since then, the use of reanalyses for quantifying global and regional sea level trends has widely spread, in particular to investigate the roles of steric and mass components. The global sea level mass component from reanalyses is somehow less exploited, likely due to the fact that river and land-ice discharge are usually taken from climatological dataset, thus being less reliable in terms of inter-annual variations.

The use of reanalyses in capturing the steric sea level change and variability offers some advantages compared to objective analyses, which is the statistical mapping of observations without a dynamical ocean model. Indeed, additional information coming, e.g., from atmospheric forcing, multiple observing networks, and model physics may overcome to some extent the lack of in-situ observation sampling. Obviously, the synthesis of all these different sources of information implies that errors can potentially propagate from atmospheric forcing, etc. to the reanalysis, which was found non-negligible during unobserved periods or over poorly sampled regions [55,56]. Nevertheless, the accuracy of reanalyses has greatly improved over time. For instance, Storto et al. [41] showed that new vintages of the CMEMS (Copernicus Marine Environment Monitoring Service) reanalyses exhibit improved steric sea level inter-annual variability compared to their predecessors, and their consistency increases especially in under-sampled regions, such as the Southern Ocean. Thus, state-of-the-science reanalyses have acquired the confidence of the climate community.

Within the Ocean Reanalysis Intercomparison Project (ORA-IP, Balmaseda et al. [57]), Storto et al. [40] compared steric sea level from an ensemble of ocean reanalyses and objective analyses against sea level data inferred from altimetry and gravimetry. The results highlighted the better performances of ocean reanalyses compared to objective analyses in reproducing the inter-annual variability of extratropical steric sea level changes, while the reproduction of seasonal variations is comparable between the two families of products. The ability of reanalyses was found particularly evident in the Southern Ocean, likely due to the sparse in-situ observational network and the enhanced mesoscale activity therein, whose effects are compensated in reanalyses by the dynamical information about barotropic and baroclinic forcing variability coming from the atmospheric reanalysis forcing.

Chepurin et al. [58] also showed that an ensemble of ocean reanalyses generally exhibits good consistency with tide gauge records, which represent independent sea level measurements because tide gauge data are not assimilated in reanalyses. They showed that the local mismatch between reanalyses and tide gauges is mostly due to meteorological forcing. However, during the pre-altimetry period (i.e., before 1993), small consistency is found between reanalyses and tide gauge-based sea level reconstructions, the latter generally outperforming the former in reproducing the inter-annual variability measured by tide gauges [33]. This might be partly attributed to the fact that some tide gauge data records are not representative of large-scale dynamics.

Pfeffer et al. [59] showed that multi-decadal variability linked to prominent climate indices is well explained in the ORAS4 reanalysis [44] but not in objective analyses, which in turn perform well enough only after the full Argo floats deployment. Objective analyses are only able to explain the variability of sea level associated to shorter than multi-decadal time scales, such as ENSO, while reanalyses also capture low-frequency (multi-decadal) variability. For instance, Calafat et al. [60] successfully used ocean reanalyses rather than observation-only products to analyze the baroclinic contribution to multi-decadal sea level change in the eastern North Atlantic Ocean and Mediterranean Sea.

On the other hand, ocean reanalyses prove significantly superior to free-running ocean model simulations in reproducing altimetry-derived trends. For instance, Köhl et al. [45] showed that unrealistic trends in ocean simulations, especially at high latitudes, are rectified in the ocean reanalysis, where data assimilation modifies the steric sea level variations and, to a lesser extent, the mass-induced variations.

In summary, reanalyses have shown accuracy in capturing sea level variability at least comparable (if not superior) to observation-only or model-only products, while their partial dynamical consistency allows for process-oriented and research-driven investigations to some extent.

1.4. Using Reanalyses to Partition Causes of Sea Level Variability

Beside the good performance in reproducing sea level variability for several regions and periods, ocean reanalyses offer the possibility of investigating in detail and attributing inter-annual sea level changes to different mechanisms, partitioning mass versus steric contributions, and investigating individual ocean process contributions. This is otherwise non-trivial by using observation-only products, which are not physically consistent by construction.

Several studies exist that have used ocean state estimates, such as ECCO [61] to attribute sea level changes. At the global scale, Piecuch and Ponte [62] showed that low-frequency changes in the global mean steric sea level are due to small imbalances between the compensating effects of atmospheric forcing and ocean transports, which can be investigated and quantified through the evaluation of reanalysis output fields. Piecuch and Ponte [63] exploited the dynamical consistency of the ECCO reanalysis to show how ocean transports are responsible for steric sea level variations rather than surface buoyancy fluxes at the global scale. Forget and Ponte [36] used the ECCO reanalysis to investigate the partitioning of steric and barystatic sea level and inferred their complementary impact on sea level changes, with the former dominating low-frequency variability (>1 year) and the latter more prominent at high frequencies (<3 months). At the regional scale, advection and diffusion also play an important role, and surface heating may induce significant variations as well, especially in the tropical regions [64]. The redistribution of water mass may be further identified as the predominant process for regional sea level decadal changes [65], although regional exceptions (e.g., in the western tropical Pacific and eastern tropical Indian oceans) exist. By separating several sea level components in an ocean state estimate, Köhl [66] suggests that the contribution from isopycnal motion to the thermosteric sea level below the mixed layer dominates the steric sea level inter-annual variability. These latter studies suggest the intrinsic low predictability of sea level changes on the multi-decadal timescale, because of the significant uncertainty of water mass and property redistribution through oceanic transports.

Compensation between temperature and salinity through adiabatic processes is also quantified through the use of ocean reanalyses: Wunsch et al. [67] found that, at the global scale, about 25% of the temperature contribution to the steric sea level is compensated by salinity changes. Locally, e.g., in the North Atlantic Ocean, the thermal and haline compensation becomes even larger and predominant [20,68], and the correct partitioning between temperature and salinity in reanalyses greatly improves after full deployment of Argo floats in the North Atlantic basin [41].

The goal of the present work was to assess the ability of state-of-the-science ocean reanalyses in representing sea level variability and trend signals, while other works in the present Special Issue focus on estimates from objective analyses of observations only. This article is structured as follows: after this introduction, the material and methods are presented in Section 2. Section 3 summarizes sea level change results from centennial (20th century) ocean reanalyses. Section 4 presents results on contemporary sea level changes from state-of-the-science reanalyses, focusing on trends, inter-annual variability, and partitioning the trends between mass and steric components, at both a global and regional scale. Section 5 summarizes and concludes.

2. Material and Methods

2.1. Formulation of the Sea Level Balance

The sea level anomaly evolution can be formalized from the hydrostatic relation [69] as Equation (1), meaning that changes in sea level anomaly ($\eta$) can be expressed as the sum of three components: (i) The atmospheric pressure term, ${\eta }_a$ (or inverse barometer effect, IB), which depends on the fluctuations of the atmospheric pressure, $p_a$, around its mean $\overline{p_a}$; (ii) the steric component, ${\eta }_s$ (local or global), which is approximated by the vertical integral of the negative density anomaly, $\delta \rho$, accounting for the expansion and contraction of the water column, where the constant volume approximation is adopted (i.e., vertical integration from $\eta =-h$ to $\eta =0$, see, e.g., Yin et al. [47] and Griffies and Greatbatch [70]. Note that the “non-Boussinesq steric effect” is assumed negligible [70] and is however not possible to calculate with current ocean models where the Boussinesq approximation is applied; and (iii) a mass component term, ${\eta }_b$ (or “barystatic” when globally averaged), proportional to the change of the ocean bottom pressure, $p_b$, which responds to the changes in water mass (e.g., continental ice melting/growing and, locally, barotropic motions):

(1)$\eta ={\eta }_a+{\eta }_s+{\eta }_b=\frac{\overline{p_a}-p_a}{g{\rho }_0}-\frac{1}{{\rho }_0}{{\displaystyle \int }}_{\eta =-h}^{\eta =0}\delta \rho \left(T,S\right) dz+\frac{p_b}{g{\rho }_0},$

where $g$ and ${\rho }_0$ are the gravitational acceleration and the reference seawater density, respectively. The inverse barometer effect has practically no contribution at the global scale and it is not further analyzed, although locally it may exhibit a non-negligible contribution to the sea level trends [71,72]. The steric sea level component can be further approximated to the sum of the thermo (${\eta }_t$) and halo (${\eta }_h$) steric components by computing the density anomaly at constant time-averaged salinity ($S^{*}$) or temperature ($T^{*}$), respectively, i.e.,:

(2)${\eta }_s\simeq {\eta }_t+{\eta }_h= -\frac{1}{{\rho }_0}{{\displaystyle \int }}_{\eta =-h}^{\eta =0}\delta \rho \left(T,S^{*}\right) dz-\frac{1}{{\rho }_0}{{\displaystyle \int }}_{\eta =-h}^{\eta =0}\delta \rho \left(T^{*},S\right)dz.$

Although the seawater equation of state, $\delta \rho \left(T,S\right)$, depends on the pressure as well (pressure-steric sea level contribution), such a contribution to the steric sea level anomaly is assumed negligible and is not considered in Equations (1) and (2) [73]. Furthermore, differences between ${\eta }_s$ and the sum of ${\eta }_t+{\eta }_h$ can be assumed negligible [44], that is, we consider ${\eta }_s\simeq {\eta }_t+{\eta }_h$. The global average of the previous equations provides the formulation of GMSL. Note that all the models used in the present study adopt Boussinesq’s approximation and are volume-conserving models. This implies that the global steric effect is calculated as a posterior correction as in Greatbatch [74], Griffies and Greatbatch [70], and Griffies et al. [73] and eventually added to the model output of sea surface height. Depending on the specific reanalysis, the global steric effect is already included in the sea surface height (e.g., Estimating the Circulation and Climate of the Ocean version 4 release 3, ECCOv4r3) or has been added during the post processing of the data.

2.2. Datasets Used in the Present Work

We devote this section to the introduction of the reanalyses used throughout the paper and the observation-based products employed for validation. For investigating centennial trends, we chose several reanalysis products available for the study period. In particular, four reanalysis products were used (see Table 1, first part); for each reanalysis product, we extracted 4 members, for a total of 16 ensemble members used hereafter. In practice, the first four members were extracted from each dataset. Note that, while some of these products can provide more than four members, we chose this number to equally weigh the multi-model members, with four being the minimum common ensemble size across the different products (from CMCC Historical Ocean Reanalysis, CHOR, see Table 1 first part).

The four families of reanalyses generally differ for the atmospheric forcing, the spatial resolution, the ocean model version and configuration, the data assimilation scheme, and observational datasets assimilated. Thus, the multi-model ensemble implicitly spans several sources of uncertainty, although a detailed discussion on the sources of uncertainty represented therein is out of the scope of the present work. In order to assess the steric sea level trend from historical reanalyses, we analyzed the ensemble mean and ensemble standard deviation (or spread) coming from this 16-member multi-model ensemble.

Tide gauge data used in Section 3.1 for the validation of historical reanalyses were obtained as referenced yearly means from the Permanent Service for Mean Sea Level (PSMSL, Holgate et al. [86]) and only the 39 stations with valid data for the period 1920–1990 were chosen for the comparison.

Regarding the contemporary reanalyses (Section 4.1), we used the combination of altimetry and gravimetry data to validate the steric sea level from reanalyses and measure the increase of accuracy of reanalyses across subsequent versions of reanalysis products. To this end, we evaluated the global steric sea level evolution by subtracting the barystatic signal taken from the gravimetry datasets (GRACE RL05, Chambers and Bonin, [87]) to the global sea level signal provided by the reprocessed altimetry dataset from the ESA CCI project [34]. Such timeseries (called SAT) was compared with the latest combination of reanalysis ensemble mean (called REA) and with the ensemble of reanalyses produced before 2013 and published in Storto et al. [40], called ORAIP. Results also report the ensemble mean from the subsampled members that are present in both REA and ORAIP (called “REA subsampled” and “ORAIP subsampled”) for sake of comparing the same family of products, in terms of ensemble size and combinations of ocean models and data assimilation methods. A similar comparison was also performed by Storto et al. [41], but with a smaller subsample of products that included only the CMEMS reanalyses.

The steric component of sea level in contemporary reanalyses (Section 4.2 to Section 4.3) was then estimated and inter-compared from five reanalysis products (ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2), which were chosen because they are public datasets providing both temperature and salinity, and sea surface height. The second part of Table 1 summarizes the main characteristics of these reanalyses. Observed GMSL from satellite altimetry (AVISO, Ablain et al. [9]) was also used to verify the sea surface height (SSH) of ocean reanalyses. The glacial isostatic adjustment (GIA, e.g., Peltier [88]) is corrected in such an altimetry record of GMSL.

Finally, we introduced the reanalyses for the European Seas: Baltic Sea (BAL), North Western Shelf (NWS), Iberian-Biscay-Ireland (IBI), and Mediterranean Sea (MED) in the third part of Table 1. Each dataset considered here represents a state-of-the-science ocean reanalysis, which takes into account all available and quality flagged observational data (in-situ and remote-sensing) during the reanalysis integration. Data are freely available through the Copernicus Marine Environment and Monitoring Service (CMEMS) web portal [89] and their accuracy is systematically assessed against independent observations, as documented in the CMEMS quality information documents (QUIDs; [90]). The specifics for each re-analysis considered are listed in Table 1 (third part). A 24-year period, from 1993 to 2016, was selected in our analyses to overlap the time-window spanned by the different reanalyses and to cover almost the entire satellite era (1993–ongoing; [91]).

3. Steric Sea Level Trends in Centennial Reanalyses

3.1. Preliminary Validation

In this section, we assess and estimate the steric sea level trends using centennial reanalyses that span the period from 1900 onwards. Given the large uncertainties in early periods already mentioned, using a multi-model ensemble of historical reanalyses appears a meaningful choice in order to mitigate drifts and systematic errors.

Although validating steric sea level from centennial reanalyses is obviously difficult due to the lack of global observations of sea level before the altimetry era, we briefly provide an assessment of the ensemble mean steric sea level, against yearly mean sea level values from long-record tide gauge measurements. As we compared steric sea level from reanalyses against full sea level signal from tide gauges, time series in both observations and reanalyses were detrended in order to exclude the contribution of the mass component from the comparison. Doing this, we assumed that only the steric component has a significant effect on the spatially-varying total sea level variability captured by the tide gauges, which may not hold for all regions (e.g., in the Southern Ocean). Removing the trend implies that the centennial linear thermal expansion was also removed, namely only the consistency of the inter-annual fluctuations between these two datasets were assessed in this exercise, while no assessment can be performed on the overall trends. Furthermore, steric sea level from reanalyses does not include coastal dynamics, which may further reduce the consistency between tide gauge data and steric sea level from reanalyses.

Figure 1 summarizes the comparison, showing green circles in the correspondence of tide gauge stations that exhibit statistically significant correlations with the ensemble mean steric sea level. Values from reanalyses are taken from the grid-point nearest to the tide gauge locations. Red and black circles correspond to the non-significant correlation case, where at least one or no ensemble member has a significant correlation, respectively. There is indeed one station (south west of the British Island) where no member exhibits a significant correlation, likely meaning that such a station is representative of, e.g., coastal processes that are not resolved in the reanalyses. In general, the comparison shows that the ensemble mean is able to capture the sea level variability spanned by tide gauges during the 20th century, except for a few stations off the Netherlands and Germany (North Sea) and in the Baltic Sea. For most stations located in Europe, North America, Indian, and Southern oceans, correlations are significant and offer a qualitative justification for the use of centennial reanalyses in sea level trend investigations. A few stations in Scandinavia, the eastern coast of the US, and Oceania have correlation values exceeding 0.6.

3.2. Centennial Trends

We now evaluate the sea level trends from the centennial reanalyses. We focused our investigation on the steric sea level, because the barystatic signal is potentially less reliable due to the use of climatological (rivers) or potentially biased (precipitation) datasets for the freshwater input. Figure 2 shows the globally averaged steric and thermosteric sea level (Equation (2)) for the entire common period of the 20th century reanalyses (1900–2009), separately for different vertical ranges (0–300 m, 0–700 m, 0–2000 m, and 0 m–bottom). Black lines are from the ensemble mean of the 16 reanalysis members, while dark or light shaded areas correspond to 1 or 3 ensemble standard deviations, respectively. Centennial trends for the four depths are 0.14 ± 0.02 (0–300 m), 0.37 ± 0.02 (0–700 m), 0.58 ± 0.04 (0–2000 m), and 0.47 ± 0.04 (0 m–bottom) mm year⁻¹, where trend uncertainty was calculated with the 95% confidence level estimated by means of ordinary bootstrapping with 1000 replicates [92]. Note that for the total column (full depth) and the periods 1900–1990 and 1971–2010, trends are equal to 0.38 ± 0.04 and 0.86 ± 0.13 mm year⁻¹, respectively, which are in close agreement with the estimates provided by the IPCC AR5 report (Table 13.1 in Church et al. [24]) using CMIP5 (Coupled Model Intercomparison Project) coupled models (equal to 0.37 ± 0.30 and 0.96 ± 0.45 mm year⁻¹, respectively). For the latter period of 1971 to 2010, the IPCC AR5 report indicates that the observed trend of the steric sea level is equal to 0.80 ± 0.45 mm year⁻¹, in close agreement with the 0.86 ± 0.13 mm year⁻¹ trend from the historical reanalysis ensemble. This indicates the close agreement of reanalysis-based estimates with multi-model ensemble, and, during the well-observed period, with observation-based estimates.

All vertical ranges exhibit an increase of steric sea level during the study period, except during the first two decades of the 20th century. This latter feature may be partly due to the signal of the Santa Maria volcanic eruption (1902) and partly to suboptimal reanalysis initialization [56]. The top 300-m sea level shows sea level fluctuations around the years of World War II, due to the lack of observations during this period. It can also be noticed that sub-decadal variability is enhanced in the 0–300 m layer, especially during the period 1940–1965 due to the effect of an intermittent observing network as a consequence of World War II and the development of the XBT-based global observing network afterwards. However, such variability is smoothed away for deeper ocean.

The total column steric sea level shows a nearly constant increase of sea level during the whole reanalyzed period after the first two decades. The comparison between total steric (thin dotted lines) and thermosteric sea level (thick solid lines) indicates the impact of the halosteric component. It turns out that for all vertical regions except the top 300 m, the contribution of the halosteric sea level is negligible at the global scale, and the steric and thermosteric sea level rises tend to coincide. However, for the 0–300 m vertical range, the halosteric presents a weak decrease over the century, leading to a steric sea level rise slightly smaller than the thermosteric one. With the steric sea level calculated for each product as an anomaly with respect to its long-term mean, the ensemble spread (shaded areas in Figure 2) is close to 0 by construction near the middle of the time series and increases towards both the beginning and the end of the timeseries. Thus, it can only be used to detect some year-by-year spread change rather than estimating decadal uncertainty.

In order to deepen the investigation about the steric sea level trends from historical reanalyses, Figure 3 shows in particular the 20-year running trends calculated from the reanalysis ensemble mean. In particular, for each year between 1910 and 2000, the trend centered over 20 years is shown as a timeseries (running trend), with the aim of highlighting periods of increasing or decreasing steric sea level trend. Both thermosteric and total steric trends are shown for the same vertical regions presented in Figure 2. Moreover, red lines represent the running trends computed over the global ocean bounded between 60° S and 60° N, in order to visualize the possible impact of the high latitude areas on the results. The bottom area of each of the four panels of Figure 3 shows the statistical significance, as explained in the caption of Figure 3. We evaluated both the robustness of the ensemble mean trend and the dispersion of the trends from the individual members.

In general, there are three periods of maximum steric sea level increase, which correspond to the 1925–1945 period, the 1970s, and the uninterrupted period from the beginning of the 1990s onward. Considering the total column panel (Figure 3, bottom right panel), the trends are always positive after World War I, almost always statistically significant based on the Mann–Kendall test, and with the ensemble mean trend exceeding the spread during the periods of maximum steric expansion mentioned above. Slowdowns of trends occur after notable volcanic eruptions, in particular following climate-impacting eruptions in 1902, 1963, 1982, and 1991 [93]. The maximum value of the trend is the latest trend (1990–2009), equal to 1.5 mm year⁻¹; the increasing trend from 1990 indicates that the historical reanalysis ensemble exhibits acceleration in the trend from the 1980s. From approximately 1995, the total steric sea level rise reaches an unprecedented rate (about 1 mm year⁻¹). Total column trends also suggest marginal differences between thermosteric and total steric trends, namely the halosteric contribution is negligible. Small differences between the global ocean trends and those for the region bounded between 60° S and 60° N are visible. The 60° S to 60° N trends have values systematically larger than the global ocean ones, indicating that in the centennial reanalyses, the high latitudes slightly compensate the low- and mid-latitude sea level rise.

The region 0–2000 m behaves closely to the total column trends and indicates the small contribution of the water masses below 2000 m. Due to the better constraining effect of the in-situ observing network, the layers 0–300 and 0–700 m are characterized by longer periods of significant trends based on the ensemble spread test (orange to purple years). In the 0–300 m, the slightly negative trend during the period 1945–1970 suggests a vertical redistribution of heat during these years, with positive trends below 300 m of depth compensating for the upper ocean contraction. In the top 300 m, the 20-year periods centered between 1930 and 1935 exhibit a steric sea level increase as high as the latest one centered in 2000. The negative contribution of the halosteric appears also significant, combining indeed with this large thermosteric expansion during the period 1920–1945.

The regional distribution of steric and thermosteric sea level trends is shown in Figure 4. The total column trend map highlights the expansion at low and mid latitudes, with peaks in the North Atlantic (south off the Gulf Stream separation), western Pacific Ocean, and Southern Ocean north of the Antarctic Circumpolar Current (ACC) front, up to 3 mm year⁻¹ during the 1900–2010 period. Especially in the tropical Pacific Ocean, the thermosteric sea level displays a positive trend in the western and a negative trend in the central eastern tropical Pacific Ocean. This signal of thermosteric sea level from the ensemble mean historical reanalyses is consistent with a previous study [95] showing a positive temperature (surface and subsurface) trend and a negative temperature (surface and subsurface) trend in the central eastern tropical Pacific mainly due to enhanced tropical Pacific circulation associated with strengthening of the tropical trade winds. Negative trends are found at high latitudes in both hemispheres. In the Arctic Ocean, the trends are driven by the contribution of the halosteric component (deductible by the difference when comparing the total steric and the thermosteric); this is found to be dominated by an outlier product, and the corresponding uncertainty and spread is large, therefore there is no confidence about this specific result. By comparing the thermosteric and total steric trend maps, it turns out that the reanalysis ensemble mean reproduces the compensation between the thermo and halo steric components prominently in the Atlantic and Indian oceans [67,68]. In most areas except the Indian Ocean, trends of the ensemble mean are statistically significant with respect to the ensemble spread.

This sketch appears almost identical for the 0–2000 m region. Upper ocean regional trends exhibit qualitatively similar results for the total steric trends, except, e.g., in the Indian Ocean, where the waters below 700 m of depth dominate the total column thermosteric sea level rise.

4. Contemporary Sea Level Change from Ocean Reanalyses

4.1. Assessment of Reanalyses

In this section, we analyze the contemporary sea level change and variability. Before investigating the sea level trends and inter-annual variability, we start the discussion by validating the contemporary steric sea level trends and variability reproduced by global ocean reanalyses. Contemporary steric sea level rise is here defined as the rise occurring during the last two to three decades. Thanks to the availability of satellite data records, independent estimates of steric sea level are available since approximately 2003, when the first accurate records of mass contribution to sea level from gravimetric missions were available (see also Section 1.2 and Section 2.1).

Results of the validation exercise are summarized in Table 2, in terms of correlation against the satellite-derived dataset and inter-annual trends (see Section 2.1 for the nomenclature). The table indicates an improvement of the correlation skill score from ORAIP to REA (from 0.63 to 0.86), which is statistically significant according to Steiger’s Z test [96], during the overlapping period 2003–2010. In this period, REA trends are overestimated with respect to the steric SAT dataset (1.57 ± 0.16 versus 1.05 ± 0.18), although closer to the observed trends than the underestimated ORAIP trend (0.33 ± 0.09).

During the period 2003–2014, REA provides similar correlation (0.87), but with the linear trends as in the satellite dataset. According to the steric sea level variability (Figure 5), the mismatch during the period 2003–2010 seems to be related to a faster rise during the years right after 2003, likely due to the adjustment after Argo floats’ deployment, which is somehow absorbed in the 2003 to 2014 increase rate. Looking at the subsampled ensemble mean products, the REA subsampled confirms the significantly higher correlation compared to ORAIP subsampled (0.84 versus 0.44), and the correctly larger trend. ORAIP subsampled presents an almost neutral rise and Figure 5 indicates that the sea level rise is clearly under-estimated especially in the first years of the timeseries (2003–2007). These results suggest that for the period where independent data are available (from 2003 onward), the latest vintage of reanalyses provides an accurate reconstruction of the steric sea level evolution, and clearly improves its representation with respect to the previous vintage of reanalyses.

4.2. Contribution of Steric Sea Level to Sea Surface Height Trends

In this section, the steric component of the sea level in ocean reanalyses is estimated and inter-compared between five different products (ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2), which were chosen because they are public datasets providing both temperature and salinity, and sea surface height. Table 1 (second part) includes some general characteristics of the reanalyses. The multi-product ensemble mean was also used to evaluate the overall performance of those products. Here, the period was focused on the satellite era (1993–2015) when altimetry data are usually assimilated in ocean reanalyses.

GMSL between 70° S–70° N was first calculated (Figure 6). During 1993–2015, the multi-product ensemble mean of global SSH has a positive trend of 2.5mm year⁻¹, which is 79% of the observed trend for GMSL (3.2mm year⁻¹) for the same period (Figure 6a). This suggests that ocean reanalyses slightly underestimate the GMSL rise rate. In the meantime, there is a divergence for the GMSL rise rate between the ocean reanalyses (Figure 6b–f). GLORYS4 produces the fastest GMSL rise rate (3.8 mm year⁻¹), while SODA reanalyses have the smallest rates (2.1 and 1.9 mm year⁻¹), with ECCO products having a moderate rate (2.5 mm year⁻¹). Part of the difference might be due to altimetry data that are not consistently assimilated in the products.

We then estimated the contributions of steric height to GMSL rise in ocean reanalyses. For the multi-product ensemble mean, the steric component of the sea level has an increasing trend of 1.1 mm year⁻¹ globally (Figure 6a), which accounts for 44% of the total sea level trend. Steric heights at different depths (300 m, 700 m, 2000 m, full depth) were also calculated, with the 300-m depth having the largest relative contribution (0.6 mm year⁻¹). The contribution of steric height is varying within different ocean basins (Figure 7a). In the Atlantic and the Southern Ocean, steric height contributes to 50% to 60% of the total SSH change rate; in the tropical Pacific and Indian Ocean, steric height has a lesser contribution (20%–40% of SSH change rate).

The effect of steric height to total SSH was also evaluated in individual reanalyses. Global steric height (full depth) is increasing during the period in all products (Figure 6b–f), at a rate varying from 0.4 to 1.8 mm year⁻¹. Steric height has a stronger effect on the global SSH trend in SODA reanalyses, while it has a weaker effect in the ECCO reanalyses. In SODA, ocean is thermally expanding at a fast rate close to the total SSH increase rate. In ECCO, the thermal expanding rate of the ocean is less than 20% of the total SSH increase rate. This can be further demonstrated in the basin scale (Figure 7b–f). In SODA products, steric height accounts for more than 80% of the total SSH change rate in the vast majority of the ocean. However, in ECCO products, the largest contribution of steric is reduced to 30% to 50%, which is in the western Pacific and eastern Indian Ocean. Moreover, in the eastern Pacific, the steric component of these products has negative trends, which are opposite to SSH trends. Thus, the contribution of steric to total sea level trend is region-dependent and product-dependent.

We found that the variety of steric contribution in reanalyses is associated with the deep ocean (>2000 m depth) (Figure 6b–f), which is greatly warming in SODA but less warming or even slightly cooling in ECCO and GLORYS4, during the period. This is particularly true in the Pacific basin (not shown). The upper ocean (0–300 m depth) shows similar steric trends between these five products (not shown), especially at low latitudes. This suggests that assimilating ocean profile observations, which are widely observed in the upper ocean, reduces the differences of the ocean state produced by the models. However, in the deep ocean, due to the limited source of observations, ocean analysis makes little change for the steric sea level, which is in turn more likely to retain the variety of ocean models.

Figure 6 also shows that the contribution of the mass component (SSH minus steric) to the SSH trend varies between ocean reanalyses. For example, the ocean mass in ECCO and GLORYS4 is increasing at a higher rate than in SODA. Although the steric effect accounts more for the SSH trend in SODA, negligible contributions of ocean mass cause the SSH trend tot still be far less than other products and observations. More efforts are needed in the future to reconcile the representation of ocean mass in reanalyses.

4.3. Inter-Annual Variability

Inter-annual variability of observed GMSL was reasonably captured by the SSH of multi-product ensemble mean (Figure 8), with a correlation equal to 0.40 (significant at 95% confidence level). The inter-annual variability is best re-produced in ECCOv4r3, with a correlation equal to 0.91. However, when not assimilating altimetry data, SODA products fail to reproduce the inter-annual variability of global SSH. In ECCO-KFS and GLORYS4, the correlation of SSH inter-annual variability is 0.60 and 0.54, respectively.

At the global scale, the contribution of steric sea level to the inter-annual variability of total SSH is essential in ocean reanalyses (Figure 8). Regional SSH variations are provided in the Supplementary Material Figure S1. The correlation between steric height at full depth and global SSH is 0.76 (significant at 95% confidence level) for the multi-product ensemble mean. For individual products, a correlation between steric and SSH anomaly is above 0.8, except ECCOv4r3, where r = 0.14 (not significant at 95% confidence level). The departure of steric height from SSH indicates other effects, e.g., the ocean mass. Like the SSH trend (Figure 6), these effects have more impacts on SSH variability in ECCO products than in SODA products (Figure 8).

It is worth noticing that the steric height inter-annual variability in the upper ocean (0–300 m) is consistently produced in those five ocean products, with an inter-product correlation greater than 0.5. This suggests that data assimilation corrects the modelled background of the ocean, and this correcting effect tends to make the analyzed ocean states from different ocean models resemble the observations. In the upper ocean, where the density of observations is high, steric height converges between products. However, in the ocean interior, where observations are rare, steric height is diverging between products due to the limited constraints of the models. This divergence is larger when ocean models are substantially different, such as those used in ECCO and SODA products.

We also evaluated the effect of steric height (full depth) on inter-annual variability of SSH at regional scales (Figure 9). For the multi-product ensemble mean, the steric component of the sea level explains more than 60% of SSH variability in the vast majority of regions. In the tropical Pacific and Indian Ocean, over 90% of SSH variability is introduced by the local steric expansion. This feature is consistently produced in the reanalyses except GLORYS4 We further confirmed that the steric height above the 300-m depth is mostly responsible for SSH variability in those basins (Supplementary Material Figure S2), and it is significantly modulated by ENSO at inter-annual timescales (Supplementary Material Figure S3).

However, differences between the reanalyses still exist in some regions (Figure 9). In the Southern Ocean, the steric component accounts for less than 50% of sea level variance in GLORYS4, but it does more than 70% to 80% in SODA products. Steric height has an intermediate effect in ECCO products (50%–70%). In the Atlantic, the discrepancy between GLORYS4 and SODA products becomes even larger. Some of those discrepancies are related to thermal expansion in the deep ocean. For example, in the Southern Ocean and the western boundary currents (i.e., the Gulf Stream and Kuroshio Extension), SSH variance is large in SODA products (Supplementary Material Figure S1), and it is largely associated with the steric height below 300-m depths (Supplementary Material Figure S3). However, in ECCO products, steric height lacks variability below the 300-m depth. This might be related to the ocean mixed layer depth, which is not deep enough to resolve the variability below the subsurface.

It is also worth noticing that GLORYS4 produces similar steric heights as SODA products across the basins, both in the subsurface and deep ocean (Supplementary Material Figures S2 and S3). However, SSH variability in GLORYS4 is not explained by steric as much as in SODA (Figure 9). This suggests that other effects, such as ocean mass, are also playing an important role in SSH variability in GLORYS4. This, again, indicates that efforts are needed in the future to vigorously evaluate those effects in ocean reanalyses.

4.4. Sea-Level Trends in Regional Reanalyses for the European Seas

In this section, we investigate the variability of the sea-level trend at the regional scales, starting from multi-resolution reanalyses for the European Seas: Baltic Sea (BAL), North Western Shelf (NWS), Iberian-Biscay-Ireland (IBI), and Mediterranean Sea (MED). In general, sea level trends obtained from multi-resolution reanalyses (Table 3) show positive values, which were in agreement with the most recent estimates obtained for the European Seas from satellite altimetry [1,2].

The sea level trend estimates obtained in the IBI (3.7 ± 0.15 mm year⁻¹) and NWS (2.12 ± 0.38 mm year⁻¹) regions were in the order (and in the uncertainty range) of those obtained in the literature [2] from satellite altimetry data over the same regions and during the same period (1993–2016).

In the Baltic Sea, sea level trends can be strongly influenced by fresh water flux variations [97]. A positive sea level rate was observed also in this basin (2.66 mm year⁻¹), where the halosteric component has a leading contribution compared to the thermosteric sea level trend. The largest uncertainty was observed in this region (>1 mm year⁻¹), which was in line with the findings of von Schukmann et al. [2], who underlined the influence of salt inflow at the entrance of the basin on the mean sea level evolution in the Baltic Sea (Major Baltic Inflow; e.g., Hordoir et al. [98]; Mohrholz et al. [99]). Departures from the sea level trend estimates obtained from satellite altimetry [2] can be explained through GIA-induced sea level rates typically applied to observation-based estimates, which can have a large contribution in the northern part of the basin [100], and which are not represented in the OGCMs physics.

In the Mediterranean Sea, the sea level trend estimation (2.8 ± 0.26 mm year⁻¹) was consistent with the results of previous studies based on in-situ and remote-sensing observations [34,101] when the resultant of the mass and thermosteric components was considered. This was not the case when the halosteric component was included, whose trends have high negative values that can lead to discrepancies when compared with the sea level trends obtained from observations (e.g., remote sensing), as underlined by Legeais et al. [34].

In order to further investigate the sea-level trend variability in the European Seas, a merged high-resolution data-set was obtained interpolating each multi-resolution reanalysis over the spatial grid (1/12˚; 7 km) used by the CMEMS Global Monitoring and Forecasting Center global ocean operational analysis and forecast system (GLO-MFC; Lellouche et al. [102]). During the merging, the sea level mass component in each (multi-resolution) reanalysis was reduced to a common reference considering the offset from a global mean, obtained as the temporal average of the CMEMS GLO-MFC fields over a 10-year period (2007–2017). A 10-year period was selected to consider a temporal window spanning both satellite altimetry [91] and ARGO floats [103] eras.

Figure 10 shows the marked sea level (mass and steric components) trend spatial variability in the European Seas, as obtained from the merged dataset. Considering the steric contribution up to the 700-m depth, sea level trends during 1993 to 2016 ranged between −8 and +8 mm year⁻¹, showing how sea level trends at the regional scale can largely deviate from the global mean sea level trend (3.35 mm year⁻¹; e.g., Ablain et al. [9]). This is in agreement with the findings of Legeais et al. [34], who starting from satellite altimetry data, observed similar sea level trend ranges in the European Seas over the same period.

Considering the steric component up to the 300-m depth, positive homogeneous trends (2–3 mm year⁻¹) were observed in the IBI and NWS regions, mostly due to thermosteric contributions (not shown). It is interesting to notice how sea level trends tended to increase (up to 8 mm year⁻¹) in occurrence of the North Atlantic Drift and Azores Current when the steric component was considered up to a 700-m depth, showing the influence of the main features of the circulation in the upper ocean on the sea level trend spatial variability. Larger positive values, up to 6 mm year⁻¹, were observed in the Baltic Sea considering the steric component (full depth), associated with large uncertainties (up to 3 mm year⁻¹). In this basin, negative trends were also observed in the Gotland basin, due to halosteric contributions (not shown) in the deepest areas of the Baltic Sea.

In the Mediterranean Sea, sea level trends (up to 700-m depth) clearly showed the signature of the most prominent features of the ocean circulation in the basin (e.g., gyres and eddies; Pinardi et al. [104]). Positive trends were found in the Adriatic Sea and Aegean Sea (up to 5 mm year⁻¹) and in the Levantine basin (up to 8 mm year⁻¹) where Pelops and Mersa-Matruh gyres occur. Smaller positive peaks were observed along the Lybian coasts, linked to the Syrte gyre. In agreement with the relevant studies acknowledged in this section, negative trends were observed in the Ionian Sea, as a consequence of an important change in the circulation observed in this basin since the beginning of the 1990s [105], and south-east of Crete associated with the Ierapetra gyre.

Considering the steric height at full depth, sea level trends showed the contribution of the halosteric component in the Mediterranean Sea, which has high negative trends and that, as already mentioned considering multi-resolution reanalyses, can lead to unrealistic sea-level trend estimates. Results from the global reanalysis GLORYS12V1 performed at 1/12° degree of horizontal resolution [85] (see also Table 1) are shown in Figure S5 in order to provide a comparison for the merged regional reanalysis trends. Such maps indicate very large trends in the Baltic Sea, dominated by the halosteric contribution (not shown). On the contrary, sea level rise in the European north-western shelf area is attenuated compared to the regional reanalyses, and the Mediterranean Sea sea level fall, due again to the halosteric sea level, is also mitigated. The comparison confirms that diversity from reanalyses is mostly due to the freshwater budget representation. Considering as a reference the recent sea level trend estimates in the European Seas (1993–2016) obtained from satellite altimetry [2], results from regional reanalyses appear more realistic in the Baltic Sea and less in the Mediterranean Sea in comparison with the global product.

Looking at the mean sea level trend, the results showed that during the period 1993–2016, the mean sea level in the European Seas rose, with rates that ranged between 1.6 and 2.6 mm year⁻¹ (±0.3 mm year⁻¹), according to the depth ranges considered to integrate the steric contributions (Table 4). In the upper ocean (up to 700 m), thermosteric sea level trends (positive) showed larger values than those due to the halosteric component (negative), while the opposite was observed considering steric contributions at full depth. The global reanalysis GLORYS12V1 (Table 5) provides comparable estimates for the layer 0–300 m and 0–700 m, within the respective error bars. Below these layers, the halosteric contributions significantly differ between the merged and the global products, making the total sea level trend estimates larger in GLORYS12V1 by more than 1 mm year⁻¹ in the 0–2000 m layer and the total column.

In the last part of this section, we focus on the European seas’ mean sea level variability and trend, considering steric contributions in the upper ocean (0–700 m). The top panel in Figure 11 shows the temporal evolution of total (red lines) and steric (black lines) sea level obtained from the merged reanalysis dataset during the period 1993–2016. In order to investigate the non-linear sea level variability and long-term trends of residual [101,106,107] in the European Seas, empirical mode decompositions (EMDs) [108,109,110] were used. EMD applied to sea level data decompose the signals in oscillatory modes (intrinsic mode functions, IMFs), representing different oceanic processes from the highest frequency to the lowest frequency oscillating mode. The remaining non-oscillating mode is the long-term trend of sea level residual [106]. Considering 24 years of monthly data, EMD gave seven IMFs both for total and steric sea level signals. Total and steric sea level signals were significantly correlated (0.9) and had a similar temporal evolution and amplitude when the combination of IMFs that explains the low-frequency modes of oscillation (~10 years) was considered, underlining the contribution of the steric component to the sea level temporal evolution at these periodicities (not shown). EMD trends (Figure 11, bottom panel) were also highly correlated (0.8), and both showed a non-linear trend with an increased slope in recent years. Interestingly, sea level residual trends were significantly correlated with those obtained considering the most prominent climate patterns in the North Atlantic (Figure 11, dashed and dotted lines). The Atlantic Multidecadal Oscillation (AMO; e.g., Enfield et al. [111]), which has its principal expression in the SST variability in the North Atlantic, showed a high positive correlation (0.9) with both sea level and steric height residuals. On the other hand, the North Atlantic Oscillation (NAO; e.g., Hurrell and Deser [112]), which is mainly an atmospheric variation mode, showed a negative correlation (−0.6) with sea level residual, due to opposite temporal evolution patterns over the period 1993–2000. These results were in line with those observed by other authors in the Southern European Seas, starting from in-situ [113] and remote-sensing [114] sea-level estimates, and underlined the role of the North Atlantic climate as one of the main drivers of the long-term trends in the European Seas. Sea level residual trends from GLORYS12V1 (Figure S6) match closely and are significantly correlated with the ones from the merged reanalysis product (with a correlation higher than 0.9), indicating the consistency between the two products in capturing the long-term sea level variability.

5. Summary and Discussion

In this work, we reviewed the use of ocean reanalyses for sea level variability investigations and assessed the trends from state-of-the-science reanalyses, including two sets of ensemble ocean reanalyses in which one ensemble reanalysis covers the entire 20th century (centennial) and the other covers the altimetry era only (from 1993 onwards, contemporary). The results show that ocean reanalyses are a powerful tool for steric sea level investigations, as they provide performances comparable to more traditional methods (e.g., objective analyses, sea level reconstructions) and allow for process-oriented studies because of their consistent multi-variate ocean state reconstruction.

Sea level trends are a fundamental metric to monitor the global warming and its potential menace on populations living near coastal areas. GMSL variability is in general dominated by thermosteric and barystatic components, while locally, sea level variability may also be affected by the halosteric component [6,65,115,116,117] and the correct partitioning of the causes of rise is crucial to understand the impact of global climate change on the maritime environment.

In spite of the large uncertainties in the early period (before the 1960s) because of poor observational sampling, historical reanalyses can still be used with a multi-system ensemble approach to shed light on the recent sea level rise in comparison with early century rise rates. In particular, a 16-member ensemble was used to assess the century-long trends, which are taken in turn from four different reanalysis systems. The ensemble mean validation against tide gauge data suggests that its inter-annual variability is found significantly correlated with the independent sea level data over most of the stations with long records. The centennial reanalyses unambiguously identify the latest two decades as the ones affected by the largest steric sea level increase. Unprecedented total column steric sea level rise was found since about 1995 (using 20-year running trends), when the rise starts to exceed 1 mm year⁻¹ and continues to increase afterwards. The centennial linear trend (1900–2010) was estimated as equal to 0.47 ± 0.04 mm year⁻¹, with a clear acceleration since the beginning of the 1990s. These trend estimates were found in close agreement with the estimates provided by the IPCC AR5 report [24] and compiled either from CMIP5 multi-model ensemble (period 1900–1990) or observations (1993–2010), motivating further use of ensemble historical reanalyses to complement coupled model simulations and statistical mapping of observations. Regional trends indicate that the steric expansion is maximum at low- and mid-latitudes, with peaks up to 3 mm year⁻¹ during the entire 20th century in the North Atlantic (south off the Gulf Stream separation) and western Pacific Oceans, and the Southern Ocean north of the ACC front.

The dynamical consistency of ocean reanalyses may be further exploited in centennial sea level reconstructions, where reanalyses provide the training data through which long records of tide gauge data are used to reconstruct the spatial patterns of sea level variability [33]. With the systematic production of multi-model ocean reanalyses, this application of reanalysis is expected to further develop and replace early attempts of using climate model simulations to draw error statistics for climate reconstructions [118]. For instance, Zanna et al. [119] successfully used an ocean reanalysis to estimate time-mean Green’s functions that were in turn exploited to reconstruct historical ocean heat content and transport changes from sea surface temperature anomaly data. Atmospheric reanalyses starting around the middle of the 19th century will allow even longer climate records for centennial investigations.

Contemporary sea level trends were obtained through analyzing an ensemble of five reanalyses. Preliminary validation against independent steric time series obtained from altimetry and gravimetry data indicates that the ensemble of reanalyses captures the recent steric sea level variability and outperforms previous vintages of reanalyses. To some extent, total sea level trends are underestimated in the comparison versus altimetry data, due to the fact that the barystatic component in reanalyses is not reproduced correctly. This is linked to the suboptimal application of freshwater input data, coming from atmospheric reanalyses (precipitation) and climatological datasets (river and land-ice discharge). More efforts are needed in the future to allow the barystatic sea level from reanalyses to show realistic variability, e.g., through sophisticated input bias correction procedures and use of time-varying freshwater datasets. The recent beginning of the production of coupled (Earth system) reanalyses [75] may in principle alleviate the misrepresentation of the barystatic component if the freshwater input is taken by the land model component, which has not yet been tested so far. However, having a closed and realistic freshwater budget in Earth system models is still challenging and shall be tackled with sophisticated methods [120].

The steric sea level trend and variability below 2000 m of depth is also less consistent between the products, as a consequence of the very limited sampling of in-situ observations therein and in agreement with previous studies [40]. Forthcoming programs, such as the deep Argo network [121], are expected to fill the gap and provide precise estimates of the steric sea level rise contribution from deep and abyssal waters.

The production of eddy resolving global reanalyses has just begun (see e.g., the new GLORYS12V1 reanalysis from Mercator Ocean [122]) and the increase of horizontal and vertical resolution obviously allows for more processes to be represented in reanalyses. However, data assimilation methods still need to be optimized for such configurations [48]. However, at the regional scale where regional reanalyses in the European seas already reach the eddy-permitting resolution, sea level trend estimates are consistent with altimetry data. In particular, considering steric contributions in the upper ocean (0–700 m) the sea level trend shows a marked spatial variability, associated with the most prominent feature of the ocean circulation (e.g., IBI and NWS regions, Mediterranean Sea). EMD analysis performed to investigate non-linear sea-level variability showed that total and steric sea-level are significantly correlated considering long-term residual trends, which also show the influence of the NAO and AMO on the long-term sea level trend in the European Seas, and in agreement with the GLORYS12 eddy-resolving global reanalysis. Considering the steric height at full depth, high negative values are observed in the Mediterranean Sea due to halosteric contributions, which lead to large departures from sea level trends obtained from observations (e.g., altimetry) and global ocean synthesis (i.e., GLORYS12V1). At such a scale, the halosteric sea level variability becomes a locally significant contributor of sea level rise (e.g., in the Baltic and Mediterranean Seas), calling again for efforts in improving the freshwater budget variability in reanalyses.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/10/1987/s1, Figure S1: SSH inter-annual variations (cm) from multi-product ensemble mean during 1993–2015. SSH inter-annual variations (cm) from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2 and SODA 3.4.2, respectively, during 1993–2015; Figure S2: Inter-annual variations (cm) of steric height above 300 m depth from multi-product ensemble mean during 1993–2015. Inter-annual variations (cm) of steric height above 300m depth from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2 and SODA 3.4.2, respectively, during 1993–2015; Figure S3: Same as Figure S2, but for steric height below 300 m of depth (300-bottom); Figure S4: Correlation of steric height above 300 m depth from multi-product ensemble mean with ENSO during 1993–2015. Correlation of steric height above 300 m depth from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2 and SODA 3.4.2, respectively, with ENSO during 1993–2015. Blank areas show the correlations that are not significant at 95% confidence level; Figure S5: As Figure 10, but using data from the GLORYS12V1 global reanalysis; Figure S6: As Figure 11, but using data from the GLORYS12V1 global reanalysis.

Author Contributions

A.S. has prepared Section 1, Section 2 and Section 5 and contributed to Section 3 and Section 4; C.Y. has prepared Section 3; X.F. and A.B. have prepared Section 4. All authors have revised the manuscript.

Funding

X.F. was partly funded by the Belt and Road Special Foundation of the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering in China. A.B. was supported by the Initiative and Networking Fund of the Helmholtz Association through the project “Advanced Earth System Modelling Capacity (ESM)”.

Acknowledgments

The authors want to acknowledge all the producers of the dataset used in this work, in particular E. de Boisseson and M. Balmaseda (ECMWF), B. Giese (Texas A&M University), the ECCO consortium (https://ecco.jpl.nasa.gov/), the SODA reanalyses (http://www.soda.umd.edu), and all global and regional reanalyses included in the Copernicus Marine Environment Monitoring Service (CMEMS). The Permanent Service for Mean Sea Level (PSMSL, https://www.psmsl.org) and the ESA Sea Level CCI (http://www.esa-sealevel-cci.org) provided tide gauge and altimetry data, respectively. GRACE ocean data were processed by Don P. Chambers, supported by the NASA MEaSUREs Program, and are available at http://grace.jpl.nasa.gov.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Validation of the centennial reanalysis steric sea level against long-record tide gauge stations: Green bullets are located where the stations exhibit statistically significant correlations (at 95% confidence level); red bullets when the stations show non-significant correlation but at least one ensemble member shows significant correlation; black bullet where neither the ensemble mean nor any individual reanalysis show significant correlation. The value of the correlation is indicated by the size of the green bullets, as in the legend. Zoom over north European coasts is shown in the bottom panel.

Figure 2. Global mean thermosteric sea level for the four vertical layers 0–300 m, 0–700 m, 0–2000 m, and 0–bottom from the 16-member ensemble mean centennial reanalysis. Total steric sea level from the ensemble mean is shown by the thin dotted lines. Dark and light orange shaded areas correspond to one and three ensemble standard deviations, respectively.

Figure 3. Twenty-year running steric sea level trends for the four vertical layers 0–300 m, 0–700 m, 0–2000 m, and 0–bottom. Black lines: Thermosteric global sea level trends; red lines: Thermosteric sea level trends in the 60° S to 60° N region; dashed lines: Global steric sea level trends. At the bottom of each panel, green bars correspond to the years when the trend is significant, based on the Mann–Kendall trend test [94] at the 95% confidence level. Orange, red, or purple years correspond to the years where the ensemble mean trends exceed one, two, or three times the ensemble spread of the 16-member trends, respectively (corresponding approximately to 70%, 95%, and 99% confidence level. Dashed and dotted horizontal lines correspond to the time-mean thermosteric and total steric trends, respectively.

Figure 4. Maps of 1900–2010 thermosteric (top four panels) and steric sea level trends (bottom four panels) in mm year−1 of the ensemble mean for the four vertical layers 0–300 m, 0–700 m, 0–2000 m, and 0–bottom. Gray crosses correspond to statistically significant trends, based on the comparison with the ensemble spread, at the 95% confidence level.

Figure 4. Maps of 1900–2010 thermosteric (top four panels) and steric sea level trends (bottom four panels) in mm year−1 of the ensemble mean for the four vertical layers 0–300 m, 0–700 m, 0–2000 m, and 0–bottom. Gray crosses correspond to statistically significant trends, based on the comparison with the ensemble spread, at the 95% confidence level.

Figure 5. Global mean steric sea level (2003–2015) from different ensemble reanalysis products as described in the text and Table 1 and the satellite steric sea level estimates (OBS) derived from altimetry and gravimetry data. Data are monthly means with the seasonal cycle removed.

Figure 6. (a) Annual global SSH and steric height at different depths from multi-product ensemble mean during 1993–2015. (b–f) Annual global SSH and steric height at different depths from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2, respectively, during 1993–2015. The global MSL with GIA corrected from AVISO is also presented (“obs” dashed line), which provides a validating dataset for the sea surface height (“ssh”) timeseries.

Figure 7. (a) Ratio of SSH trend over steric height trend at full depth from multi-product ensemble mean during 1993–2015. (b–f) Ratio of SSH trend over steric height trend at full depth from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2, respectively, during 1993–2015.

Figure 8. (a) Inter-annual variability of global SSH and steric height at different depths from multi-product ensemble mean during 1993 to 2015. (b–f) Inter-annual variability of global SSH and steric height at different depths from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2, respectively, during 1993 to 2015. The global MSL with GIA corrected from AVISO is also presented.

Figure 9. (a) Percentage of SSH inter-annual variability explained by steric height at full depth from the multi-product ensemble mean during 1993 to 2015. (b–f) Percentage of SSH inter-annual variability explained by steric height at full depth from ECCO-KFS, ECCOv4r3, GLORYS4, SODA3.12.2, and SODA 3.4.2, respectively, during 1993 to 2015.

Figure 10. Sea level (mass and steric contributions) trend from merged high-resolution regional reanalyses in European seas. In the panels, sea-level rates are shown as obtained considering the sea level as the resultant of the mass component and steric height integrated from the surface to (a) 300 m, (b) 700 m, (c) 0–2000 m, and (d) 0–bottom.

Figure 11. Empirical mode decomposition analysis (EMD) performed considering sea level (red lines) and steric height (0–700m, black lines). Top panel shows the sea level and steric height temporal evolution (cm). The bottom panel shows the EMD residual signal (centered and scaled) considering sea level, steric height, NAO (North Atlantic Oscillation, green dotted line), and AMO (Atlantic Multidecadal Oscillation, blue dashed line) climate indexes. X-axis shows the years considered (1993–2016).

Table 1

List of all reanalyses used in this study, along with ensemble size, period, ocean model and resolution, data assimilation method and observations assimilated (either directly or through relaxation), atmospheric forcing, and reference. The table is divided in three parts: the first part lists the centennial reanalyses (Section 3), the second part the global ocean contemporary reanalyses (Section 4.1, Section 4.2 and Section 4.3), and the third part the European regional reanalyses (Section 4.4).

Reanalysis	Native Ensemble Size	Period	Ocean Model (Resolution)	Data Assimilation Method (Observations)	Atmospheric Forcing	Reference
CERA-20C	10	1900–2010	NEMO (1° × 1° L42)	3DVAR (SST, in-situ)	Coupled	Laloyaux et al. [75]
CHOR	4	1900–2010	NEMO (0.5° × 0.5° L75)	3DVAR (SST, in-situ)	20CRv2c; ERA-20C	Yang et al. [55]; Yang et al. [56]
ORA-20C	10	1900–2010	NEMO (1° × 1° L42)	3DVAR (SST, in-situ)	ERA-20C	de Boisseson et al. [76]
SODA	8	1870–2013	POP (0.4° × 0.25° L40)	OI (SST)	20CRv2c	Giese et al. [77]
GREP	4	1993–2017	NEMO(0.25° × 0.25° L75)	3DVAR; SEEK(SLA, SST, in-situ, SIC)	ERA-Interim	Storto et al. [41]
GLORYS4	1	1993–2017	NEMO(0.25° × 0.25° L75)	SEEK(SLA, SST, in-situ, SIC)	ERA-Interim	Garric et al. [78])
SODA3.4.2	1	1980–2017	GFDL CM2.5/MOM5(0.25° × 0.25° L50)	EnKF(SST, in-situ, SIC)	ERA-Interim	Carton et al. [79]
SODA3.12.2	1	1980–2016	GFDL CM2.5/MOM5(0.25° × 0.25° L50)	EnKF(SST, in-situ, SIC)	JRA-55DO	Carton et al. [79]
ECCO-KFS	1	1993–2018	MITgcm(1° × 1° L46)	KFS (SLA, SST, in-situ)	NCEP	Fukumori [80]
ECCOv4r3	1	1992–2015	MITgcm(1° × 1° L50)	4DVAR (SLA, SST, in-situ, SIC, SSS, MDT)	ERA-Interim	Forget et al. [61]
BAL	1	1993–2017	NEMO(4 km L56)	LSEIK(SST, in-situ)	HIRLAM-UERRA	CMEMS-BAL-QUID [81]
NWS	1	1992–2018	NEMO (7 km L24)	3DVAR(SST, in-situ)	Era-Interim	CMEMS-NWS-QUID [82]
IBI	1	1992–2018	NEMO (5–6 km L56)	SEEK(SLA, SST, in-situ)	Era-Interim	Sotillo et al. [83]
MED	1	1987–2017	NEMO (6–7 km L72)	3DVAR(SLA, SST, in-situ)	Era-Interim	CMEMS-MED-QUID [84]
GLORYS12V1	1	1993–2017	NEMO(1/12° × 1/12° L50)	SEEK(SLA, SST, in-situ, SIC)	Era-Interim	Gasparin et al. [85]

Abbreviations: Coupled Earth system ReAnalysis of the 20th century (CERA-20C); Ocean ReAnalysis 20th Century (ORA-20C); Simple Ocean Data Assimilation (SODA); Global Reanalysis Ensemble Product (GREP); Global Ocean ReanalYsis System (GLORYS); Kalman Filter/Smoothr (KFS); Baltic Sea (BAL); European North-Western Shelf area (NWS); Ireland Biscay Iberia region (IBI); Mediterranean Sea (MED).

Table 2

Correlation between SAT (satellite steric sea level as altimetry minus gravimetry) and reanalysis sea level, and their linear trends for the different ensemble mean steric sea level datasets explained in the text.

Period	2003–2010		2003–2014
Dataset	Correlation	Trend (mm year−1)	Correlation	Trend (mm year−1)
SAT		1.05 ± 0.18		1.08 ± 0.10
REA	0.86	1.57 ± 0.16	0.87	1.08 ± 0.09
ORAIP	0.63	0.33 ± 0.09
REA subsampled	0.84	1.64 ± 0.16
ORAIP subsampled	0.44	0.15 ± 0.12

Table 3

Sea level trend (mm year⁻¹) in the European Seas (multiple resolution reanalysis). Rows: sea level trends in the Baltic Sea (BAL), North Western Shelf (NWS), Iberian-Biscay-Ireland (IBI) region, and Mediterranean Sea (MED). Columns 1–3: Thermosteric, halosteric, and steric (full depth) sea-level trends. Column 4: Sea level trend estimates obtained from the resultant of mass and steric (full depth) components. Column 5: As in column 4 but considering the resultant of mass and thermosteric components.

Region	Thermosteric	Halosteric	Steric	Sea Level	Sea Level (Thermosteric)
BAL	0.14 ± 0.10	−0.40 ± 0.14	−0.26 ± 0.10	2.66 ± 1.15	2.98 ± 1.15
NWS	−0.37 ± 0.14	0.53 ± 0.10	0.16 ± 0.10	2.12 ± 0.38	1.60 ± 0.40
IBI	1.80 ± 0.23	−0.82 ± 0.20	0.97 ± 0.10	3.70 ± 0.15	4.57 ± 0.24
MED	2.56 ± 0.10	−6.35 ± 0.10	−3.6 ± 0.37	−3.50 ± 0.30	2.85 ± 0.26

Table 4

Sea level trends (mm year⁻¹) from merged reanalysis. Rows: depth ranges considered to compute steric contributions. Columns 1–3: Thermosteric, halosteric and steric sea-level trends. Column4: Sea level trend obtained as a result of the mass component and steric contributions (integrated at different depths). Column 5: As in column 4 but considering the resultant of the mass and thermosteric components.

Layer	Thermosteric	Halosteric	Steric	Sea Level	Sea Level (Thermosteric)
0–300 m	0.73 ± 0.10	−0.61 ± 0.10	0.13 ± 0.06	2.22 ± 0.30	2.82 ± 0.16
0–700 m	1.58 ± 0.20	−1.03 ± 0.10	0.54 ± 0.15	2.63 ± 0.32	3.16 ± 0.17
0–2000 m	1.77 ± 0.22	−1.87 ± 0.21	−0.10 ± 0.05	1.98 ± 0.30	3.86 ± 0.20
0–bottom	1.82 ± 0.21	−2.31 ± 0.20	−0.50 ± 0.10	1.60 ± 0.35	3.91 ± 0.18

Table 5

As in Table 4, but for data from the GLORYS12V1 global reanalysis.

Layer	Thermosteric	Halosteric	Steric	Sea Level	Sea Level (Thermosteric)
0–300 m	0.95 ± 0.24	−0.47 ± 0.04	0.48 ± 0.24	2.73 ± 0.55	3.20 ± 0.55
0–700 m	1.63 ± 0.25	−0.87 ± 0.05	0.76 ± 0.25	3.01 ± 0.56	3.88 ± 0.56
0–2000 m	1.82 ± 0.26	−1.04 ± 0.07	0.78 ± 0.25	3.03 ± 0.57	4.07 ± 0.57
0–bottom	2.01 ± 0.26	−1.25 ± 0.08	0.76 ± 0.25	3.01 ± 0.57	4.26 ± 0.57