Open AccessArticle

Sea Level Rise Estimation on the Pacific Coast from Southern California to Vancouver Island

Xiaoxing He

Jean-Philippe Montillet

^2,3

Rui Fernandes

Timothy I. Melbourne

⁴,

Weiping Jiang

⁵ and

Zhengkai Huang

^1,6,*

School of Civil and Surveying & Mapping Engineering, Jiangxi University of Science and Technology, Ganzhou 341000, China

Physikalisch-Meteorologisches Observatorium Davos/World Radiation Center (PMOD/WRC), CH-7260 Davos, Switzerland

Institute Dom Luiz, University of Beira Interior, 6201-001 Covilhã, Portugal

⁴

Pacific Northwest Geodetic Array, Department of Geological Sciences, Central Washington University, Ellensburg, WA 98926-7417, USA

⁵

GNSS Research Center, Wuhan University, Wuhan 430079, China

⁶

School of Transportation Engineering, East China Jiao Tong University, Nanchang 330013, China

Author to whom correspondence should be addressed.

Remote Sens. 2022, 14(17), 4339; https://doi.org/10.3390/rs14174339

Submission received: 18 June 2022 / Revised: 25 August 2022 / Accepted: 30 August 2022 / Published: 1 September 2022

(This article belongs to the Special Issue Modelling Geodetic Time Series and Applications for Earth Science and Environmental Monitoring)

Download

Browse Figures

Graphical abstract
"> Figure 1
Spatial distribution over the entire west coast of North America showing the 405 GNSS stations and 31 TG sites used in this study. Note that <a href="#remotesensing-14-04339-f0A1" class="html-fig">Figure A1</a> (in the <a href="#app1-remotesensing-14-04339" class="html-app">Appendix A</a>) displays the western coast in more detail (the black lines are fault boundaries). "> Figure 2
Histogram of the percentage of data gaps for all GNSS time series analyzed in this study. "> Figure 3
Up component of the GNSS daily position time series for the SEAT station with functional model on top (red line) of the observations. "> Figure 4
Coastal latitudinal profile of VLM for British Columbia [BC] (Vancouver Island, Canada) and Washington State [WA] (inland Puget Sound, USA) based on the PANGA solution (with the BIC_tp information criterion). The stars are the location of the TGs. The red is the uncertainty associated with the VLM estimates, the green lines are the mean values of the VLM for each region. "> Figure 5
VLM estimated for each region (Vancouver Island (Canada), Salish Sea + Inland Puget S. (USA), Olympic Peninsula (USA), South WA (USA) + Oregon plains (USA), Cape Blanco Triple junction + Point Arena (USA), Southern California (USA), whole Pacific coast). The VLM is estimated using either the PANGA or the NMT solution. Note that we have also included the results using the various ICs (AIC, BIC, BIC_tp). "> Figure 6
Coastal profile of the vertical land motion for British Columbia [BC] (Vancouver Island, Canada), Washington state [WA] (Olympic Peninsula, south of WA, USA), Oregon [OR] plains (USA), California [CA] (north and south, USA). The stars are the location of the TG. The red is the uncertainty associated with the VLM estimates. The green lines are the mean values of the vertical land motion. "> Figure 7
Interpolation of the VLM uncertainties for both NMT and PANGA products based on BIC_tp (the fault zone/line is in red). "> Figure 8
Relative Sea Level Rise Rate and associated uncertainties estimated at various TGs along the Pacific coast from Vancouver Island (Canada) to Southern California (USA) using different stochastic noise models. The abbreviations of the TG names refer to <a href="#remotesensing-14-04339-t001" class="html-table">Table 1</a> and <a href="#remotesensing-14-04339-f001" class="html-fig">Figure 1</a>. "> Figure 9
Sea Level Rise Rate (blue) and ASLR (red) for the TG stations in the Pacific coast. Note that the RSLR is produced using the BIC_tp. The VLM estimates used to process the ASLR are based on the BIC_tp and the GGM model. Here, only the results with the PANGA solution are displayed. The flooding risk is added (center to 0 mm/yr) as the vertical bar from low (blue) to high (yellow). The grey line is the GIA estimated from an ensemble of models. "> Figure 10
Histogram of the absolute SLR for the various regions along the Pacific coast from Vancouver Island (Canada) to Southern California (USA). "> Figure 11
Absolute Sea Level Rise Rate and the mean sea level with satellite altimetry (SSH) estimated at the same location as the TG of PANGA solution. "> Figure A1
Zoom of <a href="#remotesensing-14-04339-f001" class="html-fig">Figure 1</a>. "> Figure A2
Absolute SLR (ASLR) estimated using TG and VLM (NMT) and the satellite altimetry mean sea level (SSH) estimated at the same location of the TG. ">

Versions Notes

Abstract

Previous studies have estimated the sea level rise (SLR) at various locations on the west coast of the USA and Vancouver Island in Canada. Here, we construct an entire SLR profile from Vancouver Island in the Pacific Northwest to San Diego in Southern California. First, we process global navigation satellite system (GNSS) measurements at 405 stations blanketing the whole coast to generate a profile of vertical land motion (VLM) known to bias century-long tide gauge (TG) measurements recording relative SLR (RSLR). We are then able to estimate the absolute SLR (ASLR) by correcting the SLR with the VLM. Our study emphasizes the relationship between the various tectonic movements (i.e., the Cascadia subduction zone, the San Andreas strike-slip fault system) along the Pacific coast which renders it difficult to accurately estimate the SLR. That is why we precisely model the stochastic noise of both GNSS and tide gauge time series using a combination of various models and information criterions (ICs). We also use the latest altimetry products and sea surface height (SSH) to compare it with ASLR at the same location as the TGs. This study supports previous analysis that the power law + white noise and generalized Gauss–Markov + white noise models are the best stochastic noise models for the GNSS time series. The new coastal profile confirms the large variability of VLM estimates in the Pacific Northwest around the Cascadia subduction zone in agreement with previous studies, and a similar result when the San Andreas fault comes onshore in Central California (San Francisco Bay). Negative RSLR values are mostly located in the Pacific Northwest (Vancouver Island and Olympic Peninsula). We also observe a much bigger variation (about 90–150%) of the ASLR in the Pacific Northwest which is predominantly due to glacial isostatic adjustment (GIA). Moreover, the comparison between the ASLR and the SSH estimates shows similarities in the center of the studied area (South Washington, Oregon planes, and some parts of Southern California) where the tectonic activity does not significantly influence the TG measurements. Finally, the twentieth-century satellite geocentric ocean height rates show a global mean of 1.5 to 1.9 mm/yr. Our estimates based on ASLR and SSH are within this interval.

Keywords:

stochastic noise model; geodetic time series; GNSS; tide gauge; sea level rise; Pacific coast

Graphical Abstract

1. Introduction

Present-day crustal uplift on the west coast of the USA and Vancouver Island (Canada) is controlled by two main solid Earth processes: (1) the interseismic loading of the Cascadia subduction fault on the Northwest Pacific coast and the San Andreas strike-slip fault system in Central and Southern California; (2) the post-glacial rebound mostly in the Northwest Pacific caused by the cyclic growth and decay of ice sheets during the quaternary until the end of the last ice age [1,2,3,4,5].

Large global-scale sea level variations (around 100–200 mm) occurred on very long scales (geological timescales of ≈100 million years) due to the tectonic processes which shape the Earth’s crust [5,6,7,8,9]. One of the lowest sea levels was observed during the last glacial maximum 26,500 years ago. With so much of the planet’s water tied up in ice, the global sea level was more than 122 m (≈400 feet) lower than the present-day estimates [10]. In recent times (i.e., the last 200 years) since the beginning of the industrial era, scientists have observed that the rapid variations of mean sea level are directly related to anthropogenic global warming [11,12]. It is one of the indexes monitoring the response of our environment to this phenomenon. There is a very high probability that sea waters expand and sea levels rise (SLR) due to the warming of the oceans. The melting of mountain glaciers due to the increase in air temperature contributes to SLR with the input of fresh water into the oceans.

Using geodetic observations, previous studies [13,14,15,16] have estimated the relative sea level rise using tide gauges (TGs) [17,18,19,20]. However, TGs cannot measure the absolute sea level change, but the height of the sea surface relative to crustal reference points that may move with tectonic activity or local subsidence. In other words, the TG observations are biased by local and regional processes that are linear or non-linear over a multi-decade timescale. Linear processes include glacial isostatic adjustment (GIA) and interseismic tectonic strain accumulation, whereas the non-linear ones include, for example, earthquakes; the non-linearity of earthquakes generally consists of all the transient signals such as the post-seismic relaxation recorded in the time series [21,22,23,24]. Therefore, the SLR estimated from TGs must be corrected from the vertical land motion (VLM) in order to get a precise estimate of the absolute SLR (ASLR) [25,26,27,28]. When dealing with century-long TG records, the estimation of the SLR and associated uncertainties is a source of error due to the inherent stochastic noise [29]. Therefore, one must carefully model the various processes and the temporally correlated noises in the TG measurements in order to estimate accurately the rate and the associated uncertainty, which is called the relative SLR (RSLR) [24]. Temporally correlated noises affect different types of time series including geodetic time series [30,31,32]. It results that each observation can be correlated with previous ones. Various models have been developed [33,34,35] in geodetic time series analysis. However, previous works (e.g., [31]) show that the power spectral density (PSD) of the TG measurements can be described by a power-law stochastic model, which is a colored noise. It is a 1/f^α noise, with α varying between [0, 2]. When the exponent of the colored noise is set to 0, the noise is called white noise, at 1 it is defined as flicker noise, and at 2 it corresponds to random walk [36].

For comparison with ASLR estimated at a TG, we also use the sea surface height (SSH) product, estimating the location of each TG modulo the resolution of the satellite altimetry product. Satellite altimetry measures the SSH above a benchmark or datum, whereas the TG benchmark is on the land close to the instrument. TG thus observes the relative sea level, respective to the elevation of the benchmark. Sea level altimetry measures the sea level with the reference to the geoid. The SSH is the height of the sea surface above a reference ellipsoid [37]. This is the direct product recorded by the satellite altimetry. The SSH values are provided along the satellites’ ground tracks or at regular grids interpolated from the values determined along the satellite tracks, e.g., the Copernicus Marine Environment Monitoring Service (CMEMS) provides regular and systematic reference information (data products) on the physical and biogeochemical ocean and sea ice state for the global ocean and the European regional seas. The CMEMS is implemented and operated by Mercator Ocean, which provides oceanographic products and services for maritime safety, coastal and marine environment, climate and weather forecasting, and marine resource users [38,39].

Besides, GNSS daily position time series contain unmodeled signals (e.g., small offsets, episodic tremors, and slips) which affect the precise estimation of the geophysical signals [36]. In this research, we continue the work of He et al. [40,41] where a methodology based on the ICs was developed to select the optimum noise model together with a functional model to describe the various geophysical signals contained in the time series. It allows estimating robustly the VLM from the GNSS daily position time series (daily average solution) [36,40,41]. Using 405 permanent GNSS stations blanketing the Pacific Northwest down to the west coast of the USA, from Vancouver Island to the bottom of California (San Diego), with a concentration of stations around the Los Angeles area monitoring the San Andrea fault, it allows us to build a coastal profile of the VLM. In addition, we compare our estimates using the GNSS solutions from two processing centers, namely the Pacific Northwest Geodetic Array (PANGA) and New Mexico Tech (NMT). The PANGA and NMT products are available freely at https://www.unavco.org (accessed on 10 January 2022).

The stochastic noise model selection is further applied to estimate the SLR from TG records and SSH measurements. We focus on investigating the impact of stochastic noise selection on the estimated RSLR and ASLR. Intrinsically, we need to evaluate the impact of the complex geodynamics characterizing the west coast of the USA. Note that in this work, we are correcting the TG measurements with the collocated GNSS stations following previous studies [24,42,43,44,45,46].

2. Materials and Methods

2.1. Data Processing

GNSS Data Processing: We used observations from continuously operating GNSS receivers distributed over the west coast of the USA (see red circle points in Figure 1).

The daily position time series (daily average solution) results from the processing released by the PANGA and the NMT [41] were computed within the International Terrestrial Reference Frame 2008 [47]. The data are available online, and Herring et al. [48], Montillet et al. [24], and He et al. [41] described comprehensively the GNSS data processing strategy for the two processing centers. In summary, the raw GNSS phase and code observations were point-positioned with ambiguity resolution using GIPSY software developed and supported by the NASA Jet Propulsion Laboratory (JPL, [49]). Satellite ephemerides, clock corrections, and wide-lane phase bias estimates were provided by JPL [50]. In our approach, station positions were loosely constrained during initial estimation and subsequently transformed into the ITRF08 using only the translation and rotation, but not scale, components of the JPL-provided Helmert transformations. The NMT processing was carried out using the software GAMIT/GLOBK [51,52] where additional global stations were also included for the stability of the reference frame. The difference with the processing carried out at NMT is mainly due to how the scale parameter is handled. Montillet et al. [24] emphasized that the choice of whether to include a radial scaling degree of freedom during daily reference frame realization primarily impacts the average network radial height and produces apparent height anomalies in excess of 5 mm that persist for months. In order to reduce the impact of missing data on the estimated results, we selected the data period from 2008 to 2018, ensuring that the data missing rate of all stations is less than 5.5%. Figure 2 shows the data gap of the analyzed 402 GNSS sites; it can be seen from Figure 2 that the analyzed 402 GNSS sites have low data gaps. Three GNSS stations have a larger data gap (i.e., P426, COUP, and GLDR are about 24.0%,12.0%, and 14.0%, respectively) but they were included to compare the results with previous studies. The GNSS time series were analyzed with FN+WN, RW+FN+WN, GGM+WN (generalized Gauss–Markov + white noise), and PL+WN (power law + white noise) stochastic noise properties [41].

Tide Gauges Processing: monthly data were downloaded from the Permanent Service for Mean Sea Level website. The monthly files list the date (year–month in decimal form) and the average monthly mean sea level value [53,54]. Figure 1 shows the spatial distribution of the 31 TG stations used in this study. The TG ID, latitude, and longitude are also displayed in Table 1. Previous studies show that the stochastic properties of TG time series can be modeled with an autoregressive moving average model (ARMA (p,q)), autoregressive fractionally integrated moving average (ARFIMA (p,d,q)), and GGM [55]. Therefore, in this research, we adopted the same models to analyze the noise and estimate the sea level rate from the TG and SSH time series. Note that p and q are the lags in the model and d is the fractional parameter [24,55]. Here, we set

p \in (0, 5)

q \in (0, 5)

SSH product: The SLR is estimated as the SSH above a benchmark or datum. For a TG, this benchmark is some adjacent land point. Therefore, a TG measures RSLR (relative to the elevation of the benchmark) as opposed to the height of the sea surface above a reference ellipsoid (e.g., ASLR). This is the direct product of satellite altimetry. Additionally, the SSH values are provided along the satellites’ ground tracks or at regular grids interpolated from the values determined along the satellite tracks, e.g., CMEMS provides regular and systematic reference information on the physical and biogeochemical ocean and sea ice state for the global ocean and the European regional seas [38,39]. Here, we used a “GLOBAL_REANALYSIS_PHY_001_030” reanalysis data product to obtain global sea surface high monthly grid data with a data resolution of 5’ covering the period 1993 to 2019 (26 years in total). It cannot properly cover all areas due to insufficient resolution of the gridded data for the SSH product. For this purpose, the coast closest to the tide gauge station was used as the interpolation point for the search within 15 km [27], especially for the TG located closest to inland narrow channels (e.g., SEATTLE, PORT CHICAGO). For the TG stations closest to inland narrow channels, we interpolated to the nearest coast which is in the resolution of our product, but the value can be highly spatially correlated with other nearby TG stations.

2.2. Stochastic and Functional Model Estimation

The GNSS position time series is a sum of stochastic processes and geophysical signals. Previous studies [40,41,56,57,58] have shown that long GNSS time series (i.e., over 8 years) are more suitable to extract accurately geophysical signals by jointly fitting a stochastic and functional model. The time-dependent model is generally a linear sum of a seasonal signal, tectonic rate, co-seismic offsets, and random stochastic processes [36,57].

Here, we followed the methodology developed in Montillet and Bos [36], He et al., 2019, and He et al., 2021 [40,41]. We processed the time series before estimating jointly the stochastic and functional models using a maximum log-likelihood estimator via the Hector software [59]. This methodology includes first removing outliers and correcting known offsets retrieved from the station’s information. Note that spurious observations were considered outliers if these points were larger than 3 times the interquartile range of the residual time series following Langbein and Bock [60]. The residual time series was obtained by coarsely removing the tectonic rate from the original observations. Finally, the automatic offset detection algorithm developed by Fernandes and Bos [61] was applied to detect undocumented offsets. Following He et al. [40], the trajectory model was a function of a linear trend, an annual and semi-annual periodic signal, together with the offsets found in the previous step.

To determine the optimal stochastic model of the GNSS and TG time series, we used the Ics described in [36,40]:

A I C = - 2 \log (L) + 2 v

(1)

B I C = - 2 \log (L) + \log (n) v

(2)

B I C_t p = - 2 \log (L) + \log (\frac{n}{2 π}) v

(3)

where n is the number of data points, a larger number of parameters v increases the AIC and BIC value and thus serves as a penalty term, and

\log (L)

is the log-likelihood function. For more details see He et al., 2019 [31,36].

However, some time series from stations mostly located in the Pacific Northwest and some in Southern California are subjected to specific events due to local geodynamics such as episodic tremors and slips (ETS) in the Cascadia subduction zone [62,63]. We used the arctangent function [64] in order to model the post-seismic relaxation characteristic of ETS. Note that the time of the slow slip event and the delay of the post-seismic deformation are required as input parameters for the estimation of ETS using Hector [41]. The time of the slow slip events can be requested from the Pacific Northwest Geodetic Array website or with a careful analysis of the time series (see the additional data material). The start of a slow slip event was evaluated via the correlation of seismic data together with a careful check of each time series [65]. In the remainder of this work, we use only a delay equal to 100 days for the post-seismic relaxation because it is difficult to precisely estimate the duration of the crustal relaxation due to the noise. This delay is a conservative number following the study of [41], knowing that the repetition of the ETS events is ≈14 months as evaluated by previous geophysical studies of Cascadia [65]. This value represents a tradeoff in not modeling enough the relaxation phenomenon, and oppositely, absorbing other geophysical phenomena due to an overestimation of the relaxation time [24,41,66]. Figure 3 displays an example of the functional model including slow slip events superimposed on the observations at the up component of the SEAT GNSS station from the PANGA solution.

3. Results

3.1. Overview of the Tectonics of the Pacific Coast and Coastal Uplift Profile

Broadly speaking, the Pacific coast is divided by two different geodynamic areas: the Pacific Northwest and the center to Southern California. The Pacific Northwest down to the triple junction at Cape Mendocino is a subduction zone defined by the Cascadia mountains where the Pacific plate subducts under the North American plate. It is a very long, sloping subduction zone where the Explorer, Juan de Fuca, and Gorda plates move to the east and slide below the much larger mostly continental North American Plate [67]. The triple junction is the location of a change in the broad plate motions which dominate the west coast of North America, linking convergence of the northern Cascadia subduction zone and translation of the southern San Andreas fault system. The fault passes along the towns of Gorman, Tejon Pass, and Frazier Park, and begins to bend northward, forming the “Big Bend”. This restraining bend is thought to be where the fault locks up in Southern California. The southern segment, which stretches from Parkfield in Monterey County all the way to the Salton Sea, and 56.0 km close to the northeast of Los Angeles, can generate large earthquakes (wit moment magnitude > 8.0). The fault is a right-lateral strike-slip with various degrees of compression forces. Along the fault, the locking zone varies due to the geodynamic complexity created from:

(1) The network of faults running parallel to the San Andreas fault in Northern and Central California.

(2) The change of geometry (e.g., upturned around the Salton Sea area).

(3) The various degrees of compressional forces between the Northern American plate towards the Pacific Plate, forcing the San Andreas fault to jog westward. The complexity of this fault system, which includes the sister faults (e.g., Rodgers Creek Fault, Hayward Fault, Calaveras Fault) can be witnessed by the creation of various mountain ranges along the fault, e.g., the Transverse Ranges in Southern California and Santa Cruz Mountains [68,69,70].

3.2. Vancouver Island, the Olympic Peninsula, and Puget Sound

Previous studies [1,24,71] show that the Cascadia fore arc is divided in three areas. The Cascadia subduction zone generates a large uplift rate observed on the northern part of Vancouver Island and the Olympic peninsula with an order of magnitude 2.0 mm/yr on average, resulting from the combination of the postglacial rebound and the subduction interseismic strain, Figure 4 supports previous results, where the VLM values around Vancouver island are gradually increasing from the Olympic peninsula (1.40 mm/yr) reaching around 2.3 mm/yr on average (see green line in Figure 4), and large uplift at the north of the island (e.g., Quadra island-QUAD 4.2 ± 0.5 mm/yr (PANGA); 3.5 ± 0.2 mm/yr (NMT), Port Alberni (PTAL) 3.5 ± 0.3 mm/yr (PANGA), 2.5 ± 0.1 mm/yr (NMT)). We can underline that there is a discrepancy between our results and Montillet et al., 2018 at some stations. For example, ELIZ experiences a 2.8 ± 0.2 mm/yr (PANGA) whereas in the previous study it was reported an uplift of 2.5 ± 0.2 mm/yr. In addition, we estimate for KTBW with an uplift of −0.1 ± 0.3 mm/yr (PANGA), whereas in the previous study the estimate was about −0.1 ± 0.2 mm/yr. These discrepancies are due to the specific modelling of the geophysical signals (i.e., considering the ETS) and the optimum stochastic noise model selection.

On the opposite, the inland waterways and Puget Sound have a uniform small amplitude, with some places with negative values generated by an oblique fault with shallow fault systems in the central Puget Sound lowland [72,73]. Overall, Figure 4 supports previous studies [1,24]. Inland Puget sound is characterized by small positive or negative VLM values but with associated large uncertainties, e.g., SEAT at 0.5 ± 0.7 mm/yr (PANGA), −1.4 ± 0.8 mm/yr (NMT); P423 −0.2 ± 0.4 mm/yr (PANGA); −0.8 ± 0.2 mm/yr (NMT)). Thus, the new profile confirms the large variability of VLM estimates in the Pacific Northwest around the Cascadia subduction zone in agreement with previous studies [1,24]. Looking at Table A1 (in the Appendix A), some stations are experiencing large discrepancies with previous studies, namely Montillet et al., (2018) and Mazzotti et al., (2007) as discussed above for the uplift on Vancouver Island. Many of the Cascadia slow slip events are affecting the time series for the stations west of the Cascades range, and are hard to model away without a full slip model for each event, knowing that GNSS stations installed in the late 1990s may have recorded more than hundreds of them [70]. Our approach is based on estimating all the slow slip events together with a stochastic noise model selected using different information criterions (Ics, e.g., AIC/BIC/BIC_tp) [34]. Therefore, the parameter estimation suffers from under-/overestimation depending on factors such as the stochastic noise characteristics of the data and the type of estimator used [40,41]. In our case, we are using the maximum likelihood estimator [36].

Looking at the differences between the various Ics in Table A2 (see Appendix A), the stochastic properties of VLM estimates do not vary for about 98.0% of the stations for both the PANGA and NMT solutions. This is due to the low percentage of stochastic noise model change. Table A2 (in the Appendix A) shows the difference in model selection for each station and for both the PANGA and NMT of the analyzed GNSS stations. The results show that the difference of the AIC, BIC, and BIC_tp selection is marginal for most of the sites, with only a few sites displaying a significant difference. Therefore, only the BIC_tp is used to select the optimal stochastic noise model in the following analysis.

Table 2 shows that the PL+WN model still appears to be the best noise model, accounting for 81.0% and 61.0% of the PANGA and NMT solutions, and the GGM+WN accounts for 14.0% and 34.0%, respectively, consistent with He et al. [40,41]. Besides, we made a statistical analysis of the ratio of velocity and velocity uncertainty with the selected optimal noise model for the PANGA and NMT solutions. The NMT is noisier than the PANGA solution when we estimate the ratio between velocity and velocity uncertainty (cf. Table 3).

3.3. South of Washington State and Oregon

Figure 3 shows that VLM decreases towards the south of Washington State (WA), reaching almost zero at the top and central parts of Oregon (Coos Bay). The VLM estimates are below 1 mm/yr (at P396 0.8 ± 0.2 mm/yr (BIC_tp–PANGA), −0.2 ± 0.4 mm/yr (BIC_tp–NMT)); at CHZZ −0.2 ± 0.3 mm/yr (PANGA), and −0.8 ± 0.4 mm/yr (NMT–BIC_tp)). This result is further emphasized by looking at the mean of the VLM (green line) which is around 0.0 mm/yr between central Oregon (plains) and the south of Washington state.

Figure 5 compares the mean and variance region-by-region where all the 405 stations processed are used. We can observe that the two solutions are comparable in terms of variance for all regions of the Pacific coast. However, there is an offset in terms of mean value ≈±1.2 mm/yr depending on the selected region. For example, for the Oregon–Southern Washington region, the mean is centered at 0.0 for the PANGA solution and ≈−1.2 mm/yr for the NMT solution. This result underlines the difference in the strategy of these two processing centers.

3.4. Cape Blanco–Cape Mendocino and Point Arena in Northern California

Towards the south of Oregon, at Cape Blanco and down to Cape Mendocino, the VLM increases progressively (CABL 1.5 ± 0.1 mm/yr (PANGA); 0.8 ± 0.2 mm/yr (NMT)) to reach more than 2mm/yr for the PANGA solution at PTSG 2.87 ± 0.34 mm/yr, and ≈2.0 ± 0.5 mm/yr for the NMT product. However, the uplift estimated from the NMT solution increases up to 1.4 ± 0.3 mm/yr when not estimating any post-seismic relaxation or transient signal (just co-seismic offsets). These uplift increases are due to the geophysical/geodynamic activities at the Mendocino triple junction, where the Cascadia interseismic strain accumulation locks the subducting Gorda slab underneath the North American plate [74].

Figure 6 also shows that the VLM coastal profile between Cape Blanco and south of Cape Mendocino (Northern California) is not smooth with drops (at TRND −1.0 ± 0.6 mm/yr (PANGA); −1.2 ± 0.9 mm/yr (NMT)). This instability is due to the numerous co-seismic offsets and geophysical signals necessary for estimation caused by the intense geodynamic activity in this particular area. The San Andreas fault originates offshore from the triple junction bending to follow the Northern California coastline underwater up to Point Arena, then with various on- and offshore areas (e.g., Point Reyes), and continues completely inland along Southern California from south of the San Francisco Bay [75]. We can observe in Figure 6 that the uplift is not large when the fault is offshore (e.g., near Punta Gorda at P157 0.5 ± 0.4 (PANGA) and 0.6 ± 0.5 mm/yr (NMT)), whereas the uplift is larger when the fault is onshore (at Point Arena—P187 3.3 ± 0.8 mm/yr (PANGA); 2.8 ± 0.6 mm/yr (NMT)). Therefore, the non-smoothness of the coastal VLM profile translates into the intense geophysical activity intrinsic to this area. Looking at Figure 5, the PANGA uplift estimates in this region are on average close to zero mm/yr, but with a standard deviation of 2.0 mm/yr, which is twice the standard deviation of the southern Washington and Central Oregon area. Note that the NMT values are similar, with a slightly smaller mean value (≈−0.5 mm/yr).

3.5. Central and Southern California

The southern San Andreas Fault system is onshore passed San Francisco Bay. The center section runs near Parkfield and Hollister, experiencing a phenomenon called aseismic creep [76,77], where the fault slips continuously without causing any earthquakes. However, the southern part, running through Parkfield, has experienced several large earthquakes, e.g., the 1989 Loma Prieta earthquake with a moment magnitude of 6.9 [78]. The southern segment is generally defined as the part of the fault running from Parkfield in Monterey County to Bombay Beach on the border of the Salton Sea, passing through the base of the San Gabriel and San Bernardino mountains. These mountain ranges are the testimony of San Andreas fault activity [79].

Figure 6 shows that the coastal VLM profile decreases from Point Arena to a minimum at the south of the San Francisco Bay (P178: −2.1 ± 0.4 mm/yr, PANGA), however, the NMT product provides an estimate of the same order of magnitude (i.e., −2.7 ± 0.6 (NMT)) of P187 (Point Arena). This shows the difficulty of estimating reliably the VLM closest to an active fault where the two strategies have a noticeable impact on the uplift estimates. However, the VLM profile decreases further from the fault and from the shoreline. For example, after Point Reyes and the back of L.A., at the CIRX station we estimate 0.7 ± 0.3 mm/yr (PANGA) and −0.2 ± 0.5 mm/yr (NMT). Towards the south of L.A. and San Diego the VLM decreases (P478 −1.5 ± 0.4 mm/yr (NMT); −0.6 ± 0.3 mm/yr PANGA). In Figure 5 the scatter of the uplift estimates for the whole Southern California region is ≈3.0 mm/yr, bigger than for Vancouver Island and the Olympic Peninsula and slightly smaller than in the Inland Puget Sound area for the PANGA values. On average, the NMT product for the whole Southern California region provides comparable values estimated for Vancouver Island and the Olympic Peninsula combined together. The homogeneity of the VLM coastal profile is difficult to realize in Southern California due to the nature of the San Andreas Fault (right-lateral strike-slip) with various degrees of compression along the fault, compared with the Cascadia subduction mechanism and the VLM pattern in the Pacific Northwest [80]. However, we can notice that Southern California and the Inland Puget Sound are two areas with a complex geophysical mechanism, therefore generating a non-smooth VLM coastal profile [81].

Finally, Figure 7 displays the linearly interpolated uncertainties of the VLM around the whole west coast for both the PANGA and NMT products. It clearly shows the active geodynamic areas along the Pacific coast which agree with the VLM coastal profile derived in Figure 4 and Figure 6. For example, the subduction zone in Inland Puget Sound (WA) produces large uncertainties, together with the activity linked to the triple junction in Northern California. One can also underline the difference between the two solutions, where the large uncertainties seem to decrease in amplitude when the San Andreas fault comes onshore for the PANGA estimates, but it extends further south with the NMT solution. Moreover, the NMT product shows more smoothness around geodynamic active regions (Vancouver Island + Puget Sound, Southern California).

4. Discussion

4.1. RSLR and ASLR Estimation along the Pacific Coast

We estimate the RSLR [82] using 31 TGs located along the Pacific coast between 32°N and 51°N (see Figure 1 and Figure 8, Table 1 and Table A3 together with RSLR using SSH product). To recall Section 2, the same methodology is used to select the stochastic noise model with the ICs when estimating the relative SLR, and to a certain extent the mean sea level with the SSH product. Negative RSLR values are located in the Pacific Northwest—Vancouver Island and Olympic Peninsula—at stations such as Campbell River (CAMP) (−1.8 ± 0.4 mm/yr) and Neah Bay (Ne.B) (−1.8 ± 0.1 mm/yr). The average over the whole of Vancouver Island is −1.0 ± 0. 6 mm/yr. Two factors are here biasing these values: (i) the large uplift estimated from the GNSS time series due to the geodynamics of the Cascadia subduction zone; (ii) the GIA activity [24,63]. For the Puget Sound area, the mean RSLR value is around 1.3 ± 0.3 mm/yr. For example, at the TG of Seat, the value is ≈2.1 ± 0.1 mm/yr (see Table A3). Towards the south of Washington and the center of Oregon, the RSLR values are on average 1.1 ± 0.8 mm/yr, except at Astoria (ASTO) (−0.2 ± 0.2 mm/yr) due to the location of the TG inside the estuary of the Columbia River as discussed in [1,24,81].

Towards the triple junction, the RSLR decreases abruptly at Crescent (CRES see Table 1 and Figure 1) (−0.8 ± 0.2 mm/yr) and increases at its maximum along the whole Pacific coast at Humboldt Bay (HUMB) (5.7 ± 0.9 mm/yr) which is located near Cape Mendocino. From Humboldt Bay to Southern California, the RSLR values are large, with an average of 1.8 ± 0.4 mm/yr (La Jolla—LaJo). The large increase at Humboldt is probably due to a combination of various factors. First, the TG is located relatively near the triple junction tectonic activity, where the VLM is negative (P162 −1.6 ± 0.8 mm/yr PANGA) which may be the result of the tectonic activity. Secondly, we cannot discard the effect of ocean eddies on TG measurements due to the geophysical location of this TG. The authors in [83] have discussed the potential bias due to strong ocean eddies on the estimation of mean sea level rise. Central and Southern California show a large RSLR, which can be associated with a high risk of flooding in coastal areas and cities as discussed in NRC (2014) and in [82,84].

Note that in terms of comparing the differences of RSLR estimates between the stochastic noise models, we can only observe a marginal difference, especially between the GGM and ARMA models with a level of difference below 0.1 mm/yr on average, i.e., nothing significant. The ARFIMA can generate much larger uncertainties such as Port Chicago, CA (P. Ch. see Table 1) with a RSLR ≈2.4 ± 1.3 mm/yr compared with the GGM estimate of ≈2.0 ± 0.9 mm/yr and the ARMA estimate ≈2.0 ± 0.8 mm/yr. One may also observe similar results between the three ICs, especially BIC and BICtp. The RSLR estimates with the BIC_tp, can vary slightly. For example, at Patricia Bay (P.Ba), the RSLR is ≈0.4 ± 0.2 mm/yr (ARFIMA–BIC_tp) and ≈1.1 ± 1.1 mm/yr (ARFIMA–BIC and BIC_tp).

4.2. Absolute SLR along the Pacific Coast of the USA Using SSH and Corrected TG Measurements

We correct the RSLR with VLM interpolated from our coastal profile in order to estimate the ASLR [85]. Figure 9 displays the RSLR (BIC_tp), ASLR (BIC_tp), glacial isostatic adjustment (GIA, ensemble of models from NRC 2015), and the risk of coastal flooding (see line center on 0 mm/yr). The results displayed are only for the PANGA solution, but the associated values for both the PANGA and NMT are listed in Appendix A (see Table A4 and Table A5). To simplify this study, we select here only the BIC_tp criterion (for there are not many differences between the ICs for the VLM and the RSLR. However, AIC and BIC have known biases, which is why He et al., 2019 and 2021 [40,41] have shown that BIC_tp is more reliable) in the optimal noise model selection in order to estimate the VLM profile used in the interpolation at each TG. This simplification is based on the results in Section 3 where we expose the marginal differences between the criterions. In addition, we have only chosen the GGM model following our previous results on the RSLR. The GIA model (gray line) is produced by an ensemble of models provided by the previous study by the NRC (2014). In addition, the indices of flooding are provided by the work of [86].

Figure 10 clearly shows that when using the PANGA product the average ASLR for the Pacific Northwest (Vancouver Island and Puget Sound) is ≈3.0 mm/yr (PANGA–GGM), whereas the value is shifted to ≈1.2 mm/yr in Central and Southern California. Compared with the NMT product, the mean ASLR values between the Pacific Northwest and Central/Southern California are ≈3.1 mm/yr and ≈0.8 mm/yr, respectively. Therefore, we observe a much bigger variation (about 90.0–150.0%) after correcting from VLM the RSLR in the Pacific Northwest which is predominantly due to the GIA. The VLM correction to the RSLR accounts for both GIA and tectonic strain accumulation. Isolating the two sources of uplift is a difficult exercise that requires precise knowledge of mantle viscosity, recent glaciation history, and interseismic coupling of the faults around the coast [24]. It was shown by Montillet et al. [24] that this correction of the RSLR with VLM in the Pacific Northwest indicates that half of the VLM on Vancouver Island arises from GIA and the rest from the subduction-related strain. The GIA correction is much smaller in Central and Southern California. It is also observed that the estimates are close to the profile estimated from satellite altimetry, especially in the Los Angeles and San Diego areas for both the PANGA and NMT. Therefore, we can conclude that in this area, the geodynamics of the continent are not biasing the RSLR as much as in other places along the Pacific coast such as in the Pacific Northwest–Vancouver Island, Olympic Peninsula, and inland Puget Sound.

Moreover, we can compare the estimation of the SLR with either the SSH or the ASLR (cf. Figure 11 for the PANGA solution and Figure A2 for the NMT solution). Overall, the SLR estimates with SSH are all positive values across the entire coast. This result is expected because the satellite altimetry is not affected by the underlying geodynamical movements due to the vertical land motion affecting the TG measurements. Because satellite altimeters have a limited ability to measure sea level height in coastal regions where the TGs are installed, the SSH measurements are not affected by the VLM [87].

The second result is the general trend. The velocity uncertainty estimated with the SSH are small, on average 1.3 ± 0.2 mm/yr (ARMA), compared with the ASLR derived from the TG with a value to 1.3 ± 0.3 mm/yr (ARMA–PANGA) or 0.5 ± 0.4 mm/yr (ARMA–NMT). In the areas with underlying complex geodynamics, we have a greater SLR. The SLR shows a large average value of 2.9 ± 1.1 mm/yr, whereas the ASLR (TG+VLM) is around 2.0 ± 1.3 mm/yr (ARMA–PANGA). In Southern California, starting at Point Reyes, we have a steady increase with an average SSH of 1.4+/1.7 mm/yr. The ASLR estimated in the same region is around 1.6 ± 1.3 mm/yr (PANGA–ARMA) or 0.8 ± 1.1 mm/yr (NMT–ARMA). Overall, the SSH and ASLR products are generally comparable, especially with the PANGA solution. Furthermore, one factor generating the differences between the two products is the ocean eddies, which have been known to impact the monitoring of SLR [83]. Several methods have been developed that mix the TG and SSH products [88], but there is no consensus on how to reconcile both the TG and SSH measurements due to many underlying geodynamical phenomena intrinsic to each product.

Finally, twentieth-century rates for the Western Pacific are more similar to twenty-first-century global rates of 3.2 ([2.8–3.6]) mm/yr [89], whereas the Eastern Pacific satellite rates are significantly lower than the global mean value. However, twentieth-century satellite geocentric ocean height rates show a global mean of 1.5 to 1.9 mm/yr, or roughly half that of twenty-first-century rates (IPCC, 2013). We estimate the mean SLR over the whole Pacific coast as shown in Table 4.

Table 4 shows that the uncertainties of the ASLR using the NMT product to estimate the VLM are much higher than the ASLR value for all three models. The VLM estimates from the NMT product have been shown to differ substantially from the PANGA product. This result was already investigated in [24] where the analysis of the processing from the two centers reveals that this discrepancy is due to a processing strategy and how the scale from the Helmert transformation is handled as discussed in Section 2.1. This issue was also highlighted in [41] when comparing both processing strategies (PANGA, NMT) and emphasizing the largest variations for the NMT product in the estimated tectonic rates and vertical land motion values from an analysis of hundreds of GNSS stations located in North America. Furthermore, we have used stations from the whole Pacific coast which includes two regions with high tectonic activity (i.e., the Cascadia subduction zone and the San Andreas fault). Even though we have carefully modeled each GNSS time series, multidecade transients will impact rate estimates differently based on both the time series duration and the structure of any known long-term transients. Compared with several studies, e.g., [24], we have here used criteria to optimally select the stochastic noise models. Our analysis emphasizes that the GNSS product [e.g., PANGA, NMT] should be carefully chosen when estimating the ASLR. The second result from the analysis of Table 4 is that the ARFIMA model tends to produce larger uncertainties. This result should be further investigated in the next study.

Finally, our estimates with the PANGA and SSH are within [1.5–1.9] mm/yr. There is an ongoing discussion in the scientific community about the discrepancy between the GMSL and the ocean height which focuses on the estimation of the acceleration in sea level rise rates (IPCC, 2013). However, this topic was left for a future study.

5. Conclusions

This paper investigates the selection of the stochastic noise models of geodetic time series (i.e., GNSS daily position, TG, and SSH) and the estimation of sea level rise on the Pacific coast from Vancouver Island (CA) to Southern California (USA). First, we analyze the stochastic noise properties of 405 permanent GNSS stations distributed over the west coast of the USA using the NMT and PANGA products. The selection of the optimum stochastic noise (i.e., FN+WN, RW+FN+WN, GGM+WN, and PL+WN) model is based on the ICs (i.e., AIC/BIC/BIC_tp). The crustal uplift is estimated by a joint inversion of functional (i.e., geophysical signals) and stochastic noise models. We then provide a coastal profile for the VLM.

Secondly, we estimate the RSLR from 31 century-long TG records. The estimates are then corrected with the VLM interpolated at each TG location from our coastal profile in order to obtain the ASLR. Besides, for comparison with our ASLR estimates, we use the SSH product estimated at the location of each TG modulo the resolution of the selected satellite altimetry product (with grid data with a data resolution of 5′).

The following conclusions are drawn from our study:

(1): For the 405 analyzed GNSS daily position time series, the PL+WN model still appears to be the best noise model, i.e., about 81.0% and 61.0% of the PANGA and NMT solutions. The GGM+WN accounts for about 14.0% and 34.0% of the PANGA and NMT, respectively. Overall, the values for the NMT product are noisier than the PANGA solution, which is consistent with [40,41]. Besides, the stochastic properties of VLM estimates are not varying using the various ICs for about 98.0% of the stations for both the PANGA and NMT solutions.
(2): The Cascadia forearc is divided into three areas: Vancouver Island, the Olympic Peninsula, and Puget Sound, among them the Cascadia subduction zone generates a large uplift rate observed on the northern part of Vancouver Island and the Olympic Peninsula with an order of magnitude about 2.0 mm/yr on average, which is caused by the combination of the postglacial rebound and the subduction interseismic strain, whereas the inland Puget Sound is characterized by small positive or negative values VLM values. This result supports previous studies (e.g., Mazzotti et al., 2007; Montillet et al., 2018) that the VLM values around Vancouver Island are gradually increasing from the Olympic Peninsula [1,24]. We also underline that some stations do not experience as much uplift as reported in the previous work of Montillet et al., (2018) [24]. These discrepancies are due to the specific modeling of the geophysical signals (e.g., the ETS events) and the optimum stochastic noise model selection. In addition, the new profile confirms the large variability of VLM estimates in the Pacific Northwest around the Cascadia subduction zone in agreement with previous studies.
(3): The VLM decreases towards the south of WA and the Oregon region. We also conclude that the PANGA and NMT processing are comparable in terms of variance for all regions of the Pacific coast. From Cape Blanco down to Cape Mendocino the VLM increases progressively, which is due to the geophysical activities at the Mendocino triple junction. For Central and Southern California, the NMT product for the whole Southern California region provides comparable values estimated for Vancouver Island and the Olympic Peninsula combined.
(4): We estimate the RSLR and ASLR along the Pacific coast. The negative RSLR values are mostly located in the Pacific Northwest—Vancouver Island and the Olympic Peninsula—with stations such as Campbell River (Camp).
(5): We observe a much bigger variation (about 90.0–150.0%) of the ASLR in the Pacific Northwest which is predominantly due to the GIA. Moreover, we compared the estimation of the ASLR with the SSH. The SLR estimated with the SSH product are all positive values across the entire coast. This result is expected because the satellite altimetry is not affected by the underlying geodynamical movements due to the VLM affecting the TG measurements. They are comparable for the center of the coast (Southern WA, Oregon planes, and some parts of Southern California) where the tectonic activity does not influence the TG measurements. However, the discrepancy between the SLR and the SSH is still discussed within the scientific community due to many factors such as the underlying geodynamics and ocean eddies. Our analysis also emphasizes the need to carefully chose the GNSS product that can introduce different variations of the VLM and then influences the estimated ASLR.
(6): Finally, we compare our various estimates with the twentieth-century satellite geocentric ocean height rates, which are between 1.5 and 1.9 mm/yr. Our estimates with the PANGA and SSH are consistent with the previous studies.

As Hammarklint [89] pointed out, global ocean levels have always fluctuated with changes in our climate. The underlying factors causing sea level change and the relationship between sea level change and climate need to be further studied.

Author Contributions

X.H. and J.-P.M., writing—original draft preparation; R.F., T.I.M. and W.J., methodology, review and editing; Z.H. and X.H., data processing and figure plotting. All authors have read and agreed to the published version of the manuscript.

Funding

This work was sponsored by National Natural Science Foundation of China (42104023), Jiangxi University of Science and Technology High-level Talent Research Startup Project (205200100564, 205200100588), Youth Talent Plan of the Science and Technology Think Tank of China Association for Science and Technology in 2022 (No.207), 2021 Jiangxi Province Higher Education Teaching Reform Research Project (JXJG-21-7-37), and Jiangxi Provincial Natural Science Foundation (20212BAB204030).

Data Availability Statement

The data for both the PANGA and NMT processing centers are available freely at https://data.unavco.org/ (cwu.final_nam08.pos.tar and nmt.final_nam08.pos.tar) (acceesed on 10 January 2022). The TG measurements are available at https://www.psmsl.org/ (acceesed on 10 January 2022). The SSH product is downloaded from https://resources.marine.copernicus.eu/ (acceesed on 10 January 2022).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1. Comparing the velocity (u) and uncertainty (sigma) with previous studies of Montillet et al. (2018) and Mazzotti et al. (2007).

Site	This Work				Montillet et al. 2018 [24]				Mazzotti et al. 2007 [1]
	PANGA		NMT		PANGA		NMT		Mazzotti et al. 2007 [1]
	u	Sigma	u	Sigma	u	Sigma	u	Sigma	u	Sigma
ALBH	1.2	0.1	0.3	0.1	0.7	0.2	0.8	0.3	1.1	0.9
DRAO	0.4	0.2	−0.1	0.2	1.0	0.2	1.2	0.3	1.2	0.7
NANO	1.7	0.5	1.7	0.6	2.2	0.3	1.8	0.4	2.5	0.9
NEAH	3.4	0.1	2.2	0.1	3.2	0.2	3.2	0.3	3.5	1.0
PGC5	1.0	0.2	0.4	0.1	0.8	0.2	0.1	0.5	1.8	1.0
SEAT	0.5	0.7	−1.4	0.8	0.1	0.3	−0.2	0.3	−0.6	0.9
UCLU	1.8	1.0	−0.7	1.1	2.5	0.2	1.9	0.3	2.7	0.9
BAMF	2.1	0. 4	1.9	0.8	2.7	0.4	1.8	0.4
BCOV	2.7	0.2	1.8	0.3	2.8	0.2	3.6	0.7
CABL	1.5	0.1	0.8	0.2	1.2	0.2	1.4	0.2
CHZZ	−0.2	0.3	−0.8	0.4	0.2	0.4	0.8	0.2
ELIZ	2.8	0.2	1.7	0.1	2.5	0.2	2.6	0.4
HOLB	1.9	0.2	1.3	0.2	2.4	0.2	0.9	1.0
KTBW	−0.1	0.3	−1.0	0.5	−0.5	0.2	−0.4	0.3
NTKA	2.7	0.2	2.0	0.3	3.6	0.2	4.3	0.4
P159	−1.2	0.4	−1.7	0.6	−0.8	0.3	−1.6	0.3
P161	−1.5	0.5	−2.1	0.8	−1.0	0.2	−1.5	0.3
P162	−1.6	0.8	−2.2	0.7	−1.2	0.2	−1.6	0.3
P316	−2.1	0.8	−2.3	0.7	−2.2	0.5	−2.1	0.6
P362	2.0	0.3	1.6	0.3	2.8	0.3	2.1	0.4
P364	1.9	0.2	1.2	0.2	2.3	0.3	1.7	0.4
P365	0.5	0.2	−0.1	0.2	1.0	0.3	0.0	0.4
P366	0.4	0.2	−0.7	0.2	0.7	0.3	−0.6	0.3
P367	−0.3	0.3	−1.1	0.2	−0.2	0.3	−0.8	0.4
P395	0.6	0.4	0.2	0.7	0.2	0.4	−0.2	0.3
P396	0.8	0.2	−0.2	0.4	1.1	0.5	0.2	0.4
P398	0.7	0.4	−0.1	0.1	1.5	0.3	0.6	0.4
P402	2.1	0.2	1.4	0.2	2.5	0.2	1.7	0.5
P423	−0.2	0.4	−0.8	0.2	−0.4	0.2	−0.9	0.3
P435	−0.2	0.3	−0.6	0.3	0.6	0.4	0.1	0.4
P437	0.6	0.8	−1.4	0.8	−0.4	0.3	−1.4	0.7
P439	1.0	0.6	−0.6	0.7	0.0	0.2	−0.3	0.4
P734	2.4	0.4	1.7	0.5	3.2	0.3	2.0	0.4
PABH	0.5	0.2	−0.3	0.1	0.2	0.2	0.2	0.3
PCOL	−0.0	0.9	−2.4	0.9	−0.6	0.3	−0.6	0.3
PTAL	3.5	0.3	2.5	0.1	3.5	0.1	0.0	0.6
PTSG	2.9	0.3	2.0	0.5	3.6	0.2	3.0	0.3
QUAD	4.2	0.5	3.5	0.2	4.3	0.4	3.9	0.4
SC04	1.4	0.2	0.7	0.3	1.2	0.2	1.0	0.2
TPW2	0.5	0.1	−0.3	0.1	0.2	0.2	0.5	0.2
TRND	−1.0	0.6	−1.2	0.9	−0.9	0.3	−0.7	0.3

Table A2. Optimal noise model differences between the various ICs for the analyzed 405 GNSS sites.

Solution	Site	AIC	BIC	BIC_tp
PANGA	CHWK	GGMWN	PLWN	PLWN
	MIDA	PLWN	FNWN	FNWN
	P283	PLWN	FNWN	FNWN
	P315	PLWN	FNWN	FNWN
	P316	PLWN	FNWN	FNWN
	SHLD	GGMWN	PLWN	PLWN
	KTBW	GGMWN	PLWN	PLWN
NMT	P156	PLWN	FNWN	FNWN
	P178	PLWN	FNWN	FNWN
	P188	PLWN	FNWN	FNWN
	P267	PLWN	FNWN	FNWN
	P273	PLWN	FNWN	FNWN
	P312	PLWN	FNWN	FNWN
	PVRS	PLWN	FNWN	FNWN
	KTBW	GGMWN	PLWN	PLWN

Table A3. Relative sea level rise (RSLR) derived from TG and SSH time series.

RSLR TG	GGM		ARFIMA BIC_tp		ARMA BIC_tp		ARFIMA BIC		ARMA BIC		ARFIMA AIC		ARMA AIC
RSLR TG	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma
0010	1.5	0.1	1.5	0.2	1.5	0.1	1.5	0.2	1.5	0.1	1.5	0.2	1.5	0.1
0127	2.1	0.1	2.1	0.2	2.1	0.1	2.1	0.2	2.1	0.1	2.1	0.2	2.1	0.1
0158	2.3	0.2	2.2	0.1	2.3	0.1	2.2	0.1	2.3	0.1	2.2	0.1	2.3	0.1
0165	−1.4	0.3	−1.3	0.2	−1.4	0.3	−1.3	0.2	−1.4	0.3	−1.3	0.2	−1.4	0.3
0166	0.8	0.1	0.8	0.1	0.8	0.1	0.8	0.1	0.8	0.1	0.8	0.1	0.8	0.1
0245	1.1	0.2	1.1	0.1	1.1	0.2	1.1	0.1	1.1	0.2	1.1	0.1	1.1	0.2
0256	2.2	0.2	2.1	0.1	2.1	0.1	2.1	0.1	2.1	0.1	2.1	0.1	2.1	0.1
0265	−0.2	0.2	−0.2	0.2	−0.2	0.2	−0.2	0.2	−0.2	0.2	−0.2	0.2	−0.2	0.2
0377	1.8	0.2	1.8	0.2	1.8	0.2	1.8	0.2	1.8	0.2	1.8	0.2	1.8	0.2
0378	−0.8	0.2	−0.8	0.1	−0.8	0.2	−0.8	0.2	−0.8	0.2	−0.8	0.1	−0.8	0.2
0384	1.2	0.2	1.2	0.1	1.2	0.2	1.2	0.1	1.2	0.2	1.2	0.1	1.2	0.2
0385	−1.8	0.2	−1.8	0.1	−1.8	0.2	−1.8	0.1	−1.8	0.2	−1.8	0.1	−1.8	0.2
0437	0.9	0.2	0.9	0.2	0.9	0.2	0.9	0.2	0.9	0.2	0.9	0.2	0.9	0.2
0508	1.0	0.3	1.0	0.2	1.0	0.2	1.0	0.2	1.0	0.2	1.0	0.2	1.0	0.2
0527	−0.6	1.3	−0.7	1.6	−0.6	1.0	−0.7	1.6	−0.6	1.0	−0.7	1.6	−0.6	1.0
1152	0.9	0.7	1.0	1.0	0.9	0.7	1.0	1.0	0.9	0.7	0.7	0.4	0.9	0.7
1196	1.7	0.4	1.7	0.6	1.7	0.5	1.7	0.6	1.7	0.5	1.8	0.3	1.7	0.5
1242	0.0	0.4	0.0	0.4	0.0	0.5	0.0	0.4	0.0	0.5	0.0	0.4	0.0	0.5
1269	1.1	0.5	1.0	0.4	1.1	0.5	1.0	0.4	1.1	0.5	1.0	0.4	1.1	0.5
1323	−1.8	0.5	−1.8	0.4	−1.8	0.5	−1.8	0.4	−1.8	0.5	−1.8	0.4	−1.8	0.5
1325	1.8	0.5	1.8	0.4	1.7	0.5	1.8	0.4	1.7	0.5	1.8	0.4	1.7	0.5
1352	1.6	0.5	1.5	0.4	1.6	0.5	1.5	0.4	1.6	0.5	1.5	0.4	1.6	0.5
1354	0.2	0.9	0.1	0.4	0.2	0.9	0.1	0.4	0.2	0.9	0.1	0.4	0.2	0.9
1394	2.2	0.5	2.2	0.5	2.2	0.5	2.2	0.5	2.2	0.5	2.2	0.5	2.2	0.5
1639	5.6	0.8	5.7	1.2	5.6	0.7	5.7	1.2	5.6	0.7	5.7	1.2	5.6	0.7
1640	1.0	1.0	1.0	0.9	1.0	1.0	1.0	0.9	1.0	1.0	1.0	0.9	1.0	1.0
1799	−0.1	0.9	−0.2	0.5	−0.1	0.9	−0.2	0.5	−0.1	0.9	−0.2	0.4	−0.1	1.0
2125	1.1	0.8	1.2	1.1	1.1	0.7	1.2	1.1	1.1	0.7	1.2	1.1	1.1	0.7
2126	2.3	1.1	2.2	0.7	2.3	1.0	2.2	0.7	2.3	1.0	2.2	0.7	2.3	1.0
2127	0.4	0.5	0.4	0.7	0.3	0.4	0.4	0.7	0.3	0.4	0.4	0.7	0.4	0.5
2330	2.0	0.7	2.2	1.1	2.0	0.7	2.2	1.1	2.0	0.7	2.2	1.1	2.0	0.7
RSLR SSH	GGM		ARFIMA BIC_tp		ARMA BIC_tp		ARFIMA BIC		ARMA BIC		ARFIMA AIC		ARMA AIC
RSLR SSH	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma
0010	1.3	0.3	0.9	1.6	1.2	0.8	0.9	1.6	1.2	0.8	0.9	1.7	1.2	0.8
0127	1.2	0.4	1.2	1.2	1.2	0.7	1.2	1.2	1.2	0.7	1.2	1.1	1.2	0.8
0158	2.0	0.4	1.6	2.0	1.9	0.7	1.6	2.0	1.9	0.7	1.6	1.9	1.9	0.7
0165	0.7	0.4	0.7	2.1	0.8	0.9	0.7	2.1	0.8	0.9	0.7	1.4	0.8	0.9
0166	1.0	0.4	1.0	1.8	1.0	0.8	1.0	1.8	1.0	0.8	1.0	1.6	1.0	0.8
0245	1.4	0.9	1.5	1.6	1.8	0.7	1.5	1.6	1.8	0.7	1.5	1.7	1.8	0.7
0256	1.9	0.4	1.5	1.8	1.8	0.7	1.5	1.8	1.8	0.7	1.5	1.8	1.8	0.7
0265	0.7	0.4	0.6	1.1	0.7	0.9	0.6	1.1	0.7	0.9	0.6	1.1	0.7	0.9
0377	1.4	0.9	1.5	1.7	1.9	0.5	1.5	1.7	1.9	0.6	1.5	1.7	1.9	0.6
0378	0.8	0.4	0.5	1.5	0.8	0.6	0.5	1.5	0.8	0.8	0.6	1.3	0.8	0.6
0384	1.0	0.4	1.0	1.1	1.0	0.5	1.0	1.1	1.0	0.5	1.0	1.2	1.1	0.7
0385	1.0	0.4	1.0	1.2	1.0	0.8	1.0	1.4	1.0	0.8	1.0	1.2	1.0	0.8
0437	1.3	0.3	0.9	1.6	1.2	0.8	0.9	1.6	1.2	0.8	0.9	1.6	1.2	0.8
0508	1.8	0.4	1.4	1.5	1.8	0.6	1.4	1.5	1.8	0.6	1.4	1.5	1.8	0.5
0527	0.9	0.4	1.0	2.3	1.0	0.9	1.0	2.3	1.0	0.9	0.9	1.4	1.0	0.9
1152	1.0	0.4	1.0	1.4	1.0	0.8	1.0	1.4	1.0	0.8	1.0	1.7	1.0	0.8
1196	0.7	0.3	0.6	1.3	0.7	0.8	0.6	1.3	0.7	0.8	0.6	1.3	0.7	0.8
1242	0.9	0.3	0.9	1.9	0.9	0.8	0.9	1.9	0.9	0.8	0.8	1.4	0.9	0.9
1269	0.8	0.4	0.6	1.4	0.8	0.8	0.6	1.4	0.8	0.8	0.6	1.4	0.8	0.8
1323	0.8	0.4	0.9	2.3	0.9	0.8	0.9	2.2	0.9	0.7	0.9	2.3	0.9	0.8
1325	1.2	0.4	1.2	1.2	1.2	0.7	1.2	1.2	1.2	0.7	1.2	1.1	1.2	0.8
1352	1.5	0.4	1.2	1.4	1.5	0.5	1.2	1.4	1.5	0.5	1.2	1.4	1.5	0.5
1354	0.7	0.4	0.6	1.1	0.7	0.9	0.6	1.1	0.7	0.9	0.6	1.2	0.7	0.9
1394	1.3	0.3	0.8	1.6	1.2	0.7	0.8	1.6	1.2	0.7	0.8	1.6	1.2	0.7
1639	0.9	0.3	0.6	1.5	0.8	0.8	0.6	1.5	0.8	0.8	0.6	1.5	0.8	0.8
1640	0.8	0.4	0.6	1.4	0.8	0.8	0.6	1.4	0.8	0.8	0.6	1.3	0.8	0.7
1799	0.5	0.4	0.4	1.5	0.5	0.9	0.4	1.5	0.5	0.9	0.4	1.5	0.5	0.9
2125	1.2	0.3	0.8	1.3	1.1	0.5	0.8	1.3	1.1	0.5	0.7	1.5	1.1	0.5
2126	1.9	0.4	1.6	1.5	1.9	0.6	1.6	1.5	1.9	0.6	1.6	1.5	1.9	0.6
2127	1.3	1.8	1.2	1.3	1.2	0.8	1.2	1.3	1.2	0.8	1.2	1.2	1.2	0.8
2330	0.6	1.8	0.7	2.0	1.1	0.8	0.7	2.0	1.1	0.8	0.8	1.6	1.1	0.8

Table A4. Absolute SLR estimation derived from PANGA.

TG Site	GGM		ARFIMA BIC_tp		ARMA BIC_tp
TG Site	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma
0010	0.1	0.5	0.1	0.5	0.1	0.5
0127	2.5	0.7	2.4	0.7	2.5	0.7
0158	1.3	0.3	1.3	0.3	1.3	0.3
0165	0.9	0.8	0.9	0.7	0.9	0.8
0166	3.6	0.6	3.6	0.6	3.6	0.6
0245	1.3	0.3	1.3	0.3	1.3	0.3
0256	1.3	0.3	1.3	0.3	1.3	0.3
0265	0.5	0.3	0.5	0.3	0.5	0.3
0377	3.1	0.6	3.1	0.6	3.1	0.6
0378	0.4	0.5	0.4	0.4	0.4	0.5
0384	2.1	0.6	2.1	0.6	2.1	0.6
0385	−0.7	1.2	−0.7	1.2	−0.7	1.2
0437	−0.3	0.5	−0.2	0.5	−0.3	0.5
0508	2.8	0.5	2.8	0.5	2.8	0.5
0527	2.6	1.4	2.5	1.6	2.6	1.0
1152	3.7	0.8	3.9	1.1	3.7	0.8
1196	1.6	0.6	1.5	0.7	1.6	0.6
1242	2.3	0.6	2.2	0.6	2.3	0.6
1269	2.0	0.5	2.0	0.4	2.0	0.5
1323	2.2	0.7	2.2	0.6	2.2	0.7
1325	2.4	0.7	2.4	0.7	2.4	0.7
1352	1.3	0.7	1.1	0.6	1.3	0.6
1354	1.0	0.9	0.9	0.5	1.0	0.9
1394	1.5	0.6	1.5	0.6	1.5	0.6
1639	4.4	1.1	4.5	1.4	4.3	1.0
1640	2.6	1.0	2.6	0.9	2.6	1.0
1799	1.7	1.0	1.7	0.6	1.7	0.9
2125	0.2	0.8	0.3	1.1	0.3	0.8
2126	3.7	1.1	3.6	0.8	3.7	1.0
2127	0.9	0.6	0.9	0.8	0.8	0.5
2330	1.3	0.8	1.5	1.2	1.3	0.8

Table A5. Absolute SLR estimation derived from NMT.

TG Site	GGM		ARFIMA BIC_tp		ARMA BIC_tp
TG Site	Velocity	Sigma	Velocity	Sigma	Velocity	Sigma
0010	−0.6	0.7	−0.6	0.7	−0.6	0.7
0127	0.7	0.8	0.7	0.8	0.7	0.8
0158	0.4	0.4	0.4	0.4	0.4	0.4
0165	−1.0	0.8	−1.0	0.8	−1.0	0.8
0166	1.4	0.7	1.4	0.7	1.4	0.7
0245	0.3	0.3	0.3	0.3	0.3	0.3
0256	0.4	0.4	0.4	0.4	0.4	0.4
0265	−0.3	0.3	−0.3	0.3	−0.3	0.3
0377	2.0	0.6	2.0	0.6	2.0	0.6
0378	−0.2	0.6	−0.2	0.6	−0.2	0.6
0384	0.7	0.6	0.7	0.5	0.7	0.6
0385	−3.5	1.1	−3.5	1.1	−3.5	1.1
0437	−0.9	0.7	−0.9	0.7	−0.9	0.7
0508	2.0	0.6	2.0	0.6	2.0	0.6
0527	1.6	1.3	1.5	1.6	1.6	1.0
1152	2.1	0.9	2.2	1.2	2.1	0.9
1196	0.8	0.5	0.8	0.7	0.8	0.5
1242	1.9	0.8	1.9	0.8	1.9	0.8
1269	1.3	0.6	1.3	0.5	1.3	0.6
1323	1.5	0.5	1.5	0.4	1.5	0.5
1325	0.9	0.7	0.9	0.6	0.9	0.7
1352	0.9	0.8	0.7	0.7	0.9	0.7
1354	0.2	0.9	0.1	0.5	0.2	0.9
1394	1.1	0.7	1.1	0.7	1.1	0.7
1639	3.8	1.1	4.0	1.4	3.8	1.0
1640	2.0	1.0	2.0	0.9	2.0	1.0
1799	1.0	1.0	1.0	0.6	1.0	0.9
2125	−0.1	1.0	0.0	1.2	−0.1	0.9
2126	2.7	1.2	2.6	0.9	2.7	1.1
2127	0.0	0.6	0.1	0.8	0.0	0.5
2330	0.9	0.9	1.1	1.3	0.9	0.9

Figure A1. Zoom of Figure 1.

Figure A2. Absolute SLR (ASLR) estimated using TG and VLM (NMT) and the satellite altimetry mean sea level (SSH) estimated at the same location of the TG.

References

Mazzotti, S.; Lambert, A.; Courtier, N.; Nykolaishen, L.; Dragert, H. Crustal uplift and sea level rise in northern Cascadia from GPS, absolute gravity, and tide gauge data. Geophys. Res. Lett. 2007, 34. [Google Scholar] [CrossRef]
Yousefi, M.; Milne, G.; Li, S.; Wang, K.; Bartholet, A. Constraining Interseismic Deformation of the Cascadia Subduction Zone: New Insights from Estimates of Vertical Land Motion Over Different Timescales. J. Geophys. Res. Solid Earth 2020, 125, e2019JB018248. [Google Scholar] [CrossRef]
Gray, H.J.; Shobe, C.M.; Hobley, D.E.; Tucker, G.E.; Duvall, A.R.; Harbert, S.A.; Owen, L.A. Off-fault deformation rate along the southern San Andreas fault at Mecca Hills, southern California, inferred from landscape modeling of curved drainages. Geology 2018, 46, 59–62. [Google Scholar] [CrossRef]
Clark, J.; Mitrovica, J.X.; Latychev, K. Glacial isostatic adjustment in central Cascadia: Insights from three-dimensional Earth modeling. Geology 2019, 47, 295–298. [Google Scholar] [CrossRef]
Mey, J.; Scherler, D.; Wickert, A.D.; Egholm, D.L.; Tesauro, M.; Schildgen, T.F.; Strecker, M.R. Glacial isostatic uplift of the European Alps. Nat. Commun. 2016, 7, 13382. [Google Scholar] [CrossRef] [PubMed]
Haq, B.U.; Schutter, S.R. A chronology of Paleozoic sea-level changes. Science 2008, 322, 64–68. [Google Scholar] [CrossRef] [PubMed]
Müller, R.D.; Sdrolias, M.; Gaina, C.; Steinberger, B.; Heine, C. Long-term sea-level fluctuations driven by ocean basin dynamics. Science 2008, 319, 1357–1362. [Google Scholar] [CrossRef]
Miller, K.G.; Mountain, G.S.; Wright, J.D.; Browning, J. A 180-million-year record of sea level and ice volume variations from continental margin and deep-sea isotopic records. Oceanography 2011, 24, 40–53. [Google Scholar] [CrossRef]
Young, A.; Flament, N.; Williams, S.E.; Merdith, A.; Cao, X.; Müller, R.D. Long-term Phanerozoic sea level change from solid Earth processes. Earth Planet. Sci. Lett. 2022, 584, 117451. [Google Scholar] [CrossRef]
Lambeck, K. Of Moon and Land, Ice and Strand: Sea Level during Glacial Cycles. Of Moon and Land, Ice and Strand; Leo S. Olschki: Florence, Italy, 2014; pp. 1–86. [Google Scholar]
Sutton, R.T. ESD Ideas: A simple proposal to improve the contribution of IPCC WGI to the assessment and communication of climate change risks. Earth Syst. Dyn. 2018, 9, 1155–1158. [Google Scholar] [CrossRef] [Green Version]
Sutton, R.T.; Hawkins, E. ESD Ideas: Global climate response scenarios for IPCC assessments. Earth Syst. Dyn. 2020, 11, 751–754. [Google Scholar] [CrossRef]
Gregory, J.M.; Church, J.A.; Boer, G.J.; Dixon, K.W.; Flato, G.M.; Jackett, D.R.; Lowe, D.R.; O’Fallerll, S.P.; Roeckner, E.; Stouffer, R.J. Comparison of results from several AOGCMs for global and regional sea-level change 1900–2100. Clim. Dyn. 2001, 18, 225–240. [Google Scholar] [CrossRef]
Holgate, S.J.; Woodworth, P.L. Evidence for enhanced coastal sea level rise during the 1990s. Geophys. Res. Lett. 2004, 31. [Google Scholar] [CrossRef]
Hannah, J.; Bell, R.G. Regional sea level trends in New Zealand. J. Geophys. Res. Ocean. 2012, 117, 1004. [Google Scholar] [CrossRef]
Raj, N.; Gharineiat, Z.; Ahmed, A.A.M.; Stepanyants, Y. Assessment and Prediction of Sea Level Trend in the South Pacific Region. Remote Sens. 2022, 14, 986. [Google Scholar] [CrossRef]
Camuffo, D.; Bertolin, C.; Schenal, P. A novel proxy and the sea level rise in Venice, Italy, from 1350 to 2014. Clim. Chang. 2017, 143, 73–86. [Google Scholar] [CrossRef]
Marcos, M.; Tsimplis, M.N. Coastal sea level trends in Southern Europe. Geophys. J. Int. 2008, 175, 70–82. [Google Scholar] [CrossRef]
Bornman, T.G.; Schmidt, J.; Adams, J.B.; Mfikili, A.N.; Farre, R.E.; Smit, A.J. Relative sea-level rise and the potential for subsidence of the Swartkops Estuary intertidal salt marshes, South Africa. S. Afr. J. Bot. 2016, 107, 91–100. [Google Scholar] [CrossRef]
Qu, Y.; Jevrejeva, S.; Jackson, L.P.; Moore, J.C. Coastal Sea level rise around the China Seas. Glob. Planet. Chang. 2019, 172, 454–463. [Google Scholar] [CrossRef]
Baker, T.F. Absolute sea level measurements, climate change and vertical crustal movements. Glob. Planet. Chang. 1993, 8, 149–159. [Google Scholar] [CrossRef]
Cazenave, A.; Nerem, R.S. Present-day sea level change: Observations and causes. Rev. Geophys. 2004, 42. [Google Scholar] [CrossRef]
Emery, K.O.; Aubrey, D.G. Sea Levels, Land Levels, and Tide Gauges; Springer Science & Business Media: New York, NY, USA, 2012. [Google Scholar]
Montillet, J.P.; Melbourne, T.I.; Szeliga, W.M. GPS vertical land motion corrections to sea-level rise estimates in the Pacific Northwest. J. Geophys. Res. Ocean. 2018, 123, 1196–1212. [Google Scholar] [CrossRef]
Wöppelmann, G.; Marcos, M. Coastal sea level rise in southern Europe and the nonclimate contribution of vertical land motion. J. Geophys. Res. Ocean. 2012, 117. [Google Scholar] [CrossRef]
Wöppelmann, G.; Marcos, M. Vertical land motion as a key to understanding sea level change and variability. Rev. Geophys. 2016, 54, 64–92. [Google Scholar] [CrossRef]
Bitharis, S.; Ampatzidis, D.; Pikridas, C.; Fotiou, A.; Rossikopoulos, D.; Schuh, H. The role of GNSS vertical velocities to correct estimates of sea level rise from tide gauge measurements in Greece. Mar. Geod. 2017, 40, 297–314. [Google Scholar] [CrossRef]
Poitevin, C.; Wöppelmann, G.; Raucoules, D.; Le Cozannet, G.; Marcos, M.; Testut, L. Vertical land motion and relative sea level changes along the coastline of Brest (France) from combined space-borne geodetic methods. Remote Sens. Environ. 2019, 222, 275–285. [Google Scholar] [CrossRef]
Christiansen, B.; Schmith, T.; Thejll, P. A surrogate ensemble study of sea level reconstructions. J. Clim. 2010, 23, 4306–4326. [Google Scholar] [CrossRef]
Press, W.H. Flicker noises in astronomy and elsewhere. Comments Astrophys. 1978, 7, 103–119. [Google Scholar]
Agnew, D.C. The time-domain behavior of power-law noises. Geophys. Res. Lett. 1992, 19, 333–336. [Google Scholar] [CrossRef]
Beran, J. Statistical methods for data with long-range dependence. Stat. Sci. 1992, 7, 404–416. [Google Scholar]
Nerem, R.S.; Chambers, D.P.; Choe, C.; Mitchum, G.T. Estimating mean sea level change from the TOPEX and Jason altimeter missions. Mar. Geod. 2010, 33, 435–446. [Google Scholar] [CrossRef]
Church, J.A.; White, N.J. Sea-level rise from the late 19th to the early 21st century. Surv. Geophys. 2011, 32, 585–602. [Google Scholar] [CrossRef] [Green Version]
Hughes, C.W.; Williams, S.D. The color of sea level: Importance of spatial variations in spectral shape for assessing the significance of trends. J. Geophys. Res. Ocean. 2010, 115. [Google Scholar] [CrossRef]
Montillet, J.P.; Bos, M. Geodetic Time Series Analysis in Earth Sciences; Springer Geophysics; Springer: Berlin/Heidelberg, Germany, 2019. [Google Scholar] [CrossRef]
Armitage, T.W.; Bacon, S.; Ridout, A.L.; Thomas, S.F.; Aksenov, Y.; Wingham, D.J. Arctic sea surface height variability and change from satellite radar altimetry and GRACE, 2003–2014. J. Geophys. Res. Ocean. 2016, 121, 4303–4322. [Google Scholar] [CrossRef]
Le Traon, P.Y.; Reppucci, A.; Alvarez Fanjul, E.; Aouf, L.; Behrens, A.; Belmonte, M.; Bentamy, A.; Bertino, L.; Brando, V.E.; Kreiner, M.B.; et al. From observation to information and users: The Copernicus Marine Service perspective. Front. Marine Sci. 2019, 6, 234. [Google Scholar] [CrossRef]
Liibusk, A.; Kall, T.; Rikka, S.; Uiboupin, R.; Suursaar, Ü.; Tseng, K.H. Validation of copernicus sea level altimetry products in the baltic sea and estonian lakes. Remote Sens. 2020, 12, 4062. [Google Scholar] [CrossRef]
He, X.; Bos, M.S.; Montillet, J.P.; Fernandes, R.M.S. Investigation of the noise properties at low frequencies in long GNSS time series. J. Geod. 2019, 93, 1271–1282. [Google Scholar] [CrossRef]
He, X.; Bos, M.S.; Montillet, J.P.; Fernandes, R.; Melbourne, T.; Jiang, W.; Li, W. Spatial Variations of Stochastic Noise Properties in GPS Time Series. Remote Sens. 2021, 13, 4534. [Google Scholar] [CrossRef]
Santamaría-Gómez, A.; Bouin, M.N.; Collilieux, X.; Wöppelmann, G. Correlated errors in GPS position time series: Implications for velocity estimates. J. Geophys. Res. Solid Earth 2011, 116, B01405. [Google Scholar] [CrossRef]
Klos, A.; Kusche, J.; Fenoglio-Marc, L.; Bos, M.S.; Bogusz, J. Introducing a vertical land motion model for improving estimates of sea level rates derived from tide gauge records affected by earthquakes. GPS Solut. 2019, 23, 1–12. [Google Scholar] [CrossRef]
Hawkins, R.; Husson, L.; Choblet, G.; Bodin, T.; Pfeffer, J. Virtual tide gauges for predicting relative sea level rise. J. Geophys. Res. Solid Earth 2019, 124, 13367–13391. [Google Scholar] [CrossRef]
Varbla, S.; Ågren, J.; Ellmann, A.; Poutanen, M. Treatment of Tide Gauge Time Series and Marine GNSS Measurements for Vertical Land Motion with Relevance to the Implementation of the Baltic Sea Chart Datum 2000. Remote Sens. 2022, 14, 920. [Google Scholar] [CrossRef]
Miller, M.M.; Dragert, H.; Endo, E.; Freymueller, J.T.; Goldfinger, C.; Kelsey, H.M.; Humphreys, E.D.; Johnson, D.J.; McCaffrey, R.; Oldow, J.S.; et al. Precise measurements help gauge Pacific Northwest’s earthquake potential. Eos Trans. Am. Geophys. Union 1998, 79, 269–275. [Google Scholar] [CrossRef]
Altamimi, Z.; Collilieux, X.; Métivier, L. ITRF2008: An improved solution of the international terrestrial reference frame. J. Geod. 2011, 85, 457–473. [Google Scholar] [CrossRef]
Herring, T.A.; Melbourne, T.I.; Murray, M.H.; Floyd, M.A.; Szeliga, W.M.; King, R.W.; Phillips, D.A.; Puskas, C.M.; Santillan, M.; Wang, L. Plate Boundary Observatory and related networks: GPS data analysis methods and geodetic products. Rev. Geophys. 2016, 54, 759–808. [Google Scholar] [CrossRef]
Zumberge, J.F.; Heflin, M.B.; Jefferson, D.C.; Watkins, M.M.; Webb, F.H. Precise point positioning for the efficient and robust analysis of GPS data from large networks. J. Geophys. Res. Solid Earth 1997, 102, 5005–5017. [Google Scholar] [CrossRef]
Bertiger, W.; Desai, S.D.; Haines, B.; Harvey, N.; Moore, A.W.; Owen, S.; Weiss, J.P. Single receiver phase ambiguity resolution with GPS data. J. Geod. 2010, 84, 327–337. [Google Scholar] [CrossRef]
Herring, T.A.; King, R.W.; McClusky, S.C. GAMIT Reference Manual. GPS Analysis at MIT. Release 10.4; Massachussetts Institute Technology: Cambridge, MA, USA, 2010. [Google Scholar]
Herring, T.A.; King, R.W.; Mc Clusky, S.C. GLOBK: Global Kalman filter VLBI and GPS Analysis Program, Release 10.4; Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology: Cambridge, MA, USA, 2010; 91p. [Google Scholar]
Woodworth, P.L.; Player, R. The permanent service for mean sea level: An update to the 21stCentury. J. Coast. Res. 2003, 19, 287–295. [Google Scholar]
Holgate, S.J.; Matthews, A.; Woodworth, P.L.; Rickards, L.J.; Tamisiea, M.E.; Bradshaw, E.; Foden, P.R.; Gordon, K.M.; Jevrejeva, S.; Pugh, J. New data systems and products at the permanent service for mean sea level. J. Coast. Res. 2013, 29, 493–504. [Google Scholar]
Bos, M.S.; Williams, S.D.P.; Araújo, I.B.; Bastos, L. The effect of temporal correlated noise on the sea level rate and acceleration uncertainty. Geophys. J. Int. 2014, 196, 1423–1430. [Google Scholar] [CrossRef]
Blewitt, G.; Lavallée, D. Effect of annual signals on geodetic velocity. J. Geophys. Res. Solid Earth 2002, 107, ETG-9. [Google Scholar] [CrossRef]
Bevis, M.; Brown, A. Trajectory models and reference frames for crustal motion geodesy. J. Geod. 2014, 88, 283–311. [Google Scholar] [CrossRef]
Montillet, J.P.; Williams, S.D.P.; Koulali, A.; Mc Clusky, S.C. Estimation of offsets in GPS time-series and application to the detection of earthquake deformation in the far-field. Geophys. J. Int. 2015, 200, 1207–1221. [Google Scholar] [CrossRef] [Green Version]
Bos, M.S.; Fernandes, R.M.S.; Williams, S.D.P.; Bastos, L. Fast error analysis of continuous GNSS observations with missing data. J. Geod. 2013, 87, 351–360. [Google Scholar] [CrossRef]
Langbein, J.; Bock, Y. High-rate real-time GPS network at Parkfield: Utility for detecting fault slip and seismic displacements. Geophys. Res. Lett. 2004, 31. [Google Scholar] [CrossRef]
Fernandes, R.M.S.; Bos, M.S. Applied Automatic Offset Detection Using HECTOR within EPOS-IP; Time- series, AGU Fall Meeting, G51A-1084; American Geophysical Union: San Francisco, CA, USA, 2016. [Google Scholar]
Rogers, G.; Dragert, H. Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. Science 2003, 300, 1942–1943. [Google Scholar] [CrossRef] [PubMed]
Melbourne, T.I.; Webb, F.H. Slow but not quite silent. Science 2003, 300, 1886–1887. [Google Scholar] [CrossRef]
Viesca, R.C.; Dublanchet, P. The slow slip of viscous faults. J. Geophys. Res. Solid Earth 2019, 124, 4959–4983. [Google Scholar] [CrossRef]
Szeliga, W.; Melbourne, T.; Santillan, M.; Miller, M. GPS constraints on 34 slow slip events within the Cascadia subduction zone, 1997–2005. J. Geophys. Res. Solid Earth 2008, 113. [Google Scholar] [CrossRef]
Miller, M.M.; Melbourne, T.; Johnson, D.J.; Sumner, W.Q. Periodic slow earthquakes from the Cascadia subduction zone. Science 2002, 295, 2423. [Google Scholar] [CrossRef]
Gulick, S.P.; Meltzer, A.S.; Henstock, T.J.; Levander, A. Internal deformation of the southern Gorda plate: Fragmentation of a weak plate near the Mendocino triple junction. Geology 2001, 29, 691–694. [Google Scholar] [CrossRef]
Prescott, W.H.; Lisowski, M.; Savage, J.C. Geodetic measurement of crustal deformation on the San Andreas, Hayward, and Calaveras faults near San Francisco, California. J. Geophys. Res. Solid Earth 1981, 86, 10853–10869. [Google Scholar] [CrossRef]
Galehouse, J.S.; Lienkaemper, J.J. Inferences drawn from two decades of alinement array measurements of creep on faults in the San Francisco Bay region. Bull. Seismol. Soc. Am. 2003, 93, 2415–2433. [Google Scholar] [CrossRef]
Streig, A.R.; Weldon, R.J.; Biasi, G.; Dawson, T.E.; Gavin, D.G.; Guilderson, T.P. New Insights into Paleoseismic Age Models on the Northern San Andreas Fault: Charcoal Inbuilt Ages and Updated Earthquake Correlations. Bull. Seismol. Soc. Am. 2020, 110, 1077–1089. [Google Scholar] [CrossRef]
Melbourne, T.I.; Szeliga, W.M.; Miller, M.M.; Santillan, V.M. Extent and duration of the 2003 Cascadia slow earthquake. Geophys. Res. Lett. 2005, 32. [Google Scholar] [CrossRef]
Keranen, K.M.; Mace, C. Oblique Fault Systems Crossing the Seattle Basin: Seismic and Aeromagnetic Evidence for Additional Shallow Fault Systems in the Central Puget Lowland. J. Geophys. Res. Solid Earth 2011, 117, GP33B-04. [Google Scholar]
Michel, S.; Gualandi, A.; Avouac, J.P. Similar scaling laws for earthquakes and Cascadia slow-slip events. Nature 2019, 574, 522–526. [Google Scholar] [CrossRef]
Hammond, W.C.; Blewitt, G.; Kreemer, C. GPS Imaging of vertical land motion in California and Nevada: Implications for Sierra Nevada uplift. J. Geophys. Res. Solid Earth 2016, 121, 7681–7703. [Google Scholar] [CrossRef]
Langenheim, V.E.; Graymer, R.W.; Jachens, R.C.; McLaughlin, R.J.; Wagner, D.L.; Sweetkind, D.S. Geophysical framework of the northern San Francisco Bay region, California. Geosphere 2010, 6, 594–620. [Google Scholar] [CrossRef]
Rosen, P.; Werner, C.; Fielding, E.; Hensley, S.; Buckley, S.; Vincent, P. Aseismic creep along the San Andreas Fault northwest of Parkfield, CA measured by radar interferometry. Geophys. Res. Lett. 1998, 25, 825–828. [Google Scholar] [CrossRef]
Carpenter, B.M.; Marone, C.; Saffer, D.M. Weakness of the San Andreas Fault revealed by samples from the active fault zone. Nat. Geosci. 2011, 4, 251–254. [Google Scholar] [CrossRef]
Bürgmann, R.; Segall, P.; Lisowski, M.; Svarc, J. Postseismic strain following the 1989 Loma Prieta earthquake from GPS and leveling measurements. J. Geophys. Res. Solid Earth 1997, 102, 4933–4955. [Google Scholar] [CrossRef]
Fumal, T.E. Timing of large earthquakes during the past 500 years along the Santa Cruz mountains segment of the San Andreas fault at Mill Canyon, near Watsonville, California. Bull. Seismol. Soc. Am. 2012, 102, 1099–1119. [Google Scholar] [CrossRef]
Lindsey, E.O.; Fialko, Y.; Bock, Y.; Sandwell, D.T.; Bilham, R. Localized and distributed creep along the southern San Andreas Fault. J. Geophys. Res. Solid Earth 2014, 119, 7909–7922. [Google Scholar] [CrossRef]
Sweet, W.V.; Park, J. From the extreme to the mean: Acceleration and tipping points of coastal inundation from sea level rise. Earths Future 2014, 2, 579–600. [Google Scholar] [CrossRef]
Boruff, B.J.; Emrich, C.; Cutter, S.L. Erosion hazard vulnerability of US coastal counties. J. Coast. Res. 2005, 21, 932–942. [Google Scholar] [CrossRef]
van Westen, R.M.; Dijkstra, H.A. Ocean eddies strongly affect global mean sea-level projections. Sci. Adv. 2021, 7, eabf1674. [Google Scholar] [CrossRef]
Gall, M.; Boruff, B.J.; Cutter, S.L. Assessing flood hazard zones in the absence of digital floodplain maps: Comparison of alternative approaches. Nat. Hazards Rev. 2007, 8, 1–12. [Google Scholar] [CrossRef]
De Biasio, F.; Baldin, G.; Vignudelli, S. Revisiting vertical land motion and sea level trends in the Northeastern Adriatic Sea using satellite altimetry and tide gauge data. J. Mar. Sci. Eng. 2020, 8, 949. [Google Scholar] [CrossRef]
Oelsmann, J.; Passaro, M.; Dettmering, D.; Schwatke, C.; Sánchez, L.; Seitz, F. The zone of influence: Matching sea level variability from coastal altimetry and tide gauges for vertical land motion estimation. Ocean Sci. 2021, 17, 35–57. [Google Scholar] [CrossRef]
Yang, L.; Jin, T.; Gao, X.; Wen, H.; Schöne, T.; Xiao, M.; Huang, H. Sea Level Fusion of Satellite Altimetry and Tide Gauge Data by Deep Learning in the Mediterranean Sea. Remote Sens. 2021, 13, 908. [Google Scholar] [CrossRef]
Church, J.A.; White, N.J. A 20th century acceleration in global sea-level rise. Geophys. Res. Lett. 2006, 33. [Google Scholar] [CrossRef]
Jevrejeva, S.; Moore, J.C.; Grinsted, A.; Matthews, A.P.; Spada, G. Trends and acceleration in global and regional sea levels since 1807. Glob. Planet. Chang. 2014, 113, 11–22. [Google Scholar] [CrossRef] [Green Version]

Figure 1. Spatial distribution over the entire west coast of North America showing the 405 GNSS stations and 31 TG sites used in this study. Note that Figure A1 (in the Appendix A) displays the western coast in more detail (the black lines are fault boundaries).

Figure 2. Histogram of the percentage of data gaps for all GNSS time series analyzed in this study.

Figure 3. Up component of the GNSS daily position time series for the SEAT station with functional model on top (red line) of the observations.

Figure 4. Coastal latitudinal profile of VLM for British Columbia [BC] (Vancouver Island, Canada) and Washington State [WA] (inland Puget Sound, USA) based on the PANGA solution (with the BIC_tp information criterion). The stars are the location of the TGs. The red is the uncertainty associated with the VLM estimates, the green lines are the mean values of the VLM for each region.

Figure 5. VLM estimated for each region (Vancouver Island (Canada), Salish Sea + Inland Puget S. (USA), Olympic Peninsula (USA), South WA (USA) + Oregon plains (USA), Cape Blanco Triple junction + Point Arena (USA), Southern California (USA), whole Pacific coast). The VLM is estimated using either the PANGA or the NMT solution. Note that we have also included the results using the various ICs (AIC, BIC, BIC_tp).

Figure 6. Coastal profile of the vertical land motion for British Columbia [BC] (Vancouver Island, Canada), Washington state [WA] (Olympic Peninsula, south of WA, USA), Oregon [OR] plains (USA), California [CA] (north and south, USA). The stars are the location of the TG. The red is the uncertainty associated with the VLM estimates. The green lines are the mean values of the vertical land motion.

Figure 7. Interpolation of the VLM uncertainties for both NMT and PANGA products based on BIC_tp (the fault zone/line is in red).

Figure 8. Relative Sea Level Rise Rate and associated uncertainties estimated at various TGs along the Pacific coast from Vancouver Island (Canada) to Southern California (USA) using different stochastic noise models. The abbreviations of the TG names refer to Table 1 and Figure 1.

Figure 9. Sea Level Rise Rate (blue) and ASLR (red) for the TG stations in the Pacific coast. Note that the RSLR is produced using the BIC_tp. The VLM estimates used to process the ASLR are based on the BIC_tp and the GGM model. Here, only the results with the PANGA solution are displayed. The flooding risk is added (center to 0 mm/yr) as the vertical bar from low (blue) to high (yellow). The grey line is the GIA estimated from an ensemble of models.

Figure 10. Histogram of the absolute SLR for the various regions along the Pacific coast from Vancouver Island (Canada) to Southern California (USA).

Figure 11. Absolute Sea Level Rise Rate and the mean sea level with satellite altimetry (SSH) estimated at the same location as the TG of PANGA solution.

Table 1. Station information of the processed 31 TG sites, the short name for each TG is indicated between (.).

Station	ID	Lat	Lon	Station	ID	Lat	Lon
SAN FRANCISCO (Sanf)	0010	37.81	−122.47	SOUTH BEACH (S.Be)	1196	44.63	−124.04
SEATTLE (Seat)	0127	47.60	−122.34	BAMFIELD (Bamf)	1242	48.85	−125.13
SAN DIEGO (SanD)	0158	32.71	−117.17	CHARLESTON II (Char)	1269	43.35	−124.32
TOFINO (Tofi)	0165	49.15	−125.92	CAMPBELL RIVER (Camp)	1323	50.02	−125.23
VICTORIA (Vict)	0166	48.42	−123.37	PORT TOWNSEND (P.To)	1325	48.11	−122.76
LOS ANGELES (LosA)	0245	33.72	−118.27	MONTEREY (Mont)	1352	36.61	−121.89
LA JOLLA (LaJo)	0256	32.87	−117.26	TOKE POINT (Toke)	1354	46.71	−123.97
ASTORIA (Asto)	0265	46.21	−123.77	POINT REYES (P. Re)	1394	38.00	−122.98
SANTA MONICA (S.Mo)	0377	34.01	−118.50	N. SPIT (Humb)	1639	40.77	−124.22
CRESCENT CITY (Cres)	0378	41.75	−124.18	PORT ORFORD (P. Or)	1640	42.74	−124.50
FRIDAY HARBOR (Fr.H)	0384	48.55	−123.01	WINTER HARBOUR (WinH)	1799	50.52	−128.03
NEAH BAY (Ne.B)	0385	48.37	−124.61	ARENA COVE (Aren)	2125	38.91	−123.71
ALAMEDA (Alam)	0437	37.77	−122.30	SANTA BARBARA (S.Ba)	2126	34.41	−119.69
PORT SAN LUIS (P.Sa)	0508	35.18	−120.76	PORT ANGELES (P. An)	2127	48.13	−123.44
PORT ALBERNI (P.Al)	0527	49.23	−124.82	PORT CHICAGO (P. Ch)	2330	38.06	−122.04
PATRICIA BAY (Pa.B)	1152	48.65	−123.45	Note: details see www.psmsl.org (accessed on 10 January 2022)

Table 2. Statistics of optimal noise model of the 405 sites (vertical).

Model	PANGA			NMT
Model	AIC	BIC	BIC_tp	AIC	BIC	BIC_tp
FN+RW+WN	0	0	0	0	0	0
FN+WN	17	21	21	11	18	18
GGM+WN	58	55	55	138	137	137
PL+WN	330	329	329	256	250	250
∑	405	405	405	405	405	405

Table 3. The ratio of velocity/velocity uncertainty of the analyzed 405 sites.

Solution	AIC	BIC	BIC_tp
PANGA	4.8	4.8	4.8
NMT	3.8	3.7	3.7

Table 4. Mean sea level estimated from the SRL and ASLR over the whole pacific coast from Vancouver Island (BC, Canada) to Southern California (USA). For the ASLR, we use the ARMA, ARFIMA, and GGM models to estimate the trend at the TGs and the interpolated VLM as described in Section 3. The SLR is estimated from the SSH described in Section 2. Mu (u) and Sig. are the mean value and standard deviation, respectively.

Mean Sea level Pacific Coast (mm/yr)	Model
	ARMA		ARFIMA		GGM
	u	Sigma	u	Sigma	u	Sigma
SLR (SSH)	1.9	1.8	1.8	2.0	1.9	1.7
ASLR (NMT)	0.8	1.7	0.8	1.7	0.8	1.7
ASLR (PANGA)	1.8	1.5	1.8	1.5	1.8	1.5

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

He, X.; Montillet, J.-P.; Fernandes, R.; Melbourne, T.I.; Jiang, W.; Huang, Z. Sea Level Rise Estimation on the Pacific Coast from Southern California to Vancouver Island. Remote Sens. 2022, 14, 4339. https://doi.org/10.3390/rs14174339

AMA Style

He X, Montillet J-P, Fernandes R, Melbourne TI, Jiang W, Huang Z. Sea Level Rise Estimation on the Pacific Coast from Southern California to Vancouver Island. Remote Sensing. 2022; 14(17):4339. https://doi.org/10.3390/rs14174339

Chicago/Turabian Style

He, Xiaoxing, Jean-Philippe Montillet, Rui Fernandes, Timothy I. Melbourne, Weiping Jiang, and Zhengkai Huang. 2022. "Sea Level Rise Estimation on the Pacific Coast from Southern California to Vancouver Island" Remote Sensing 14, no. 17: 4339. https://doi.org/10.3390/rs14174339

APA Style

He, X., Montillet, J. -P., Fernandes, R., Melbourne, T. I., Jiang, W., & Huang, Z. (2022). Sea Level Rise Estimation on the Pacific Coast from Southern California to Vancouver Island. Remote Sensing, 14(17), 4339. https://doi.org/10.3390/rs14174339

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu