Multidecadal Sea Level Variability in the Baltic Sea and Its Impact on Acceleration Estimations

Multidecadal sea level variation in the Baltic Sea is investigated from 1900 to 2020 deploying satellite and in situ datasets. As a part of this investigation, nearly 30 years of satellite altimetry data are used to compare with tide gauge data in terms of linear trend. This, in turn, leads to validation of the regional uplift model developed for the Fennoscandia. The role of North Atlantic Oscillation (NAO) in multidecadal variations of the Baltic Sea is also analyzed. Although NAO impacts the Baltic Sea level on seasonal to decadal time scales according to previous studies, it is not a pronounced factor in the multidecadal variations. The acceleration in the sea level rise of the basin is reported as statistically insignificant in recent studies or even decelerating in an investigation of the early 1990s. It is shown that the reason for these results relates to the global warming hiatus in the 1950s−1970s, which can be seen in all eight tide gauges used for this study. To account for the slowdown period, the acceleration in the basin is investigated by fitting linear trends to time spans of six to seven decades, which include the hiatus. These results imply that the sea level rise is accelerated in the Baltic Sea during the period 1900–2020.


INTRODUCTION
Sea level rise as a result of climate change is one of the significant environmental threats which requires a deep understanding not only globally but also on regional scales. This will lead to having more realistic future scenarios and, thus, helping decision-makers in planning a safer future for coastal regions (Dangendorf et al., 2014;Frederikse et al., 2016;Prandi et al., 2021). Although, because of the crustal land uplift, a substantial part of the Baltic Sea coast is safe from future sea level rise, understanding the mechanism of sea level variability in the region is important since several reliable sea level measurements stretching back to the 18th century are located in this region. To be able to use these measurements in global-scale analyses, the geophysical and climate-related processes need to be understood to separate the global effects from regional drivers. This is the reason that there have been numerous studies addressing these variations in the Baltic Sea (Hünicke and Zorita, 2008;Hünicke et al., 2015).
The Baltic Sea is a semi-enclosed basin, under the constant impact of external and internal forcings . Water discharge of the surrounding catchments, precipitation, and water inflow/outflow from the North Sea are the external variability factors. On the other hand, air pressure, wind, and density variations account for the internal forcing. Freshwater input acts as internal and external forcings as it changes both the volume and salinity of the water in the basin. The classification of contributors to external and internal forcings are according to Samuelsson and Stigebrandt (1996) that shows the external forces are responsible for 50-80% of sea level variations in the Baltic Sea for periods longer than 1 month, from 1977 to 1987. In the other two complementary studies Meier and Kauker, 2003), half of decadal sea salinity variation in the basin are attributed to freshwater inflow in the twentieth century. Additionally, sea level pressure (SLP) is found as a substantial factor on this time scale. Most recently, Gräwe et al. (2019) conducted thorough analyses on decomposing the role of different drivers in Baltic Sea level rise for the period 1950-2015 using different reanalyzes and in situ measurements. They found that 75% of the sea level rise in the basin is due to the water inflow from the North Sea.
High energy North Atlantic climate has a substantial influence on sea level variability of the Baltic Sea and consequently, it has been subject to different studies. In this respect, North Atlantic Oscillation (NAO) is a major driving force for the Baltic Sea on seasonal to decadal time scales (Andersson, 2002;Yan et al., 2004;Hünicke and Zorita, 2006;Hünicke et al., 2015). Gräwe et al. (2019) reasoned that NAO impacts the Baltic Sea by weakening or strengthening westerlies. This, in turn, leads to variations in inflow/outflow rate from/to the North Sea. Even though this is in line with the findings of Dangendorf et al. (2014) for the North Sea, they expressed a paucity of a thorough investigation on this topic. Other studies tie NAO to sea level variability in the Baltic Sea by thermal expansion due to temperature increase (Hünicke and Zorita, 2006;Dangendorf et al., 2012) or freshwater input by rising the river discharge to the basin and also direct contribution because of higher precipitation Zorita, 2006, 2008;Lehmann et al., 2011). Karabil et al. (2017) investigated the decadal variability of the Baltic Sea by analyzing a wide range of climate data (i.e., SLP, air temperature, precipitation, and climate indices) and concluded that the response of sea level to these factors is not homogeneous in the basin where they weaken toward the southern coasts. Pajak and Blaszczak-Bak (2019) reflected the asymmetry in terms of linear change using altimetry data between 1993 and 2017. Similar asymmetric variabilities were found by Barbosa (2008) on a centennial time scale. Karabil et al. (2018) introduced a new climate index to the basin, which represents the atmospheric circulation better than NAO. Passaro et al. (2015) found that the main driver for the annual cycle of sea level variations around the Baltic coasts is wind stress. They also found that steric parameters influence the inner basin variabilities. Passaro et al. (2021) investigated contemporary sea level variations, local drivers, and spatial patterns in the Baltic Sea. They optimized the satellite altimetry measurements for the basin to minimize the impact of rugged coastlines and sea ice cover in wintertime. They also found that wind is a major driver for the gradient of the sea level rise in the altimetry era. Pajak and Kowalczyk (2019) used tide gauges along the Polish coast and satellite altimetry data to analyze the seasonal variations and reported a substantial agreement between these two datasets in this respect. In another study on the seasonal scale, Hünicke and Zorita (2008) investigated the annual signal changes in the basin and reported centennial variations in this signal. They speculated that precipitation drives these low-frequency variations. Xu et al. (2015) reported a fall in the correlation between the wintertime NAO and sea level variation in the last decade at the time of the study using altimetry data. Hünicke and Zorita (2008) showed that precipitation and temperature variations contribute to sea level changes in the Baltic Sea seasonally, but they are not as substantial as SLP.
The Baltic Sea basin is heterogeneously subject to postglacial rebound resulting from the last glacial age. The range of crustal uplift ranges from −1 to +12 mm/year in the basin although this range changes slightly in different models. This crustal movement requires accurate modeling to enable the use of sea levels recorded in the basin in a geocentric system. In addition to global glacial isostatic adjustment (GIA) models (e.g., Steffen and Wu, 2011;Roy and Peltier, 2018), the Nordic Commission of Geodesy (NKG) has released a semi-empirical land uplift model in which a GIA model is combined with geodetic observations (Vestøl et al., 2019). Madsen et al. (2019) used this model to reconstruct sea level for the Baltic Sea, and they found an average of 1.3 ± 0.3 mm/year sea level rise in the basin for the twentienth century. Richter et al. (2012) estimated a similar trend for the same period but used GNSS observations to correct the land uplift signal. GNSS-derived land uplift might not be enough to completely explain the GIA effect since it does not take the gravitational attraction of the accumulated mass on sea level into account (Tamisiea and Mitrovica, 2011). In the Baltic Sea, however, the ratio of sea level rise to uplift (because of geoid change) is negligible, 1-20. Thus, this ratio does not affect the agreement between the two aforementioned studies. Gräwe et al. There were also attempts to determine a possible acceleration in the sea level rise of the basin but their outcome are irreconcilable. Woodworth (1990) estimated decelerations for the majority of tide gauges located on the Baltic coasts using varying lengths of tide gauges from 1870 to 1981 [See Table  1(a) of the study for further information]. Donner et al. (2012) reported that there are decadal-scale accelerations and decelerations, using tide gauge data of varying lengths (1849-2006) but did not draw any solid conclusion. Stramska and Chudziak (2013) estimated the acceleration in the altimetry era and concluded that there is no statistically significant acceleration in the basin. Spada et al. (2014), on the other hand, reported a secular acceleration in the Baltic Sea and related it to isostatic adjustment. Hünicke and Zorita (2016) carried out a thorough investigation of the acceleration determination in the basin using multiple tide gauges and methods from 1900-2012. They found that although all-station-average-record yielded "almost statistically significant" acceleration, the individual estimations were statistically insignificant but positive for most of the tide gauges.
All aforementioned studies provide valuable information about the Baltic Sea level variability and the role of contributors on these variations in different periods. The analyses of all these studies were conducted on seasonal to decadal time scales.
In contrast, the main focus of this study is the multidecadal variability of the Baltic Sea using tide gauge data and altimetry data. Since the length of the datasets may not lead to a reliable outcome when spectral methods are used, most of the analyses are carried out by comparing linear trends fitted to the multidecadal periods. As such, we first use three decades of altimetry data to evaluate the consistency between VLMcorrected tide gauge records and VLM-free satellite data by comparing sea level trends. The role of NAO on sea level rise and variation on the multidecadal time scale is also investigated. Furthermore, we analyze the multidecadal sea level variations in the tide gauge data and show that the global warming hiatus can be seen in the basin between the 1950s and 1970s. The link between this hiatus and sea level acceleration estimation is discussed in the last section of the study.

DATA
Five datasets are used in this study, and they are briefly explained in this section.

Altimetry
Altimetry data are acquired from Radar Altimetry Database System (RADS; Scharroo et al., 2013) with standard geophysical and range corrections applied on them, except for wet tropospheric correction. The dataset is comprised of collinear measurements of Topex/Poseidon, Jason-1, Jason-2, and Jason-3 missions spanning the period from 1993 to 2020. The collinear data are clustered according to the equator pass time of the satellite for each pass and, thus, provide ready-to-use sea level anomaly (SLA) time series. The onboard radiometer correction for wet tropospheric correction is replaced by the model correction provided by European Centre for Mediumrange Weather Forecasts (Xu et al., 2015). The poor performance of the radiometer is due to the large footprint of the instrument and hence land contamination (Brown, 2010). Although this replacement increases the number of altimetry measurements, it does not improve the number of time series which are long enough to use in the analyses. Time series with at least 75% of valid measurements are involved in the analyses; out of 1,286 time series in the basin, only 311 fulfill this condition (Figure 1). In addition to being relatively shallow, the basin has several islands and narrow subbasins which increase the land contamination on the backscattered waveforms of altimetry (Deng and Featherstone, 2006).
The main reason for the absence of any valid time series above 63 • in latitude is the presence of sea ice in wintertime (Omstedt et al., 2004). This is flagged in the surface class of each measurement by RADS and these measurements are discarded from the dataset. Soomere (2016, 2017) reported that in RADS datasets the surface flag of an altimetry measurement is set to sea ice when the ice concentration is 50% and showed that this threshold is suitable for the Baltic Sea. The same studies also recommended discarding measurements with proximities of closer than 0.2 • to the nearest coast. The dataset used in this study follows this condition as measurements closer than 20 km to the nearest coast are excluded. Radar Altimetry Database System uses Dynamic Atmospheric Correction (DAC) as an alternative to Inverted Barometer (IB) correction because it represents the atmospheric impact better than IB, which is a statistic model (Scharroo et al., 2016). To compare the altimetry measurements to the tide gauge records, this correction is only added back to the SLA data when the role of IB on trend estimation is analyzed in Altimetry Era section. The altimetry data in this study are also corrected for the geoid change of the GIA effect (further information about the GIA effect on geoid is provided in Land Uplift Model Section).

Tide Gauge
Monthly tide gauge records are provided by Permanent Service for Mean Sea Level (PSMSL; Holgate et al., 2012). One of the objectives of this study is detecting low-frequency variations in the Baltic Sea, which are generally of low amplitudes. Detection of these signals requires a certain level of accuracy in the observed data. For this reason, only datasets of the twentieth century are included in the dataset. Moreover, tide gauges with long gaps (longer than 2 years) are excluded. These criteria narrow down the number of tide gauges to eight (Figure 1 and Supplementary Table 1). Tidal modeling is not applied on the tide gauges records as these variations are almost negligible due to (1) the effect of the Danish Straits (Gräwe et al., 2019) and (2) the basin size, which is not considerably big to develop substantial tidal effects.

Land Uplift Model
To transform the relative sea level measurements to a geocentric system, an uplift model is required for the Baltic Sea because of the substantial rebound due to the last deglaciation. Following previous studies (Gräwe et al., 2019;Madsen et al., 2019), the uplift model developed by NKG is adopted to apply the crustal uplift correction. This model is released in two versions, NKG2016LU_abs and NKG2016LU_lev (Vestøl et al., 2019). For the latter version, the geoid variations due to GIA is not considered. However, the difference between the two models does not exceed 10% which is slightly higher than the previous studies estimated the geoid to land variation rate to be 5% (Tamisiea and Mitrovica, 2011). Tide gauge measurements are not deployed in developing this model (c.f., Hill et al., 2010) so that it can be used to correct the tide gauge measurements as an independent model. The model geographically extends from 0 • to 50 • in longitude and 49 • -75 • in latitude. The maximum uplift in the region is 9.6 mm/year, close to the Ratan tide gauge station. GIA renders uplift on all the tide gauges used in this study except for Warnemünde, which is located on the bulge and experiences a slight depression. The uncertainty of the estimated uplift varies along with the geographic extent of the model depending on the availability of geodetic measurements. However, along the Baltic coasts, it ranges from 0.2 to 0.3 mm/year (Vestøl et al., 2019). This study uses the NKG2016LU_abs version to correct the vertical land motion at the tide gauge locations.

Sea Level Pressure
Sea level pressure is one of the main drivers of sea level in the Baltic Sea and it is governed by NAO, particularly during wintertime. For SLP, the results of twentieth-century Reanalysis Project Version 3 (Compo et al., 2011;Giese et al., 2016;Slivinski et al., 2019), provided by the National Oceanic and Atmospheric Administration, is used. The dataset covers the whole twentieth century but extends only until 2016. Therefore, the dataset is augmented with reanalysis data of the National Centers for Environmental Prediction and the National Center for Atmospheric Research (Kalnay et al., 1996) to cover the gap between 2016 and 2020. To estimate the inverted barometric effect of SLP, the following equation is used (Ponte, 2006;Piecuch et al., 2016): where p a is the SLP in tide gauge location, p a is mean global SLP over oceans, ρ = 1029 kg/m 3 reference ocean water density, and g = 9.81 m/s 2 is the gravitational acceleration. Mean global SLP is estimated by masking the land in the models. It should be noted that this equation is based on the fact that the sea surface follows the variations in SLP. This is not completely valid for the Baltic Sea due to the impact of the Danish Straits which confines the water exchange with the North Sea. Hence, the SLP and the SLA are not synchronized. On the other hand, due to the gradient in the water density, the sea level in the north of the basin stands 30-35 cm higher than that of the Danish Straits. Nevertheless, this does not affect the SLA variations (Ekman and Mäkinen, 1996).

North Atlantic Oscillation
North Atlantic Oscillation is provided by the Climate Data Guide initiative, which is a project led by the University Corporation for Atmospheric Research and National Corporation for Atmospheric Research (Schneider et al., 2013). NAO is an atmospheric pattern that influences the sea level variation in the basin by intensifying the westerlies when it is positive and, thus, causes higher water inflow from the North Sea (Hurrell and Deser, 2009). NAO is estimated in two forms, station-based and principal component-based, and covers the period of 1899-2020. The latter form is used in this study as it is less prone to measurement noise than the former one and it represents the spatial pattern better.

RESULTS
Sea level rise is analyzed in different multidecadal periods. It starts with the altimetry era and comparing the sea level rise estimated by satellite altimetry and tide gauge data. The impact of NAO on sea level rise on multidecadal periods is investigated. Multidecadal sea level variations since the beginning of the twentieth century and how they influence the sea level acceleration estimations are also scrutinized.

Altimetry Era
Satellite altimetry enabled the estimation of the absolute sea level rise since 1993 in nearly global coverage. This technique is more valuable in regions with substantial crustal motions such as the Baltic Sea. To this end, the sea level trend in the basin is estimated by a regression model and a proper noise model (Figure 2).
where η is dynamic sea level, a, b, and c are coefficients to be estimated, t is time, ω i is the frequency corresponding to annual and semiannual signals, and ε is the noise model. The noise model is autoregressive fractionally integrated with order one, which is also denoted as ARFIMA (1, d, 0) in which d and 0 show fractional memory parameter and order of moving average, respectively. To decide on a proper noise model, five noise models are tested and scored according to the Akaike information criterion (AIC) and Bayesian information criterion (BIC). ARFIMA (1, d, 0) was identified as the most suitable noise model for the majority of the time series. Details about the noise model analyses can be found in Supplementary Materials. Sea level trend for the same period is also estimated in the tide gauges after applying the IB and land uplift corrections. The mean sea level rise in the region estimated by altimetry and tide gauge data are 4.8 ± 0.7 and 3.8 ± 1 mm/year, respectively. The difference between the mean trends can be attributed to the heterogeneous distribution of both datasets throughout the area and the underestimation of the uplift in all or some of the stations. There is a north-south pattern in the estimated trend in which the rise in the Gulf of Bothnia is 1.5 times larger than in southern coasts (the Baltic Proper). The uncertainty map One of the geodetic applications of satellite altimetry observations is the estimation of crustal land uplift by comparing them to tide gauges (e.g., Nerem and Mitchum, 2002;Kuo et al., 2004Kuo et al., , 2008. The results of the tide gauges are compared with the altimetry wherever there are altimetry measurements inside the 2 • × 2 • box around the station. This narrows down the number of tide gauges to three. Considering the error budgets of observations and land uplift model, the agreement between the trends estimated by altimetry and tide gauges are comparatively high in two of these stations, Stockholm and Kungsholmsfort. Passaro et al. (2021) reached the same results for the trend difference between altimetry and tide gauge although a greater set of tide gauge and altimetry pairs was deployed. However, the estimated trend in Ölands Norra Udde station is not in line with the results of the other two stations with a 1.5 mm/year discrepancy. This does not necessarily mean that the uplift model is inaccurate as the sea level variations in island stations can be different from those on the mainland due to hydrodynamic effects.
To investigate the effect of SLP on the sea level rise in the altimetry era, the trend map in Figure 2 is reestimated by adding DAC to the altimetry observations (results are not shown). Although the variability of DAC and IB corrections is not exactly coherent in terms of the correlation coefficient, their contributions to the estimated trend are analogous. Removing these corrections does not change the overall outcome; however, an approximately constant fall of −0.2 mm/year can be seen in both altimetry and tide gauge results. The effect of SLP, therefore, is <10% and it slightly decreases the sea level rise in the region in the altimetry era.

NAO Related Sea Level Variation and Rise
North Atlantic Oscillation (NAO) is a climate mode that impacts the sea level variability in the North Atlantic by having imprints in buoyancy flux and surface wind stress with interannual to decadal variability (Curry and Mccartney, 2001). It has been clearly shown in previous studies (Andersson, 2002;Yan et al., 2004;Hünicke and Zorita, 2006;Hünicke et al., 2015;Xu et al., 2015;Karabil et al., 2017) that NAO is one of the governing factors in sea level variability of the Baltic Sea on seasonal to decadal periods. However, the role of NAO on multidecadal sea level variability and how it affects the sea level rise estimations on this time scale has yet to be studied.
The long-term records of tide gauges in the Baltic Sea offer an opportunity of investigating the long-term coherence between the sea level variation and NAO. Figure 3 shows the wavelet coherence (Grinsted et al., 2004) of the tide gauges and NAO. The effect of NAO can also be analyzed spatially by partitioning the basin into three subregions of north, middle, and south represented by Ratan, Stockholm, and Warnemünde, respectively. As it has been reported, NAO has a substantial impact on sea level variation of the basin not only in the altimetry era but also throughout the twentieth century in interannual and decadal scales. Additionally, this impact is stronger in the north and it fades out toward the south. The main question left, however, is the effect that NAO has on sea-level rise estimation in the Baltic Sea.
One of the most frequent methods to model internal variability is to augment the linear regression model with climate indices which govern the climate of the study region (Zhang and Church, 2012;White et al., 2014;Frankcombe et al., 2015). In the case of this study, this method can be formulated as The effect of NAO on sea level variation of the Baltic Sea is not consistent in time, and it shows strong variations on seasonal to multidecadal time scales. Karabil et al. (2018) estimated correlations between two long-term tide gauges (Warnemünde and Stockholm) and NAO separately for winter and summer in 21-year gliding segments. They showed that correlation for different seasons and stations can range from −0.4 to 0.8. Additionally, the varying coherencies between NAO and the SLA time series (Figure 3) demonstrate that the influence of NAO depends not only on the period but also on the time scale. Due to the time variability of NAO impact, Equation (3) does not eliminate the NAO-related variations. To clarify this, another analysis is carried out. For periods longer than 16 years, the results of wavelet coherence are not reliable considering the cone of influence.
To see a clear picture of the variations in decadal and multidecadal scales, a low-pass bandwidth filter is applied on both SLA and NAO time series. Continuous wavelet transform is used to filter out signals with periods shorter than 4 years. This threshold is set according to visual analyses of the spectrograms which show that signals with periods shorter than this threshold are mostly nonstationary. The filtering is applied on the SLA of each tide gauge and then they are averaged to represent the basin-wide sea level variation (Figure 3, bottom panel). The main remark is the time variable impact of NAO on sea level variability of the basin. This leads to complexity in modeling its footprint on sea level rise estimations. For instance, the NAO influence on the mean SLA of the basin is more pronounced in the 1940s while it comparatively becomes less effective in the 1950s. A similar conclusion can be extended to sea level variations in the altimetry era in the middle and southern subregions by analyzing the coherence in the Stockholm and Kungsholmsfort stations.
Additionally, the impact of NAO is not identical for different seasons of the year. This causes that the fitting process in Equation (3) depends on the smoothing applied to both the SLA and NAO time series. This has been quantified by the coefficient of determination, which is a measure to quantify the fitting performance of a functional model. It is defined as 1 − (SS res /SS tot ) where SS res represent the sum of the square of residuals and SS tot shows the variance of the data (the SLA). For instance, the coefficient of determination for basin-average sea level trend accounts for only 0.29, while when a 13 months smoothing is applied and seasonal terms are removed from the equation, the coefficient of determination rises to 0.45. Low-pass filtering of time series, in turn, may lead to manipulations in the temporal correlation of the observations in time series, hence underestimation of uncertainty. Therefore, there is a trade-off between smoothing and uncertainty estimations.
Despite all these obstacles in modeling NAO-related variations of sea level, the silver lining is the quasi-periodic nature of these variations with seasonal to decadal components. These types of variations have generally insignificant contributions to trend estimation after three or more cycles (Karimi and Deng, 2020). Additionally, according to the coherence figures, the NAO signal can be considered stationary in the altimetry era for the northern subregion (Vasa, Ratan, and Oulu stations) where the coherence is above 0.5 and in phase. Thus, a multivariate regression model can provide some insights into the role of NAO in SLA trend estimations. Figure 4 shows the difference in trend estimation using Equations (2) and (3). The maximum contribution of NAO to SLA trend is less than 0.4 mm/year. Therefore, the internal variability imposed by NAO accounts for a maximum of ∼10% of sea-level rise in the basin in the altimetry era.

Multidecadal Sea Level Variations and Acceleration Estimations
The Baltic Sea is one of the best-monitored basins where the sea level has been recorded almost regularly since the beginning of the twentieth century or even earlier. This provides the possibility of analyzing low-frequency variations and investigate their occurrences and periodicities. We assume that the land uplift model completely removes the GIA effect and that the basin is not subject to any nonlinear solid earth signal which might leave an imprint on sea level variation of the twentieth century. By analyzing the multidecadal variations, we try to reconcile the contradicting acceleration results presented in the early 1990s and most recently.
The tide gauge records are low passed filtered using the same approach in the previous section but with a wider bandwidth which removes all decadal or shorter periods (Figure 5). The sea level rise in the basin can be divided into three phases where the breakpoints are 1945 and 1978. These breakpoints are chosen according to the results of spectral analyses in Figure 3. After a substantial fluctuation between 1940 and 1950, the low-frequency signal of the SLA follows a downward trend until it reaches another low in 1978, despite some dips and peaks in the interim. A linear regression model is also fitted to the monthly records of each phase to show the general tendency of the sea level rise in these phases. Moving the breakpoints a few years changes the estimated trends but the overall outcome for the three phases does not change. The estimated trends imply a multidecadal variation in almost all tide gauges. These findings can be linked to the results of other studies on a global scale. Dangendorf et al. (2019) and Frederikse et al. (2020) reported that the rate of global sea level rise experiences a pronounced drop in the period of 1950s−1970s. Frederikse et al. (2020) reported that this is a consequence of lower mass contribution to sea level rise with respect to before and after this period. The presence of such multidecadal variations can be a substantial obstacle to estimate the sea level acceleration in the basin (Calafat and Chambers, 2013). Depending on the length of the dataset used for acceleration analyses, the outcome can be negative, positive, or insignificant when the range and magnitude of acceleration estimated for sea-level on a global scale is considered, ∼0.01-∼0.025 mm/year 2 (Church and White, 2006;Calafat and Chambers, 2013). Overlooking the physical processes with multidecadal periods can make the mathematical modeling results sensitive to the length of the dataset and the period in which it is applied on.
In the case of this study, if the dataset is confined, for instance, to the period of 1944-2020, the acceleration estimations will yield considerable figures comparing with any other period (as shown in Supplementary Table 2). These interpretations, however, are apart from the statistical significance of the estimated accelerations when a proper noise model is deployed.
The results of the previous studies on the sea level acceleration of the Baltic Sea can be interpreted by these multidecadal variations. For instance, Woodworth (1990) estimated negative accelerations for most of the tide gauges in the Baltic Sea. Considering the slowdown period between the 1950s and 1970s, this is an expected outcome if one used the data until 1990. In another comprehensive study, Hünicke and Zorita (2016) found that most of the accelerations estimated for the individual tide gauges were statistically insignificant in all three methods used in the study. The first method in their study is a regression model with a quadratic term that represents the acceleration/deceleration, similar to Woodworth (1990). The results of this method resemble those presented in Supplementary Table 2 although more tide gauges are used, including those of this study. When this model is used to estimate the acceleration, it is very critical to include the low-frequency variations in the regression model. Otherwise, the estimated acceleration might be statistically insignificant due to the misfit raised from these variations.    There is a uniform pattern for all tide gauges when the trends are compared with different periods. The higher trend in the most previous period demonstrates the acceleration in the rate of sea level rise in the basin.
In a case study on the Australian coasts, Agha Karimi (2021) reported how multidecadal variations can manipulate the magnitude of the estimated acceleration and its uncertainty when a quadratic term is used to represent the acceleration. Thus, if an additional term representing this low-frequency signal is added to the functional model of the acceleration estimation, a statistically significant figure may be estimated for individual tide gauges. Additionally, attributing a single figure for the rate of sea level rise for the whole period of analyses might be misleading as we are not aware of the exact point at which acceleration starts and also if it is subject to any change over the time series. The other two methods used by Hünicke and Zorita (2016) are based on a gliding linear trend (11-year segment) and fitting linear trend to the annual increment of the SLA. The latter suffers the same problem as the regression model with the quadratic term. The gliding linear trend, on the other hand, needs to have a longer segment length to filter out the impact of decadal and multidecadal signals. Due to the masking effect of multidecadal variations, Hünicke and Zorita (2016) could not reach a reliable figure for the acceleration of the Baltic Sea level. Visser et al. (2015) listed 30 different methods to analyze the sea level trend and acceleration and reported that there is not a perfect approach to quantify these parameters. They showed that every method has shortcomings and recommended using a combination of methods to draw a solid conclusion. However, this combination should be complementary in a way that the deployed methods cover the shortcomings of others. In the case of this study, the wavelet transform lacks uncertainty measures. However, it gives a general tendency of the SLA variation in the basin which helps in forming up the phases and breakpoints. To overcome the uncertainty measure, we use the linear trend fitted to each tide gauge in three different periods covering two or three of the phases ( Table 1). Considering the length of each period (at least 70 years), it can be assumed that the internal variability has only a trivial effect on the estimated trends. The hiatus is also intentionally involved in all three periods to account for its effect. All of the tide gauges show higher rates of rise in the period of 1944-2020 indicating that the sea level rise in the basin is accelerated over the last 120 years. If the rates from the altimetry era are taken into account, the acceleration is even more pronounced. Therefore, when this hiatus is considered, the sea level rise in the Baltic Sea is accelerating according to the trends fitted to the different multidecadal periods ( Table 1). It should be noted that the tide gauge locations are confined to the northern and eastern coasts of the Baltic Sea and this conclusion may not be extended to Latvian and Polish coasts due to some local drivers.

CONCLUSION
Multidecadal variations of the Baltic Sea level are analyzed using tide gauge data and altimetry measurements. Due to the paucity of reliable datasets with multicentennial length, multidecadal variations are investigated by fitting a linear regression model to multidecadal periods of the datasets and interpretation of lowpassed time series. As an independent source of data, satellite altimetry is used to validate the uplift model of NKG in two tide gauge stations by comparing the linear trends. NAO influences the sea level variations in the basin on seasonal to decadal time scales. It is shown that this impact is not substantial on the multidecadal scale by fitting a multivariate linear regression model which accounts for the NAO variations and also applying bandwidth filtering.
Previous studies could not detect any statistically significant sea level acceleration in the Baltic Sea level. It is demonstrated that the presence of multidecadal variations can lead to this outcome. As such, the global warming hiatus, in the 1950-1970s, causes a slowdown in the rate of sea level rise in the basin which prevents quantitative determination of the acceleration. To overcome this, the mechanism (or sea level contributor) which leads to this slowdown needs to be identified, modeled, and added to the functional model of the acceleration estimations. Despite the lack of this information, we tried to reveal the acceleration in the basin by comparing the sea level trend in periods double of this hiatus (6-7 decades) since other methods such as gliding segment with decadal length are also sensitive to multidecadal variations. The results used in this study imply the sea level acceleration in all eight tide gauges.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.