Sea Level Rise in New Zealand: The Effect of Vertical Land Motion on Century‐Long Tide Gauge Records in a Tectonically Active Region

Historically tide gauge (TG) data have been used to estimate global sea level rise. Critical to the analysis of TG records is the assumption that the TG sites are stable and not affected by vertical land motion (VLM). We analyze century‐long TG records from New Zealand that have been affected by VLM due to both major and transient earthquake events at a regional level as well as local instabilities. Using combined GPS and precise leveling, we estimate relative VLM between the GPS and TG of up to 1 mm/year. Based on 15–20 years of GPS data, the effect of seismic activity and slow slip events has uplifted sites by up to 0.8 mm/year on average in Wellington and Dunedin. Updated estimates of long‐term relative sea level (RSL) at four New Zealand TGs, as well as an estimate of RSL at a new fifth TG (New Plymouth), are determined using data through 2013—an additional 13 years compared to the previous study. The VLM corrected RSL rates gives our best estimate absolute sea level of +1.45 ± 0.28 mm/year (1891−2013). It is important that absolute sea level derived from RSL rates include realistic estimates of VLM, especially at TG sites that are close to plate boundaries, located in seismically active regions, affected by glacial isostatic adjustment, land ice loss, or gas/oil/water extraction.


Introduction
The long-term trend in global mean sea level (GMSL) is recognized as an important indicator of climate change. A rise in sea levels over the long term has the potential to adversely impact societies-especially those in vulnerable low lying coastal areas where population density tends to be high. For this reason, the change in GMSL has received considerable attention during the past several years, gaining particular prominence in publications such as those of the Intergovernmental Panel on Climate Change (e.g., IPCC, 2007IPCC, , 2013. For example, the Church and White (2006) analysis, based on long record tide gauge (TG) records, suggests a twentieth century rate of 1.7 mm/year, while Church and White (2011) show a small acceleration in GMSL.
Satellite altimetry is becoming increasingly useful in estimating GMSL (e.g., Cazenave & Llovel, 2010;Nerem et al., 2010;Church & White, 2011). However, the rate of GMSL derived from altimetry records relies upon a relatively short time series of data that is affected by short-term climate variability (e.g., El Niño-Southern Oscillation [ENSO]; Cazenave et al., 2014) and has been shown by Watson et al. (2015) to be affected significantly by the altimeter calibration used. TG records, on the other hand, contain a long time series of historical information, dating in some cases back to the nineteenth century, and can provide calibration data for the altimetry. Their importance is especially true in the Southern Hemisphere where long-term (~100 years) TG records are sparse (Woodworth & Player, 2003).
At interannual and decadal time scales TG records typically show variations in annual means on the order of 5-10 cm [e.g., Holgate & Woodworth, 2004;Hannah & Bell, 2012] and a possible 60-year oscillation of up to 3 cm (Chambers et al., 2012). Although to estimate the long-term trend in relative sea level (RSL), five to six decades of reliable and essentially complete TG records are typically required (e.g., Douglas, 1997), techniques exist that can use sparse data sets (e.g., Hannah & Bell, 2012). A TG measures the RSL trend or the change in sea level relative to the adjacent land. Independent estimates of the vertical land motion (VLM) are required to estimate the absolute sea level (ASL). Such motion may occur due to tectonic movements, ice mass loading, glacial isostatic adjustment (GIA), or local site instabilities due to, for example, compaction, sedimentation, or water, oil, or gas extraction (Carter et al., 1989).
The stability of a TG site is a critical issue that requires verification (Bevis et al., 2002) and needs to be considered at both local and regional scales. Historically, many TGs are located in shipping ports where harbor infrastructure is often developed on reclaimed or modified land. Occasionally, port operations may require a TG to be relocated. Local TG stability is traditionally monitored by establishing a local benchmark network that is connected to the TG using spirit leveling. While at least one benchmark is assumed to be stable, it is preferable that the stability of several marks be demonstrated with respect to each other. In addition, largerscale regional VLM such as tectonic activity also needs to be accounted for; otherwise, the resulting estimate of the mean sea level trends will be biased.
It is only in the last two decades that space geodetic techniques such as GPS have allowed VLM to be measured to a precision better than 1 mm/year (e.g., Lidberg et al., 2010). Nevertheless, this is a technically challenging task given the systematic errors that contaminate the GPS vertical component. Such errors include, but are not limited to, antenna phase center variations (Rothacher, 2001;Ortiz de Galisteo et al., 2010), atmospheric path delay biases, (Steigenberger et al., 2009), ionospheric signal bending (Petrie et al., 2010a(Petrie et al., , 2010b, hydrological loading (Santamaría-Gómez & Mémin, 2015), present-day ice mass loss (Riva et al., 2017), and atmospheric loading (van Dam et al., 2010) as well as the effect that reference frame stability has on vertical velocity rates (e.g., Schöne et al., 2009;Collilieux & Wöppelmann, 2010).
Recent studies in which vertical crustal motion rates have been estimated at a global set of TG sites [e.g., Wöppelmann et al., 2009;Santamaría-Gómez et al., 2012], the "absolute" vertical rates, while determined with respect to the mass center of the Earth, are actually relative to the chosen reference frame such as ITRF2008 (e.g., Wöppelmann et al., 2007;Collilieux & Wöppelmann, 2010). Although absolute vertical rates will be subject to some error due to the realization of the reference frame, recent realizations of the ITRF have resulted in stable velocities within a small fraction of a 1 mm/year. New Zealand has an important part to play in these global studies. Not only does it have four long-term (over 120 years) TG records, but it also borders the Southland Front/Southland Current, which in turn affects sea level variability at local and regional scales (Wunsch et al., 2007). The global distribution of century-long TG records is biased toward the Northern Hemisphere, making New Zealand TG sites valuable for both regional and global studies. However, New Zealand straddles the plate boundary between the Pacific and Australian tectonic plates, thus giving particular importance to the task of estimating and monitoring the VLM at the TG sites. Seismic events and interseismic tectonic deformation have the potential to deform land where the TG recorder is located, thus biasing the sea level estimates. While there is no indication of significant seismic events occurring in the Auckland region during the period of the measurements reported here, Wellington, Lyttelton, and Dunedin have all been subjected to both major earthquakes and transient events.
The most significant (and damaging) earthquake events affecting a New Zealand TG were the 4 September 2010, M w 7.1 Darfield earthquake, the M w 6.3 Christchurch earthquake on 22 February 2011, the M w 5.9 Godley Head earthquake on 13 June 2011, and a sequence of earthquakes on 23 December 2011 (Kaiser et al., 2012), all causing significant VLM at the Lyttelton TG site (Table 1 and Figure 1). Historically, Canterbury is a low seismic region for which paleoseismology suggests that the penultimate rupture on the Greendale fault (Darfield) occurred between 20 and 30 kyr ago (Hornblow et al., 2014) and that the Christchurch events are a once in 10-kyr occurrence (Mackey & Quigley, 2014). All of these events caused ground motion at the Lyttelton TG, the largest single vertical displacement was >50 mm (cumulative displace-ment~110 mm). Also of concern for the Wellington region, and only recently detectable through dense continuous GPS (cGPS) networks, are slow slip events (SSEs) that are transient in nature and cause nonlinear ground motion at sites overlying the Hikurangi subduction zone along the east coast of the North Island and the north of the South Island. (e.g., Beavan et al., 2007;Wallace & Beavan, 2010;McCaffrey, 2014;Wallace et al., 2017). The south west tip of New Zealand is also an active seismic region with the Australian plate subducting beneath the Pacific plate along the Psuysegur trench. Even though Dunedin is some 200-300 km from the region, it has been affected by recent significant large earthquakes in Fiordland such as the George Sound 2007 (Petersen, et al., 2009) and the Dusky Sound 2009 (Beavan, Samsonov, et al., 2010). In addition, the M w 8.1 Macquarie Island 2004 event (Watson et al., 2010) affected the whole of New Zealand.
Although no GPS measurements are available to quantify VLM at the TG sites prior to 2000, Table 1 (bottom half) lists all earthquakes with M > 7 that have occurred between 1855 and 2003, which may potentially have contributed to VLM in the New Zealand region. These data are compiled from the New Zealand GeoNet Quake Search webpage (quakesearch.geonet.org.nz) and include only earthquakes with depths of less than 35 km. The depth cutoff was included to avoid earthquakes that represent readjustment within the subducting Pacific and Australian plates since these earthquakes do not normally produce significant surface displacements. A starting date of 1855 was selected as this coincides with the major Wairarapa earthquake of 1855, which causes significant VLM at the meter level in the Wellington region.
While the Wairarapa earthquake occurred 36 years before the start of the Wellington TG record, its size (M w 8.2) and proximity to Wellington means that postseismic relaxation may affect the initial TG measurements. Past experience with large subduction zone earthquakes (Suito & Freymueller, 2009) shows that long-lived postseismic deformation can play an important role in subduction zone earthquakes. For example, Suito and Freymueller (2009) show that significant postseismic uplift associated with longer-term viscoelastic transients persisted for over 35 years following the M w 9.2 1964 Alaska earthquake, while Denys et al. (2016) show that velocity changes are associated with viscoelastic transients after the M l 7.8 2009 Dusky Sound earthquake. However, even for the 1964 Alaska earthquake, the rate of postseismic relaxation has decreased to less than 1 mm/year by 2000.
In the New Zealand case, of greater concern is the coseismic displacement from the earthquakes listed in Table 1. In particularly the 9 March 1909 Murchison earthquake and the 12 February 1930 Napier earthquake. Given the large size (M l 7.8), proximity (200-250 km from Lyttelton and Wellington gauges, respectively) and combined with the fact that both earthquakes were predominantly thrust events, significant VLM is possible. Indeed, our experience with the Kaikoura 2016 earthquake (M w 7.8) recorded maximum vertical shifts of 4 cm at 200 km and 3 cm at 250 km . Similarly, Dunedin is potentially affected

Methodology
In this study, the ASL trends at New Zealand's longest-running TGs are estimated. At four of the TG sites, colocated (<500 m) cGPS station was installed in the early 2000s enabling the estimation of VLM associated with the TG. We thus determine the ASL trend as where RSL is the RSL estimated from the TG and VLM is the VLM estimated from the cGPS data.
Where possible, vertical trends are also determined at GPS sites situated within 50 km of the TG to test the consistency of VLM at the regional scale compared to the TG site. We model these biases (as well as antenna offsets) in the GPS position time series analysis. Fault lines (red) are from the GNS Science Active Fault Database (Langridge et al., 2016). In the bathymetry image, the orange color indicates continental crust while green through blue shows progressively deeper oceanic crust.
Particular care has been taken in monitoring the stability of the TG benchmark (TG BM ) with respect to a network of benchmarks. Within the constraints of each TG site, the intent has been to have a network of benchmarks colocated or close to the tide pole, TG and cGPS antenna site (TG GPS ), and some that are further away from the immediate harbor infrastructure. Some networks have more than 10 benchmarks with every attempt being made to have at least one benchmark set in bedrock.
The New Zealand TGs are located at an active plate boundary with significant earthquakes, interseismic, and postseismic deformation and/or SSEs affecting the sites. In determining VLM, we recognize that the cGPS data span is limited to approximately 15 years, while the TG data span is close to~100 years. Given the level of tectonic activity experienced in New Zealand, this value must be interpreted with care since it assumes that the rate measured over the last 15 years applies equally to the past 100 years.

The Data Set
We have analyzed and combined heterogeneous data from TG recorders, precise leveling, GPS leveling, and continuous GPS to determine the stability of sea level change at each TG site. A "site" is defined as the TG, the local benchmarks, the TG GPS station, and nearby (regional) GPS sites. The range of available data is shown in Table 2.
The TG records from the Ports of Auckland, Wellington, Lyttelton, and Dunedin ( Figure 1) provide the longest time series of sea level data in New Zealand. This data record (at places incomplete), and its construction, has previously been discussed in Hannah (1990Hannah ( , 2004. The TG site at New Plymouth (located at Port Taranaki) is a new addition to the series being the only long-term TG record available for the west coast of New Zealand. The data of yearly means derived from the TG record for the five sites are given as supporting information (Data Sets S1-S5).
When this project was started, the key objective was to colocate cGPS at each of the four long record TG sites in order to measure VLM. To establish the stability of each TG/cGPS, precise leveling is undertaken to a network of benchmarks, typically located within 2-3 km of the TG/cGPS site. These networks have included as many as possible historical benchmarks as well as new and reliable benchmarks.
The port TGs are typically located on reclaimed land surrounded by harbor infrastructure that may affect the performance of the cGPS. While the GPS quality control statistics do not indicate adverse antenna/receiver performance, such issues cannot be completely ruled out. All four TG sites have potential sky visibility issues: The Auckland and Wellington sites are adjacent to the central business districts, while the Lyttelton and Dunedin harbors are within extinct volcanoes. The Auckland and Lyttelton TG GPS sites are located on wharves directly beside the TGs, but this was not possible at Wellington. The Dunedin TG GPS site is located on a wharf beside a secondary TG some 7 km away from the primary TG. The two TGs have been linked by precise leveling, tidal datum transfer, and by GPS leveling.
Owing to the lack of a suitable GPS location at the Wellington TG site, the TG GPS station (WGTT) was installed on top of a large building with deep foundations (Te Papa, the National Museum of New Zealand), located some 500 m away. Due to the impracticability of leveling to the top of a building, it was decided to monitor the vertical tie using a short-baseline 24-hr GPS survey between WGTT and a benchmark 50 m from the TG. The final 50 m of the tie is completed by precise leveling. Initially, the leveling connection was measured several times during the 24-hr period in order to monitor any tidally induced tilting of the wharf, but this was found to be unnecessary. As no significant annual cycle variation was detected when the tie was initially observed every 3 months, the frequency was reduced to 6 months or longer.

TG Record
Two different procedures have been used in the analysis of the TG records reported here. In the first instance, and so as to provide consistency with earlier analyses, the procedures outlined by Hannah (1990) were again adopted and were applied to each of the five TG data sets. With the exception of New Plymouth, these records are typically of about 100 years in length. While the mathematical model can be used to estimate any unknown datum offsets, this feature was only used at the Wellington site where a new TG was installed at a new location in 1944 (but with an unknown vertical offset). This is described in greater detail in Hannah (1990). The model not only calculates the linear trend of the TG record, but additional terms are included to account for annual variations in mean sea level due to pressure and temperature variations as well as the influence of the two long-term lunar tides (8.85 and 18.6 years). The model does not attempt to determine parameters associated with other periodic phenomena such as ENSO and Interdecadal Pacific Oscillation (IPO) effects, which do not have strictly regular periods. In this approach each annual mean sea level is treated as an independent observation, albeit weighted according to the completeness of the given annual data set.
In the second instance, the HECTOR software (Bos et al., 2013) has been used to estimate the sea level rates. The noise model used is ARMA(1,0) (autoregression-moving average with 1 AR and 0 MA coefficients), which assumes (weakly) stationary noise. While it is the results from this second methodology that are shown in Table 4, the results from the first methodology have been included as supporting information. For three of the four longest records (Auckland, Wellington, and Dunedin), the differences between the two approaches were 0.03 mm/year or less. At Lyttelton this difference was 0.1 mm/year. However, the derived standard deviations increased by a factor of nearly 2 (range 1.1-1.9) (see also Table 4). This was an expected outcome given the different noise models used in the two processing approaches. The uncertainties as derived from the HECTOR approach are to be preferred.
A critical step in deriving the RSL rates has been the assessment of the reliability of the TG records. This issue has previously been discussed in Hannah (1990Hannah ( , 2010. As part of this project, however, a complete reassessment of all the datum information related to each of the existing long-term TG records was undertaken. In addition, a first complete assessment of the New Plymouth gauge was undertaken. With the exception of Lyttelton, the results of these assessments, in the form of changes that were made to the observed sea level record prior to final data analysis, are shown in Table 3. The Auckland data from 2001-2013 were sourced from Land Information New Zealand (LINZ). On 1 January 2001 a new gauge commenced operation but in May 2003 its zero level was found to be incorrectly set. LINZ applied a +0.0347-m datum correction for the period 1 January 2001 to 8 May 2003 prior to delivery of the data (Table 3). We follow LINZ in assuming that this correction is applicable for the entire 2.35year period. Note. In each case the time interval between the start and finish of data collection are shown, although different data spans have sometimes been used in the various analyses (see text and Tables 5 and 6). a Because of the harbor configuration, an additional GPS leveling line was needed in Wellington (see section 3 for explanation).
The reconstruction of the data from New Plymouth for the period 1956-1973 was a particularly difficult task with a number of datum shifts (typically of 1, 2 or 3 feet; equivalently 0.3048, 0.6096 or 0.9134 m) being encountered. There were also periods of time when this gauge did not function satisfactorily.
At Lyttelton the reassessment of the gauge history resulted in three small corrections (Table 3)  At Dunedin, and with the advantage of new leveling, a previous assumption that the TG data should be corrected for wharf subsidence of −0.3 mm/year since 1964 was found to be incorrect.

TG Trends
The estimates of RSL rise for both the historical and the new analyses at each TG site are shown in Table 4. The data used in the most recent analysis are also shown in Figure 2. The weighted mean for New Zealand as a whole which, with the exception of New Plymouth, generally spans the period 1900-2013, falls in the range of 1.6-1.7 mm/year. This is in excellent agreement with other estimates derived from global data sets that span a similar period [Church & White, 2011;IPCC, 2013. Since 1990, the RSL trend at Auckland TG has increased by 0.23 ± 0.19 mm/year, a time interval that covers several ENSO cycles and perhaps one IPO cycle [Hannah & Bell, 2012]. When the entire data set is broken into two parts (i.e., 1899-1960 and 1961 to present day), their respective linear trends are significantly different with the trend as determined from the latter data set being ≈+0.6 mm/year higher. However, the IPO index over the period 1990-2013 does show a substantial trend, which may affect the New Zealand sea levels and could result in the apparent increase in trends.
The sea level trend at the New Plymouth TG has not previously been determined from a long-term series of data. The record, particularly between 1955 and 1976, is influenced significantly by datum shifts of one type or another. A great deal of effort was spent in attempting to resolve these issues with the information available. The RSL trend is strongly influenced by the mean value for RSL for the period 1918-1921. If the value is removed, the trend is not significantly different (95% confidence level) compared to the full data set. This single value, which was found in old correspondence files, is known to have been derived from the original tide charts. Unfortunately, these charts were discarded  Journal of Geophysical Research: Solid Earth decades ago, making verification impossible. The RSL trend, while lower than at most other ports, has a standard deviation that is 3 times larger, reflecting both the shorter span of data and its inherent uncertainties.
In previous analyses, the linear trend derived from RSL data at the Lyttelton gauge has been higher than the New Zealand average. It has been hypothesized in the past that the primary reason for this has been the lack of data prior to 1924 [Hannah, 1990]. While this is still likely, the revised assessment of the datum history (described in section 4.1) has had by far the greatest influence in reducing the previous trend value.
The linear RSL trend in Dunedin has shown the greatest change since the 2004 analysis. This is due to new evidence that points to the stability of the wharf pile to which the gauge is attached. In the 1990 analysis the wharf pile was assumed to be stable, but in the 2004 analysis the pile was assumed to be sinking along with nearby local benchmarks. However, new leveling collected over the last decade has confirmed both the stability of the wharf pile and the instability of the nearby local benchmarks. The trend thus reverts back to a similar figure as derived in the 1990 analysis when the gauge was assumed stable.
The variation between the five sites in the most recent estimate is~0.8 mm/year. Differences in RSL rates between the five TGs are most likely due to one or more of the following sources: errors in the TG data (including datums and referencing to the local benchmarks), long-term changes in salinity or temperature at the TG site and VLM.

Precise Leveling Data
For each site, the precise leveling data are adjusted relative to a single benchmark that is considered stable. We check this, in part, by analyzing its stability relative to other marks in the network. This provides a time series of benchmark heights, which in turn are used to estimate the vertical stability of each benchmark. Since the precise leveling data are observed using similar procedures and instrumentation, all individual observations are given the same weight. The root-mean-square errors in the precise leveling, which are thus a function of the length of the leveling run, are between 0.5 and 1.0 mm.

Precise Leveling Trends
Long-term TGs, such as those established in New Zealand, are often located in working ports. Inevitably, some TGs or TG structures have been moved, which makes the monitoring of these sites difficult. The Auckland TG site is a good example where both the TG and cGPS site had to be relocated due to wharf redevelopment. In addition, the harbor is close to the Auckland central business district where there have been frequent redevelopments that may affect the benchmark network. Prior to the TG relocation in 2007, all benchmarks related to the TG site are regarded as stable. For the new TG site, the network has not been observed sufficiently frequently to determine the benchmark trends reliably.
Similar to Auckland, the Wellington TG is located close to the central business district, which makes maintaining a level network problematic.
Relative to the leveling network, the TG benchmark, TG BM , is considered stable at −0.06 ± 0.08 mm/year (1 sigma uncertainty). However, as described in section 3.2, it has been necessary to combine precise leveling with GPS leveling. The combined precise and GPS leveling height difference shows a clear reduction of the height difference (−1.0 ± 0.1 mm/year) between the TG GPS (WGTT) and the TG BM (see Figure 3).
Port Lyttelton was substantially uplifted due to the earthquake events in 2010 and 2011 (DR10 and CH11, Table 1). The level network was observed at approximately 4-month intervals during the 2 years following the earthquakes, so although the benchmark trends include more noise due to seismic relaxation following the events, we do have a reliable set of measurements that clearly show the degree of movement. Prior to 2004, the TG BM had been leveled periodically for 20 years and was shown to be stable. In 2004, following the establishment of a new tide pole, it was found to be subsiding at −0.5 mm/year level from 2004 to 2012. Following the earthquakes in 2012, another TG BM site was established and the current rate of subsidence is −1.3 mm/year (based on 2 years of leveling data). The TG GPS benchmark, which is colocated with the TG structure, shows that the structure is uplifting at the rate of +0.29 ± 0.08 mm/year (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).
The current Dunedin TG was relocated in the mid-1990s and is now located in an area where there has been substantial land reclamation. Since 1998, the TG BM has been subsiding at the rate of −0.42 ± 0.31 mm/year.

GPS Data and Analysis
The raw cGPS data have been processed using the Bernese GNSS software (v5.2) (Dach et al., 2015) to determine daily positions in ITRF2008 reference frame (a detailed description is given in supporting information Text S1 and Table S1). For this analysis, we have kept the ITRF2008 solution rather than reprocessing with ITRF2014. We do not expect the impact on the velocities to be significant in the region of interest (<0.2 mm/year, Ballu et al., 2019) The data set spans 5.9-20.6 years with a median data availability of 96% (Tables 5 and S3). In 2009, the Auckland cGPS site (TAKL) was replaced (AUKT) resulting in position time series spanning approximately 6 years for each sites (total data span 14 years). The TG cGPS sites have had various data outages that have resulted in data availability of between 69% and 94%. We use position time series analysis of the daily coordinates for each site. This is carried out using all three dimensions (east, north, and height), but only the height component is reported here. There are two key challenges that need to be addressed, namely, the effect of land deformation caused by earthquake events and the time-correlated (colored) noise of the position time series.
It is well established that GPS time series are time-correlated (colored) noise (Santamaría-Gómez et al., 2011). The largest contribution to the colored noise is likely to be due to unmodeled errors caused by the satellite (orbit and clock) and atmosphere (troposphere and ionosphere) but can also be due to the local environment where the monument stability and VLM contributes to the time series noise. The key A linear velocity term is used to parameterize secular site motion. Offset terms are included to account for equipment changes (e.g., antenna) and environmental changes are accounted for using annual and semiannual cyclic terms. To deal with the various seismic events, we include an offset parameter to account for the coseismic displacement and use standard functions to account for SSEs and postseismic decay. The additional station parameters estimated for each cGPS site are given in Table S2.
Earthquake and transient deformation events require careful analysis especially with respect to the time scales. For example, Jiang et al. (2012) used cubic splines to model SSE in Costa Rica. For our analysis, we follow the procedure outlined in Denys and Pearson (2015), which uses the error function to model the SSE transient deformation, Δ sse , with where A e is the amplitude of the event, t and t e are the time and the mean event time, respectively, and 2τ e is the duration of the event. The error function erf x ð Þ ¼ 1 ffiffi ffi π p ∫ x −x e −t 2 dt represents a sigmoid curve (S-shaped curve) that models position offset caused by the SSE motion as a transient ramp.
The pattern of subduction along the East Coast of the North Island and top of the South Island varies spatially (Wallace et al., 2004, Wallace & Beavan, 2006, Wallace & Beavan, 2010, Wallace et al., 2017. The coupling of the subduction interface changes from the north (freely slipping where the creep is occurring steadily over time) to the south (locked where slip is released periodically on the subduction interface). The nature of the SSEs and therefore the subsidence rate changes to reflect the changing subduction zone, but we find that over small regions (e.g., 50 × 100 km) the displacements are spatial time correlated. To overcome the difficulty of identifying the start and end times of an event, we stack all three components of the time series for multiple sites in close proximity (i.e., <25 km) and use nonlinear least squares estimation to determine the common start and end times and hence the duration and mean time of the event. We estimate different amplitudes for each site and hence different total displacements. Unlike a SSE, the time of an earthquake event is precisely known, but the decay time scale and decay function have to be determined. The type of function, for example, logarithmic, exponential, power law, is dependent upon the type of seismic process, for example, after-slip or viscoelastic response. If the earthquake event is sufficiently large, it may require months to years of data to determine the appropriate function(s). For the earthquake events that have occurred during this project, we find that the logarithmic function is the most appropriate and use this function to estimate the coseismic and postseismic displacement, Δ ps , with where O k is the coseismic offset, A k is the decay amplitude, t and t k are the time and the event time, respectively, and τ k is the decay time scale. Similar to SSE, postseismic deformation has a well-defined regional signal and we use nonlinear least squares to estimate the decay time scale by stacking multiple sites and estimating a common decay time scale but different coseismic displacements and decay amplitudes for each site. Zhang et al. (1997) and Mao et al. (1999) showed that formal errors based solely on white noise can result in the underestimation of the velocity errors by up to an order of magnitude. It is important to determine realistic uncertainties and most analysis describe cGPS position time series noise using temporally correlated (colored) noise. The influence of time-correlated errors in position time series and hence on the estimated parameters, for example, site velocity, needs to be accounted for in the velocity precision estimates (Hackl et al., 2011;Bos et al., 2013). Common methods include the spectral analysis and maximum likelihood estimation (MLE) (Williams et al., 2004) but can be computational intensive. A widely used approach is the CATS software [Williams, 2008] where the time correlated noise in position time series can be approximated with one or more independent components of power law noise.
To overcome the computation burden, Hackl et al. [2011Hackl et al. [ , 2013 introduced the robust Allan variance rate (AVR) method. This approach is based on the Allan variance, originally developed by Allan [1966] to analyze atomic frequency standards, which defines the Allan variance as half the average of the sum squared differences between consecutive measurements sampled using a bin size of τ. To evaluate the noise characteristics of the position time series, the mean of each bin is replaced by its slope.
We use two approaches that are computationally efficient to determine the velocity precision, namely, empirical formulae and a reformulation of the MLE algorithm by Bos et al. [2013]. (Table S3). The approximate formulae [Zhang et al., 1997;Mao et al., 1999;Williams, 2003] incorporates the variances of random walk noise, σ 2 RWN , and white noise, σ 2 WN , to estimate the velocity precision with where A RWN and A WN are the amplitudes of the random walk noise and white noise respectively, T is the time series duration, and N is the number of observations. Due to the short time series of 2-5 years the random walk noise determined by Beavan [2005] is probably too high A RWN ≈2 mm= ffiffiffiffi ffi yr p À Á and we assume amplitude of A RWN ¼ 1 mm= ffiffiffiffi ffi yr p , consistent with Dmitrieva et al. [2015] and Shen et al. [2011]. The amplitude of white noise is determined from the individual time series. Table S3 gives the estimated rates based on a standard least squares solution for comparison purposes. Using the described formula, the range of the velocity uncertainty is ±0.11-0.26 mm/year. By assuming that the noise process is stationary, Bos et al. [2013] implemented the software, HECTOR, which reformulated the standard MLE method to improve the efficiency of the algorithm. The reformulation takes advantage of faster matrix operations and compared to CATS reduces the computational time by up to 2 orders of magnitude. As most GNSS time series exhibit power law noise [Williams et al., 2004;Santamaría-Gómez et al., 2011] we assume a white plus power law noise model to estimate the linear trend and known offsets (e.g., earthquakes and antenna) using HECTOR (Table 5). For all sites (Table S3), the mean power law index is −0.6 and the range of the velocity precision is between 0.10 and 1.06 mm/year (mean 0.33 mm/year).

cGPS Trends
The estimated velocity and one sigma uncertainty using a white noise and power law model are tabulated in Table 5. In addition to the velocity and seasonal terms, the number of offset terms (changes in equipment (i.e., antenna) or earthquake events (coseismic displacement)) are given (Table 5) as are the overall data availability and time span of the data. Also tabulated, where relevant, are the earthquake events that are modeled by including SSE (error function) or postseismic decay (logarithmic) terms. To demonstrate the (in)consistency of the rate determined for some sites, alternative solutions are given in the Table S3. These additional solutions compare the complete or full time series model with alternative site parameterization or a simplified/reduced model. An example of an alternative parameterization is the inclusion of a logarithmic term and a simplified model determines the rate using preearthquake data only. The velocity trend lines are shown in Figure 5. Auckland The Auckland sites, AUCK, WARK, TAKL, and AUKT, are located on the Australian plate and are regarded as being tectonically stable, although this region has been affected by volcanic activity in the past. In 2007, the Auckland TG was relocated due to wharf reconstruction, which required the establishment of a new site AUKT in 2009 to replace TAKL. As the two sites are 50 m apart and on the same wharf structure, they are regarded as being colocated, the combined data set covering a time period of over 14 years.
Overall, there is good agreement of the rates in the Auckland region (Table 5 and Figure 5). The AUCK rate of −0.56 ± 0.17 mm/year, based on nearly 20 years of data, is slightly less than the other three sites. TAKL, AUKT, and WARK have all been operating individually for approximately six years (Table S3). The combined rate for TAKL and AUKT agrees with the individual site rate at the 0.1-mm/year level. Wellington The Wellington sites, WGTN, WGTT, and AVLN, have been affected by SSE along the Hikurangi subduction zone, while two strike-slip earthquakes in Cook Strait in 2013 had no apparent effect on the vertical component. The average interseismic subsidence rates are high at~3 mm/year, but during this time period of nearly 20 years, there are three SSEs starting in 1999SSEs starting in , 2007SSEs starting in , and 2013, each with a duration of approximately 1 year. Although the three SSEs are clearly seen in the horizontal components, the vertical displacement of the first event in 1999 is marginally disenable at <1 mm.
In the absence of SSEs the vertical velocity of the Wellington site is dominated by strain accumulation on the locked subduction interface, which causes subsidence due to the sinking of the subducted plate. The effect of a SSE is a decoupling of part of the plate interface allowing the overriding plate to rebound and the surface to move up. This results in uplift with vertical displacements of +15 mm (WGTN, 19 years), +13 mm (WGTT, 15 years), and +19 mm (AVLN, 9 years). The uplift effectively counteracts the subduction induced subsidence and reduces (less negative) the average subsidence rate by~0.8 mm/year (WGTN and WGTT) and 2 mm/year (AVLN, due to the shorter time period). At all three sites, the subsidence is not steady state and to understand the VLM, it is the average rate for the duration of the TG data (subduction plus uplift) rather than the actual subduction rate that we need to determine. Lyttelton The two sites, MQZG and LYTT, are located on the south side of Banks Peninsular (MQZG) and a wooden wharf in the Lyttelton Port (LYTT). Prior to the Darfield (2010) earthquake, the region was not considered to be within an active fault zone. However, both sites were significantly affected by the Darfield and Christchurch (DR10 and CH11) earthquake events and aftershocks. Vertically, MQZG was uplifted by approximately 10 mm and LYTT by over 100 mm. The position time series has been parameterized by including offsets (earthquake and coseismic) and a postseismic logarithmic decay term (Solution F, Tables 4 and S2). Since the Darfield 2010 event, the combined effect of the coseismic displacements of the Christchurch 2011 earthquake events and postseismic deformation has resulted in subsidence at MQZG that, after 5 years, equates to the original coseismic uplift. Conversely, at LYTT (Figure 5), the postseismic displacement has continued to uplift the site (10 mm over 5 years).
The consistency of the estimated rates at each site has been estimated using different parameter sets. These include (see Table S3), only using the time series data prior to the Darfield 2010 event (Solution B); including a only coseismic offsets for each earthquake event (Solution D); and by excluding the GPS data between the Darfield 2010 and Christchurch December 2011 events (Solution E), that is, only one offset is included to account for the vertical displacement caused by all earthquake events. For both sites, all the estimated rates for the alternative solutions are consistent and agree to better than 0.3 mm/year. Dunedin For the four sites close to Dunedin, OUSD, OUS2, DUND, and DUNT, the estimated rates are consistent at the 1.2 mm/year level (mean rate −1.1 mm/year) when an offset term for the coseismic displacements for each earthquake event is included (Solution D, Tables 4 and S2). It is evident that the earthquakes with epicenters in the Fiordland region affect the cGPS sites. The Dusky Sound 2009 event (DS09), the largest Fiordland earthquake since the TG records began, displaced sites over the whole of the lower South Island (south of Christchurch) and there is ongoing postseismic decay that can be clearly seen in the hori- Similar to the Lyttelton sites, the rates have been estimated using data prior to the Dusky Sound 2009 event only (Solution C, Table S3) and including a postseismic decay (logarithmic) term (Solution F). However, the rates derived from these alternative parameterizations are not consistent. For example, there is a decrease in the rates using the pre-Dusky Sound data (Solution C). A difference between the preearthquake and postearthquake horizontal velocities has been reported by Denys et al. [2016]. The largest change was 10.1029/2019JB018055 Journal of Geophysical Research: Solid Earth observed at DUND (1.3 mm/year), which may, in part, be attributed due to the short time span of the pre earthquake time series (3.9 years).
While it is clear that uplift has occurred, estimating reliable offsets can be difficult due to the small vertical displacements of each event, multiple terms (both equipment/antenna changes) and the less precise nature of the vertical component of the time series. Other factors that affect the rate estimate are periods of noisy data due to receiver/equipment (OUS2) and data outages especially at the time of the Dusky Sound 2009 event (DUNT).
Including a logarithmic term to model the postseismic deformation decreases the rates by 0.2-0.8 mm/year (Solution F). Even though the postseismic deformation is clear in the horizontal component, it is likely there is insufficient signal in the vertical component for the postseismic terms to be estimated accurately. For the Dunedin sites, a conservative approach is adopted and only offset terms are included (Solution D). The uplift due to the four earthquake events at OUSD is 8 mm, which is equivalent to reducing the rate by +0.5 mm/year. Multiple earthquake events affect the site, which account for +4-8 mm of uplift. Clearly, the effect will be different for each site dependent upon the length of the time series data. An average rate does not have a lot of meaning although clearly each site is affected.

Sea Level Rise and VLM
For each site we estimate the secular (long) term vertical velocity, v GPS (Table 6). If applicable, we add corrections for a local trend, v LOCAL , and an earthquake (tectonic process) trend v EQ. The local trend accounts for a vertical velocity between the TG GPS and the TG BM and the earthquake trend is a time average of the offsets caused by earthquake events. Of the two corrections, the earthquake trend is less well defined as it is dependent upon the time interval over which the events have been measured (typically a maximum of 15-20 years for this project) and how regularly the events reoccur.

Tectonic and Local Signals
At the Auckland TG, the local network precise leveling is stable and the vertical rate of the TG BM close to the TG is close to 0, indicating no localized settling of the wharf close to the TG site. The cGPS rate, v GPS (Table 6 ), is derived from the combined TG cGPS sites (TAKL+AUKT). Given that there is little differential vertical motion between the regional cGPS sites (WARK, AUCK, and TAKL+AUKT), this velocity correction should be representative of the entire region.
The Wellington benchmark network and TG BM is stable over the period of measurement. However, the combined GPS and precise leveling height difference between the TG GPS and TG has been reducing consistently over 15 years at the rate of v LOCAL = −1 mm/year. The long-term subsidence rate of WGTT (v GPS = −2.8 ± 0.18 mm/year) reflects the ongoing subduction of the Pacific plate under the Australian plate. However, the total subsidence is significantly less than the subsidence rates would indicate due to the periodic uplift effect of the SSEs. Over the period of time of this project there have been three events that account for uplift of approximately ±15 mm resulting in a rate v EQ of ≅+0.8 mm/year. The v EQ uncertainty (Table 6) is determined from the uncertainty of the offsets caused by each event (the amplitude parameter, A e , in equation (2)) with σ 2 vEQ ¼ ∑σ 2 Ae . The combined VLM trend from the TG GPS , local and SSEs (v GPS + v LOCAL + v EQ ) at the TG is −1.06 mm/year. As the cGPS could not be directly colocated with the TG due to reduced elevation masks, it demonstrates the importance of accurately determining the trend of the height difference between the cGPS site and TG site. This is an example where it cannot be assumed that the VLM at the cGPS and TG sites is the same even though the sites are only 500 m apart.
There are two related tectonic effects causing vertical motion at Wellington. One is the coupling between the Pacific and Australian plates on the subduction interface that underlies Wellington [Darby & Beavan, 2001;Wallace et al., 2004] and the other is SSEs [Wallace & Beavan, 2010]. The coupling causes the Wellington region to be dragged westward relative to the interior of the Australian plate, and it also causes uplift or subsidence of the region depending on the exact distribution of the coupling between the plates. Wellington is located close to the hinge line between subsidence and uplift, so that relatively minor changes in the coupling distribution with time can change not only the magnitude but also the sense of tectonic induced VLM (uplift or subsidence) near Wellington. The M w 8.2 Wairarapa earthquake in 1855 probably involved slip on the lower part of the subduction interface [Darby & Beanland, 1992], and changes in coupling (and possibly postseismic effects) following this event probably persisted at least into the early part of the Wellington TG record. As well as any long-term changes in coupling, there are short-term changes of coupling as a result of SSEs on the deeper part of the subduction interface west of Wellington [Wallace & Beavan, 2010]. The interface between about 40-and 70-km depth remains coupled together for years to decades and causes the Wellington region to subside. The coupling at this depth range then releases and the two sides of the subduction interface slide past each other over a period of~1 year.
Three SSE events have been recorded in the Wellington region since 1997 when the first cGPS instrument (WGTN) capable of detecting these events was installed. The secular trend plus SSE trend (v GPS + v EQ ) therefore represents the average subsidence over the time of the cGPS record. It is not possible to quantify the difference between the 2000-2015 and 1891-2013 VLM rates because of the lack of any VLM or SSE records prior to 1997, but the effect (v LOCAL + v EQ ) is in the correct sense (Table 6) in that it reduces the Wellington TG subsidence rate.
At Lyttelton, the TG and cGPS are colocated on the same structure. The secular cGPS trend at Lyttelton is v GPS = −0.5 mm/year (Table 6) but the leveling measurements since 2001 show that the TG structure (and hence wharf) is uplifting at +0.3 mm/year. As the wooden wharf shows a seasonal motion (horizontal amplitude~4 mm), it is supposed that the structure is being slowly wobbled upward. Since the leveling network reference benchmarks and the TG benchmark are stable (V BM = 0 mm/year), this implies that the level network is subsiding at v LOCAL = −0.3 mm/year. Hence, the TG VLM correction is v GPS + v LOCAL = −0.5-0.3 = −0.8 mm/year. The Christchurch (DR10 and CH11) earthquake sequence uplifted the Lyttelton TG site by over 100 mm. In addition, with 3-4 years of postearthquake data, it is now apparent that there is a small, but significant, postseismic signal that has resulted in >10 mm of uplift since the earthquakes. As the coseismic and postseismic displacements only affect the very last of the TG record, neither displacements have been accounted for in the sea level record but will need to be in the future.
The Dunedin cGPS sites show that the Dunedin region is subsiding at between −1 and −3.2 mm/year (Table  S3). The variation is largely due the cGPS data and the time series parameterization used. Several Fiordland earthquake events give rise to localized tectonic motion, at least over the last decade of a combined displacement of up to +8 mm of uplift. In Table 6, the time-averaged rate has been applied to OUSD as indicative of a change in the trend for Dunedin. Over the time of this project, this equates to a positive correction of v EQ = +0.5 mm/year, which reduces the average subsidence.
At Dunedin, the v GPS is based on the OUSD site, which ironically, is closer to the TG than the TG GPS receiver (DUNT). (This anomaly is something of an historical accident that acts as a constraint on the quality of this result.) Using the OUSD rate based on the <20 years of data to correct the RSL data from the last 100 years gives an apparent mean sea level rise in Dunedin of 0 mm/year. While this rate cannot be discounted, it is likely that periodic uplift due to earthquake events could contribute to this discrepancy. The 100-year average RSL rate at Dunedin may be low due to the absence of significant recent tectonic induced uplift that 10.1029/2019JB018055

Journal of Geophysical Research: Solid Earth
would produce an average of v EQ = +1 mm/year. Uplift could be small and occur reasonably frequently as a result of the Fiordland subduction zone (as demonstrated in the recent position time series record of Pearson et al., 2015), or from nearby regional structures, for example, the Akatore fault [Litchfield & Norris, 2000].

Solid Earth Deformation and GIA Signals
In many, particularly global studies, the observed RSL rates are routinely corrected for GIA, for example, ICE-6Gc model (VM5a) [Argus et al., 2014;Peltier et al., 2015]. Although a viscoelastic GIA correction may be appropriate for formerly glaciated high-latitude regions (e.g., Scandinavia, Canada, and Antarctica), as shown by Houston and Dean (2012), outside these regions the GIA vertical rates are a small fraction of the GPS VLM. In TG sites in New Zealand, the main effect is to introduce a roughly uniform 0.3 mm/year correction in our sea level estimates due to a change in the geoid, v GIA, [Peltier et al., 2015, Figure 23], which is subtracted from the RSL trend, i.e., v TG − v GIA .
The Riva et al. [2017] model of elastic solid earth deformation (SED) shows a consistent but changeable rate of land ice mass throughout New Zealand over the last century ( Figure 4). The total elastic SED displacement for the period 1902 to 2012 (113 years) is approximately −30 mm for the South Island sites (−28.6 and −31.0 mm for Lytteton and Dunedin, respectively) and −36 mm for the North Island sites (−37.0, −36.1, and −35.7 mm for Auckland, New Plymouth, and Wellington, respectively). Compared to the North Island sites, the smaller South Island elastic SED subsidence is due to uplift centered along the South Island's Southern Alps (maximum +16 mm) where the greatest land ice mass has been lost.
Also clearly shown in Figure 4 is the variation in the elastic SED rate over the last century. The maximum rates are >−0.1 mm/year in the mid-1920s, decreasing to <−0.5 mm/year in 1940, a gradual increase for 45 years to 1985, followed by a rapid decrease in rates approaching −0.5 mm/year in 2016. The average over the last 15 years till 2016 is −0.50 mm/year (South Island sites) and −0.47 mm/year for the North Island sites. To account for the changing land ice loss over the twentieth century, we apply v SED correction to each site, which is approximately v SED = −0.3-(−0.5) = +0.2 mm/year. The −0.3 mm/year accounts for the average SED displacement over all 113 years of the SED model, while the −0.5 mm/year is the current day SED, which must be removed as it is included in the GPS estimated VLM rate (v GPS ).

ASL
If VLM is measured and applied correctly, the ASL rates should be consistent at the regional level even if the RSL rates are different. For the New Zealand TGs, correcting the RSL rates using the estimated GPS, local, earthquake SED and GIA rates (v GPS + v LOCAL + v EQ + v SED − v GIA ) results in a mean ASL rate of +1.45 ± 0.36 mm/year. The regional consistency of Sea Level Rise (SLR) has improved from a range of 0.8 to 0.4 mm/year (Table 6). Improved regional SLR consistency can be considered an indicator of better (relative) VLM corrections, which in the New Zealand case includes estimates of local and earthquake corrections. Provided the TG sites and leveling networks continue to be monitored, it is anticipated that the trends will become more reliable as the height time series lengthens. This estimate is on the low side compared to most GMSL rise estimates and is most probably indicative of regional variability [e.g., Hay et. al., 2015] but is consistent with the difference between the Northern and Southern Hemispheres estimates of GMSL rise as shown by Wöppelmann et al. [2014]. There are also specific physical processes that might account for this difference, such as ocean dynamics and biases caused by reference frame issues.

Conclusions
New Zealand's long-term TGs have records that extend back to the late nineteenth and early twentieth centuries. The TG measurements record the change in level between the mean sea level surface and the adjacent land. It is clear that the TG records can be affected by plate tectonic driven seismic events, both instantaneous earthquakes and SSEs and solid earth elastic deformation due to land ice wastage through the twentieth century. To understand the VLM at these sites, we established cGPS at four TG sites and created a leveling benchmark network to monitor the local VLM. We have also included in the analysis, where available, nearby (<50 km) regional cGPS sites to establish the consistency of regional trends. The project has demonstrated that it is critical to measure and understand VLM at TG sites, which can then be applied to TG trends.
The updated analyses of the RSL TG data alone continue to support previously published information regarding RSLs in New Zealand. However, this study indicates that regional tectonic signals have affected the Wellington, Christchurch (Lyttelton), and Dunedin TG records over at least the last two decades. The Auckland TG site has an excellent TG record and VLM appears to be consistent at the regional level. This site should be regarded as the most reliable of the New Zealand TG sites.
The East Coast subduction zone results in a large subsidence rate at the Wellington TG. This rate is offset by a differential rate between the TG and TG GPS (+1 mm/year), determined by a combination of GPS and precise leveling and SSEs (~+1 mm/year) averaged over the GPS measurement period. Because there is clear nonlinearity in the VLM history at Wellington, and because this cannot be reliably documented prior to 1997, we can only determine the cumulative effect of the SSEs as indicative of the long-term trend.
The Christchurch earthquake events resulted in over 10 cm of uplift of the TG structure. If the Christchurch earthquakes are assumed to be rare events with long recurrence intervals, then the secular vertical velocity should be a reliable estimate of the long-term vertical trend. There is a small postseismic signal, which needs to be monitored to confirm that there is no significant long-term change in the vertical rate.
The current analysis gives the RSL rate for Dunedin as being the lowest and this discrepancy may be due to the effect of earthquakes in the Fiordland region over the long term. All the Dunedin sites have been affected by a number of earthquakes that have probably uplifted the Dunedin region. The effect of seismic events during the cGPS record are well constrained, but there have been possibly two earthquakes SW of Fiordland (1945 and1979) that may have occurred on the Puysegur Subduction Zone. Such events are of a size and sufficiently close to Dunedin to cause land movement, but there is insufficient information for the VLM to be quantitatively estimated. The vertical contribution of these earthquakes is the biggest uncertainty in the VLM rates.
Predicted GIA, being a far-field signal in New Zealand, is uniformly similar for all TG sites used in this study. The geoid component of GIA, which is not measured by GPS, increases sea level by~+0.3 mm/year and is in the opposite sense compared to the GPS derived VLM. In contrast, current day land ice loss is measured by GPS but has only been possible to do so for two decades. Since the 1980s, the rate of solid Earth deformation 10.1029/2019JB018055

Journal of Geophysical Research: Solid Earth
has accelerated and is currently~−0.5 mm/year, compared to an average (but variable) rate of~−0.3 mm/year over the twentieth century.
Applying VLM estimates for tectonic, local, GIA, and SED gives a mean ASL of +1.45 ± 0.36 mm/year, which is consistent with global estimates. In applying VLM to each TG site, the regional consistency of the SLR at the regional level improves to within 0.4 mm/year (range +1.2 to +1.6 mm/year). We will always have to make assumptions concerning the appropriateness of applying the VLM rates determined over short time periods to the TG record extending back 100 years and longer. Clearly, VLM contributions from tectonic, local, and SED are significant and must be taken into consideration with the view that the corrections will improve the consistency of regional TG rates when we sample more of these tectonic related events over time. It also implies that unknown and unmodeled VLM is likely to occur at other global long record TG sites and is currently not being taken into account and therefore leading to bias GMSL estimates. Thus, more work is needed to understand the physical processes responsible for lower rates of VLM-corrected sea level change in NZ when compared with GMSL.