Sea Level Acceleration in the China Seas

Abstract

While global mean sea level rise (SLR) and acceleration (SLA) are indicators of climate change and are informative regarding the current state of the climate, assessments of regional and local SLR are essential for policy makers responding to, and preparing for, changes in sea level. In this work, three acceleration detection techniques are used to demonstrate the robust SLA in the China Seas. Interannual to multidecadal sea level variations (periods >2 years), which are mainly related to natural internal climate variability and significantly affect estimation of sea level acceleration, are removed with empirical mode decomposition (EMD) analysis prior to the acceleration determination. Consistent SLAs of 0.085 ± 0.020 mm·yr⁻² (1950–2013) and 0.074 ± 0.032 mm·yr⁻² (1959–2013) in regional tide gauge records are shown to result from the three applied approaches in the Bohai Sea (BS) and East China Sea (ECS), respectively. The SLAs can be detected in records as short as 20 years if long-term sea level variability is adequately removed. The spatial distribution of SLA derived from a sea level reconstruction shows significant SLA in the BS, Yellow Sea (YS) and Yangtze River Estuary.

1. Introduction

The Chinese coasts have been classified as one of the most vulnerable areas under climate change scenarios [1,2,3,4]. The rates of sea level rise (SLR) exhibit a large spatial variability, ranging from −1 to 10 mm·yr⁻¹ since 1950 [5]. During the satellite altimetry era, the Bohai Sea, Yellow Sea and East China Sea (BYECS, Figure 1a) experienced SLR rates of 4.44 mm·yr⁻¹, 2.34 mm·yr⁻¹ and 3.02 mm·yr⁻¹, respectively, (without glacial isostatic adjustment correction, GIA) in the time span of 1993–2012 [6]. In the last decade (2002–2012), ice melting from Greenland, Antarctic, and Glaciers are the main contributors to the coastal SLR (Figure 2 in [7]). Understanding local and regional sea level variability in China Seas is fundamental for implementing appropriate adaptation and mitigation measures to local SLR along the “wealth-belt” of China and low-lying coastal communities (e.g., Yangtze River Estuary, [8,9,10]), which are facing an increase in the frequency of flooding [11], salty water intrusions [12], coastline erosion [13] and ecosystem changes [14].

Knowledge regarding sea level acceleration (SLA, note that hereafter “A” in SLA is referred to “acceleration”, not the often-used “anomaly”) is crucial for future sea level projections. If the SLR rate is not constant (i.e., accelerates, decelerates, or changes due to long-term oscillations), it is difficult to accurately project the future SLR based on past sea level statistics, which is important for coastal protection planning and climate change-related measures. Recently, many efforts have been devoted to studying the regional or global SLA (e.g., [15,16,17,18,19,20,21]). Spada et al. [21] reviewed global mean SLR estimates published in scholarly journals since the 1940s for a heuristic assessment of global sea level acceleration. Due to the spatial variability of SLR, however, global mean SLA is not particularly useful for the public and policy makers [20]. Indeed, regional and local SLR and SLA can be quite different from global values [17,19], with some areas experiencing rates of SLR that are four times the global average during the altimeter era. In China Seas, the tide gauge data reflects a regional SLR rate of close to 3 mm·yr⁻¹ during the period of 1980–2014 [22], but questions remain regarding the influence of variability that can impact the estimation of the rate of SLR.

Long tide gauge records have been widely used to investigate the regional and global sea level accelerations [23]. Woodworth et al. [24] reviewed the evidence for accelerations in regional and global average sea level on several decades and longer timescales. Jevrejeva et al. [25] suggested an acceleration of 0.02 ± 0.01 mm·yr⁻² in tide gauge data based on global sea level reconstruction (1807–2009), which might have started at the end of 18th century [26]. There is paleo-evidence for the acceleration at the end of the 19th century [27,28]. Calafat and Chambers [29] suggest an acceleration of 0.022 ± 0.015 mm·yr⁻² between 1952 and 2011 after removing the contribution of internal climate variability.

Detecting and understanding the SLA along coastal regions is often complex with spatial variation relative to the global SLA [16,19,30]. Perhaps the best approach to allow for the earliest possible detection of a significant sea level acceleration lies in improved understanding (and subsequent removal) of interannual to multidecadal variability in sea level records [20]. On a regional scale over short periods of time, the uncertainty on acceleration estimates is dominated by decadal and multidecadal variability, which can largely mask a possible underlying acceleration due to global warming [29].

In addition to the tide gauge data, sea level reconstructions are also used to identify accelerations by focusing on global averages [25,26,31] and spatial distribution of the accelerations [24]. In this study, we explore the regional SLA based on tide gauge data, satellite altimetry data, and sea level reconstruction in the BYECS regions. The relative sea level changes in China Seas differ considerably from place to place [22,32]. As a result, the accelerations are estimated for the Bohai Sea (BS), Yellow Sea (YS), and East China Sea (ECS) separately. Visser et al. [33] reviewed 30 trend models applied in the field of sea level research and concluded that varying trend approaches may lead to contradictory acceleration-deceleration inferences. Therefore, the present study compares three different methods, including the linear least squares fitting method, evolution of SLR rates method [20], and sea level rate difference method [29,34] to determine the SLA in the past 55–65 years over the China Seas. The impact of long-term oscillations on the detection of SLA is also evaluated.

2. Data and Methods

2.1. Data

Monthly sea level records from nine tide gauge stations (Figure 1a) were obtained from the Permanent Service for Mean Sea Level [3]. The data include five, two, and two tide gauge stations along the BS, YS, and ECS coasts, respectively (LAOHUTAN and DALIAN have the same position). Information on the tide gauges can be found in

Table 1. Information of the selected tide gauges from PSMSL (* denote the data are quality control flagged). In the Bohai Sea, the correlations are calculated between the tide gauge records (IB and GIA applied) at DALIAN and other sites. In the Yellow Sea and East China Sea, the correlations are calculated from the selected two gauges, respectively.

ID	GIA	Time Span	Data Completeness	Correlations
Bohai Sea (1950–2013)
TANGGU	−0.36	1975–1994	100%	0.85
LAOHUTAN	−0.31	1980–1997	100%	0.97
QINHUANGDAO	−0.33	1950–1994	99%	0.86
DALIAN	−0.31	1954–2013	79%	–
YANTAI *	−0.28	1954–1994	100%	0.89
Yellow Sea (1975–1994)
SHIJIUSUO	−0.37	1975–1994	100%	1.00
LIANYUNGANG	−0.39	1975–1994	100%	–
East China Sea (1959–2013)
KANMEN	−0.37	1959–2013	99%	0.73
LUSI *	−0.41	1961–2013	83%	–

. The classical inverted barometer (IB) correction is applied at each tide gauge site using NCEP/NCAR reanalyzed (2.5° × 2.5° grids) monthly sea level pressure data [35] and the GIA correction is based on ICE-5G (VM2) model (Table 1) [36]. Flagged values with quality issues are removed at YANTAI (earlier than 1970). The data prior to 1965 are removed at LUSI, which are also data quality control flagged in PSMSL Revised Local Reference time series, for sparse available measurements.

Figure 1b shows the sea level time series at each selected site with IB and GIA applied. Consistent standard deviations of 11–15 cm, 11–12 cm, and 12 cm are shown in the BS, YS, and ECS, respectively (see Table 1 in [37]). In the BS, the high correlations (Table 1, on 95% confidence level) between DALIAN and other sites imply similar sea level variability at the sites. High correlations are also shown in the YS and ECS. In order to fill the data gaps in the tide gauge data and provide regional long sea level records for computing SLA, we averaged the tide gauge data in each regional sea based on the time of the measurements at each site (for details see [38]). Over 50 years of high quality data were available in the BS (1950–2013) and ECS (1959–2013), respectively. However, only 20 years of data can be used in the YS (1975–1994). Regarding the fluctuation of seasonal cycle over the study area [37], the signal in monthly tide gauge data are removed with the cyclostationary empirical orthogonal function analysis (CSEOF) to reduce the uncertainties in estimation of SLA [39,40].

The weekly sea level reconstructions used in the present study are on a spatial and temporal resolution of half a degree and span a period from January 1950 to June 2009 [40,41,42,43]. The dataset was developed based on the CSEOF analysis on combined sea level measurements from tidal gauge stations (collected by the PSMSL) and satellite altimeters (from Ssalto/DUACS). The sea level fields are reconstructed by fitting satellite altimetry-derived basis functions to historical tidal gauge data. Since the sea level reconstruction does not contain the seasonal signal and to allow for a direct comparison between the different sources of sea level data, the seasonal signal is calculated and removed from each dataset with the least squares fitting method.

2.2. EMD Analysis

In the present study, the sea level variability with periods >2 years are removed based on empirical mode decomposition (EMD) to further reduce the uncertainties of long term sea level variability on the estimation of SLA [26]. Due to the wide range of time scales involved and the non-stationary nature of climate change, the empirical mode decomposition and Hilbert-Huang transformation (EMD/HHT, [44,45]) was implemented for SLR and variability studies [16,17,18,46]. EMD/HHT is a non-parametric analysis for non-stationary time series. EMD decomposes a time series into a finite number, N, of intrinsic oscillatory modes $c_i\left(t\right)$ and a residual “trend” $r\left(t\right)$, so that an EMD of a sea level record from a given location (or from regional mean) is represented by

$$\mathrm{h}\left(t\right)={\displaystyle \sum }_{i=1}^N{c}_i\left(t\right)+r\left(t\right)$$

(1)

Particular modes do not necessarily represent specific processes, but the EMD analysis allows for the evaluation of relations between different data that may depend on different time scales. The EMD analysis is used to calculate and remove the long-term sea level variability in the tide gauge data and sea level reconstruction. For more details about the method refer to the studies mentioned above.

It is acknowledged that some studies indicate uncertainties in the ability of EMD to detect long-term variations and acceleration in sea level [30,47], especially if the records are not long enough, though this problem is not unique to EMD and according to [33], it is not clear which of the many possible analysis methods is preferred for detection of sea level acceleration. Furthermore, a recent study [18] evaluated the statistics of ensemble EMD calculations using very long European sea level records and shows that for records long enough (even with some gaps) EMD can detect acceleration that is almost identical to that obtained by the standard least squares fitting methods. The EMD calculations in those tests also demonstrate how to detect robust long-term oscillations and their statistical confidence levels.

Comparing the EMD to the least squares fitting method, the latter has the constraint that data is fitted to a particular shape (linear or polynomial, for example), so there is no assurance that the chosen shape is really the best choice. The EMD is more general than such fitting methods and can systematically filter out oscillating modes with unknown and variable frequencies.

2.3. Least Squares Fitting Method

We apply the least squares fitting method [19] to the tide gauge data and sea level reconstruction after removing the seasonal signal by:

$$h\left(t\right)=a+bt+\frac{c}{2}{t}^2$$

(2)

for the acceleration, or:

$$h\left(t\right)=a+bt$$

(3)

for linear regression. Here t is the time relative to the center of the selected time period and the coefficients a and b in Equations (2) and (3) denote the initial mean sea level and linear trend, respectively, and c in Equation (3) is acceleration. In regions that do not experience major, irregular changes in land level because of tectonic processes, one can make an assumption that vertical land movements have a linear time dependence over a century or so, suggesting quadratic or higher-order time dependence can be extracted from the records separately from vertical land motion effects [24]. Adopting the approach of Haigh et al. [20], the present work systematically estimates the accelerations (Equation (2)) or trends (Equation (3)) for all possible start dates and all possible date lengths in the BS and ECS. For example, when considering 50-year data lengths starting in 1950 in the BS, the subsets 1950–2009, 1951–2010, and so on, up to 1954–2013 were analyzed. The minimum subset length considered here is 10 years.

To calculate the uncertainties in the curve fitting, the confidence bounds (Figure 2) for fitted coefficients (e.g., trend or acceleration) are given by $C=b \pm r\sqrt{S}$, where b are the coefficients produced by the fit, r depends on the confidence level (on a 95% confidence interval), and S is a vector of the diagonal elements from the estimated covariance matrix of the coefficient estimates, (X^TX)⁻¹s². In a linear fit (Equation (2)), X is the design matrix, while for a nonlinear fit (Equation (3)), X is the Jacobian of the fitted values with respect to the coefficients. X^T is the transpose of X, and s² is the mean squared error. The uncertainties in the determined trend and accelerations are given as half of the interval width, i.e., (± (upper bound − lower bound)/2) (Figure 2).

Another method to estimate SLA is using sea level rate differences (SLRD, [25,30]). The given window size $\mathsf{\tau }$ is first split into two halves, and the SLRDs between the two half-window series can be calculated as:

$$SLRD (t_c,\mathsf{\tau })={\mathsf{\beta }}_{2 }(\frac{\mathsf{\tau }}{2}) - {\mathsf{\beta }}_{1 }(\frac{\mathsf{\tau }}{2})$$

(4)

where $t_c$ denotes the central year of the window. For the 60 year window over the 1950–2009 time period, the two half windows (30 years) of 1950–1979 and 1980–2009 are used to calculate the linear trends (using Equation (3)) ${\mathsf{\beta }}_1$ and ${\mathsf{\beta }}_2$ in Equation (4). For a more detailed (e.g., the calculation of the uncertainties of the method) description we refer to [25,30]. In this work, all uncertainties represent the 95% confidence interval.

3. Results

3.1. Effects of Lunar Nodal Cycle on Sea Level Acceleration Estimation

The effects of 18.61 years of the lunar nodal cycle on the SLA estimation have been investigated in the past for its considerable amplitude in the study area [8,48,49,50,51] (i.e., by adding item $\text{Acos}\frac{2\mathsf{\pi }\text{t}}{{\mathrm{T}}_{18.61}}+\mathrm{Bsin}\frac{2\mathsf{\pi }\text{t}}{{\mathrm{T}}_{18.61}}$ in Equation (1), and $\sqrt{{\mathrm{A}}^2+{\mathrm{B}}^2}$ denotes the amplitude of the lunar nodal cycle).

Table 2. Sea level acceleration (mm·yr⁻²) estimated from tide gauge data on the basis of the least squares fitting method in the selected time period. (* in brackets denotes the lunar nodal cycle is kept for SLR acceleration estimation.).

Time Span	1950–2013	1975–1994	1959–2013
Bohai Sea	0.086 ± 0.024 (0.100 ± 0.023 *)	0.305 ± 1.116 (0.719 ± 0.381 *)	0.098 ± 0.031 (0.097 ± 0.031 *)
East China Sea	–	−0.410 ± 1.103 (−0.059 ± 0.372 *)	0.044 ± 0.036 (0.046 ± 0.036 *)
Yellow Sea	–	0.364 ± 1.049 (0.577 ± 0.384 *)	–

lists the SLA calculated from least squares fitting (Equation (1)) in local tide gauge records. In the BS and ECS, the SLAs of 0.086 ± 0.024 mm·yr⁻² and 0.044 ± 0.036 mm·yr⁻² are obtained over a time span of 1950–2013 and 1959–2013, respectively (Table 2). The SLA of 0.098 ± 0.031 mm·yr⁻² is obtained in the BS in the period of 1959–2013, which tentatively demonstrates the impact of length of sea level records on local SLA estimation. The rates estimated from 20 years of tide gauge records are not statistically different from zero in the BS and ECS.

Note that in long records, the lunar cycle has little influence on the calculated acceleration, but in shorter records the results can be strongly affected by this cycle if it is not removed, as demonstrated in Table 2. Figure 2a,b show the scatterplots of tide gauge data with the fitting curve (sum of linear trend and acceleration in Equation (1)) overlaid in the BS and ECS, respectively. The grey and black curves denote the upper bound and lower bound (± (upper bound − lower bound)/2 denotes the uncertainty of the curve fitting) of the computed trend, respectively. In the BS, SLA of 0.100 ± 0.023 mm·yr⁻² (1950–2013) is obtained and the acceleration is shown from ~1980. In the ECS, the acceleration is observed from the start of the tide gauge records with a rate of 0.046 ± 0.036 mm·yr⁻² (1959–2013). The values are similar to those calculated with lunar nodal cycle removed. On the contrary, the computed SLA of 0.719 ± 0.381 mm·yr⁻² implies that the nodal cycle significantly affects the determined SLA in 20 years sea level records in the BS and YS. The accelerations are not statistically different from zero in the ECS over the selected 20 year time window.

3.2. Effects of Long-Term Sea Level Variability

3.2.1. On Sea Level Acceleration Estimation

Douglas [52] pointed out that one of the main problems in the detection of acceleration in observational records is an inadequate method for accounting for interannual and decadal variability in sea level time series [15,53]. The decadal variation of sea level dominates the acceleration for records shorter than ~50 years, so some suggest this record length as a minimum required for calculations of SLA from tide gauge data alone [15]. The accuracy of acceleration in the shorter records is not good and it is likely affected by unresolved long-term oscillations. Hence, the accelerations in the BS and ECS are estimated in the time spans of 1950–2013 and 1959–2013, respectively.

Figure 3 shows the decadal and multidecadal sea level variability (defined as sum of modes with T > ~5 years) in the tide gauge records estimated from EMD analysis without the residual ‘trend’. Significant long-term sea level variability is observed in the BS and ECS, which is mainly related to changes in atmospheric forcing (wind and pressure) and ocean mass or heat redistribution arising—for instance—from regional climate influence such as the Pacific Decadal Oscillation (PDO), El Niño–Southern Oscillation (ENSO), and North Pacific Gyre Oscillation [5,20,37,54,55]. The long-term sea level variations in the YS are consistent with that in the BS and ECS in the period of 1975–1994.

In Figure 3, the opposite phase of sea level variations between BS and ECS before 1980 compared to a coherent pattern after ~1980 might relate to the observed significant increases of sea level that started around 1980 (Figure 2), near the time when the PDO phase had shifted from negative to positive (see Figure 7 in [24]). The correlation coefficients of 0.72 (−0.81) and 0.85 (−0.91) were obtained between sea level and long-term Southern Oscillation Index (SOI) (Pacific Decadal Oscillation, PDO) variations (e.g., calculated with EMD analysis) in the BS and ECS over 1980–2013, respectively. Before 1980, the correlation coefficients of 0.55 (−0.61) and 0.77 (0.48) were obtained between sea level variations and SOI (PDO) in the BS and ECS, respectively. Therefore, in recent years, long-term sea level variations are more strongly coupled with climate indices than they were in the earlier period.

Table 3. Determined sea level acceleration (mm·yr⁻²) using the least squares fitting method with long term sea level variations in tide gauge data removed on the basis of EMD analysis. The number in the brackets denote the acceleration calculated from (1) evolution of sea level rise rates and (2) sea level rate differences.

Time Span	1950–2013	1959–2013
Bohai Sea	0.085 ± 0.020 (0.081 ± 0.021 1 0.085 ± 0.008 2)	0.075 ± 0.028
East China Sea	–	0.074 ± 0.032 (0.070 ± 0.038 1 0.073 ± 0.010 2)
Yellow Sea (1975–1994)	−0.022 ± 0.337	–

lists the sea level accelerations with the interannual, decadal, and interdecadal sea level variations were estimated and removed with EMD analysis (Equation (1)). The SLA of 0.085 ± 0.020 mm·yr⁻² and 0.074 ± 0.032 mm·yr⁻² are obtained in the BS and ECS, respectively (Table 3). Consistent SLA is obtained in the BS and ECS over the same selected time period of 1959–2013. Compared with the values listed in Table 2, the differences indicate the significant effects of long-term sea level variability on estimation of SLA in the study area, particularly in the ECS. In the YS, the SLA is not statistically different from zero because of the short (20 years) tide gauge records.

The atmospheric pressure forcing significantly affects coastal sea level variability in the BS [5,37]. Using the same method, we estimated the SLA in the BS and ECS, but without IB correction applied. The SLAs of 0.112 ± 0.024 mm·yr⁻² and 0.082 ± 0.034 mm·yr⁻² are obtained in the BS and ECS, respectively. Compared with the results in Table 1, the differences are attributed to the contributions of atmospheric pressure on regional SLA estimation should be considered in the study area, particularly in the BS [37].

3.2.2. On Sea Level Rise Rate Estimation

Recent attempts to detect a significant acceleration in regional SLR might underestimate the impact of natural variability [56]. In the BS, SLR rates of 1.6 ± 0.2 mm·yr⁻¹ and 1.4 ± 0.2 mm·yr⁻¹ are obtained over the periods of 1950–2013 with and without the long-term sea level variability removed based on the EMD analysis. The local SLR rates could be higher if the land vertical motion linear rates of tide gauge stations are considered [57]. In the ECS, the rates of 3.4 ± 0.3 mm·yr⁻¹ and 4.0 ± 0.2 mm·yr⁻¹ are obtained from the sea level time series with and without long-term sea level variability removed over the time span of 1959–2013, respectively. For the selected time windows, the removal of long term signals significant decreases and increases the determined SLR rates in the BS and ECS, respectively.

3.3. Sea Level Acceleration Estimated with Evolution of Sea Level Rise Trend Approach

Another commonly used approach for acceleration detection is the estimation of linear rates for consecutive overlapping periods. The method is applied to the tide gauge records with long term sea level variability (>2 years) in the entire available time series removed using the EMD analysis. The linear regression (Equation (2)) over overlapping periods (10 years, 20 years, …, 50 years and 60 years) are calculated to determine the evolution of the SLR trend (Figure 4).

In Figure 4, the time denotes the central year of the selected consecutive overlapping periods (when considering 10-year data lengths on a sea level record starting in 1950, the subsets 1950–1959, 1951–1960, and so on, up to 2004–2013 were analyzed; when considering 50-year data lengths, the subsets 1950–1999, 1951–2000, and so on, up to 1964–2013 were analyzed). The sample size chosen influences one’s findings on accelerations. With the length of the window increasing, the linear trends are more statically significantly and different from zero. When the time window of >20 years (Figure 4b–f) is used, stable increases of the SLR trend presented from the start of the selected period in the BS. That is different with the trend pattern with the long term signal remaining in Figure 2a, when the SLA shown started at ~1980.

In the BS, the SLR trend varies from 1.2 ± 0.2 mm·yr⁻¹ to 1.6 ± 0.2 mm·yr⁻¹, corresponding to the periods 1950–1979 and 1954–2013, respectively (Figure 4e). In the ECS, the SLR trend changes from 3.9 ± 0.3 mm·yr⁻¹ to 4.2 ± 0.3 mm·yr⁻¹ over the time periods 1959–1983 and 1989–2013, respectively (Figure 4e). The patterns also confirm the robust positive SLA that exists in the study area, since a linear increase in SLR as shown in Figure 4 indicates a constant acceleration. The rise rate of the SLR trend of 0.081 ± 0.021 mm·yr⁻² and 0.070 ± 0.038 mm·yr⁻² are in agreement with that estimated by SLA of the least squares fitting method in Table 3.

On a local scale, it is difficult to assess whether observed changes are due to natural or anthropogenic causes because the presence of considerable natural and anthropogenic variability means that the observed sea level curve is not necessarily represented by Equations (1) or (2). Hence, it can lead to estimates of rates of change (if real on a local scale) that mostly reflect internal natural variability, rather than providing evidence for, or against, externally forced accelerations [20]. In Figure 4, the same method is used to compare SLR obtained with and without the long-term sea level variability, and the impact of the long-term variability on the detection of SLA is very clear. Significant SLR rate variability over ~30 years are shown with selected window of 20 years (Figure 4b). That implies the removal of long-term sea level variability allows for detection of such accelerations. Furthermore, the SLAs of 0.087 ± 0.014 mm·yr⁻² and 0.069 ± 0.005 mm·yr⁻² are obtained from linear rates for consecutive overlapping period of 20 years in the BS and ECS, respectively (Figure 4b). The results are close to the estimations from 60 year and 50 year overlapping periods in the BS and ECS, which implies the removal of long-term sea level variability shortens the detection time of significant SLA in the study area.

To further demonstrate the removal of long-term signals on the evolution of SLR trends, the blue and red dashed curves with circles in Figure 4d,e show the evolution of SLR trends, but with the EMD analysis applied to each 50-year and 40-year data segments in the BS and ECS, respectively. Compared with the EMD analysis applied to all available data prior to SLR trend computation, strong variations are shown in the determined SLR trend time series. This implies the differences of the long-term signal in each selected segment may affect the determination of SLA. In the BS (ECS), the SLAs of 0.086 ± 0.003 mm·yr⁻² (0.067 ± 0.006 mm·yr⁻²) and 0.114 ± 0.034 mm·yr⁻² (0.031 ± 0.063 mm·yr⁻²) are obtained from the sea level time series with EMD applied to all available data and each selected 50 (40)-year segment, respectively.

3.4. Sea Level Acceleration Estimated with the SLRD Approach

The third commonly-used method for SLA estimation is the SLRD approach applied to the sea level time series. Figure 5 shows the evolution of the SLRD using 50-years (e.g., half-window of 25 years and the SLRD difference between the subsets 1975–1999/1950–1974, 1976–2000/1951–1975 and so on up to 1989–2013/1964–1988 were analyzed) and 60-years (half window of 30-years) windows in the BS and ECS (Equations (3) and (4)). The calculations start from 1950 and 1959 for the data in the BS and ECS, respectively. Note that, unlike Figure 4 where an upward trend indicates positive acceleration, in Figure 5, SLRD > 0 indicates positive acceleration. The small variations in SLRD demonstrates that SLA is almost constant in the tide gauge data (Figure 5). The SLRD is calculated from all of the available data for the BS and ECS, respectively. The obtained SLA of 0.085 ± 0.008 mm·yr⁻² and 0.073 ± 0.010 mm·yr⁻² are consistent with the results obtained from the other approaches (Table 3). If the same SLR rate continues over the next 100-years it would result in SLR of ~42 cm and ~36 cm in addition to the rise due to the linear rate in the BS and ECS coastal regions.

In Figure 5, the SLRDs, with the long-term sea level variability included, are overlaid for comparison. The sea level deceleration is seen in the ECS and significant variations shown in the 50-year SLRD evolution without the removal of long-term signals. Overall, robust SLAs are obtained from three approaches with long-term sea level variability removed prior to the estimation. The removal of long-term sea level variability allows the acceleration to be detected earlier.

3.5. Spatial Pattern of Sea Level Acceleration

In studies of global sea level change, reconstruction techniques [40,41,42,43] are often used to combine tide gauge and altimeter data and account for internal variability. Sea level reconstructions are particularly useful to study spatial distribution of sea level accelerations that is not provided by tide gauge data alone [24]. In this section, the long term sea level variability is left in the tide gauge data and the spatial pattern of SLA is determined with two approaches.

The geographical distribution of SLA and uncertainties estimated from weekly sea level reconstruction (1950–2008) using the linear least squares fitting method (Equation (2)) is shown in Figure 6a,b, respectively. In Figure 6a, significant SLA of >0.06 mm·yr⁻² is observed in the BS, YS, and coastal regions of the ECS. Higher SLA (>0.08 mm·yr⁻²) is presented at DALIAN (Figure 1), the western YS and Yangtze River Estuary with uncertainties <0.015 mm·yr⁻². The SLA decrease to 0.04 mm·yr⁻² in the open ocean of the BS and ECS with an uncertainty of ~0.005 mm·yr⁻². The distribution patterns estimated using the SLRD method (Figure 6c, Equation (4)), without the removal of long-term signals, are similar to Figure 6a. The similarity between the patterns in Figure 6a,c confirms again that the two methods are consistent with each other when applied to reconstruction data, as they were for the tide gauge data.

The estimated linear SLR rates (Equation (3)) for the same time period are also shown in Figure 6d. Compared with Figure 6a, the regions with a trend higher than 2 mm·yr⁻¹ are accompanied with the SLA > 0.08 mm·yr⁻². The similarity between the pattern of the linear and non-linear sea level trends in Figure 6 (with both showing the highest rates in the YS and lower in the ECS) is quite interesting and may indicate local effects from coastal dynamics, such as wind-driven currents, river flow, and thermostatic effects.

4. Discussion

Detecting sea-level acceleration is especially difficult because of the need to separate between oscillatory changes and long-term trends [18]. Woodworth [24] suggested that a possible contribution to regional sea level accelerations may come from the ENSO pattern. Changes in the frequency of El Nino events are known to have occurred during the 20th century, most notably and recently during the climate shift around 1976 related to the phase change of PDO. The associated changes in heat content and wind field could have resulted in corresponding sea level changes. Along the western coasts of the North Atlantic, ocean dynamics, and variations in the Gulf Stream in particular, show considerable impact on SLA, whereas a climate-related slowdown of the Gulf Stream seems to lead to an acceleration along the U.S. East Coast [17,34]. However, in the China Seas the influence on coastal sea level from the Kuroshio is not yet fully understood and requires further investigation in the future and, instead, impacts from local coastal dynamics and Pacific variations, such as PDO, may contribute to variations in SLR rates.

Past studies exhibit apparent regional sea level variability on interannual and longer time scales in the China Seas [54]. The sea level variations are highly related to atmospheric forcing over the continental shelf in the BS and YS [5,37,58]. The regional sea level variation in the ECS was influenced not only by local factors but also by remote wind from the adjoining ocean with the oceanic connectivity influenced by upper-ocean circulation [59]. ENSO could influence the local wind field and, thus, the interannual sea level variation indirectly, and that impact is more apparent in the southern area than in the northern area in the ECS [60]. In the China Seas, the decadal sea level variability is mainly modulated by the wind stress curl variability related to the PDO [55], which is strongly correlated with vertically-integrated transport carried by the Kuroshio through the ECS [61].

The effects of long-term sea level variability on the detection of SLA is, thus, addressed in the present work, using three different approaches. Standard curve-fitting and filtering methods would have difficulties separating the sea level trend from long-term variations, such as the global 60-year oscillation cycle [26,62]. This difficulty led to the development of a new analysis method to remove high-frequency oscillations, as well as cycles with long periods that are comparable to the record length itself [17,18].

The accelerations might be detected much earlier on local scales if the processes behind internal variability are adequately understood, which would subsequently allow the removal of these components from the records [19]. At the regional scale, the required period is much longer than the global sea level trend and acceleration determination and decadal variability determines the regions with larger detection times. The needed detecting time (minimum record length) is larger when/where the actual trend or acceleration is smaller or when/where the underlying interdecadal variability is larger [63].

5. Conclusions

The calculated SLA is dependent on the length of sea level records and the particular method that is used. In the present study, the long-term sea level variability, including the 60-year oscillation cycle mode in front of the last residual “trend” mode in EMD analysis, are removed prior to determination of SLA. Based on three commonly-used approaches for SLA estimation, consistent SLA of 0.085 ± 0.020 mm·yr⁻² (1950–2013) and 0.074 ± 0.032 mm·yr⁻² (1959–2013) are observed in the BS and ECS, respectively. The acceleration is not statistically different from zero in the YS, where the record length (20 years) is relatively short.

In the China Seas, the sea level variability on time scales of 20–30 years are shown in the tide gauge data (Figure 3 and Figure 4a–c). With the long-term sea level variability removed, the evolution of SLR is stable and 20-year tide gauge data may be sufficient for the determination of SLA in the BS and ECS if care is given to remove long-term oscillations (Figure 4b,c).

We also estimated the distribution of SLA on the basis of sea level reconstruction in the study area. Although the sea level constructions demonstrate weaker long-term sea level variability, the dataset still shows significant temporal viability. Coherent SLA spatial patterns are shown between the linear least square fitting and SLRD approaches. Significant SLA is shown in the BS, YS and Yangtze River Estuary with relatively lower SLA presented in the open ocean of the BS and ECS. The results presented here demonstrate ways to detect acceleration in sea level studies, which will eventually improve estimates of future SLR trends. This, in turn, will help the development of mitigation and adaptation strategies for shallow depths offshore and low lying regions, which are under risk of increasing flooding by tropical cyclones as a result of accelerated sea-level rise in the BS and ECS [64].