<figure><img src="image_1.png" /></figure>

HAL
open science

# Observations of ocean tidal load response in South America from subdaily GPS positions

Hilary R. Martens, Mark Simons, Susan Owen, Luis Rivera

► To cite this version:

Hilary R. Martens, Mark Simons, Susan Owen, Luis Rivera. Observations of ocean tidal load response in South America from subdaily GPS positions. Geophysical Journal International, 2016, 205 (3), pp.1637-1664. 10.1093/gji/ggw087. hal-02497900

HAL Id: hal-02497900

https://hal.science/hal-02497900

Submitted on 27 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Observations of ocean tidal load response in South America from subdaily GPS positions

Hilary R. Martens,<fnref n="1" /> Mark Simons,<fnref n="1" /> Susan Owen<fnref n="2" /> and Luis Rivera<fnref n="3" />

Accepted 2016 February 29. Received 2016 February 28; in original form 2015 November 14

## SUMMARY

We explore Earth's elastic deformation response to ocean tidal loading (OTL) using kinematic Global Positioning System (GPS) observations and forward-modelled predictions across South America. Harmonic coefficients are extracted from up to 14 yr of GPS-inferred receiver locations, which we estimate at 5 min intervals using precise point positioning. We compare the observed OTL-induced surface displacements against predictions derived from spherically symmetric, non-rotating, elastic and isotropic (SNREI) Earth models. We also compare sets of modelled predictions directly for various ocean-tide and Earth-model combinations. The vector differences between predicted displacements computed using separate ocean-tide models reveal uniform-displacement components common to all stations in the South America network. Removal of the network-mean OTL-induced displacements from each site substantially reduces the vector differences between observed and predicted displacements. We focus on the dominant astronomical tidal harmonics from three distinct frequency bands: semidiurnal ($M_2$), diurnal ($O_1$) and fortnightly ($M_f$). In each band, the observed OTL-induced surface displacements strongly resemble the modelled displacement-response patterns, and the residuals agree to about 0.3 mm or better. Even with the submillimetre correspondence between observations and predictions, we detect regional-scale spatial coherency in the final set of residuals, most notably for the $M_2$ harmonic. The spatial coherency appears relatively insensitive to the specific choice of ocean-tide or SNREI-Earth model. Varying the load model or 1-D elastic structure yields predicted OTL-induced displacement differences of order 0.1 mm or less for the network. Furthermore, estimates of the observational uncertainty place the noise level below the magnitude of the residual displacements for most stations, supporting our interpretation that random errors cannot account for the entire misfit. Therefore, the spatially coherent residuals may reveal deficiencies in the *a priori* SNREI Earth models. In particular, the residuals may indicate sensitivity to regional deviations from standard globally averaged Earth structure due to the presence of the South American craton.

**Key words:** Time-series analysis; Satellite geodesy; Tides and planetary waves; South America.

## 1 INTRODUCTION

Tidal forces, generated primarily by gravitational interactions with the moon and sun, deform the Earth both directly through the gravitational potential (body tides) and indirectly through the periodic redistribution of fluid mass loading Earth's surface (e.g. oceanic and atmospheric load tides). The density and elastic structure of Earth's interior controls the spatiotemporal characteristics of the deformation response; thus, observations of surface displacements caused by the tidal potential and tidal loading may potentially be

exploited to study the material properties of the solid Earth (e.g. Love 1909; Melchior 1983; Baker 1984; Bos et al. 2015).

Whereas body tides are long-wavelength phenomena that sample a large-scale average of Earth structure (e.g. Latychev et al. 2009), ocean tidal loads are shorter-wavelength features that probe Earth's material properties at finer spatial scales (e.g. Baker 1984; Jentzsch 1997; Agnew 2015). Moreover, whereas the spatial distribution of amplitude and phase for body tides generally follows that of the equilibrium tide derived directly from the gravitational potential, ocean tides exhibit a complicated spatial pattern due to

<footnote n="1">¹Seismological Laboratory, Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA.
E-mail: hmartens@caltech.edu</footnote>

<footnote n="2">²Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA</footnote>

<footnote n="3">³Institut de Physique du Globe de Strasbourg (UMR7516), Université de Strasbourg/CNRS, Strasbourg, France</footnote>
interactions with continental boundaries and bathymetry (e.g. Zahel 1997; Pugh & Woodworth 2014), thereby generating a rich spectrum of deformation responses.

Seismology remains a primary means for investigating Earth's interior structure, but is also limited in its ability to distinguish between variations in the two elastic moduli and density inside the Earth. Separating the three unknown parameters contained within measurements of body wave velocities necessarily involves scaling laws and assumptions about Earth's interior state (e.g. Dziewonski & Anderson 1981), or additional information, such as normal mode or surface wave data (Ishii & Tromp 1999; Lin et al. 2012), which are limited in spatial resolution and/or depth extent. Alternatively, ocean tidal loading (OTL) excites both elastic and gravitational deformation responses in the Earth at a variety of wavelengths that sample the crust and the upper mantle. Combinations of OTL-response and seismic observations should allow the two elastic moduli and the density to be independently constrained as a function of depth.

Although Earth's deformation response to OTL may be monitored using a variety of techniques, including very long baseline interferometry (VLBI), gravity, tilt, strain, and Global Positioning System (GPS)-inferred displacements (e.g. Baker 1980, 1984; Baker et al. 1991, 1996; Petrov & Ma 2003; Pugh et al. 2011; Yuan et al. 2013; Penna et al. 2015), the load Green's functions (LGFs) for tilt and strain decrease as $r^{-2}$ from the load point, where r is distance, and gravity and displacement LGFs decrease as $r^{-1}$ from the load point. Thus, the gravitational and displacement responses are sensitive to a wider expanse of ocean loads and regional solid Earth structure than tilt and strain (Baker 1984). In addition, local inhomogeneities in mechanical properties limit the effectiveness of using strain and tilt observations to investigate OTL response (Baker 1984; Jentzsch 1997). VLBI installations measure tidal response with high accuracy (Petrov & Ma 2003; Thomas et al. 2007), but only a small number exist worldwide (King 2006). GPS stations, in contrast, are now deployed extensively around the globe and also measure tidal deformation very accurately. Therefore, GPS has emerged as a preferred method for OTL-response investigations (e.g. Schenewerk et al. 2001; King 2006; Thomas et al. 2007; Ito & Simons 2011; Yuan & Chao 2012; Bos et al. 2015).

Recent studies have explored various GPS processing techniques used to isolate tidal harmonic signals (e.g. King 2006; King et al. 2008; Penna et al. 2015). In particular, two main strategies have been proposed for extracting tidal harmonic information from GPS data: static and kinematic. In the static approach, also referred to as the harmonic parameter estimation approach (Penna et al. 2015), tidal harmonic coefficients are estimated along with daily station positions as part of a static GPS solution, generally using precise point positioning (PPP) or network estimation techniques (e.g. Schenewerk et al. 2001; Allinson et al. 2004; King 2006; Thomas et al. 2007; Yuan & Chao 2012; Yuan et al. 2013). Series of static solutions, along with their full variance-covariance matrices, are subsequently combined to derive the final estimates for the harmonic coefficients, typically using a Kalman filter. For the kinematic approach, station positions are estimated from the GPS data at sub-daily intervals, without simultaneously estimating the OTL-induced displacements (e.g. Khan & Tscherning 2001; King & Aoki 2003; King 2006; Penna et al. 2015). A post-processing analysis is then performed to extract harmonic coefficients from the kinematic time-series. We adopt the kinematic approach in this study.

Since deformation responses induced by OTL depend on the material properties of Earth's interior, precise observations may

potentially be exploited to constrain solid-Earth structure. Motivated by recent advancements in ocean-tide models and deployments of dense GPS networks, Ito & Simons (2011) were the first to attempt to invert GPS-inferred observations of OTL-generated surface displacements for one-dimensional profiles of the elastic moduli and density through the crust and upper mantle beneath the western United States. Their forward modelling procedure, however, neglected to account for the motion of the geocentre induced by variations in the surface mass loads, which resulted in an inconsistent mapping between observations and predictions (Blewitt 2003; Fu et al. 2012; Wu et al. 2012; Desai & Ray 2014). The use of inconsistent reference frames introduced long-wavelength artefacts into the displacement-response residuals, which critically biased the analysis and results.

Furthermore, ocean-tide models have long been considered a dominant source of error in predicting Earth's response to OTL, particularly near the coast (e.g. Francis & Mazzega 1990; Agnew 1997; Khan & Scherneck 2003; Bos & Baker 2005; Penna et al. 2008; Yuan & Chao 2012; Yuan et al. 2013). Citing the large coastal errors, Yuan & Chao (2012) and Yuan et al. (2013) opted to examine OTL and body-tide response residuals only at GPS stations located more than 150–200 km inland of the coast. Operating under the assumption that, at the accuracy of their GPS observations, the OTL prediction error could be neglected beyond 150–200 km of the coast, they attributed continental-scale spatial coherency with non-diminishing amplitudes exhibited by the inland residuals to possible elastic and anelsatic deficiencies in the *a priori* body tide model. They made no attempt to invert the residuals for perturbations to their pre-assumed Earth model, however, claiming that sensitivity to structure appears to be small compared with sources of error and requires further investigation.

Penna et al. (2015) and Bos et al. (2015) recently completed a study examining OTL-induced surface displacements in western Europe using kinematic GPS processing methods. For the $M_2$ harmonic, they found statistically significant residual displacements that exceeded the observational uncertainties of 0.2–0.4 mm, and suggested that the discrepancy could be explained in large part by accounting for mantle anelasticity and anisotropy.

Improving the ability to extract tidal harmonics from GPS data and to forward-model Earth's response to OTL not only enhances the possibility of constraining material properties from observations of Earth's OTL-induced deformation, but also of better accounting for the effects of loading response when examining other geodetic signals of interest, such as aseismic transients. Here, we investigate the precision and spatiotemporal characteristics of observed and predicted OTL-induced surface displacements across a regional network of GPS stations in Brazil, Argentina, and Uruguay. Since the tidal harmonics within a given tidal species exhibit similar patterns (e.g. Pugh 1987), we focus our analysis on the dominant tides from three distinct frequency bands: the principal lunar semidiurnal tide ($M_2$), which has a period of 12.42 hr; the principal lunar diurnal tide ($O_1$), which has a period of 25.82 hr; and the principal lunar fortnightly tide ($M_f$), which has a period of 13.66 d. The global distributions of the tide amplitudes for the three harmonics are shown in Fig. 1.

We begin by introducing our procedure for modelling OTL-induced surface displacements and exploring the sensitivity of the predicted deformation to a selection of ocean-tide and elastic Earth models. We then transition into a discussion of the observational methods and results. Finally, we compare the observed and predicted OTL-induced surface displacements and consider the implications of the residual displacements. Additional details of our GPS
<figure><img src="image_4.png" /><figcaption>Figure 1. Spatial distribution of tide amplitudes based on the FES2012 ocean-tide model (Lyard et al. 2006; Carrère et al. 2012) for the (a) principal lunar semidiurnal harmonic, M₂; (b) principal lunar diurnal harmonic, O₁; and (c) principal lunar fortnightly harmonic, M₃. The black box in panel (a) outlines our study area. Note also that the tides remain active beneath the floating ice shelves in Antarctica.</figcaption></figure>

processing and harmonic analysis techniques are provided in the Appendices.

## 2 PREDICTIONS

To predict Earth's elastic displacement response to surface mass loading, we convolve an ocean-tide model with displacement LGFs that represent the response of a spherically symmetric, non-rotating, elastic and isotropic (SNREI) Earth to a point load of unit mass. The predicted surface displacements induced by OTL are given by

$$L(r, S, Z, \rho_{\text{sea}}) = \int_{\Omega} \rho_{\text{sea}}(r') G(|r - r'|, S) Z(r') \, d\Omega, \quad (1)$$

where L is the complex-valued response of Earth at the observation point r, $\rho_{\text{sea}}$ is the density of sea water at the load point r', G is the Green's function per kg load, and Z is the complex-valued height of the ocean tide at the load point. The integral is taken over the surface area of the oceans, $\Omega$. The LGFs depend on the angular distance to the load as well as Earth structure, S, where the structure is assumed SNREI (e.g. PREM). The LGFs are formed by spherical harmonic combinations of load Love numbers, which are derived from integrating the equations of motion through the layered Earth structure with normal-traction boundary conditions applied at the

surface (Farrell 1973). We evaluate eq. (1) using software developed in-house, LoadDef, which is parallelized and written in Python.

Initially, we compute the load Love numbers with respect to the centre of mass of the solid Earth, known as the CE reference frame (Blewitt 2003). The GPS orbit and clock products from our observational analysis, however, are provided in the CM reference frame, which is referenced to the centre of mass of the entire Earth system (including the solid Earth, oceans, and atmosphere) (Desai & Ray 2014). For investigations into Earth's OTL response, predictions and observations must be computed in the same reference frame; otherwise, positioning errors of order 1 mm or more may arise (Fu et al. 2012; Wu et al. 2012; Desai & Ray 2014). Thus, we transform the degree-1 load Love numbers into the CM reference frame prior to computation of the LGFs (Blewitt 2003; Agnew 2012; Wang et al. 2012).

For the convolution (eq. 1), we adopt a station-centred template grid, which simplifies the integration of the LGFs across individual cells, easily facilitates the inclusion of multiple loading models, and naturally allows us to refine the resolution of the integration mesh near the station (Goad 1980; Agnew 1997, 2012). To develop the template grid, we place a station at the pole of a spherical coordinate system, where $\theta$ represents the polar angle and $\alpha$ represents the azimuthal angle, and vary the resolution in polar angle as a function
of distance to the station: $Δθ = 0.001°$, or about 100 m, within $θ = 1°$; $Δθ = 0.01°$ from $θ = 1–10°$; $Δθ = 0.1°$ from $θ = 10–90°$; and $Δθ = 1.0°$ beyond $θ = 90°$. For each $θ$ in the mesh, $Δα = 0.1°$.

To also refine the integration mesh around the coastline, we first extrapolate the ocean model inland by one grid cell and then apply a land-sea mask. In the far-field, we define the land-sea mask using ETOPO1, which provides global topographic and bathymetric relief information at 1 arc-minute resolution (Amante & Eakins 2009). In the near-field (within 1.5° of a station), we use the Global Self-consistent, Hierarchical, High-resolution Shoreline (GSHHS) database at full resolution (Wessel & Smith 1996). Since ETOPO1 registers floating ice shelves as landmasses, but the tides remain active in those regions (Fig. 6), we use the Scientific Committee on Antarctic Research (SCAR) Antarctic Digital Database version 6.0 to define the Antarctic coastline. The value for seawater density in eq. (1) should represent the density of ocean water at the seafloor (Ray 2013); here, we adopt a uniform value of 1035 kg m⁻³.

LoadDef differs from most other OTL-response modelling software in that it derives LGFs on-the-fly from an input Earth model, rather than uses pre-computed LGFs from published tables. The convolution portion of LoadDef, however, is modelled after, and therefore not significantly different algorithmically from, the SPOTL package (Agnew 2012) and both codes yield very similar results (to

within about 1–2 per cent for SPOTL version 3.3.0.2). Minor differences arise from details of the coastline refinement, resolution of the integration mesh, and values adopted for seawater density.

In addition to an SNREI Earth model, the forward model (eq. 1) also requires the input of an ocean-tide model. One class of modern ocean-tide models assimilates satellite altimetry and tide gauge data into global hydrodynamic simulations. A second class of ocean-tide models involves purely empirical developments, relying primarily on satellite altimetry constraints. The models are often made available on regularly spaced latitude–longitude grids, with amplitude and phase values supplied for up to 30 tidal harmonics or more. The ocean-basin and Earth-rotational effects produce tidal circulation systems, centred on points of zero tidal amplitude called amphidromes. Complicated ocean-land interactions present challenges for global hydrodynamic modelling. Furthermore, altimetry satellites, which provide some of the primary data for ocean-tide models, have difficulty sampling at high latitudes and near coastlines. Thus, the ocean tides are notoriously difficult to model, particularly in the polar regions and shallow seas, as depicted in Fig. 2. The resolution and accuracy of the ocean-tide models, however, have improved considerably over the last decade (e.g. Ray 2013; Stammer et al. 2014).

To investigate the sensitivity of the predicted OTL-induced surface displacements to changes in the ocean-tide model, Z, or SNREI

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Differences in tide amplitude between the (a) FES2012 and GOT4.10, (b) FES2012 and TPX08-Atlas and (c) FES2012 and EOT11A ocean-tide models for the $M_2$ harmonic. The differences, denoted by the colour bar, represent the magnitude of vector differences between the complex-valued tide amplitudes at each grid cell.</figcaption></figure>
Earth structure, S, we compare the vector differences, $|L_1 - L_2|$, between pairs of forward models (eq. 1). For an individual tidal harmonic, the OTL-induced displacements at a particular geographic location may be represented by a closed ellipse in 3-D space and traced out completely during each tidal period. To illustrate the harmonic displacements graphically, we depict the combined east and north displacement by a horizontal particle motion ellipse (PME) and denote the vertical motion by the colour of the ellipse.

## 2.1 Ocean-tide model comparisons

Acknowledging the availability of a large number of ocean-tide models, we opt for a representative sampling in the interest of clarity. We consider FES2012 (Lyard et al. 2006; Carrère et al. 2012), TPXO8-Atlas (Egbert & Erofeeva 2002; Egbert et al. 2010), EOT11A (Savcenko & Bosch 2012) and GOT4.10 (Ray 1999, 2013). FES2012 and TPXO8-Atlas were generated from global hydrodynamic simulations that assimilated tide gauge and satellite

altimetry data, whereas EOT11A and GOT4.10 rely primarily on empirical altimetry observations. Moreover, the models FES2012, TPXO8-Atlas and EOT11A include estimates for the fortnightly astronomical tide, $M_f$. TPXO8-Atlas also has local tide models incorporated into its final solution.

As discussed in Desai & Ray (2014), a majority of altimetry-based ocean-tide models do not yet account for the effects of tide-induced geocentre variations on altimetric determinations of ocean-tide heights. To the best of our knowledge, the four models considered here are no exception. Recently, however, the altimetry observations used to constrain GOT4.10 were adjusted for tidal geocentre variations, culminating in an updated model: GOT4.10c. To remain internally consistent in our comparisons of ocean-tide models, we focus here on the four models that are presumably uncorrected for geocentre variations (i.e. FES2012, TPXO8-Atlas, EOT11A and GOT4.10), but include a basic assessment of GOT4.10c with the Discussion.

Fig. 3 shows the vector differences between predicted OTL-induced surface displacements as PMEs for the $M_2$ tidal harmonic.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Particle motion ellipses (PMEs) depicting the vector differences between pairs of predicted OTL-induced surface displacements for the $M_2$ tidal harmonic computed using selected ocean-tide models: (a) FES2012 and TPXO8-Atlas, (b) FES2012 and GOT4.10, (c) FES2012 and EOT11A, (d) TPXO8-Atlas and GOT4.10, (e) TPXO8-Atlas and EOT11A, and (f) GOT4.10 and EOT11A. In each case, we adopted PREM as the input SNREI Earth model. The size and orientation of each ellipse represent the displacement differences for the horizontal components, with a reference ellipse shown in the lower right corner of each panel; the colour of each ellipse represents the displacement difference for the vertical component (upper colour bar). The lower colour bar depicts the $M_2$ tide amplitude difference between each model pair (Fig. 2). Ellipses outlined in black denote stations that recorded more than 1000 d of data during our study period; ellipses outlined in white indicate stations that recorded fewer than 1000 d of data during our study period. Although the forward-modelled results shown here do not depend on data length, we identify the short-record stations for later comparison with the observational results.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Fig. 3, but with the common-mode component (network-mean OTL-induced displacement) removed across the network. The common-mode ellipse that was subtracted from all stations is shown in the black-box inset of each panel. Note the large residuals remaining near the Amazon river delta and Patagonian shelf, which are notoriously difficult regions to constrain in the development of the ocean-tide models. Due to the large uncertainties, we remove stations immediately adjacent to Patagonia and the mouth of the Amazon from our analysis.</figcaption></figure>

Note the spatial coherency in size and orientation of many of the ellipses, suggesting the presence of a uniform OTL-induced displacement, or ‘common-mode’, component that is constant across the entire network. In the context of this manuscript, a common-mode component refers to a constant OTL-response amplitude and a constant OTL-response phase-lag that are common to all stations in the network for a particular tidal harmonic (not to be confused with a network-averaged displacement removed from a geodetic time-series). For example, a comparison of FES2012 with TPXO8-Atlas (Fig. 3, panel a) reveals ellipses oriented primarily in the north–south direction, with non-diminishing amplitudes inland of the coast. Removing the network-mean OTL-induced displacement from each station significantly reduces the magnitude of the differences, as depicted in Fig. 4.

With the common-mode component removed, regions of enhanced ocean-tide model uncertainties appear prominently. Unsurprisingly, the coastal areas around the Amazon river delta and Patagonian shelf, which are difficult to constrain with satellite altimetry and difficult to model hydrodynamically, exhibit substantial

inter-model discrepancies (Fig. 4). We therefore exclude nineteen of the most severely affected stations, located immediately adjacent to Patagonia and the mouth of the Amazon, from all of our subsequent analysis.

In Fig. 5, we show the root-mean-square (RMS) differences between pairs of displacement-response predictions for the $M_2$, $O_1$ and $M_f$ tidal harmonics. The hatching on the bars denotes the RMS differences after the common-mode component has been removed (see also Table 1). Note that a substantial portion of the differences between models may be explained by the common mode. For the horizontal components of the $M_2$ harmonic, in particular, removing the network-mean OTL-induced displacement reduces the sensitivity to choice of ocean-tide model by up to several folds and, in some cases, even reduces it to below the sensitivity to some choices of SNREI Earth model (cf. Fig. 9, Tables 1 and 2).

The RMS values represent the network-averaged discrepancies between predicted OTL-induced surface displacements derived using different ocean-tide models. A total of 97 stations were used in the computation of the RMS statistics, after the exclusion of the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Root-mean-square (RMS) differences between pairs of predicted OTL-induced surface displacements for the South America network made using different ocean-tide models. The hatching on the bars illustrates the RMS misfits after a common-mode component (network-mean OTL-induced displacement) was removed from all stations. The models considered include FES2012 (FES), TPXO8-Atlas (TPXO), EOT11A (EOT) and GOT4.10 (GOT). The ordering of the legend corresponds directly to the ordering of the bars in the figure. Note that the GOT4.10 model does not include the $M_f$ harmonic. In each comparison, we adopted PREM as the input SNREI Earth model. As with all RMS statistics presented in this manuscript, we exclude stations immediately adjacent to Patagonia and the mouth of the Amazon as well as stations that recorded less than 1000 d of data.</figcaption></figure>

<table><caption>**Table 1.** RMS differences between pairs of predicted OTL-induced surface displacements for various ocean-tide models after removal of the common-mode component. The values listed in the table correspond to the hatched bars in Fig. 5. In each comparison, we adopted PREM as the input SNREI Earth model. Note that the GOT4.10 model does not include the $M_f$ harmonic.</caption><thead><tr><th rowspan="3"></th><th colspan="9">Ocean-tide model comparisons</th></tr><tr><th colspan="3">East (mm)</th><th colspan="3">North (mm)</th><th colspan="3">Vertical (mm)</th></tr><tr><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th></tr></thead><tbody><tr><td>FES2012 – TPXO8-Atlas</td><td>0.038</td><td>0.017</td><td>0.008</td><td>0.064</td><td>0.023</td><td>0.010</td><td>0.190</td><td>0.068</td><td>0.023</td></tr><tr><td>FES2012 – GOT4.10</td><td>0.031</td><td>0.012</td><td>–</td><td>0.037</td><td>0.011</td><td>–</td><td>0.108</td><td>0.036</td><td>–</td></tr><tr><td>FES2012 – EOT11A</td><td>0.049</td><td>0.020</td><td>0.007</td><td>0.046</td><td>0.012</td><td>0.010</td><td>0.121</td><td>0.043</td><td>0.023</td></tr><tr><td>TPXO8-Atlas – GOT4.10</td><td>0.050</td><td>0.014</td><td>–</td><td>0.053</td><td>0.022</td><td>–</td><td>0.135</td><td>0.066</td><td>–</td></tr><tr><td>TPXO8-Atlas – EOT11A</td><td>0.052</td><td>0.025</td><td>0.008</td><td>0.060</td><td>0.024</td><td>0.007</td><td>0.178</td><td>0.065</td><td>0.014</td></tr><tr><td>GOT4.10 – EOT11A</td><td>0.070</td><td>0.019</td><td>–</td><td>0.050</td><td>0.014</td><td>–</td><td>0.125</td><td>0.046</td><td>–</td></tr></tbody></table>

Patagonian and Amazon stations as well as stations that recorded less than 1000 d of data. Although the forward models are unaffected by data length and quality, we elect to remove the short-record stations from all RMS computations in order to allow direct comparisons with observational residuals (shown later). Unless specified otherwise, we continue to show the PMEs for stations that recorded less than 1000 d of data, but we distinguish them by white outlines (e.g. Figs 3 and 4). The supporting information contains the tabulated results for all stations, regardless of the receiver location or amount of recorded data.

## 2.2 SNREI Earth model comparisons

We now explore discrepancies between predicted OTL-induced surface displacements computed using four reference Earth models: PREM (Dziewonski & Anderson 1981), STW105 (Kustowski et al. 2008), AK135f (Kennett et al. 1995; Montagner & Kennett 1996) and SNA (Grand & Helmberger 1984). PREM, STW105 and AK135f represent globally averaged structure, whereas SNA represents stable continental shield structure. Below approximately 1000 km depth, SNA assumes the structural properties of AK135f. For PREM, STW105 and AK135f, we replaced the water layer at
<table><caption>**Table 2.** RMS differences between pairs of predicted OTL-induced surface displacements for selected SNREI Earth models. For each comparison, we adopted the FES2012 ocean model. The values listed in the table correspond to the bars in Fig. 9.</caption><thead><tr><th rowspan="3"></th><th colspan="9">SNREI Earth-model comparisons</th></tr><tr><th colspan="3">East (mm)</th><th colspan="3">North (mm)</th><th colspan="3">Vertical (mm)</th></tr><tr><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th></tr></thead><tbody><tr><td>PREM – STW105</td><td>0.138</td><td>0.026</td><td>0.004</td><td>0.074</td><td>0.033</td><td>0.004</td><td>0.112</td><td>0.018</td><td>0.003</td></tr><tr><td>PREM – AK135f</td><td>0.023</td><td>0.003</td><td>0.001</td><td>0.020</td><td>0.006</td><td>0.001</td><td>0.060</td><td>0.011</td><td>0.001</td></tr><tr><td>PREM – SNA</td><td>0.032</td><td>0.010</td><td>0.001</td><td>0.040</td><td>0.013</td><td>0.002</td><td>0.045</td><td>0.009</td><td>0.001</td></tr><tr><td>STW105 – AK135f</td><td>0.125</td><td>0.026</td><td>0.004</td><td>0.069</td><td>0.029</td><td>0.003</td><td>0.118</td><td>0.018</td><td>0.003</td></tr><tr><td>STW105 – SNA</td><td>0.139</td><td>0.033</td><td>0.004</td><td>0.094</td><td>0.043</td><td>0.006</td><td>0.145</td><td>0.025</td><td>0.004</td></tr><tr><td>AK135f – SNA</td><td>0.046</td><td>0.010</td><td>0.001</td><td>0.053</td><td>0.018</td><td>0.003</td><td>0.061</td><td>0.011</td><td>0.002</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Profiles of PREM (black), STW105 (red), AK135f (green) and SNA (blue) through the crust and upper mantle. Panels (a)–(c) show the shear modulus, bulk modulus and density profiles, respectively. Panels (d)–(f) show the differences in shear modulus, bulk modulus and density, respectively, relative to PREM. Note that the elastic moduli are shown on common-log scales.</figcaption></figure>

the surface with typical values for the upper crust: $V_p = 5800$ m s⁻¹, $V_s = 3200$ m s⁻¹, and $\rho = 2600$ kg m⁻³. Profiles of the elastic moduli and density for the reference Earth models are tabulated in the supporting information and shown in Fig. 6. The displacement LGFs derived from each model are depicted in Fig. 7.

Differences between predicted OTL-induced surface displacements for the $M_2$ harmonic, derived from the SNREI Earth models depicted in Fig. 6, are shown in Fig. 8 as PMEs. Consistent with

previous studies, the largest sensitivities are observed near the coast, or at small observer-to-load angular distances (e.g. Ito & Simons 2011; Bos et al. 2015). Rather surprisingly, STW105 exhibits the largest discrepancies in OTL-induced displacement relative to the other models. Although the displacement LGFs appear to suggest strong similarities between STW105 and PREM (Fig. 7), closer inspection reveals that the STW105 LGFs are most discrepant relative to the mean of the four models. In particular, integration of the LGFs
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Displacement load Green's functions (LGFs) in the CM reference frame for the PREM (black), STW105 (red), AK135f (green) and SNA (blue) Earth models. The left panels show the horizontal-displacement component of the LGFs and the right panels show the vertical-displacement component. Panels (a) and (b) depict the LGFs over the angular distance range of 0.01°–170°. The remaining panels show a zoomed-in section from 0.1°–10°. Although the LGFs for STW105 track well the LGFs for PREM, they remain consistently negative in relation to PREM. In contrast, AK135f and SNA oscillate about PREM. Furthermore, note that the angular distances are plotted on a common-log scale. The displacement LGFs are scaled by a factor $10^{12}a\theta$, where *a* is Earth's radius (metres) and $\theta$ is the angular distance between the load point and the receiver (radians).</figcaption></figure>

shows that STW105 differs from the mean integrated-response of the four models by a factor of 0.3 per cent in the horizontal and vertical displacement components. For PREM, AK135f and SNA, the discrepancies are generally less than 0.1 per cent from the mean. Note that the LGFs in Fig. 7 have been plotted on a log-scale and that the AK135f and SNA LGFs oscillate about PREM, whereas the STW105 LGF remains consistently negative relative to PREM.

The specific structural reasons for the discrepancies between LGFs are not obvious, but the bulk modulus profile for STW105 remains relatively low within the upper 300 km, which might provide part of the explanation. Ultimately, despite the differences in elastic moduli and density, the corresponding perturbations to the LGFs, and by extension the OTL-response predictions, are very small. Fig. 9 shows the RMS differences between the selected models. The RMS differences for the $M_2$ tidal harmonic are largest, a reflection of its relatively large load amplitude. Despite the large contrast exhibited by STW105, however, the RMS differences never exceed 0.15 mm and are significantly smaller than that for both the $O_1$ and $M_f$ tides, suggesting only a subtle sensitivity to SNREI-based structural variations. Table 2 lists the RMS differences explicitly.

# 3 OBSERVATIONS

## 3.1 Kinematic GPS processing

For many geodetic studies, the OTL-response signal is an inconvenient source of noise; thus, the signal is often removed at the GPS processing stage using forward-modelled coefficients for the dominant tidal harmonics. Smaller tidal harmonics are typically modelled and removed as well by interpolation of the admittance for the dominant harmonics, where the admittance is assumed smooth across each tidal constituent band (e.g. Foreman 1977; Agnew 2012). We, however, aim to retain and isolate the OTL displacement-response signal by initially generating sub-daily time-series of site displacements and then performing harmonic analyses to extract individual tidal harmonics.

We use GIPSY version 6.2 (Zumberge et al. 1997) in PPP mode to process the GPS data at individual receiver sites without requiring interstation double-differencing. PPP implementation relies on precise satellite orbit and clock products determined from a global network of GPS satellites and permanent receivers (Zumberge et al. 1997). Our data set consists of up to 14 yr of time-series from 160
<figure><img src="image_11.png" /><figcaption>**Figure 8.** PMEs depicting the differences between pairs of predicted OTL-induced surface displacements for the M₂ tidal harmonic computed using selected SNREI Earth models: (a) PREM and STW105, (b) PREM and AK135f, (c) PREM and SNA, (d) STW105 and AK135f, (e) STW105 and SNA, and (f) AK135f and SNA. In each case, we adopted FES2012 as the input ocean-tide model. The size and orientation of each ellipse represent the differential horizontal-displacement response, with a reference ellipse shown in the lower right corner of each panel; the colour of each ellipse represents the differential vertical-displacement response (upper colour bar). The lower colour bar depicts M₂ tide amplitude in the oceans.</figcaption></figure>

stations in South America. Tables of geographic coordinates and data availability for each station are provided with the supporting information. The median time-series length is 1760 d, or nearly 5 yr.

GIPSY performs single-receiver ambiguity resolution by estimating wide-lane and phase biases for each station individually using double-differences relative to an extensive global network of other ground-based stations in view of the same satellites (Bertiger *et al. 2010*). The wide-lane and phase bias estimates are pre-computed and distributed by the Jet Propulsion Laboratory (JPL) along with the precise satellite orbit and clock products. Our standard, kinematic PPP methodology involves processing 30 s data to obtain position estimates every 5 min using a random-walk stochastic parametrization of the position estimate. We process the data in 30 hr batches (i.e. a full day plus 3 hr on either side of the day) and extract only the central 24 hr of positions in order to mitigate end effects (King & Aoki 2003; King 2006; Penna *et al. 2015*). The OTL displacement-response signals remain unmodelled at the GPS processing stage. Tidal harmonics are extracted from the position time-series using harmonic analysis (discussed later).

We use JPL's precise satellite orbit and clock products in final and fiducial form ('flinnR' format; version 'repro2.0'). The refer-

ence frame adopted for the orbit and clock products at the time of processing was IGS08 (Altamimi *et al. 2011*; Rebischung *et al. 2012*), which has its coordinate origin at the centre of mass of the total Earth system (CM; Blewitt 2003; Wu *et al. 2011*). To ensure consistency with the observations, we compute predictions of OTL-induced surface displacements in the CM frame as well (Fu *et al. 2012*). The ocean-tide model used to develop the 'repro2.0' version of the orbit and clock products was FES2004 (Lyard *et al. 2006*), which did not account for the effects of tidal geocentre variations on the altimetric observations used to constrain the ocean-tide model (Desai & Ray 2014).

Tropospheric zenith delay terms are estimated stochastically relative to nominal values provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) and mapped to lower elevation angles using the Vienna Mapping Functions (VMF1) (Boehm *et al. 2006*). We apply an elevation-angle cut-off of 7° and assume elevation-dependent weighting according to the square root of the sine of the elevation. Guided by synthetic test results discussed in Appendix A, we adopt a process noise setting for the tropospheric zenith delay of $5.0 \times 10^{-8} \text{ km s}^{-\frac{1}{2}}$, or 3 mm $hr^{-\frac{1}{2}}$ (cf. Yuan & Chao 2012; Yuan *et al. 2013*). Horizontal tropospheric gradient
<figure><img src="image_12.png" /><figcaption>**Figure 9.** RMS differences between pairs of predicted OTL-induced surface displacements for the South America network made using selected reference Earth models (Fig. 6). The input ocean-tide model, FES2012, remains consistent for each comparison. The ordering of the legend corresponds directly to the ordering of the bars in the figure.</figcaption></figure>

parameters are also estimated using a process noise value of 5.0 × 10⁻⁹ km s⁻¹, or 0.3 mm hr⁻¹ (Bar-Sever et al. 1998; Larson et al. 2010).

We apply phase-centre corrections to the receiver antennas, extrapolating the models down to 7°-elevation as needed. The effects of solid Earth and pole tides are modelled and removed according to IERS conventions (Petit & Luzum 2010). The JPL orbit and clock products did not include a second-order ionospheric correction at the time of our analysis. Since using different correction factors for a PPP analysis compared with those adopted by the orbit and clock products can lead to artefacts in the position estimates, we opted instead for first-order ionospheric corrections in our PPP computations to mitigate any potential bias.

The choice of stochastic parameters for the station-position estimates can critically control the ability to resolve an OTL-response signal. For a random-walk formulation, overly strict stochastic parameters can dampen the true signal due to excessive smoothing. In contrast, overly loose stochastic parameters, such as in the extreme case of a white-noise parametrization, might unnecessarily keep noise levels high and therefore reduce the ability to extract the signal of interest (e.g. Larson et al. 2001; King & Aoki 2003). To explore the trade-off, we performed a second series of synthetic tests to determine an appropriate coordinate process noise setting, which we ultimately set to 5.0 × 10⁻⁷ km s⁻¹. Details of the synthetic tests are provided in Appendix A.

## 3.2 Harmonic analysis

The forcing function that generates the tides, i.e. the astronomical ephemeris, may be broken down into discrete periods that can be

combined algebraically to excite responses within and on the Earth that are also periodic. A tidal harmonic, $n$, may be characterized by a harmonic expression of the form (e.g. Foreman et al. 2009):

$$A_n f_n(t) \cos(V_n(t) + u_n(t) - \phi_n), \quad (2)$$

where $A_n$ is the amplitude of tidal harmonic $n$, $V_n$ represents the astronomical argument relative to the Greenwich Meridian at time $t$, and $\phi_n$ is the phase lag in degrees measured relative to the equilibrium tide observed at Greenwich. We adopt the convention of phase lags positive. The time-dependent factors $f_n$ and $u_n$ correct for amplitude and phase modulations that arise due to the presence of subsidiary peaks in the frequency domain, which alter the complex-valued amplitude of a primary harmonic over time (e.g. Foreman 1977; Foreman et al. 2009). The subsidiary harmonics are most commonly separated from a main harmonic by cycles of the lunar perigee (8.85 yr) and lunar ascending node (18.6 yr). Some harmonics are also separated by the cycles of perihelion, but since the period of perihelion is so long (>20 000 yr), they are often neglected. Following Foreman et al. (2009), we update the astronomical argument as well as the amplitude- and phase-modulation factors at every epoch in the time-series.

For the selection of primary harmonics, we consider the Rayleigh criterion (with a cut-off factor of 1.1 cycles) to determine the ability to separate two constituents in the frequency domain over a given time span of observations. To prevent aliasing, we only consider harmonics that have frequencies less than half the sampling rate. We perform the Rayleigh comparison in hierarchical fashion, beginning with the largest-amplitude harmonics, based on the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The left panels show time-series of GPS-inferred receiver positions (top panel in each pair) and residuals (bottom panel in each pair) for station RIO2 during one week in 2010. Each pair of panels depicts the data (blue dots), model fit (solid black line) and residuals (black dots) for each spatial component of the displacement. The fits to the data for each component, which were made to the full 7 yr time-series, include tidal harmonics and a linear trend term; the mean has been subtracted from the position estimates. The right panels show histograms of residuals from the full 7 yr time-series for the east- (top), north- (centre) and vertical-displacement (bottom) components.</figcaption></figure>

Cartwright–Taylor–Edden (CTE) equilibrium tide catalogue (Cartwright & Taylor 1971; Cartwright & Edden 1973).

In the pelagic ocean, where the tide amplitudes are much smaller than the water depth and the tide wavelengths are much longer than the water depth, the tides are well described by the astronomical harmonics. In shallow seas and estuaries, however, nonlinear effects become important and produce higher-order harmonics, including overtides and compound tides (e.g. Doodson & Warburg 1941; Doodson 1957; Parker 2007; Pugh & Woodworth 2014). We account for contributions to the time-series by a selection of nonlinear, or shallow-water, harmonics that arise from distortions and interactions of the semi-diurnal and diurnal astronomical harmonics. We follow the suggestions of Godin (1972) and Foreman (1977) to guide our selection and ordering of shallow-water harmonics to test using the Rayleigh criterion. We also refrain from including a particular shallow-water harmonic until all of its contributing astronomical tides have also been included. A complete list of the nonlinear harmonics used in our study is provided with the supporting information.

Fig. 10 shows an example of our model fit to a kinematic GPS time-series. In this particular case, the fit was made to 7 yr of data recorded at coastal station RIO2, from which we extracted a random snapshot of one week. To clean the data prior to the

harmonic analysis, we removed isolated segments that spanned less than 30 d and that were separated from other data in the time-series by more than 60 d. We also removed large outliers prior to the harmonic analysis based on a running median absolute deviation, with a cut-off criterion of three standard deviations. Further details of our estimation procedure are discussed in Appendix B. After deriving an initial solution for the complex-valued amplitudes of the tidal harmonics, we applied a sidereal filter (20-day window) to the residuals to estimate the contribution of repeating multipath signals to the original time-series. An estimate of the multipath was made for every 5 min epoch in the time-series by averaging nearby position estimates, within a 10 d window on either side of each epoch, that were separated in time by integer multiples of the 23.93 hr sidereal day (rounded to the nearest epoch). We then removed the receiver-specific multipath signals from the original time-series and re-estimated the harmonic coefficients.

## 3.3 Residuals

The comparisons of predicted OTL-induced site displacements from Sections 2.1 and 2.2 provide a general appreciation for the sensitivity of Earth's elastic response to different load and structural
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Observed and predicted OTL-induced surface displacements for the M₂ (top), O₁ (centre) and Mf (bottom) tidal harmonics. The colour bars on the left denote tidal amplitudes. The colour bars on the right denote OTL-induced vertical displacements. The size and orientation of each ellipse indicate the horizontal-displacement response, with a reference ellipse shown in the lower-right corner of each panel. The predictions were computed using PREM and the FES2012 ocean-tide model. Since the PMEs do not show phase information explicitly, we include vector diagrams of the residuals between the observed and predicted OTL-induced displacements with the supporting information.</figcaption></figure>

model inputs. To ascertain the applicability of the models to the South American continent, however, we must compare the predictions against GPS-inferred observations of OTL-induced surface displacements. In particular, we are motivated to explore the suitability of different elastic Earth models with application to the South American continent, which has been shown by seismic tomography to support a deep cratonic keel (Supporting Information Fig. S5).

The observational results for the M₂, O₁ and Mf tidal harmonics are shown in Fig. 11 as PMEs at each of the GPS stations considered in our analysis. The figure also shows a corresponding set

of predictions derived from PREM and FES2012 using our convolution software LoadDef. In general, the observations and predictions show great resemblance for each harmonic, including the small-amplitude Mf harmonic. One notable discrepancy between observations and predictions, however, is the apparent northward offset of the vertical-displacement response for the Mf harmonic seen in the mid-continent (Fig. 11, panels e and f, 30°S–10°S). In particular, the smallest response amplitudes for the vertical component occur further to the north in the observational results (panel e) relative to the predicted results (panel f). Since the observations exhibit spatial coherency, the offset is unlikely to result solely from
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Residuals between observed and predicted OTL-induced surface displacements for the M₂ harmonic, shown as PMEs. The size and orientation of each ellipse represent the residual horizontal-displacement response; the colour of each ellipse represents the residual vertical-displacement response (upper colour bar). A mean particle motion ellipse (common-mode component) has been removed from the residual displacements in each panel (shown in the black box inset). The lower colour bar depicts the M₂ tide amplitude. The predicted OTL-induced displacements were computed using the following ocean and Earth model combinations: (a) PREM and FES2012; (b) STW105 and FES2012; (c) SNA and FES2012; (d) PREM and TPXO8-Atlas; (e) PREM and EOT11A; and (f) PREM and GOT4.10. Here, we have excluded stations that recorded less than 1000 d of data in addition to the stations already removed near the Amazon basin and Patagonian shelf; a total of 97 stations remain.</figcaption></figure>

random observational error. Uncertainties in the ocean-tide models, which are difficult to constrain for the small-amplitude long-period tides, might account for some of the spatial discrepancy.

Fig. 12 shows the residuals between observed and predicted OTL-induced surface displacements for the M₂ tidal harmonic based on multiple forward models. The common-mode components (network-mean OTL-induced displacements) have been removed (boxed PMEs). From Fig. 12, we note two important results. First, the residuals remain consistent across the various forward models. The consistency in the residuals, regardless of the SNREI Earth model or ocean-tide model used to generate the predictions, implies a general insensitivity to variations in the particular forward models considered here. Second, each unique set of residuals (e.g. panel a of Fig. 12) exhibits a regional spatial coherency. Random observational errors are unlikely to produce such systematic results. We therefore suggest that the spatially non-random patterns of the PMEs across South America might result from deficiencies in our

forward-model assumptions, such as spherical symmetry and pure frequency-independent elasticity.

Figs 13 and 14 show the residuals for the O₁ and Mf tidal harmonics, respectively, also with the common-mode component removed. The residuals for O₁ and Mf are smaller than for M₂, albeit not substantially. As with M₂, the residuals remain consistent regardless of the adopted forward model, corroborating the general insensitivity to choice of SNREI Earth model and modern ocean-tide model at the precision of the observations. In contrast to the residuals for M₂, the residuals for O₁ and Mf lack pronounced regional spatial coherency. However, some of the O₁ residual ellipses appear to be oriented in the direction of maximum loading, potentially indicating forward-model deficiencies (cf. Fig. 8). Moreover, some of the smaller Mf residual ellipses appear to be systematically oriented in the east-west direction. Figs 16 and 17 show the RMS misfits between the observed and predicted OTL-induced surface displacements derived from combinations of ocean-tide and
<figure><img src="image_16.png" /><figcaption>Figure 13. Same as Fig. 12, but for the O₁ tidal harmonic. Note, however, that the scales for the ellipses and ocean-tide amplitude differ from Fig. 12.</figcaption></figure>

SNREI-Earth models, respectively. Table 3 lists the RMS misfits explicitly.

## 3.4 Uncertainty estimates

To estimate the statistical significance of the residuals relative to observational noise levels, we compute errors for the derived amplitude and phase parameters using the techniques outlined in Pawlowicz et al. (2002). In particular, we determined the average power spectral density of the post-fit residuals within a frequency window of width 0.4 cycles per day centred around each tidal harmonic. The width of the frequency window is chosen such that the power spectrum may be assumed approximately flat within the window. A representative power spectrum is provided with the supporting information.

Since the average power spectral density provides a variance estimate for the harmonic coefficients (c and s in eq. B7), we must also map the error estimates into amplitude and phase values. We use a parametric bootstrap algorithm to derive distributions of several thousand amplitude and phase values for each harmonic. To derive each instance, we add noise to the original amplitude and phase values based on randomly selected samples from a normal distribution, with the variance of the normal distribution determined from the power spectra. Standard statistical analyses may then be

performed on the bootstrapped distributions. The results are shown as PMEs in Fig. 15. Complete tables of the observed OTL-induced surface displacements, including error estimates, are provided with the supporting information.

To verify the appropriateness of the error estimates, we made a second assessment of the observational uncertainties using an alternative technique. For the second method, we performed independent harmonic analyses on yearly chunks of time-series data. For thirteen stations that recorded at least 10 yr of data, we computed statistics on the distributions of harmonic coefficients derived from the yearly analyses. The 2σ standard deviations for the amplitudes and phases are listed in Table 4. Although the estimates may seem relatively high, in this case we have only performed the harmonic analyses on up to 1 yr of data at a time, rather than on multiple years of data. Performing the power-spectrum analysis (Method 1) on a single year of data also shows similarly larger errors (Table 4, centre columns) and provides an additional level of confidence in the error estimates. Moreover, the estimates of observational uncertainty are generally consistent with the 0.2 mm level uncertainties derived by Penna et al. (2015) through rigorous synthetic testing.

The accuracy of the harmonic estimates generally improves with longer time spans of data. Since good convergence for the lunar-derived harmonics may be achieved with at least ~1000 d of data (Yuan et al. 2013) and the median time-series length for our station
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Same as Fig. 12, but for the $M_f$ tidal harmonic. Note, however, that the scales for the ellipses and ocean-tide amplitude differ from Fig. 12. The GOT4.10 model does not include $M_f$; therefore, panel (f) shows the residuals between the observed and predicted OTL-induced surface displacements computed using AK135f and FES2012.</figcaption></figure>

<table><caption>**Table 3.** RMS misfits between observed and predicted OTL-induced surface displacements (from Figs 16 and 17). Solid Earth body tides were removed according to IERS conventions at the GPS processing stage. A common-mode component (network-mean OTL-induced displacement) was also removed prior to the RMS computation. The SNREI-Earth and ocean-tide models used in each comparison are noted in the first column of the table.</caption><thead><tr><th rowspan="3"></th><th colspan="9">Observational residuals</th></tr><tr><th colspan="3">East (mm)</th><th colspan="3">North (mm)</th><th colspan="3">Vertical (mm)</th></tr><tr><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th></tr></thead><tbody><tr><td>STW105 & FES2012</td><td>0.237</td><td>0.125</td><td>0.082</td><td>0.255</td><td>0.130</td><td>0.073</td><td>0.342</td><td>0.206</td><td>0.251</td></tr><tr><td>AK135f & FES2012</td><td>0.180</td><td>0.124</td><td>0.081</td><td>0.234</td><td>0.128</td><td>0.073</td><td>0.340</td><td>0.207</td><td>0.251</td></tr><tr><td>SNA & FES2012</td><td>0.188</td><td>0.124</td><td>0.081</td><td>0.239</td><td>0.129</td><td>0.072</td><td>0.330</td><td>0.205</td><td>0.251</td></tr><tr><td>PREM & FES2012</td><td>0.177</td><td>0.124</td><td>0.081</td><td>0.234</td><td>0.127</td><td>0.072</td><td>0.319</td><td>0.205</td><td>0.251</td></tr><tr><td>PREM & TPXO8-Atlas</td><td>0.167</td><td>0.129</td><td>0.078</td><td>0.253</td><td>0.130</td><td>0.074</td><td>0.298</td><td>0.214</td><td>0.248</td></tr><tr><td>PREM & EOT11A</td><td>0.163</td><td>0.118</td><td>0.079</td><td>0.248</td><td>0.127</td><td>0.072</td><td>0.340</td><td>0.208</td><td>0.249</td></tr><tr><td>PREM & GOT4.10</td><td>0.195</td><td>0.127</td><td>—</td><td>0.237</td><td>0.129</td><td>—</td><td>0.321</td><td>0.210</td><td>—</td></tr></tbody></table>

network is nearly 5 yr, we elect to exclude all stations that recorded less than 1000 d of data from the RMS computations. As stated previously, we also removed nineteen stations with demonstrable sensitivity to large ocean-tide model errors around Patagonia and the Amazon delta. Unless stated otherwise, we retain the short-record stations in the figures showing PMEs, which we distinguish from the long-record stations by white-outlined PMEs (e.g. Fig. 11).

## 4 DISCUSSION

We have derived observations and predictions of OTL-induced surface displacements across South America for the $M_2$, $O_1$ and $M_f$ tidal harmonics. For each harmonic, the observed OTL responses exhibit spatiotemporal coherency and match the predicted responses at most stations to within about 0.3 mm (Figs 11–14). Both
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Observational uncertainties computed from a windowed power spectrum and bootstrap analysis (see the text for details). The left panels show 2σ standard-deviation (SD) error ellipses for the observed OTL-induced surface displacements. Specifically, the size of each ellipse denotes the horizontal-amplitude uncertainty in both the east and north components; the colour of each ellipse denotes the vertical-amplitude uncertainty. Phase-uncertainty information is not displayed, but may be found in tables supplied with the supporting information. For comparison, the right panels (b, d, and f) show the residual surface displacements derived from PREM and FES2012, reproduced from Figs 12(a), 13(a) and 14(a). The ellipses and vertical-displacements (right colour bars) are shown on the same scale to facilitate comparison. Note that the horizontal scale for the $M_2$ residuals shown in panel b differs from the scale in Fig. 12(a) by a factor of two. The left colour bars show tidal amplitude in the oceans.</figcaption></figure>

modelling and observational uncertainties contribute to the sub-millimetre residuals, and an improved understanding of each contribution may potentially be used to refine ocean-tide models and to constrain solid Earth structure. In particular, we find large uniform-displacement components in the differences between predicted OTL-induced site displacements for various ocean-tide models. Removal of the ‘common-mode’ component (network-mean OTL-induced displacement) significantly reduces the discrepancies between predictions of OTL response made using different ocean-tide models (Figs 4 and 5). Moreover, the residuals between the ob-

served and predicted OTL-induced surface displacements are also substantially reduced by removing the common-mode component, particularly for the $O_1$ harmonic (Fig. 13).

Although a detailed investigation into the origin(s) of the common-mode component is beyond the scope of this paper, we suggest that large ocean-tide model uncertainties in the polar regions, such as under the Antarctic ice shelves, could contribute to the offset in OTL-induced displacement. Direct differences between three pairs of ocean-tide models reveal large discrepancies in the polar regions around both the Arctic and the Antarctic
<table><caption>**Table 4.** Error estimates for the observational results from 13 stations that operated for periods of at least 10 yr: BELE, BOMJ, BRAZ, CUIB, IGM1, IMPZ, POAL, RECF, RIOD, SMAR, TUCU, UNSA and VICO. All errors are quoted as $2\sigma$ standard deviations from the mean. We directly compare two methods for computing the errors based on the same set of stations. For Method 1, we computed the errors based on an analysis of the power spectrum around each tidal harmonic, following the techniques discussed in Pawlowicz *et al.* (2002). We performed Method 1 with the full data sets (>10 yr) at each station as well as just for 1 yr of data during 2013. For Method 2, we processed the data at each station in yearly chunks and then performed a statistical analysis of the resulting sets of amplitude and phase parameters.</caption><thead><tr><th rowspan="3"></th><th colspan="9">Observational error analysis</th></tr><tr><th colspan="3">Method 1 (full data set)</th><th colspan="3">Method 1 (1 yr)</th><th colspan="3">Method 2 (yearly chunks)</th></tr><tr><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th><th>M<sub>2</sub></th><th>O<sub>1</sub></th><th>M<sub>f</sub></th></tr></thead><tbody><tr><td>East amplitude (mm)</td><td>0.036</td><td>0.046</td><td>0.070</td><td>0.123</td><td>0.154</td><td>0.283</td><td>0.193</td><td>0.210</td><td>0.212</td></tr><tr><td>East phase (°)</td><td>0.913</td><td>1.316</td><td>11.666</td><td>3.090</td><td>4.656</td><td>66.316</td><td>4.625</td><td>4.928</td><td>39.738</td></tr><tr><td>North amplitude (mm)</td><td>0.036</td><td>0.041</td><td>0.069</td><td>0.125</td><td>0.140</td><td>0.273</td><td>0.209</td><td>0.267</td><td>0.301</td></tr><tr><td>North phase (°)</td><td>0.613</td><td>1.353</td><td>9.573</td><td>2.220</td><td>5.575</td><td>27.520</td><td>2.936</td><td>5.600</td><td>39.814</td></tr><tr><td>Vertical amplitude (mm)</td><td>0.147</td><td>0.161</td><td>0.206</td><td>0.469</td><td>0.515</td><td>0.674</td><td>0.600</td><td>0.588</td><td>0.720</td></tr><tr><td>Vertical phase (°)</td><td>0.931</td><td>4.396</td><td>46.217</td><td>3.055</td><td>16.031</td><td>58.333</td><td>4.500</td><td>18.170</td><td>123.031</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** RMS misfits between observed and predicted OTL-induced surface displacements made using selected ocean-tide models. In each case, we adopted the SNREI Earth structure of PREM. As with all other RMS computations, we excluded stations immediately adjacent to the Amazon river delta and the Patagonian shelf (Fig. 4) as well as stations that recorded less than 1000 d of data (Table 4). The hatching on the bars shows the RMS misfits after the removal of the common-mode component.</figcaption></figure>

(Fig. 2). In Fig. 18, we partition the direct differences for the $M_2$ harmonic into latitude bands and compute the predicted displacement response in South America due to the discrepancies between ocean-tide models. Note the strong north-south trend in the PMEs in response to the southern-most band from 90°S to 60°S (Fig. 18, panel f; cf. Fig. 3). Significant contributions, however, also come from mid-latitude bands. Similar figures showing the comparisons between FES2012 and GOT4.10 as well as FES2012 and EOT11A are included with the supporting information.

We suggest that a second contributor to the common-mode component could involve reference-frame inconsistencies. Reference-frame inconsistencies can manifest at several points in an analysis of OTL response. First, if the load Green's functions used to predict the OTL-induced surface displacements are computed in a reference frame that differs from that of the observations, then errors of order 1–2 mm may arise in the residual displacements (e.g. Fu *et al.* 2012). The JPL orbit and clock products used in our analysis (version 'repro2.0') were referenced to CM; therefore, we transformed the load Green's functions to the CM frame prior to the
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Same as Fig. 16, but showing the RMS misfits between observed and predicted OTL-induced surface displacements made using selected SNREI Earth models. In each case, we adopted the FES2012 ocean-tide model. The orange bars (PREM) therefore match exactly the orange bars (FES2012) in Fig. 16.</figcaption></figure>

convolution with a load model (eq. 1). Second, since most ocean-tide models are constrained by satellite altimetry, with orbits typically referenced to CM (Desai *et al.* 2014), the altimetric observations of sea-surface height must be adjusted for load-induced variations in the geocentre (Desai & Ray 2014). Fig. 19 shows the vector differences between predicted OTL-induced surface displacements computed using GOT4.10c and GOT4.10. The two ocean-tide models differ only in one aspect: the altimetric measurements used to constrain GOT4.10c were adjusted for tidal geocentre variations, whereas the measurements used to constrain GOT4.10 were not. Desai & Ray (2014) demonstrated that, by accounting for the effects of tidal geocentre variations on altimetry-based observations of sea-surface height, the residual variance between ocean-tide models and bottom-pressure observations could be improved by up to 30–40 per cent, with the $O_1$ and $K_1$ tidal harmonics exhibiting the largest variance reduction. For South America, the average vector differences between predicted OTL-induced surface displacements computed using GOT4.10c and GOT4.10 for the $M_2$ harmonic are 0.13, 0.29, and 0.24 mm in the east, north, and vertical components, respectively (denoted by the mean ellipse in panel b of Fig. 19). The average vector differences for the $O_1$ harmonic are 0.26, 0.40, and 0.12 mm in the east, north, and vertical components, respectively (panel d of Fig. 19). Since the geocentre-motion correction is primarily a degree-1 adjustment, the substantial diminishment of the inter-model response differences after the removal of the common-mode component is not surprising. Thus, removing the common-mode factor can effectively eliminate OTL-response discrepancies due to inconsistent reference frames. Third, the development of

CM-referenced orbit and clock products for the GPS processing requires an input ocean-tide model to account for load-induced displacements. We used JPL orbit and clock products version ‘repro2.0’ for our GPS processing, which adopted FES2004 (Lyard *et al.* 2006) as the input ocean-tide model. The altimetric measurements used to constrain FES2004 were not corrected for the effects of tidal geocentre variations (Desai & Ray 2014). To our knowledge, none of the ocean-tide models compared in Section 2.1 were adjusted to account for geocentre motion, and should therefore be internally consistent with the ‘repro2.0’ orbit and clock products. Interestingly, prior to removal of the common-mode component, the residuals between observed and predicted OTL-induced surface displacements are smaller for predictions made using GOT4.10 than for predictions made using GOT4.10c (Fig. 16). Consistent with the results of Desai & Ray (2014), the discrepancies are largest for the $O_1$ tidal harmonic. Furthermore, the discrepancies effectively vanish after removal of the common-mode component, as expected for a reference-frame inconsistency (cf. Fig. 19). Future investigations should explore these matters further and consistently use geocentre-corrected ocean-tide models for both the development of the GPS orbit and clock products (now available in version ‘repro2.1’ of the JPL orbit and clock products) as well as for the computation of the predicted OTL-induced surface displacements.

Regardless of the origin of the common-mode component, however, the removal of the network-mean OTL-induced displacement can significantly reduce the magnitude of the residuals (Figs 12–14). With the size of the residuals reduced, the remaining response appears markedly non-random for the $M_2$ harmonic (Fig. 12),
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Predicted M₂ OTL-induced surface displacements across South America due to the direct differences between the ocean-tide models FES2012 and TPXO8-Atlas, partitioned into six latitudinal bands: (a) 60°–90°N, (b) 30°–60°N, (c) 0°–30°N, (d) 30°–0°S, (e) 60°–30°S and (f) 90°–60°S. Note that the PMEs in panel (f), derived from ocean-tide model discrepancies around the Antarctic continent, are oriented strongly in the north–south direction and exhibit non-diminishing amplitudes as a function of distance from the coast (cf. Fig. 3). The colour bar to the left of panel (e) represents the magnitude of the vector differences between ocean-tide models and the colour bar to the right of panel (f) represents the predicted OTL-induced vertical displacement due to the ocean-tide differences.</figcaption></figure>

suggesting that random GPS measurement errors probably do not account for a majority of the discrepancy between observations and predictions. In contrast, the smaller amplitude harmonics, O₁ and Mf, show less obvious signs of regional spatial coherency and are closer to the level of observational uncertainty due to random errors (Fig. 15). The consistency between forward models, observed most prominently after the removal of the common-mode component, indicates that the sensitivity to errors in the ocean-tide models may no longer be a limiting factor in predicting OTL-induced surface displacements (Figs 12–14).

With the common-mode component removed, the sensitivity to choice of ocean model is reduced to a level comparable with the sensitivity to choice of SNREI Earth model. In both cases, the response differences between forward models are predominantly less than 0.1 mm. Thus, to discern the appropriateness of one model over another, the observed OTL-induced surface displacements would need to be accurate to at least within 0.1 mm and often much better. From comparisons of SNREI Earth models, the discrepancies between predicted OTL-induced displacements for our station network are generally on the order of 0.1 mm for the M₂ tidal harmonic,
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Vector differences between predicted OTL-induced surface displacements computed using the GOT4.10c and GOT4.10 ocean-tide models. The direct differences between the ocean-tide models, shown in the global maps to the left for the $M_2$ (top) and $O_1$ (bottom) tidal harmonics, reflect the influence of tidal geocentre variations on the satellite altimetry measurements of sea-surface height that are used to constrain the ocean-tide models. The elastic surface displacements generated by loading due to the direct differences between GOT4.10c and GOT4.10 are shown in panel (a) for the $M_2$ harmonic and in panel (c) for the $O_1$ harmonic. Panels (b) and (d) show the remaining elastic surface displacements after a common-mode component (mean ellipse outlined by the black box) is removed.</figcaption></figure>

0.01 mm for the $O_1$ tidal harmonic, and 0.001 mm (1 micron!) for the $M_2$ tidal harmonic (Table 2). At this level of sensitivity, supporting or rejecting SNREI Earth models based on measurements of OTL response remains tenuous, given the observational uncertainties. The large-amplitude $M_2$ and, perhaps, $O_1$ OTL responses provide the most promising outlooks at present.

The final set of residuals, between observations and predictions for our station network, exhibits RMS misfits of order 0.1–0.3 mm for each spatial component and tidal harmonic (Figs 16 and 17). Differences between the ocean-tide models or SNREI Earth models probably cannot account for all of the misfit. Indeed, swapping out the various forward models has little effect on the size of the RMS misfits after the common-mode component has been removed. Even the Earth model designed to represent stable continental shield structure, SNA, does not generate significantly better predictions of the OTL-induced surface displacements for South America than the globally averaged models. Other contributors to the misfit include observational uncertainties (e.g. GPS data acquisition, GPS data processing, harmonic analysis) and deficiencies in the forward

model (e.g. numerical integration of eq. 1, spatial variations in sea-water density, deviations from SNREI structure).

We estimated the observational error using two different techniques, which yielded similar results (Table 4). Processing multiple years of data in a single inversion clearly has the potential to improve the accuracy of observed OTL-induced surface displacements, evidenced by the reduction of error for analyses of several years of data compared to analyses of one year of data. In particular, for thirteen stations with long data records (>10 yr), we estimate that the horizontal-displacement errors for the $M_2$ and $O_1$ tidal harmonics are, remarkably, less than 0.1 mm at two standard deviations, which rivals the sensitivity to choice of ocean-tide model and SNREI Earth structure for many of the stations in the South America network.

In general, the residuals derived for the $M_2$ tidal harmonic in South America significantly outweigh the two-sigma observational uncertainties (Fig. 15), leaving open the possibility to invert for structural deficiencies at that frequency. Furthermore, the spatiotemporal characteristics of the $M_2$ residuals remain regionally coherent as well as consistent between forward models (Fig. 12),
indicating that all of the RMS misfit cannot apparently be explained by errors in the observations and adopted forward models. Rather, we suggest that a significant part of the RMS misfit may be due to deviations from our assumed SNREI Earth structure, such as laterally heterogeneous and anelastic material properties (e.g. Latychev *et al.* 2009; Yuan *et al.* 2013; Bos *et al.* 2015).

Deficiencies in the solid Earth body tide (SEBT) model, removed at the GPS processing stage, could also contribute to the residual displacements; however, solid Earth body tides operate at very long (global) wavelengths and would thus produce residuals coherent across similarly large scales. Such long-wavelength coherency in the residuals could consequently contribute to the common-mode component for a regional GPS network. The removal of a uniform-displacement component from residual OTL-induced surface displacements at each station could therefore eliminate information about the long-wavelength deficiencies in the SEBT model. Notwithstanding, a significant contribution to the common-mode component appears to be derived from inaccuracies in the ocean-tide models as well as inconsistencies in reference frames.

For the $O_1$ tidal harmonic, the residuals are more comparable to the noise levels, but nevertheless appear to exceed slightly the level of uncertainty at many of the South America stations and particularly in the horizontal components (Figs 15c and d, Tables 4 and 3). Fig. 13 shows consistency between forward models and some evidence for regional spatial coherency. In particular, some of the ellipses are oriented such that the semi-major axis points toward the maximum load amplitude. For the $M_f$ tidal harmonic, the residuals are close to or below the level of observational uncertainty at most stations and therefore less conclusive (Fig. 15). The observed OTL-induced surface displacements, however, exhibit spatial coherency as well as a coherent northward offset in the vertical component of the displacements relative to the predictions (Figs 11e and f).

## 5 SUMMARY AND CONCLUSIONS

We derived OTL-induced surface displacements from sub-daily GPS time-series using kinematic PPP and harmonic analysis for a network of stations in Brazil, Argentina, and Uruguay. After a common-mode component representing the network-mean OTL response is removed from each station, the misfits between the observed and predicted OTL-induced displacements are approximately 0.1–0.3 mm for the dominant tidal harmonics in three distinct frequency bands: $M_2$, $O_1$ and $M_f$. An assessment of the observational error suggests that, at least for the $M_2$ harmonic and sufficiently long time-series, the OTL-response residuals exceed random noise from the data processing. Therefore, OTL-response residuals may potentially be used to refine ocean-tide models and to constrain solid-Earth structure.

Comparisons of forward-modelled predictions for South America suggest that the sensitivity to choice of ocean-tide model still generally exceeds the sensitivity to choice of SNREI reference Earth model, albeit not substantially. Removal of the common-mode component across the network substantially reduces the discrepancies between ocean-tide models. We suggest that possible sources for the common-mode component might include deficiencies in the ocean-tide models at high latitude, and reference-frame inconsistencies. Based on comparisons of selected ocean-tide models and SNREI Earth models, RMS differences between predictions of OTL-induced displacements are at the sub-mm level for each harmonic considered, and often much less than 0.1 mm.

We find evidence for regional spatial coherency in the residuals between observed and predicted OTL-induced displacements

that remain consistent for a variety of ocean-tide and SNREI-Earth model combinations. We postulate that part of the spatial coherency could be attributed to deficiencies in the *a priori* Earth model, which would undoubtedly include deviations from an assumed SNREI structure (e.g. anelasticity, anisotropy, and lateral heterogeneities). As ocean-tide models, OTL-response modelling, and data processing methods continue to improve, the ability to probe Earth structure through observations of OTL-induced surface displacements becomes increasingly tractable.

## ACKNOWLEDGEMENTS

We are indebted to Duncan Agnew, Richard Ray, Matt King and an anonymous reviewer for insightful ideas and valuable critiques that greatly improved our manuscript. We also sincerely thank Shailen Desai for helpful discussions on OTL analysis and tidal geocentre variations as well as Angelyn Moore and Willy Bertiger for providing ongoing GIPSY support. Dan Bower graciously supplied scripts to extract and plot the seismic tomography data shown in Supporting Information Fig. S5. The GPS data used in our study was made available by the governments of Brazil (Instituto Brasileiro de Geografia e Estatística), Argentina (Instituto Geográfico Nacional), and Uruguay (Servicio Geográfico Militar). We used geographic information from the Scientific Committee on Antarctic Research (SCAR) Antarctic Digital Database (ADD) to develop a land-sea mask around Antarctica. We gratefully acknowledge support from the National Science Foundation Geophysics Program funding under Grant No. EAR-1417245. This manuscript is based upon work supported by the NASA Earth and Space Science Fellowship to HRM under Grant No. NNX14AO04H. Some figures were generated using Generic Mapping Tools (Wessel *et al.* 2013). Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

## REFERENCES

* Agnew, D.C., 1997. NLOADF: a program for computing ocean-tide loading, *J. geophys. Res.*, **102**(B3), 5109–5110.

* Agnew, D.C., 2012. SPOTL: Some programs for ocean-tide loading, Scripps Institute of Oceanography Technical Report.

* Agnew, D.C., 2015. Earth Tides, in *Treatise on Geophysics*, 2nd edn, vol. 3, pp. 151–178, ed. Schubert, G., Elsevier.

* Allinson, C., Clarke, P., Edwards, S., King, M., Baker, T. & Cruddace, P., 2004. Stability of direct GPS estimates of ocean tide loading, *Geophys. Res. Lett.*, **31**(15), L15603. doi:10.1029/2004GL020588.

* Altamimi, Z., Collilieux, X. & Métivier, L., 2011. ITRF2008: An improved solution of the International Terrestrial Reference Frame, *J. Geod.*, **85**(8), 457–473.

* Amante, C. & Eakins, B.W., 2009. *ETOPO1 1 Arc-minute Global Relief Model: Procedures, Data Sources and Analysis*, US Department of Commerce.

* Aster, R.C., Borchers, B. & Thurber, C.H., 2013. *Parameter Estimation and Inverse Problems*, Academic Press.

* Baker, T., 1980. Tidal tilt at Llanrwst, North Wales: tidal loading and Earth structure, *Geophys. J. Int.*, **62**(2), 269–290.

* Baker, T., 1984. Tidal deformations of the Earth, *Sci. Prog.*, **69**, 197–233.

* Baker, T., Curtis, D. & Dodson, A., 1996. A new test of Earth tide models in central Europe, *Geophys. Res. Lett.*, **23**(24), 3559–3562.

* Baker, T.F., Edge, R. & Jeffries, G., 1991. Tidal gravity and ocean tide loading in Europe, *Geophys. J. Int.*, **107**(1), 1–11.

* Bar-Sever, Y.E., Kroger, P.M. & Borjesson, J.A., 1998. Estimating horizontal gradients of tropospheric path delay with a single GPS receiver, *J. geophys. Res.*, **103**(B3), 5019–5035.
Bertiger, W., Desai, S.D., Haines, B., Harvey, N., Moore, A.W., Owen, S. & Weiss, J.P., 2010. Single receiver phase ambiguity resolution with GPS data, *J. Geod.*, **84**(5), 327–337.

Blewitt, G., 2003. Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth, *J. geophys. Res.*, **108**(B2), 2103, doi:10.1029/2002JB002082.

Blewitt, G., 2015. GPS and space-based geodetic methods, in *Treatise on Geophysics*, 2nd edn, vol. 3, pp. 307–338, ed. Schubert, G., Elsevier.

Boehm, J., Werl, B. & Schuh, H., 2006. Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data, *J. geophys. Res.*, **111**(B2), doi:10.1029/2005JB003629.

Bos, M.S. & Baker, T.F., 2005. An estimate of the errors in gravity ocean tide loading computations, *J. Geod.*, **79**(1), 50–63.

Bos, M.S., Penna, N.T., Baker, T.F. & Clarke, P.J., 2015. Ocean tide loading displacements in western Europe: Part 2. GPS-observed anelastic dispersion in the asthenosphere, *J. geophys. Res.*, **120**(9), 6540–6557.

Carrère, L., Lyard, F., Cancet, M., Guillot, A. & Roblou, L., 2012. FES2012: A new global tidal model taking advantage of nearly 20 years of altimetry, in *Proceedings of meeting "20 Years of Altimetry"*, Venice.

Cartwright, D. & Edden, A., 1973. Corrected tables of tidal harmonics, *Geophys. J. R. astr. Soc.*, **33**, 253–264.

Cartwright, D. & Taylor, R., 1971. New computations of the tide-generating potential, *Geophys. J. R. astr. Soc.*, **23**, 45–74.

Dach, R. & Dietrich, R., 2000. Influence of the ocean loading effect on GPS derived precipitable water vapor, *Geophys. Res. Lett.*, **27**(18), 2953–2956.

Desai, S.D. & Ray, R.D., 2014. Consideration of tidal variations in the geocenter on satellite altimeter observations of ocean tides, *Geophys. Res. Lett.*, **41**(7), 2454–2459.

Desai, S.D., Bertiger, W. & Haines, B.J., 2014. Self-consistent treatment of tidal variations in the geocenter for precise orbit determination, *J. Geod.*, **88**(8), 735–747.

Doodson, A., 1957. The analysis and prediction of tides in shallow water, *Int. Hydrogr. Rev.*, **41**, 85–126.

Doodson, A. & Warburg, H., 1941. *Admiralty Manual of Tides*, HMSO.

Dragert, H., James, T. & Lambert, A., 2000. Ocean loading corrections for continuous GPS: A case study at the Canadian coastal site Holberg, *Geophys. Res. Lett.*, **27**(14), 2045–2048.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**(4), 297–356.

Egbert, G.D. & Erofeeva, S.Y., 2002. Efficient inverse modeling of barotropic ocean tides, *J. Atmos. Ocean. Technol.*, **19**(2), 183–204.

Egbert, G.D., Erofeeva, S.Y. & Ray, R.D., 2010. Assimilation of altimetry data for nonlinear shallow-water tides: Quarter-diurnal tides of the Northwest European Shelf, *Cont. Shelf Res.*, **30**(6), 668–679.

Elosegui, P., Davis, J., Johansson, J. & Shapiro, I., 1996. Detection of transient motions with the Global Positioning System, *J. geophys. Res.*, **101**(B5), 11249–11261.

Farrell, W., 1973. Earth tides, ocean tides and tidal loading, *Phil. Trans. R. Soc. Lond.*, **A**, **274**(1239), 253–259.

Foreman, M., 1977. Manual for tidal heights analysis and prediction, *Pac. Mar. Sci. Rep.*, **77–10**, 1–58.

Foreman, M., Cherniawsky, J. & Ballantyne, V., 2009. Versatile harmonic tidal analysis: improvements and applications, *J. Atmos. Ocean. Technol.*, **26**(4), 806–817.

Francis, O. & Mazzega, P., 1990. Global charts of ocean tide loading effects, *J. geophys. Res.*, **95**(C7), 11411–11424.

Fu, Y., Freymueller, J.T. & van Dam, T., 2012. The effect of using inconsistent ocean tidal loading models on GPS coordinate solutions, *J. Geod.*, **86**(6), 409–421.

Goad, C.C., 1980. Gravimetric tidal loading computed from integrated Green's functions, *J. geophys. Res.*, **85**(B5), 2679–2683.

Godin, G., 1972. *The Analysis of Tides*, University of Toronto Press.

Grand, S.P. & Helmberger, D.V., 1984. Upper mantle shear structure of North America, *Geophys. J. Int.*, **76**(2), 399–438.

Ishii, M. & Tromp, J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth's mantle, *Science*, **285**(5431), 1231–1236.

Ito, T. & Simons, M., 2011. Probing asthenospheric density, temperature, and elastic moduli below the western United States, *Science*, **332**(6032), 947–951.

Jentzsch, G., 1997. Earth tides and ocean tidal loading, in *Tidal Phenomena*, pp. 145–171, eds Wilhelm, H., Zrn, W. & Wenzel, H.-G., Springer.

Kennett, B., Engdahl, E. & Buland, R., 1995. Constraints on seismic velocities in the Earth from traveltimes, *Geophys. J. Int.*, **122**(1), 108–124.

Khan, S.A. & Scherneck, H.-G., 2003. The M2 ocean tide loading wave in Alaska: vertical and horizontal displacements, modelled and observed, *J. Geod.*, **77**(3–4), 117–127.

Khan, S.A. & Tscherning, C.C., 2001. Determination of semi-diurnal ocean tide loading constituents using GPS in Alaska, *Geophys. Res. Lett.*, **28**(11), 2249–2252.

King, M., 2006. Kinematic and static GPS techniques for estimating tidal displacements with application to Antarctica, *J. Geodyn.*, **41**(1), 77–86.

King, M. & Aoki, S., 2003. Tidal observations on floating ice using a single GPS receiver, *Geophys. Res. Lett.*, **30**(3), doi:10.1029/2002GL016182.

King, M.A., Watson, C.S., Penna, N.T. & Clarke, P.J., 2008. Subdaily signals in GPS observations and their effect at semiannual and annual periods, *Geophys. Res. Lett.*, **35**(3), doi:10.1029/2007GL032252.

Kustowski, B., Ekström, G. & Dziewoński, A., 2008. Anisotropic shear-wave velocity structure of the Earth's mantle: a global model, *J. geophys. Res.*, **113**, B06306, doi:10.1029/2007JB005169.

Larson, K.M., Cervelli, P., Lisowski, M., Miklius, A., Segall, P. & Owen, S., 2001. Volcano monitoring using the Global Positioning System: filtering strategies, *J. geophys. Res.*, **106**(B9), 19453–19464.

Larson, K.M., Poland, M. & Miklius, A., 2010. Volcano monitoring using GPS: developing data analysis strategies based on the June 2007 Kīlauea Volcano intrusion and eruption, *J. geophys. Res.*, **115**(B7), doi:10.1029/2009JB007022.

Latychev, K., Mitrovica, J.X., Ishii, M., Chan, N.-H. & Davis, J.L., 2009. Body tides on a 3-D elastic earth: Toward a tidal tomography, *Earth planet. Sci. Lett.*, **277**(1–2), 86–90.

Lin, F.-C., Schmandt, B. & Tsai, V.C., 2012. Joint inversion of Rayleigh wave phase velocity and ellipticity using USArray: Constraining velocity and density structure in the upper crust, *Geophys. Res. Lett.*, **39**(12), L12303, doi:10.1029/2012GL052196.

Love, A., 1909. The yielding of the Earth to disturbing forces, *Proc. R. Soc. A*, **82**(551), 73–88.

Lyard, F., Lefevre, F., Letellier, T. & Francis, O., 2006. Modelling the global ocean tides: Modern insights from FES2004, *Ocean Dyn.*, **56**(5–6), 394–415.

Melchior, P., 1983. *The Tides of the Planet Earth*, 2nd edn, Pergamon Press.

Montagner, J.-P. & Kennett, B., 1996. How to reconcile body-wave and normal-mode reference Earth models, *Geophys. J. Int.*, **125**(1), 229–248.

Parker, B.B., 2007. *Tidal Analysis and Prediction*, vol. NOAA Special Publication: NOS CO-OPS 3, US Department of Commerce, National Oceanic and Atmospheric Administration, National Ocean Service, Center for Operational Oceanographic Products and Services.

Pawlowicz, R., Beardsley, B. & Lentz, S., 2002. Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE, *Comput. Geosci.*, **28**(8), 929–937.

Penna, N.T., Bos, M.S., Baker, T.F. & Scherneck, H.-G., 2008. Assessing the accuracy of predicted ocean tide loading displacement values, *J. Geod.*, **82**(12), 893–907.

Penna, N.T., Clarke, P.J., Bos, M.S. & Baker, T.F., 2015. Ocean tide loading displacements in western Europe: Part 1. Validation of kinematic GPS estimates, *J. geophys. Res.*, **120**(9), 6523–6539.

Petit, G. & Luzum, B., 2010. *IERS Technical Note No. 36, IERS Conventions (2010)*, International Earth Rotation and Reference Systems Service: Frankfurt, Germany.

Petrov, L. & Ma, C., 2003. Study of harmonic site position variations determined by very long baseline interferometry, *J. geophys. Res.*, **108**(B4), doi:10.1029/2002JB001801.

Pugh, D., 1987. *Tides, Surges and Mean Sea-Level*, John Wiley & Sons Ltd.

Pugh, D. & Woodworth, P., 2014. *Sea-level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-level Changes*, Cambridge Univ. Press.
Pugh, D.T., Woodworth, PL. & Bos, M.S., 2011. Lunar tides in Loch Ness, Scotland, *J. geophys. Res.*, **116**(C11), doi:10.1029/2011JC007411.

Ray, R., 2013. Precise comparisons of bottom-pressure and altimetric ocean tides, *J. geophys. Res.*, **118**(9), 4570–4584.

Ray, R.D., 1999. A global ocean tide model from TOPEX/POSEIDON altimetry: GOT99. 2., NASA Technical Memorandum 209478.

Rebischung, P., Griffiths, J., Ray, J., Schmid, R., Collilieux, X. & Garayt, B., 2012. IGS08: The IGS realization of ITRF2008, *GPS Solut.*, **16**(4), 483–494.

Savcenko, R. & Bosch, W., 2012. EOT11a–empirical ocean tide model from multi-mission satellite altimetry, DGFI Report, 89.

Schenewerk, M.S., Marshall, J. & Dillinger, W., 2001. Vertical ocean-loading deformations derived from a global GPS network., *J. Geod. Soc. Japan*, **47**(1), 237–242.

Stammer, D. et al., 2014. Accuracy assessment of global barotropic ocean tide models, *Rev. Geophys.*, **52**(3), 243–282.

Thomas, I.D., King, M.A. & Clarke, P.J., 2007. A comparison of GPS, VLBI and model estimates of ocean tide loading displacements, *J. Geod.*, **81**(5), 359–368.

Vey, S. et al., 2002. GPS measurements of ocean loading and its impact on zenith tropospheric delay estimates: a case study in Brittany, France, *J. Geod.*, **76**(8), 419–427.

Wang, H., Xiang, L., Jia, L., Jiang, L., Wang, Z., Hu, B. & Gao, P., 2012. Load Love numbers and Green's functions for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0, *Comput. Geosci.*, **49**, 190–199.

Wessel, P. & Smith, W.H., 1996. A global, self-consistent, hierarchical, high-resolution shoreline database, *J. geophys. Res.*, **101**(B4), 8741–8743.

Wessel, P., Smith, W.H., Scharroo, R., Luis, J. & Wobbe, F., 2013. Generic mapping tools: Improved version released, *EOS, Trans. Am. geophys. Un.*, **94**(45), 409–410.

Wu, X., Collilieux, X., Altamimi, Z., Vermeersen, B., Gross, R. & Fukumori, I., 2011. Accuracy of the International Terrestrial Reference Frame origin and earth expansion, *Geophys. Res. Lett.*, **38**(13), doi:10.1029/2011GL047450.

Wu, X., Ray, J. & van Dam, T., 2012. Geocenter motion and its geodetic and geophysical implications, *J. Geodyn.*, **58**, 44–61.

Yuan, L. & Chao, B.F., 2012. Analysis of tidal signals in surface displacement measured by a dense continuous GPS array, *Earth planet. Sci. Lett.*, **355**, 255–261.

Yuan, L., Chao, B.F., Ding, X. & Zhong, P., 2013. The tidal displacement field at Earth's surface determined using global GPS observations, *J. geophys. Res.*, **118**, 2618–2632.

Zahel, W., 1997. Ocean Tides, in *Tidal Phenomena*, pp. 113–143, eds Wilhelm, H., Zürn, W. & Wenzel, H.-G., Springer.

Zumberge, J., Heflin, M., Jefferson, D., Watkins, M. & Webb, F., 1997. Precise point positioning for the efficient and robust analysis of GPS data from large networks, *J. geophys. Res.*, **102**(B3), 5005–5017.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** Timeline of station activity for the first 40 receivers in the GPS network across Brazil, Argentina and Uruguay.

**Figure S2.** Continued timeline of station activity for the next 40 receivers.

**Figure S3.** Continued timeline of station activity for the next 40 receivers.

**Figure S4.** Continued timeline of station activity for the next 40 receivers.

**Figure S5.** Seismic tomography profile depicting the Amazonian Craton in South America.

**Figure S6.** Residual OTL-induced displacements shown as vectors.

**Figure S7.** Residual OTL-induced displacements shown as vectors, but with the common-mode component (network-mean OTL-induced displacement) removed.

**Figure S8.** Residuals between observed and predicted OTL-induced surface displacements for the $M_2$ tidal harmonic, without the common-mode component removed.

**Figure S9.** Residuals between observed and predicted OTL-induced surface displacements for the $O_1$ tidal harmonic, without the common-mode component removed.

**Figure S10.** Residuals between observed and predicted OTL-induced surface displacements for the $M_f$ tidal harmonic, without the common-mode component removed.

**Figure S11.** Predicted $M_2$ OTL-induced surface displacements in South America due to the direct differences between the ocean tide models FES2012 and GOT4.10, partitioned into six latitudinal bands.

**Figure S12.** Predicted $M_2$ OTL-induced surface displacements in South America due to the direct differences between the ocean tide models FES2012 and EOT11A, partitioned into six latitudinal bands.

**Figure S13.** A periodogram of time series residuals for the vertical displacement component of station RIO2.

**Data set S1.** Station names, geographic locations, number of operational days, and distance from the coast for the GPS receivers in the South America network.

**Data set S2.** Radial profiles of P-wave velocity, S-wave velocity, density, bulk modulus, and shear modulus for the isotropic and oceanless version of the Preliminary Reference Earth Model (PREM) used in our study.

**Data set S3.** Same as Data set S2, but for reference model STW105.

**Data set S4.** Same as Data set S2, but for reference model AK135f.

**Data set S5.** Same as Data set S2, but for model SNA.

**Data set S6.** Displacement load Green's functions (LGFs) for PREM in the CM reference frame.

**Data set S7.** Displacement LGFs for STW105 in the CM reference frame.

**Data set S8.** Displacement LGFs for AK135f in the CM reference frame.

**Data set S9.** Displacement LGFs for SNA in the CM reference frame.

**Data set S10.** A compressed file containing the observed OTL-induced surface displacements for GPS receivers in the South America network.

**Data set S11.** A compressed file containing the predicted OTL-induced surface displacements for multiple ocean tide and SNREI Earth model combinations.

**Data set S12.** A table of the non-linear shallow-water tides considered in our analysis.

(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggw087/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## APPENDIX A: PROCESS NOISE SETTINGS FOR GPS ANALYSIS

We performed basic synthetic tests to inform our selection of appropriate coordinate and tropospheric process noise settings for the
GPS data processing. We generated the synthetic signals within GIPSY by differencing pre-fit residuals computed using two different OTL-response models: (1) the OTL-response model for a coastal station in Brazil with a large offshore tide ('master site') and (2) the OTL-response model at a separate station in the network ('test site'). We selected six test sites for our analysis, each exhibiting small OTL-response amplitudes relative to the master site. The displacement-response signal at the master site has a vertical amplitude of 3.69 cm for the M₂ harmonic, whereas the M₂ displacement-response amplitudes at the test sites do not exceed 1 cm.

We then added the differenced pre-fit-residuals to the raw GPS data of each test site as modifications to the pseudorange and phase observables for each transmitter-receiver pair in the GIPSY quick-measurement files. The raw GPS data contain information about the actual OTL response at each test site; thus, when the differenced pre-fit residuals (i.e. master-site residuals — test-site residuals) are added to the original data, the revised data should contain only the OTL-response signal for the master site, along with any pre-existing noise and non-tidal signals. Note that we deliberately selected test sites with small OTL-response amplitudes, since errors in our prediction of the tidal response signals at the test sites can bias our attempts to recover the synthetic tidal signal (i.e. the predicted OTL response at the master site). All synthetic tests used up to one year of data from 2010 at each of the six test stations in South America.

To guide our exploration of coordinate process noise values to test, we first computed a range of theoretically suitable values. We estimated the instantaneous velocity, $v_{\text{inst}}$, experienced by any given station as:

$$v_{\text{inst}} = \frac{d}{dt}(A \cos \omega t) = -A \omega \sin \omega t, \quad (A1)$$

where A is the signal amplitude, $\omega$ is the frequency of the signal, and t is time. The absolute value of the instantaneous velocity given by eq. (A1) is maximized when the quantity $\sin \omega t = \pm 1$, which occurs at the maximum slope of the harmonic wave. For the M₂ tidal harmonic, $\omega = 0.5059$ rad hr⁻¹ and the response amplitudes in South America reach as high as ~4 cm. Thus, the maximum instantaneous velocity, $v_{\text{max}}$, expected for our network is ~20 mm hr⁻¹.

To parametrize the coordinate process noise within GIPSY, one must specify the variance per unit time of the allowed site displacement (or, more specifically, the square root of the variance per unit time). Constraints that are too strict will bias the solutions, whereas constraints that are too loose will retain large data outliers. Following Elosegui et al. (1996), the variance per unit time, $\sigma_{\text{RW}}^2$, may be related to the maximum site velocity, $v_{\text{max}}$, and the time between solution epochs, $\Delta t$, by a constant of proportionality, $\xi_{\text{RW}}$, known as the dynamic resolution parameter:

$$\xi_{\text{RW}} = \frac{\sigma_{\text{RW}}}{v_{\text{max}} \sqrt{\Delta t}}. \quad (A2)$$

For $\xi_{\text{RW}} \gg 1$, solutions will be weakly constrained; for $\xi_{\text{RW}} \ll 1$, solutions will be tightly constrained. Elosegui et al. (1996) opted for a 'standard', yet somewhat arbitrary, value of $\xi_{\text{RW}} = 10$. For $v_{\text{max}} = 20$ mm hr⁻¹ and $\Delta t = 300$ s between solution epochs, a plausible range of random-walk standard deviations to test should include $1.9 \times 10^{-7} < \sigma_{\text{RW}} < 1.9 \times 10^{-6}$ km s⁻¹, corresponding to dynamic resolution parameters in the range of $2 < \xi_{\text{RW}} < 20$.

The root-mean-square (RMS) misfits between the recovered and synthetic signals, averaged across all six test stations, are shown in

Fig. A1. We examined six settings for the coordinate process noise: $1.0 \times 10^{-8}$, $1.0 \times 10^{-7}$, $3.0 \times 10^{-7}$, $5.0 \times 10^{-7}$, $7.0 \times 10^{-7}$ and $1.0 \times 10^{-3}$ km s⁻¹. In each case, we held the process noise setting for the tropospheric zenith delay fixed at $5.0 \times 10^{-8}$ km s⁻¹, which is the GIPSY-recommended value for slow-moving platforms. We find that the tightest random-walk constraint we considered, $1.0 \times 10^{-8}$ km s⁻¹, severely dampens the amplitudes of the station-position estimates and of the recovered tidal signals. The constraint of $1.0 \times 10^{-7}$ km s⁻¹ also somewhat overdamped the recovered tidal response, particularly for the vertical component of the M₂ harmonic. The exceptionally loose random-walk constraint of $1.0 \times 10^{-3}$ km s⁻¹ produced good amplitude recovery of the synthetic tidal response, but also limited the constraint on the noise. Tests of intermediary values revealed that a coordinate process noise value of $5.0 \times 10^{-7}$ km s⁻¹ recovers both the amplitude and the phase of the synthetic tidal response to relatively high precision. Loosening the parameter does not significantly affect the ability to resolve the synthetic signals, even compared with very loose parametrizations such as $1.0 \times 10^{-3}$ km s⁻¹, yet any stricter constraints begin to bias the recovered signal towards zero. Furthermore, the preferred parametrization of $5.0 \times 10^{-7}$ km s⁻¹ is suitable for even the largest OTL-induced surface displacements observed in South America, yet strict enough to limit noise substantially.

For a random-walk coordinate setting of $\sigma_{\text{RW}} = 5.0 \times 10^{-7}$ km s⁻¹ and for $\Delta t = 300$ s between solution epochs, we derive a dynamic resolution parameter of:

$$\xi_{\text{RW}} = \frac{\sigma_{\text{RW}}}{v_{\text{max}} \sqrt{\Delta t}} = \frac{5.0 \times 10^{-7}}{\frac{20 \times 10^{-6}}{3600} \sqrt{300}} \approx 5. \quad (A3)$$

Our result of $\xi_{\text{RW}} = 5$ from the synthetic tests is very near to the value adopted by Elosegui et al. (1996), albeit with a slightly tighter constraint on the solution to limit noise contamination, thus lending support to the validity of our preferred random-walk parametrization. The selection of a relatively tight random-walk constraint also mitigates the effects of multipath (e.g. Larson et al. 2010).

In comparison, Penna et al. (2015) derived an optimal (minimum) coordinate process noise setting of $3.2 \times 10^{-6}$ km s⁻¹, which is somewhat looser than $5.0 \times 10^{-7}$ km s⁻¹. We have, however, constructed the synthetic tests to explore the ability of the GPS processing to retain the full OTL-response signals, as opposed to the residuals after OTL is removed, and we find that $5.0 \times 10^{-7}$ km s⁻¹ is sufficiently loose to recover even the largest amplitude OTL responses predicted for our station network.

In addition to the coordinate process noise setting, recovering the OTL-response signal also depends largely on the ability to account for propagation delays in the transmitted carrier wave signals through the troposphere (e.g. Bar-Sever et al. 1998; Dach & Dietrich 2000; Dragert et al. 2000; Vey et al. 2002; Khan & Scherneck 2003; Larson et al. 2010; Penna et al. 2015). The tropospheric wet delay (non-hydrostatic) arises from the interaction between the electromagnetic (EM) carrier wave signal and the static dipole moment of water molecules in the atmosphere (e.g. Blewitt 2015). The tropospheric dry delay (hydrostatic) refers to the dynamic dipole moment induced on all atmospheric molecules, including water, by the propagating EM wave. In our analysis, we account for both types of delays.

We find that the tropospheric zenith delays are best estimated stochastically along with the station coordinates in a single kinematic run. We therefore explored the effects of varying the tropospheric process noise on the ability to recover the synthetic tidal
<figure><img src="image_23.png" /><figcaption>**Figure A1.** Root-mean-square (RMS) misfits between recovered and synthetic OTL-induced surface displacements. Here, we compare coordinate process noise settings used to estimate receiver positions during the kinematic PPP GPS data processing.</figcaption></figure>

response. In each test, we supplied the kinematic runs with initial tropospheric zenith delay estimates from ECMWF. With the coordinate process noise held fixed at $5.0 \times 10^{-7}$ km s$^{-1/2}$, we examined five settings for the tropospheric process noise: $1.0 \times 10^{-8}$, $2.5 \times 10^{-8}$, $5.0 \times 10^{-8}$, $1.0 \times 10^{-7}$ and $2.0 \times 10^{-7}$ km s$^{-1/2}$. The RMS misfits between the recovered and synthetic OTL-response signals for each process noise setting are depicted in Fig. A2.

The vertical-displacement component clearly exhibits greater sensitivity to tropospheric process noise than the horizontal-displacement components. For the $M_2$ harmonic, the RMS misfits in the vertical-displacement component are minimized for a tropospheric process noise setting of $2.5 \times 10^{-8}$ km s$^{-1/2}$, followed closely by a setting of $5.0 \times 10^{-8}$ km s$^{-1/2}$. For the $O_1$ tidal harmonic, looser constraints of $1.0 \times 10^{-7}$ and $2.0 \times 10^{-7}$ km s$^{-1/2}$ seem more suitable, although a setting of $5.0 \times 10^{-8}$ km s$^{-1/2}$ yields similar results. The RMS misfits for the $M_f$ harmonic do not vary substantially with changes to the tropospheric process noise.

Consistent with the results of Penna *et al.* (2015), we do not find significant leakage of the synthetic signal between the spatial components. The synthetic signal applied to the $M_2$ harmonic band, for example, had an amplitude of 36.9, 4.4, and 7.6 mm in the vertical, east, and north components, respectively. Given that the vector differences between the recovered and synthetic signals are submillimetre in each component, we infer that nearly all of the input signal for a particular coordinate maps directly into the recovered signal for the same coordinate.

Based on the six test sites that we considered, the optimal setting for the tropospheric process noise is not sharply defined, but a

setting of $5.0 \times 10^{-8}$ km s$^{-1/2}$ performs well overall and corresponds to the value recommended by GIPSY for slow-moving objects. We therefore adopt a tropospheric process noise setting of $5.0 \times 10^{-8}$ km s$^{-1/2}$ for our analysis, recognizing that the truly optimal value will likely differ between individual tidal harmonics and various geographic locations. The optimal tropospheric process noise value found by Penna *et al.* (2015), $1.0 \times 10^{-7}$ km s$^{-1/2}$, differs from $5.0 \times 10^{-8}$ km s$^{-1/2}$ by only a factor of two, which might be due in part to different climatic settings between western Europe and South America. The two parametrizations do not, however, yield appreciably different RMS misfits in our analysis for South America (Fig. A2).

The smallest RMS misfits from Figs A1 and A2 are generally consistent with our $2\sigma$ uncertainty estimates from Table 4 for one year of data or less. It should be recalled, however, that our synthetic tests were developed by adding modelled OTL-response signals to real data in the dominant tidal bands, and therefore should not be used to make strict assertions about uncertainties in tidal-response estimation from GPS data. Notably, any errors in the predicted OTL response at a particular test site will be contained within our RMS estimates of signal recovery, and therefore compound the uncertainties derived from noise in the time-series. Although injecting the synthetic signal into a non-dominant tidal band could reduce the prediction errors considerably and thereby allow for better estimation of the uncertainties in GPS estimates of OTL response using synthetic tests (Penna *et al.* 2015), the errors should remain consistent between each synthetic test and thus should not significantly bias our selection of suitable process noise parameters.
<figure><img src="image_24.png" /><figcaption>Figure A2. RMS misfits between recovered and synthetic OTL-induced surface displacements. Here, we compare a selection of tropospheric process noise settings.</figcaption></figure>

## APPENDIX B: HARMONIC ANALYSIS PROCEDURE

The equation representing the model fit, $Z(t)$, to the GPS-inferred displacement time-series is given by

$$Z(t) = m_0 + m_1 t + \sum_{n=1}^{N} A_n f_n(t) \cos(V_n(t) + u_n(t) - \phi_n), \quad (B1)$$

where $m_0$ is a constant-offset term, $m_1$ is a linear-trend term, $f_n$ represents the harmonic-modulation correction factor for the amplitude, $u_n$ represents the harmonic-modulation correction factor for the phase, $n$ represents a particular tidal harmonic, $N$ represents the total number of tidal harmonics used in the model, and $t$ is time.

We seek the set of model terms [$m_0, m_1, A_n, \phi_n$] by minimizing the misfit between data and model using an iterative re-weighted least-squares (IRLS) approach. To perform the inversion, we first separate the harmonic portion of eq. (B1) into in-phase and quadrature components:

$$ \begin{aligned} A_n f_n(t) \cos(V_n(t) + u_n(t) - \phi_n) \\ = A_n \cos(\phi_n) f_n(t) \cos(V_n(t) + u_n(t)) \\ + A_n \sin(\phi_n) f_n(t) \sin(V_n(t) + u_n(t)). \end{aligned} \quad (B2) $$

Defining $c_n \equiv A_n \cos(\phi_n)$ and $s_n \equiv A_n \sin(\phi_n)$, eq. (B2) becomes:

$$ \begin{aligned} A_n f_n(t) \cos(V_n(t) + u_n(t) - \phi_n) \\ = c_n f_n(t) \cos(V_n(t) + u_n(t)) \\ + s_n f_n(t) \sin(V_n(t) + u_n(t)). \end{aligned} \quad (B3) $$

The amplitude and phase modulation factors, $f_n(t)$ and $u_n(t)$, as well as the astronomical argument, $V_n(t)$, are known functions derived from the astronomical ephemeris (e.g. Foreman 1977). We incorporate harmonic-modulation correction factors directly into the inversion framework (Foreman et al. 2009), rather than apply constant and thus approximate correction terms at the post-processing stage, thereby allowing multiple years of data to be processed in a single analysis.

Substituting eq. (B3) back into the full formula for a tidal signal, eq. (B1), we have:

$$ Z(t) = m_0 + m_1 t + \sum_{n=1}^{N} [c_n f_n(t) \cos(V_n(t) + u_n(t)) \\ + s_n f_n(t) \sin(V_n(t) + u_n(t))]. \quad (B4) $$

Eq. (B4) may now be used to invert real tidal data for the unknown model parameters (i.e. $m_0, m_1, c_n$ and $s_n$). Note that if $Z_0$ (mean sea level) is included as a tidal constituent, then $m_0$ must be removed from the model parameters to avoid ill-conditioning in the matrix inversion.

We develop a system of linear equations of the form:

$$ G m = d, \quad (B5) $$

where $d$ is the observed tidal data, $m$ is a vector of model parameters, and $G$ is a matrix of known quantities that interact with the model parameters. We aim to match the model, $Gm$, with the observed data, $d$, by initially minimizing the norm of the squared residuals, and then iterating on the solution residuals. The data and model vectors are given by

$$ d = [d(t_0) \ d(t_1) \ d(t_2) \ d(t_3) \ d(t_4) \ d(t_5) \ d(t_6) \ d(t_7) \ \dots]^T \quad (B6) $$
and

$$m = [m_0 \ m_1 \ c_1 \ s_1 \ c_2 \ s_2 \ c_3 \ s_3 \ c_4 \ s_4 \ \dots \ c_N \ s_N]^T, \quad (B7)$$

respectively. The subscripts for c and s represent individual tidal harmonics (e.g. M₂).

The $G^T$ matrix (transpose of G) is given by

$$G^T =$$

$$ \begin{bmatrix} 1 & 1 & 1 & \cdots \\ t_0 & t_1 & t_2 & \cdots \\ f_1(t_0) \cos[V_1(t_0) + u_1(t_0)] & f_1(t_1) \cos[V_1(t_1) + u_1(t_1)] & \cdots & \cdots \\ f_1(t_0) \sin[V_1(t_0) + u_1(t_0)] & f_1(t_1) \sin[V_1(t_1) + u_1(t_1)] & \cdots & \cdots \\ f_2(t_0) \cos[V_2(t_0) + u_2(t_0)] & f_2(t_1) \cos[V_2(t_1) + u_2(t_1)] & \cdots & \cdots \\ f_2(t_0) \sin[V_2(t_0) + u_2(t_0)] & f_2(t_1) \sin[V_2(t_1) + u_2(t_1)] & \cdots & \cdots \\ f_3(t_0) \cos[V_3(t_0) + u_3(t_0)] & f_3(t_1) \cos[V_3(t_1) + u_3(t_1)] & \cdots & \cdots \\ f_3(t_0) \sin[V_3(t_0) + u_3(t_0)] & f_3(t_1) \sin[V_3(t_1) + u_3(t_1)] & \cdots & \cdots \\ f_4(t_0) \cos[V_4(t_0) + u_4(t_0)] & f_4(t_1) \cos[V_4(t_1) + u_4(t_1)] & \cdots & \cdots \\ f_4(t_0) \sin[V_4(t_0) + u_4(t_0)] & f_4(t_1) \sin[V_4(t_1) + u_4(t_1)] & \cdots & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$

To solve for the model vector, we perform an IRLS inversion, which evaluates a series of weighted least-squares problems that converge to an L1-norm solution (Aster *et al. 2013*). The L1-norm minimizes the absolute value of the residuals and is therefore highly effective at down-weighting outliers. For the initial model vector, we compute an L2-norm solution, from which a weighting matrix may be constructed based on the residuals between the observations and the forward model. Since the weighting matrix is a nonlinear function of the model vector, the normal equations must be solved iteratively. Thus, updated model vectors are derived from subsequent L1-norm solutions to the normal equations and tested against a tolerance value. The process repeats until a suitable level of convergence is achieved.

The resulting in-phase and quadrature coefficients for each isolable tidal harmonic may be re-combined to obtain amplitude and phase values

$$A_n = \sqrt{c_n^2 + s_n^2} \qquad (B8)$$

$$\phi_n = \operatorname{atan2}(s_n, c_n). \qquad (B9)$$