An Assessment of Vertical Land Movement to Support Coastal Hazards Planning in Washington State

Abstract

The sea and land change elevation spatially and temporally from a multitude of processes, so it is necessary to constrain the movement of both to evaluate how coastlines will evolve and how those evolving coastlines will impact the natural and built environment over time. We combine land movement observations from global navigation satellite systems (GNSSs), leveling of geodetic monuments, and tide gauge records with a tectonic model of the Cascadia subduction zone to constrain absolute rates of vertical land movement in coastal Washington. We infer rates of vertical land movement in areas lacking direct observations by interpolating high-quality land movement observations and a discretely sampled interseismic locking model. Here we present a model of absolute vertical land movement that is combined with sea level rise estimates to inform local relative sea level projections on a community-scale. The most rapid vertical uplift (~3.5 mm/year) of the land is found across the northwest Olympic Peninsula, which currently outpaces sea level rise. Conversely, some areas, including a stretch of the northern Pacific Ocean coast from La Push to Kalaloch and the southern Puget Sound, are found to be subsiding at 0.5–1.0 mm/year, exacerbating the rate of relative sea level rise and thereby increasing the vulnerability of coastal communities.

1. Introduction

Coastlines support diverse ecosystems and a myriad of anthropogenic interests like communities, historical sites, and infrastructure. The safety and prosperity of coastal communities and the natural and built environment depend on our ability to plan for, predict, and adapt to changing conditions like rising mean sea level. Mean sea level change over time-frames of decades to centuries is primarily driven by thermal expansion of the ocean and mass loss from glaciers and ice sheets [1], but it can be modified locally by the vertical velocity of the land along the coast, leading to relative sea level change. Therefore, it is necessary to pair estimates of absolute sea level change with vertical land movement estimates in order to constrain relative sea level change and better assess the vulnerability of coastal resources. Vertical land movement (VLM) refers to the average long-term displacement of land surface in a geocentric reference frame that results from reorganization of the materials comprising Earth’s crust or from isostatic adjustment of the crust in response to a changing load [2]. Both anthropogenic processes (e.g., groundwater and hydrocarbon extraction, loading from water reservoirs, and anthropogenic compaction) and local natural processes (e.g., sediment compaction, tectonic deformation, and glacial isostatic adjustment) can cause vertical land movement at a variety of temporal (days to centuries) and spatial scales (millimeters to thousands of kilometers) [3,4]. Estimates of VLM are used to constrain models of these causal processes and assessments of earthquake, tsunami, and storm hazards [5,6]. The focus of this study is to constrain rates of land motion and relative sea level change to inform coastal hazards planning utilizing a method that is adaptable to other regions.

In coastal Washington state, on the Pacific coast of central North America, variable rates of coastal vertical land movement have been estimated by a variety of researchers, ranging between 4 mm/year uplift and −1 mm/year subsidence (Table S1). Regional patterns of VLM in coastal Washington are largely attributed to tectonic processes [7,8] from the Cascadia subduction zone (CSZ) to the west, and glacio-isostatic adjustment from deglaciation [9,10]. The CSZ hosts great earthquakes on a 500–530 year average recurrence interval [11] and is the predominant geodynamic process in this region generating crustal strain accumulation, in addition to clockwise forearc rotation [12]. Washington also hosts smaller crustal faults like the Seattle fault [13], as well as a multitude of volcanoes along the Cascadia volcanic arc. The geology of the coastal region is dominated by young oceanic basalt and sedimentary rock spanning the study area [14,15]. The area of interest in this study spans from the Pacific Ocean to the Cascadia volcanic arc, specifically along the coastline of Washington state.

Previous investigators have attempted to estimate rates of vertical land movement in the Pacific Northwest of the USA generally, and specifically in Washington state, motivated either by the desire to incorporate those movements into sea level assessments, or to inform seismic hazard assessments [16,17,18,19]. Most investigations have utilized one of three different methods for measuring long-term rates of vertical land movement:

(1) Leveling, which utilizes repeat measurements of geodetic monuments to derive a relative estimate of vertical land movement (i.e., vertical movement of one survey benchmark relative to another; see Burgette et al. [7]. Many early studies of vertical land movement referenced these relative rates to a geodetic reference frame by linking leveling results to estimates of vertical land movement derived from tide gauges.

(2) “Water-level differencing” methods, which difference water levels measured at tide gauges located close to each other to derive an estimate of vertical land movement at one tide gauge site relative to another. Many early studies of vertical land movement converted the relative rates derived using this approach to absolute rates using an assumed global average rate of sea level change (Table S1).

(3) Continuous global navigation satellite systems (CGNSS), which utilizes regularly-sampled measurements of position from fixed GNSS stations to derive an estimate of long-term vertical land movement rates in an absolute (i.e., geocentric) reference frame (i.e., movements relative to a center of mass of the Earth).

Mote et al. [20] conducted a review and discussion of investigations to characterize vertical land movement in Washington state. They compared three previous studies [21,22,23], all based primarily on leveling, to a preliminary CGNSS analysis. They found general agreement in regional patterns between these previous investigations, but discrepancies in the estimated rates of vertical land movement in some locations. More recent studies [16,17,20,24,25], summarized in Table S1, attempted to estimate rates and patterns of vertical land movement with the goal of resolving differences between relative sea level measured at tide gauges against vertical land movement rates at co-located CGNSS stations. These investigations provide a way to assess trends in regional absolute sea level, and compare the regional sea level trends against globally averaged rates. In general, these investigations rely on CGNSS datasets, typically from just one of several data providers, and also typically focus on a very limited number of locations in Washington state that are co-located with tide gauges. The most recent studies have also utilized GNSS glacial isostatic adjustment (GIA) corrections at tide gauges [18] and multiple CGNSS data providers [19].

In terms of resolving spatial variations in vertical land movement, most previous analyses suffered from either limited data, relying on estimates from just one data source or methodology, and/or a lack of estimates made directly on the coast. In addition, Snay et al. [16] point out discrepancies between different providers of continuous CGNSS-derived vertical land movement estimates (i.e., Table S1), which are also a consequence of the application of different methodologies, differing time-series lengths, or different analysis approaches. Finally, previous sea level rise assessments focused on Washington state [20,25] incorporated estimates of VLM into relative sea level projections, but without enough spatial resolution to resolve meaningful along-coast variations on a community scale. Partly as a consequence of these short-comings, local coastal managers in Washington state reported that they didn’t have adequate information to plan adequately for sea level rise [26]. Globally, studies have identified centimeter-scale VLM signals combining leveling and spaceborne SAR [27], however coastal Washington is heavily vegetated and regional mm/year scale VLM is outside of the observable range of standard SAR processing in our study region.

In this study, we present an approach to assimilate VLM estimates and their uncertainties, derived from multiple data types, in order to maximize the density of estimates on the coastal landscape and better account for the discrepancies between them. We also account for the primary VLM signal in areas without measurements by interpolating multiple data types and modeled vertical land velocities from interseismic tectonic locking on the CSZ, which dominates vertical velocities in instrumented areas. Our work reveals variable rates of vertical land movement in coastal Washington that may exacerbate or outpace sea level change depending on the location. We combine our model of vertical land movement with absolute sea level change projections and produce probabilistic relative sea level change projections, under two greenhouse gas scenarios and assuming no subduction zone earthquake, along the entirety of coastal Washington (first presented by Miller et al. [28]). The approach presented here is adaptable to any study area for which measurements and models of VLM exist, and represents an improvement over current methods of VLM estimation, resulting in dense and robust estimates of VLM that benefit coastal hazards planning.

2. Materials and Methods

2.1. Data

2.1.1. Continuous GPS

GNSS observations provide estimates of vertical land motion rates in a geocentric reference frame. While these estimates are relatively precise (with uncertainties on the order of several mm/year), solutions can vary depending on processing assumptions, such as regional filtering methods or how atmospheric water vapor is modeled. The observed differences in vertical rate estimates and their uncertainties across different CGNSS datasets [16] motivated us to develop an assimilative approach for deriving vertical land movement rates from CGNSS data from multiple GNSS datasets (Table S2). We based our best estimates and associated uncertainties for CGNSS sites on the distribution in estimated rates and uncertainties across datasets provided by five different providers: Nevada Geodetic Laboratory [29], NASA JPL [30], Scripps Orbit and Permanent Array Center [31], Pacific Northwest Geodetic Array Geodesy Laboratory [32], and UNAVCO [33]. All of the GNSS data are in the ITRF2008 reference frame.

Data were downloaded for all available stations in the Pacific Northwest (inclusive of western Oregon, western Washington and southern British Columbia) from the five different CGNSS data processing groups. To address the disparity in rates and uncertainties (Figure S1) across CGNSS-derived estimates we first classified CGNSS stations into high/low quality categories based on the length of the time-series over which vertical velocity solutions were estimated from a visual examination of the relationship between uncertainty and time-series length (Figure S1). Records exceeding nine years were initiallyclassified as “high-quality”. Records shorter than nine years were classified as “high-quality” only if at least three of the vertical rate estimates across the different processing groups agreed to within ±1 mm/year. Stations designated as “low-quality” via this filtering step were excluded from further analysis. This quality assessment procedure classified 117 out of 188 total CGNSS stations in our area of interest (Figure 1) as “high-quality”.

The best estimate of the vertical rate at a CGNSS site was designed to take into account differences in solutions between the five different datasets, and was arrived at through a multi-step process. First, an outlier detection algorithm was applied to the “high-quality” CGNSS sites with more than three estimates from different processing groups, in order to separate out stations with longer time-series but for which one of the processing group’s estimates was dramatically different from the others. For this assessment we applied a median absolute deviation approach to identify outliers, with an outlier cut-off value of 6. Therefore, any value with a median absolute deviation of more than 6 times the median value was eliminated from the calculation of a mean velocity across data processing centers.

Next, we calculated a weighted average of the individual vertical rate estimates, where the weighting was scaled by an adjusted uncertainty in the station vertical velocity:

$$\overline{\mathrm{X}}={{\displaystyle \sum }}_{i = 1}\frac{{\mathrm{X}}_i}{{\mathsf{\sigma }}_i^2} /{{\displaystyle \sum }}_i\frac{1}{{\mathsf{\sigma }}_i^2},$$

(1)

where x_i are the reported velocities from the different processing groups and σ_i are the adjustment uncertainties for each observation. The weighted averaging approach has the effect of making those vertical velocities with smaller adjusted uncertainties carry more weight in the vertical velocity averaging process. The magnitude of the uncertainty estimates for CGNSS-based vertical land movement observations were not consistent across data providers, with some providing a relatively narrow uncertainty range and others estimating a wide uncertainty range (Figure S1a, Table S2). However, the relative scale of the uncertainties was generally consistent within and across datasets and was assumed to be indicative of the relative scale of uncertainty of the estimates between stations. In order to utilize the weighted averaging approach described in Equation (1) the different uncertainty estimates for each station needed to be adjusted. As a result, we developed an adjustment process to rescale uncertainties and facilitate comparisons across datasets. We assume that the uncertainty in the vertical rate estimate should asymptote to a stable value as more observations become available (i.e., due to a longer observational record). We therefore corrected the uncertainty estimates by: (1) Calculating the mean uncertainty for each of the five CGNSS data sets for all high-quality stations. (2) Estimating a scaling factor for the reported uncertainties for each processing center such that the mean uncertainty for all high-quality stations is consistent between processing centers. This ensures that the variability in uncertainty is preserved between stations. (3) Multiplying all CGNSS vertical uncertainties by the appropriate scaling factor assigned to each processing center (Table S2). To estimate the uncertainty of our final, average vertical velocity rates for each station we calculated the standard deviation of the vertical rate estimates for a station, as calculated by the different processing groups (Figure 1). This approach was more conservative than propagating the error associated with each rate.

2.1.2. Tide Gauges

Water-level differencing, in our case differencing monthly-averaged water levels from tide gauges, provides an estimate of the relative difference in vertical land velocity rates between tide gauges. There are two approaches: “single-differencing”, which assumes that the stations are close enough to each other that they experience the same variations in sea level; and “double-differencing”, which does not assume that the stations share an absolute sea-level trend, but rather attempts to directly estimate the absolute sea level pattern at a tide gauge using altimetry [34]. We conducted a single-differencing analysis that employed the following steps: (1) Obtain monthly average water level data from the Permanent Service for Mean Sea Level for a set of 23 stations in Oregon, Washington and British Columbia. (2) Choose a set of “reference stations” (Table S3) that were co-located with a “high quality” CGNSS station (within 5km, though 3 of the 4 reference stations were located less than 500 m apart), and included a long (>20 years) water level record. (3) Difference each station’s monthly sea level against those from each of the four reference stations. (4) Remove the mean seasonal signal in the resulting differenced time-series, by subtracting out the long-term average for each calendar month. (5) Estimate the trend and standard error in the linear regression, accounting for auto-correlation (using the variance inflation factor approach described in Zervas [35]. (6) Adjust the estimated relative vertical velocity to an absolute rate by adding it to the best estimate vertical rate derived from the high quality CGNSS station co-located with each “reference” water level station. The uncertainties assigned to the rates derived using this approach (Figure 2) were estimated by propagating the uncertainty assigned to the co-located CGNSS station.

2.1.3. Leveling

Leveling is a surveying method that is used to measure the relative height difference between two nearby benchmarks. By repeating a survey of the benchmarks at a later time and differencing against an earlier data set, it is possible to measure the rate of change in the elevation difference between benchmarks. This method provides relative changes in height between adjacent stations (i.e., tilt) and must be tied to a nearby CGNSS or tide gauge station in order to be placed in a reference frame. We incorporate vertical rates from leveling in our analysis using the methodologies and data described in Burgette et al. [7], Krogstad et al. [36], Mitchell et al. [22], and Schmidt et al. [37], which utilize historical leveling data dating to the 1920′s (Figure 3), by adjusting earlier tide gauge-based reference frames to ours, based on CGNSS as described above. We also added two new lines from Westport and Tacoma east. These relative rates were converted to absolute rates of vertical land movement by attaching leveling lines to nearby high-quality CGNSS stations. Relative vertical rates for a series of benchmarks along a leveling line (typically located along transportation corridors) were adjusted to best fit the vertical land movement estimates of one or more nearby CGNSS stations or tide gauges. The uncertainty of the CGNSS station is assigned to the nearest leveling point, then uncertainty is propagated from that point using the measurement uncertainties of successive leveling observations. We also recalculated the uncertainties for the north-south line on the east side of the Olympic Peninsula that were too small in previous work that we incorporated.

2.2. Tectonic Model

Within our study area we found a total of 519 VLM rate estimates from CGNSS, leveling, and tide gauges (Figure 4). However, there are large stretches of our study area for which vertical land movement observations do not exist (Figure 4). The limited data that exist near the coast and the spatially-dense leveling lines (roughly perpendicular to the CSZ) suggest large gradients in vertical land movement over short distances, which could lead to a misrepresentation of rates of vertical land movement on the shoreline in areas where most observations are collected farther inland. Specifically, at locations on the north and central coast near La Push and Kalaloch, Washington, observations point to subsidence on the shoreline, and rapid uplift only a short distance inland. Further, an interpolation of the combined vertical land velocity observations yields unphysical edge effects since all data were acquired on land and there are no seafloor observations to constrain the western edge of the interpolation. To first order, the large-scale pattern of deformation revealed by vertical land movement observations in our area of interest can be explained by interseismic locking of the megathrust on the Cascadia subduction zone, and large gradients in near-coast deformation are likely [38].

We incorporated a discretely sampled velocity field derived from a forward model of a locked fault slipping at depth in the Cascadia subduction zone in our analysis to place geophysical constraints on the expected vertical signal in regions without direct measurements. The interseismic locking model uses the modeling framework of Burgette et al. [7], which was solely focused on modeling leveling data in Oregon, but expanded to Washington state and constrained by both vertical leveling and horizontal GPS velocities by Schmidt et al. [37]. This locking model was compared to and combined with other similar locking models to constrain seismic hazards in the Pacific Northwest for the United States National Seismic Hazard Maps [5] in Frankel and Petersen [8]. The interseismic locking model assumes that the slip deficit rate is equal to the plate convergence rate (i.e., fully locked) near the trench and decreases exponentially in the transition zone. The locking model is also assumed to vary smoothly along strike. The depth of the lower edge of the fully locked zone and lower edge of the transition zone as a function of latitude are optimized through a grid search. The resulting locking model is largely similar to other locking models for Cascadia [5], where fault locking is primarily confined to the megathrust at depths shallower than 20 km (coinciding with sections of the megathrust west of the coastline). In using the locking model for this paper, we use the distribution of the slip-rate deficit on the megathrust to forward predict uplift rate across our focus area. To better match the new refined data set, we have applied a uniform vertical shift to account for our different reference frame.

The locking model is intended to physically model the crustal strain associated with the convergence of the Juan de Fuca and North American plates, which is the dominant signal in the geodetic data. This strain is elastically stored in the crust and will eventually be released in a future megathrust earthquake on the subduction zone, thus our estimates of vertical land velocity are not valid after the next megathrust earthquake. Our tectonic model does not model other processes, like compaction of sediments, erosion, or the movement of shallow crustal faults that may result in the vertical movement of the ground. Vertical rate estimates were extracted from the locking model at a grid spacing of 0.2 degree (16–22 km). These were combined with our other vertical rate estimates for the purposes of interpolating a best-estimate vertical velocity field. The modeled vertical rates were arbitrarily assigned a relatively large uncertainty of ± 2 mm/year to minimize the influence of the model in areas with ample observations. The locking model is not used to replace or modify the observations, but rather is used to more accurately interpolate vertical velocity estimates at locations between the observations, particularly where large gradients exist.

2.3. Interpolation

We interpolate between all discrete observations and gridded tectonic model velocity field points to represent the vertical land velocity as a surface that is a function of the longitude and latitude, or a height field with z = f(longitude, latitude). Since we model the vertical velocity as a height field, the field may not contain vertical slopes or overhangs, neither of which are realistic representations of land movement at the scale of interest (~10 km). We first generate a scattered interpolant from the combined observed and modeled vertical velocities using Delaunay triangulation with bivariate linear interpolation [39]. Triangulation of the scattered vertical velocities in the longitude-latitude plane generates a continuous surface composed of triangular pieces joined at edges, or a triangulated irregular network, over the area of interest. The triangulated irregular network represents a piecewise interpolation scheme in which bivariate linear interpolation is applied to each triangle. We then query the interpolant at a grid spacing of 0.1 degree to yield an interpolated vertical land velocity height field (Figure 5a; Data S3).

To include a measure of uncertainty for interpolated vertical velocities we employ bootstrapping. At each location with observed or modeled vertical velocities we specify a normal distribution defined by the vertical velocity as the distribution mean and the standard deviation. We then randomly sample each distribution and generate an interpolated vertical velocity height field, repeating the process 10,000 times. The reported interpolation uncertainty is the standard deviation, with N-1 normalization, of the 10,000 interpolation iterations at each grid point, thus representing the variation in interpolated vertical velocity given hypothetical variations in the data. In essence, the uncertainty presented in Figure 5b is a measure of the stability of the interpolation results.

3. Results

We model the vertical velocity field of coastal Washington as an interpolated height field to produce the most spatially comprehensive estimate of vertical land motion for coastal Washington to date. Previous investigators (Table S1) have found variations in coastal vertical land movement ranging from ~1.0 mm/year of subsidence (in central Puget Sound) to approximately 3–4 mm/year of uplift (on the northwest tip of the Olympic Peninsula). Overall, we find the same general patterns suggested by previous investigators, and resolve many of the finer-scale variations between tide stations and leveling lines. In addition, we find the most rapid rates of coastal uplift associated with the northwest tip of the Olympic Peninsula (Figure 5). Our results also suggest:

• A strong west-to-east gradient in rates of vertical land movement across our area of interest that is particularly pronounced along a transect following the Strait of Juan de Fuca (Figure 5a), and that we attribute to subduction zone processes.
• Low rates of subsidence (<1 mm/year) in Puget Sound, but with a south-to-north gradient showing higher rates of subsidence to the south and lower rates to the north, consistent with the pattern expected for glacio-isostatic adjustment.
• Rapid uplift along the southwest Washington coast, but with a large uncertainty in estimated rates, including Westport and the Long Beach Peninsula.

We also find a broad area of coastal subsidence along the northern Pacific Ocean coast of Washington state (Figure 6). Previous studies generally do not make spatially continuous VLM estimates throughout this area, though some studies did find evidence of subsidence in the same area [23,34]. Other studies estimated uplift (e.g., [14]) along this section of coast, or observations in this area were not considered (e.g., [11]). Figure 6 shows the northwest region of study area where we find subsidence along coastal Washington, indicated by the colored grid that constitutes the interpolated vertical velocity surface, in addition to VLM observations and the gridded tectonic model that inform the interpolation. Observations of VLM from leveling, tide gauges, and CGNSS stations in this region show spatial vertical velocity gradients (0.2 mm/km) that are steeper than those present in the tectonic locking model (0.1 mm/km; Figure 6). Figure 7a shows the interpolant queried at a sub-kilometer scale along the coastline to depict smooth variations in coastal VLM, and Figure 7b shows the associated uncertainty.

Despite the density of observations utilized in this assessment, areas of large uncertainty remain, either due to uncertainties in the observations themselves (Figure S2), large gradients in the rate of vertical land movement, or a lack of observations in certain parts of the landscape. We quantify the uncertainty associated with our interpolated model of vertical land motion (Figure 5b). It is important to note, there are also areas of uplift or subsidence driven by processes other than interseismic strain accumulation and at scales smaller than the resolution of this study are not well-represented in our model or uncertainty analysis. Areas containing VLM observations with high uncertainties (observations with black outlines in Figure 4) correspond to the areas in our model with the highest uncertainty (dark blue colors in Figure 5b). The focus of this study is to provide robust estimates of vertical land movement along the coast, so our estimated rates of VLM in inland Washington, particularly near volcanoes and the eastern boundary of our model, do not fully characterize contemporary crustal movements considering interseismic strain accumulation decreases with distance from the subduction zone.

4. Discussion

Our estimates of VLM and the associated uncertainties, when combined with sea level projections, provide insights about current and future relative sea level changes for coastal communities throughout Washington that can supply local coastal managers with the information necessary to plan adequately for sea level rise (Figure 7a,b). An assumption of this study, supported by over 70 years of tide gauge and repeated leveling observations, is that the observed and modeled rates of vertical land movement are constant over the interseismic period, so we discuss the implications of a Cascadia subduction zone earthquake on vertical land motion in the region. Additionally, we further explore spatial patterns in the resulting interpolated vertical velocity field and identify a convolved north-south vertical velocity gradient, likely from glacial isostatic adjustment, that will continue whether there is a Cascadia earthquake or not. To conclude this study, we explore the effect of the assumed GIA signal on an elastic dislocation model of Cascadia subduction zone locking, and thus could impact earthquake hazard level estimates.

4.1. Relative Sea Level Change

We combine our estimate of interpolated vertical land movement with absolute sea level change projections to produce probabilistic relative sea level change projections for 171 locations throughout coastal Washington under two greenhouse gas scenarios, first presented by Miller et al. [28]. We utilize a single absolute sea level rise projection for all of Washington state, as the regional variations in absolute sea level change are not well resolved over this scale [28]. The low and high greenhouse gas scenarios (RCP 4.5 and 8.5 [40]) reflect different emissions trajectories that depend on human behavior over decades; thus, we do not speculate which scenario is more likely. The projections presented in Figure 8 incorporate uncertainties from the models and observations incorporated into the assessment, which includes the range of expected mass loss from the Antarctica and Greenland ice sheets. Our projections are presented as the likelihood that sea level change will meet or exceed a particular offset from contemporary sea level, or the probability of exceedance. We utilize an approach for estimating absolute sea level change adapted from Kopp et al. [41] for seven Washington state tide gauges, fully described by Miller et al. [28].

Because interseismic strain accumulation from tectonic locking on the CSZ is the primary VLM signal in the coastal Pacific Northwest on a regional scale, relative sea level rise patterns mirror the opposite of the tectonic uplift signal. Areas like Neah Bay that experience rapid vertical uplift (Figure 7) will, in general, experience less relative sea level rise, and sea level drop is present in cases where vertical land velocities outpace rates of absolute sea level rise. Conversely, regions experiencing land subsidence, like Tacoma (Figure 7), will have rates of relative sea level rise that are greater than the rate of absolute sea level rise. The effect of large uncertainties in vertical land velocities can be seen in the projections for sites like Long Beach peninsula, resulting in a large difference between probability of exceedance projections. Relative sea level projections for 171 locations along the Washington coast are available and we refer readers to Miller et al. [28] for additional details of our relative sea level rise assessment.

4.2. Cascadia Earthquake Scenario

We assume constant rates of vertical land movement during an interseismic period, and thus constant strain accumulation. This assumption is reasonable given that: (1) we are late in the interseismic period when deformation rates are expected to be constant in time [38], and (2) nearly constant rates of VLM are observed for the ~70 year tide gauge and leveling record. However, the uplift proximal to the coast indicates that strain energy is accumulating in the crust and will be released in the next megathrust earthquake. A large earthquake would result in coseismic land displacements that void our assumption of a constant deformation rate and produce temporally and spatially variable rates of vertical land movement. Additionally, rapid and spatially variable post-seismic deformation follows large earthquakes, resulting from viscoelastic stress transfer that continues decades after the earthquake. The Tohoku earthquake that struck central Japan on March 11 2011, for example, resulted in coastal subsidence of more than 1 m in places along the east coast of Japan [42]. Localized uplift of the coast is also possible depending on the geometry of the subduction zone and proximity of the coast to the zone of locking. Figure 9 provides simplified examples of the effect of a megathrust earthquake on relative sea level for three sites from 2020 to 2120 with and without a hypothetical megathrust earthquake at the year 2070, assuming that our estimated VLM rate is entirely due to subduction processes and full recovery of the VLM rate summed over the hypothetical interseismic period (1700–2070). The true coseismic VLM depends on the type and timing of rupture (Table S4), but previous studies that focus on constraining the complexities of coseismic displacements report a range of coseismic subsidence consistent with our first-order estimate for areas experiencing uplift [6,43,44,45]. Despite the simplifications and assumptions inherent to the hypothetical coseismic land motions depicted in Figure 9, the region of subsidence we find from La Push to Kalaloch may produce coseismic uplift rather than the coseismic subsidence expected in regions experiencing interseismic uplift, depending on the complexities of subduction processes. Thus, Figure 9 serves to reiterate the importance of constraining VLM and the driving processes to inform hazard assessments.

4.3. North-South VLM Gradient

We find two dominant signals in the data that are most evident in our study area: an east-west gradient consistent with locking along the Cascadia subduction zone (CSZ), and a smaller magnitude south-north gradient that is most apparent east of the region affected by subduction zone locking (Figure S3a). The two gradients are superimposed in the coastal region, although VLM data are routinely used to inform subduction zone models that assume the data are absent of independent sources of motion. To isolate the north-south vertical velocity gradient, we utilize the gridding of our interpolated surface and find the mean value for all rows of 0.1 degree by 0.1 degree cells east of −123.0° longitude, then we subtract the mean for each row from all 0.1 degree by 0.1 degree cells in that row. The longitude-averaged north-south vertical velocity ranges from 0.6 mm/year in the north (49.0° latitude) to −0.7 mm/year in the south (45.5° latitude). This simple estimate of the observed north-south gradient assumes that the vertical velocity gradient is strictly north-south trending, but is intended to provide only a first-order approximation of the signal.

To assess the contribution of this vertical velocity gradient, we generate simple elastic dislocation models (as in Schmidt et al. [37]) for CSZ locking with and without the north-south gradient by subtracting the observed gradient east of the Puget Sound from all of Washington, including the CSZ locking-dominated gradient in coastal Washington, west of the Cascades. We compare the predicted horizontal strain from each model to the observed horizontal strain, measured from horizontal GNSS motion. Correcting for the north-south gradient in our model reduces the RMS error by 2% and better predicts the observed strain in Washington. Additionally, the correction extends the modeled extent of tectonic locking farther downdip near La Push (Figure S3b). Since the observed vertical and horizontal strain best fit our model with the north-south regional uplift gradient removed, we hypothesize that this gradient is superimposed on the coastal region, and likely unrelated to the CSZ. We speculate that the observed north-south gradient is GIA, consistent with previous studies [4,9]; however, other mechanisms are possible, including a 3D viscoelastic response to past CSZ rupture, rotation of the forearc, and local subsidence associated with the broader Puget Sound forearc basin region. Regardless of the source, the convolved signal implies that multiple overlapping processes are driving vertical land movement in Washington state.

5. Conclusions

We generate the most spatially comprehensive estimate of VLM to date for coastal Washington and find moderate rates of subsidence throughout the Puget Sound and near La Push, separated by a broad region of rapid uplift that extends south to Oregon. Combining VLM with forecast sea level rise provides relative sea level change for coastal Washington. A north-south vertical velocity gradient, consistent with GIA, is apparent in VLM observations, and correction of an elastic dislocation model for this signal yields lower errors in the predicted horizontal strain in Washington when compared to observations. The dominant VLM pattern in this region is consistent with interseismic strain accumulation on the CSZ.

Ideally, future research will combine new techniques and spatially-comprehensive data sets to produce dense and robust estimates of VLM. Previous VLM assessments predominantly infer uplift in the region spanning La Push to Kalaloch where we find subsidence [20,25]. Conflicting VLM results pose challenges for hazards planning in coastal communities; thus, future studies that further constrain vertical uplift rates in areas with few direct observations are critical to understanding VLM and its sources throughout Washington. Since relative sea level projections utilize annually compounding rates of VLM, misestimation of vertical rates may lead to large errors in projected relative sea level and misuse of community resources leading to potential ramifications for community vulnerability. Further, VLM estimates inform coseismic deformation models that are critical for accurate estimation of tsunami hazards, requiring robust, spatially-comprehensive, and repeatable VLM estimates.