Spatial Variations of Stochastic Noise Properties in GPS Time Series

: The noise in position time series of 568 GPS (Global Position System) stations across North America with an observation span of ten years has been investigated using solutions from two processing centers, namely, the Paciﬁc Northwest Geodetic Array (PANGA) and New Mexico Tech (NMT). It is well known that in the frequency domain, the noise exhibits a power-law behavior with a spectral index of around − 1. By ﬁtting various noise models to the observations and selecting the most likely one, we demonstrate that the spectral index in some regions ﬂattens to zero at long periods while in other regions it is closer to − 2. This has a signiﬁcant impact on the estimated linear rate since ﬂattening of the power spectral density roughly halves the uncertainty of the estimated tectonic rate while random walk doubles it. Our noise model selection is based on the highest log-likelihood value, and the Akaike and Bayesian Information Criteria to reduce the probability of over selecting noise models with many parameters. Finally, the noise in position time series also depends on the stability of the monument on which the GPS antenna is installed. We corroborate previous results that deep-drilled brace monuments produce smaller uncertainties than concrete piers. However, if at each site the optimal noise model is used, the differences become smaller due to the fact that many concrete piers are located in tectonic/seismic quiet areas. Thus, for the predicted performance of a new GPS network, not only the type of monument but also the noise properties of the region need to be taken into


Introduction
Space geodetic techniques, and GNSS (Global Navigation Satellite System) in particular, are being used in the study of dynamics of the Earth's crust [1]. GPS stations have been installed worldwide to measure their changes in position over time with millimeter-level accuracy, which is associated with a number of geophysical phenomena such as plate motion [2][3][4], earthquakes (i.e., pre-, co-, and postseismic offsets and transients) [5,6], glacial isostatic adjustment [7][8][9][10], ocean tide loading [11][12][13][14][15], and atmospheric loading [16,17]. These geophysical phenomena can be described by a trajectory model (also called functional model) [18][19][20] and the remaining difference between the model and observations is noise, which is assumed to be a sum of stochastic processes. It is well known that noise is temporally correlated, for which power spectral density follows a power-law behavior with a spectral index of approximately −1 [21,22]. In addition, Langbein (2012) [23] and Dmitrieva et al. (2015) [24] found that random walk (RW) is present in some of the GPS time series. The effect of RW is almost a doubling of the uncertainty of the estimated linear motion [23]. On the other hand, He et al. (2019) described the analysis of 110 stations from the global IGS core network with more than 12 years of observations for which the power spectral density seems to flatten at low frequencies for 3-5% of horizontal components and around 12% for the vertical component [25]. The result is a reduction of the trend uncertainty by a factor of approximately two. In this research, we investigate the flattening of the power spectral density in some time series of the 568 permanent GPS stations around the USA (United States of America) and Canada, e.g., those blanketing the Pacific Northwest and the central western coast of the USA, from Vancouver Island to southern of California, with a concentration of stations around the Los Angeles area monitoring the San Andreas fault. A small amount of flattening of the power spectral density at low frequencies is also caused by the estimation of the linear trend [26,27]. Therefore, only strong flattening that starts at a period of 1-2 years will be considered. Flattening that starts at a longer period will be assumed to be artificial or undetectable.

GPS Data Processing
We used data from 568 continuously operating GPS receivers distributed over the United States, including Alaska. These data were used to compute time series consisting of position estimates with 1-day sampling by the analysis centers-PANGA/CWU (Pacific Northwest Geodetic Array/Central Washington University) and New Mexico Tech (NMT). Their solutions are given in the global International Terrestrial Reference Frame ITRF2008 reference frame [28]. The analysis center CWU computes the daily positions using the Precise Point Positioning method using the GIPSY software developed by NASA's Jet Propulsion Laboratory (JPL), which also provides the necessary satellite ephemerides [29], clock corrections, and wide-lane phase bias estimates [30]. Note that the station positions were loosely constrained during the initial estimation and subsequently transformed into the International Terrestrial Reference Frame (ITRF2008) using only the translation and rotation [28]-but not scale-components of the JPL-provided Helmert transformations. On the other hand, NMT processing was performed using the software GAMIT/GLOBK [31,32] utilizing the same stations in North America as those with the PANGA processing, but additional stations in other parts of the world were also included for the stability of the reference frame. The Vienna Mapping Function 1 (VMF1) grid was used in both processing by PANGA and NMT for handling the troposphere delay [33]. All common parameters used in both processing steps are explained in Herring et al. (2016) [34]. No common mode error filtering is used in any processing. It is important to emphasize that GAMIT double differencing does not remove common mode error. The network of GPS stations selected by NMT and PANGA is large (a quarter of the Earth's surface), which dilutes the strong common mode error that is detectable over smaller regions. The final processing of these time series described by Herring et al. (2016) [34] rotates the loosely constrained solutions provided by PANGA and NMT in the NAM08 reference frame using GLOBK [31,32].
Our study focuses only on the PANGA and the NMT solution (the original time series are cwu.final_nam08.pos and nmt.final_nam08.pos). For both the PANGA and the NMT solutions, the baseline lengths are their uncertainties do not change between the "loose" solutions submitted by PANGA and NMT and the solutions rotated/translated in the NAM08 reference frame. The difference with the processing preformed at NMT is mainly due to how the scale parameter is handled. The strategy in the latter includes the scale in the Helmert transformation, whereas the scale is not estimated directly in the former. Montillet et al.'s (2018) [35] study emphasized that the choice of including a radial scaling degree of freedom during daily reference frame realization primarily impacts the average network radial height and produces apparent height anomalies in excess of 5 mm that persist for months. A comprehensive discussion about the Helmert transformation and the scale parameter can be found in references [34,35].
In our analysis, the 568 stations have time series that began on 1 of January 2008 and ended on 1 January 2018. Our reason for choosing a fixed data time span is to reduce the differences between random models at different time scales. We also choose GPS stations with very few data gaps, less than 8%, which reduced the total number to 568 sites. In Appendix A, Table A1 shows that the percentage of the 568 permanent GPS stations are listed with less than a 3% data gap for each time series. As a result, more than 90% of the stations have more than 9.7 years of data. The average, maximum, and minimum data gaps of the 568 stations are also listed in this table to supply information (see Table A1 in the Appendix A) on the quality of the selected time series used throughout this study.
The GPS stations analyzed in this study have large diversity of monuments on which the GPS antenna has been installed. The metadata file (or log file) associated with each station provides a description of the monument, often referred to as mast, pillar, roof top, tower, or tripod [36][37][38]. In this study, we classified all monument types into four categories: concrete piers (CP), deep-drilled brace monument (DDBm), shallow-drilled brace monument (SDBm), and roof top/chimney (RTC). This classification follows previous studies [39,40]. Concrete pier (CP) is a pillar that can reach several meters attached deeply into the ground (up to 10 m below the surface). DDBm is braced monument where four or five 2.5 cm-diameter pipes are installed and cemented into inclined boreholes with the antenna attached at~1 m above the surface [40,41]. The pipes are also attached deeply below the surface (up to~10 m) using heavy motorized equipment. SDB refers to the type of equipment attached to the surface (<1 m-deep) using a hand-driller. The fourth category (RTC) gathers the antennas installed on the top of buildings sometimes using a mast attached to a wall, or with a concrete support. Note that our classification is based on the monument's description included in each log file available for each station (Table A1 in Appendix C gives more details on the monument type of the analyzed 568 sites).

GPS Time Series Analysis
First, outliers were removed from the time series. Outliers are observations that are larger than 3 times the interquartile range of the residual time series [42]. Second, the parameters of the trajectory model were estimated with weighted least-squares while the parameters of the model that describe the noise η i were estimated using maximum likelihood estimation [18,19]. For this estimation process, we used Hector software [43].
The trajectory model y(t i ) is a linear sum of the tectonic rate, seasonal signals, coseismic offsets, and random stochastic processes (e i ), see Bevis et al. [20]: where a is the initial position at the reference epoch t 0 , b is the rate, c j and d j are the periodic motion parameters (j = 1, 2 for annual and semiannual seasonal terms, respectively). The offset term g k can be caused by earthquakes, equipment (environment) changes, or human intervention, in which it is the magnitude of the change at epochs; is the total number of offsets; H is the Heaviside step function. The time of known offsets t k are retrieved from the station's metadata. Finally, the automatic offset detection algorithm developed by Fernandes and Bos (2016) is applied to detect undocumented offsets [44].
The Cascadia subduction zone is a convergent plate boundary that stretches from northern Vancouver Island in Canada to northern California in the United States. As we model and study time series of stations located in the Pacific Northwest including the Cascadia Mountains, specific events must be modeled, such as the Episodic Tremor and Slip (ETS) [45][46][47][48], for which we used the hyperbolic tangent function [49]. The amplitude of this function is described by S i . Finally, e i describes the noise/random stochastic processes. The time of the slow slip event and the delay of the postseismic deformation are required as input parameters for the estimation of the ETS using Hector [49] with a hyperbolic tangent of which the shape is prescribed by the time the ETS event occurred and its width (see in Appendix B). The times of the slow slip events can be requested from the Pacific Northwest Geodetic Array website or by a careful analysis of the time series with some training. The start of a slow slip event is evaluated via the correlation of seismic data together with a careful check of each time series [50]. In the remainder of this work, we use four delays-namely 30, 80, 100, and 130 days-for the postseismic deformation, because it is difficult to precisely estimate the duration of crustal decay. Note that these delays are conservative numbers knowing that the repetition of the ETS events is~14 months, as evaluated by previous geophysical studies of Cascadia [35,50]. These values represent a tradeoff in not modeling enough the phenomenon, and in contrast, absorbing other geophysical phenomena due to an overestimation of the decay time [50][51][52].

Stochastic Model Selection Criteria and Simulation Experiment
Power-law noise with a spectral index of −1 is called flicker noise (FN). Random walk (RW) noise has a spectral index of −2. Generalized Gauss Markov (GGM) noise is similar to power-law noise (PL) but flattens below a specified frequency. As noted in the introduction and shown in the following sections, the selection of the correct noise model has a significant influence on the trend uncertainty [21,23,53,54]. The theory of selecting the best model to describe the observation has a long history and many research areas simply use the Akaike or Bayesian Information Criterion [55,56]. However, Langbein (2004) followed a more empirical approach by performing Monte Carlo simulations using synthetic noise, and in this way, determined how much the difference between two log-likelihood values of two competing noise models must be before one can confidently choose one over the other [42].
In the report by Langbein (2004) and Santamaría-Gómez et al. (2011), the default noise model for these simulations, or the null model, was in all cases a random walk plus white noise and it was determined how much the log-likelihood value needed to be higher before one could accept another noise model as being better with 99% confidence [42,57]. The log-likelihood value will be abbreviated as MLE since it is estimated with the maximum likelihood method. Differences between two MLE values are represented as dMLE. The likelihood function represents a probability, although not normalized. Therefore, MLE and dMLE are logarithms of the probabilities and have no units.
We repeated the simulations using 5000 daily time series with a length of 10 years of synthetic random walk + white noise, each with amplitudes of 1 mm/yr 0.5 and 0.5 mm, respectively-that is, for 5000 simulations, the dMLE for which there are 50 values greater is identified as the 99% level to reject the null hypothesis. The results are shown in Table 1. The Bayesian Information Criterion (BIC) and BIC_tp are defined as follows (He et al., 2019) [25]: where MLE = ln(L), the log-likelihood value; ν is the number of parameters in the noise model; N is the number of observations. The noise model with the lowest BIC value is selected. For 10 years of daily observations, ln(N) = 8.
where v n and v b are the number of parameters and MLE n and MLE b are the MLE values of the null and new models, respectively. If this criterion is larger than zero, then the new noise model (b) is more likely than the null model (n  [25], we only consider the detection of GGM to be the most likely noise model significant if φ < 0.98; this parameter is also estimated by Hector. If this condition is not met, then the second most likely noise model is chosen. He et al. (2019) explained that this condition implies that we only detect GGM noise with flattening that already starts around a period of 1 year [25]. For the rest of this research, this extra condition of φ < 0.98 was always applied in addition to satisfying Equation (4).
The values of the parameters in the noise models used in the Monte Carlo simulations of Langbein (2004) [42] are slightly different from the noise values discussed here. To ensure the most realistic results, we determined the mean values for each estimated noise model for the horizontal and vertical components for our 5000 time series, see Tables 2 and 3. For each noise model, 5000 synthetic noise time series were generated. Each of them was analyzed using FN + WN (Flicker Noise + White noise), RW + FN + WN (Random Walk + Flicker Noise + White noise), GGM + WN (Generalized Gauss Markov + White Noise), and PL + WN (Power-law + White noise). Instead of tabulating the 99% quantile of the difference in MLE b − MLE n , we show all differences as box-whisker plots in Figures 2 and 3. Thus, if all values in the box-whisker plots were negative, one could be 100% sure that using the selected noise model (the noise model we think is correct) with the highest MLE value would indeed reflect the true underlying noise models (better than the alternative noise models). One can see that this is not always the case. Applying the BIC correction (BIC and BIC_tp), resulting in the blue box-whisker plots, reduces the MLE of noise models with more parameters than the test/null model and increases it for models with fewer parameters. Therefore, BIC helps to detect FN noise while reducing the rate of detecting GGM noise. The red box-whisker plots represent the results of using BIC_tp. BIC helps increase the number of true positives (the true noise model was selected) and reduce the number of false positives (the false noise model was selected). Overall, its performance is better than BIC; thus, it will be used for the rest of this research.  From Figures 2 and 3, it can be concluded from the RW + FN panels with positive box-whisker plots (MLE b − MLE n > 0) that in many cases, PL + WN noise is detected while in fact, the true underlying noise is RW + FN + WN [23]. On the other hand, if RW + FN + WN noise is detected, then we have high confidence that it is correct since false positives are extremely rare-see also He et al. [25]. Important for this research is the fact that for synthetic GGM, the MLE of the other noise models is always lower than that of GGM. In addition, from the other panels of Figures 2 and 3 one can also see that GGM is almost never selected when the underlying noise is not GGM. Thus, we can conclude with great confidence that any detection of GGM using BIC or BIC_tp is correct.
The random walk + flicker + white noise model was used by Langbein and Svarc [40] in the analyses of their time series. However, the analyses of our time series show that some flatten at low frequencies, which can be better described by a GGM noise model. These different conclusions might be caused by the fact that Langbein and Svarc [40] analyzed regionally filtered time series, whereas we are looking at unfiltered time series that are noisier and in which the smaller random walk signal might be hidden.
Note that Santamaría-Gómez et al. (2011) [57] concluded that using neither AIC nor BIC is recommended as a means to discriminate between models. Using similar Monte Carlo simulations, the reader should be convinced that for the time series of a length of 10 years used in this research, the performance of BIC and BIC_tp is actually similar to that of using the approach of Langbein and Svarc (2019) [40], leaving the general debate of using or not using the information criteria in the selection of the stochastic noise model for the future.

Multitaper Analysis
We have previously discussed the selection of noise models based on log-likelihood values. However, we can also illustrate the difference in stochastic properties of the noise by power spectral density (PSD) plots. Here, we used the Welch method to compute the PSD. We further used the multitaper method of Thomson [58], using the software of Prieto [59], which produces increasingly accurate estimates of the power at low frequencies [60]. Figure 4 shows three examples that illustrate the different behavior at low frequencies.
For some stations, a flattening of the power spectrum can be observed, while for other stations, the slope of the power spectra increases at low frequencies and is better described by a random walk noise model. For the left and center panels in Figure 4, the RW part of the RW + FN + WN model was very small and resulted in very similar power spectra as FN + WN.

Episodic Tremors and Slip
In this research, the focus is on the noise within GPS time series. However, the trajectory model should be accurate in order to ensure that the separation between geophysical processes and other noise sources is realistic. This is especially important for stations that experience episodic tremors and slip (ETS) events. ETS events are restricted to 20 stations located specifically in the Cascadia area. Table 4 displays the results of modeling or not the ETS events using 100 days of postseismic decay time scale for the two processing centers, PANGA and NMT. Tables A2 and A3 in Appendix A provide more details on the stochastic noise models' behavior of each East North Up (ENU) component. Note that the results using 30, 80, and 130 days of postseismic relaxation are included in the appendices (see Table A4), but do not vary significantly compared with the results shown in Table A4. It is worth emphasizing that there is a slightly higher frequency rate of selecting the RW + FN + WN model, for both PANGA and NMT processing, when including a small postseismic relaxation delay (i.e., 30 or 80 days) in the functional model, which indicates that the remaining relaxation phenomenon is modeled as RW noise. From Table 4, we can see that there is a much higher chance of selecting an RW component (i.e., the RW + FN + WN model) in the stochastic model for both processing strategies without modeling the ETS. Looking at the results using the BIC_tp, the percentages are 10.00% and 5.00% for PANGA and NMT, respectively, whereas with the ETS model, the percentages decrease to 1.67% and 0.00% for PANGA and NMT, respectively. The results are similar using the AIC. In addition, we can see that the AIC selects slightly more often the RW component in the case of PANGA processing when there is no modeling of the ETS events. This result echoes the conclusions in He et al. (2019) [25]. where it was found that the AIC criteria were slightly more sensitive to the detection of the RW component. Furthermore, with ETS modeling, the proportion of GGM + WN increases, from approximately 37% to 58% and approximately 41% to 51% for the PANGA and NMT solutions, respectively. In addition, modeling the ETS has a significant effect on the selection of the stochastic noise model of the GPS time series, as shown in Table A5 (Appendix A). As an example, the model changes when including or not including the ETS events, which is approximately 35% and 40.0% for PANGA and NMT, respectively, using the AIC, and approximately 21.7% for both PANGA and NMT with the BIC_tp criteria. In particular, for the variation in the stochastic noise models, the GGM + WN model more likely fits the time series of the Up component than the other models for the 20 sites considered here. Tables A2 and A3 show that the changes in stochastic noise models mainly occurred in the East and North components.

Spatial Variations Analysis
In this section, we use the 568 stations located in North America and displayed in Figure 5. To process the time series and select the optimal stochastic noise model using the AIC and BIC_tp criteria, we follow the same approach described in the previous sections.
Note that for the stations experiencing ETS events in the Cascadia region, the postseismic relaxation was set up to 100 days to be consistent with our results established in the previous section. In addition, we implement automatic offset detection on the GPS residual time series discussed in Fernandes and Bos (2016) [44]. Figure 5. Spatial distribution of the selected noise models using AIC (there is slight difference with BIC_tp, the main difference is that AIC is slight sensitive to the RW noise, as shown in Table 5). The selected stochastic noise models using AIC and BIC_tp over 568 stations show high consistency (approximately 98.7% and 98.6% for the PANGA and NMT solutions, respectively). The significant differences in the selection of the stochastic noise model between AIC and BIC_tp applied to the PANGA and NMT solutions are displayed in Table 5. We can see that AIC is generally more sensitive to RW noise in the selection of the RW + FN + WN and GGM + WN models, which supports the previous conclusions mentioned in He et al. (2019) [25]. Table 6 shows the percentage of stochastic noise models selected over the 568 stations for the PANGA and NMT solutions with only the AIC (for ENU and the average of the 3 components). The difference in terms of percentage is much higher in the Up component, where, for the NMT solutions, the GGM + WN model is more likely selected, and the PL + WN for the PANGA time series. The results in Table 6 also emphasize that the RW + FN + WN is approximately 8-14% for the horizontal components and only approximately 0.2% for the Up component; this is also demonstrated in Figure 5. Spatial distribution of the selected noise models using AIC (there is slight difference with BIC_tp-the main difference is that AIC is a little more sensitive to the RW noise, as shown in Table 5). The RW + FN + WN noise is much less in the Up component. In fact, we can see that all four models are selected in the horizontal components, whereas the vertical components have mostly PL + WN and GGM + WN noise models. This result holds in the analysis of the time series of both processing centers. This difference can be attributed to the fact that the vertical component is generally much noisier (~3 times greater) than the East and North components [22]. He et al. (2019) [25] concluded that GGM + WN fits slightly better than the other noise models for the vertical component, which is also supported by our results.
To further analyze the potential correlation of sites showing RW noise with geophysical phenomena, we highlight the sites for both NMT and PANGA solutions in Figure 5 (i.e., red dot for the RW + FN + WN). Most are located near the coasts (distance <10 km from the shoreline-see Figure A1 in Appendix A, which shows the spatial distribution of the RW sites for both PANGA/NMT) and in known tectonically active areas: the Cascadia subduction zone and San Andreas strike-slip fault. This result is in agreement with our results shown in the previous section, where we demonstrated that when the functional model does not completely capture all geophysical signals (due to earthquakes and/or postseismic deformation and possibly generates short-duration transient signals), the algorithm estimates higher amplitude RW noise than in much less tectonically active areas. The results also show that the detected RW amplitudes are larger than 0.5 mm/yr 0.5 . (see Figure A1 and Table A6 in Appendix A).

Relation between Type of Noise and Type of Monument
First, we analyze the velocity uncertainty of the 568 GPS stations related to type of monument. Figure 6 shows the spatial distribution of the velocity uncertainty with the optimal stochastic noise model, using AIC, for PANGA and NMT. From Figure 6, we can see that the Nevada/Idaho and Utah regions are quieter (i.e., lower uncertainty value) than the other areas, especially for the NMT solutions. In contrast, the Washington/Oregon region is noisier (i.e., higher uncertainty value). It can also be observed that there is no obvious pattern between types of monuments and velocity uncertainty. DDB monuments have the largest percentage of random walk component according to Figure 6. These monuments are mostly located along the San Andreas fault in central and southern California. Therefore, they display larger noise amplitudes due to high tectonic activities [61][62][63]. CP monuments are mostly located near Washington/Oregon states and near the Great Lakes. To recall the previous discussion, the Pacific Northwest is also a tectonic active area with the Cascadia subduction zone [64]. The Great Lakes region is not a strong tectonic active area but is known for strong postglacial rebound [65] and complex geodynamics [66]. Besides, there are other geodynamic mechanisms related to the lakes such as the lake-level variability impacted by climate change [67]. For this type of monument, GGM + WN is often the optimal model. Furthermore, the statistical analysis of Tables A7 and A8 together with Figure 7 show that there is diversity in the results depending on the types of monuments and components. Previous results established that the FN + WN model is the stochastic noise model mostly selected for East and North components. Our results support that the FN + WN and PL + WN are still the most common noise models, but in the up component, GGM + WN is more present as well for both PANGA and NMT solutions.
Note that the largest ratio of GGM + WN model mainly appears for the CP monuments (44.4% Up (PANGA) and 60% Up (NMT)). Thus, we cannot conclude on the discrimination between type of monuments and stochastic noise model. The noise properties are spatially varying so we cannot define a general 'preferred' noise model for a type of monument since it depends on the geographical distribution of the monument. The largest ratio of RW mainly appears in DDB and SDB monuments. Moreover, the largest percentage of RW component is prominent in the SDB and DDB types of monuments. Furthermore, the statistical analysis of Tables A7 and A8 together with Figure 7 show that there is diversity in the results depending on the types of monuments and components. Previous results established that the FN + WN model is the stochastic noise model mostly selected for East and North components. Our results support that the FN + WN and PL + WN are still the most common noise models, but in the up component, GGM + WN is more present as well for both PANGA and NMT solutions.
Note that the largest ratio of GGM + WN model mainly appears for the CP monuments (44.4% Up (PANGA) and 60% Up (NMT)). Thus, we cannot conclude on the discrimination between type of monuments and stochastic noise model. The noise properties are spatially varying so we cannot define a general 'preferred' noise model for a type of monument since it depends on the geographical distribution of the monument. The largest ratio of RW mainly appears in DDB and SDB monuments. Moreover, the largest percentage of RW component is prominent in the SDB and DDB types of monuments. Furthermore, the statistical analysis of Tables A7 and A8 together with Figure 7 show that there is diversity in the results depending on the types of monuments and components. Previous results established that the FN + WN model is the stochastic noise model mostly selected for East and North components. Our results support that the FN + WN and PL + WN are still the most common noise models, but in the up component, GGM + WN is more present as well for both PANGA and NMT solutions. Note that the largest ratio of GGM + WN model mainly appears for the CP monuments (44.4% Up (PANGA) and 60% Up (NMT)). Thus, we cannot conclude on the discrimination between type of monuments and stochastic noise model. The noise properties are spatially varying so we cannot define a general 'preferred' noise model for a type of monument since it depends on the geographical distribution of the monument. The largest ratio of RW mainly appears in DDB and SDB monuments. Moreover, the largest percentage of RW component is prominent in the SDB and DDB types of monuments.
From Tables A7 and A8, we can see some slight variations when comparing the results between the East and North components for each type of monument. The percentage of the selected RW + FN + WN model is the largest for the DDB monuments on the East component (~12.15% PANGA; 11.6% NMT), but not on the North component (~11.60% PANGA; 12.15% NMT). Instead, the probability of selecting this noise model is the largest for the SDB monuments on the North component (~19.10% PANGA; 16.73% NMT).
This sensitivity to RW noise is not experienced for the CP and RTC monuments. However, one needs to take into account that in our study, most of the CP monuments are installed in tectonically active areas where higher noise amplitude can mask the RW component. The RW component models the small-amplitude, short-duration, transient signals originating from regional (or local) deformations, which seem to be more present in the time series of the SDB and DDB monuments.
In addition, the stochastic noise properties of the vertical component are mostly modeled with PL + WN. However, the GGM + WN model challenges this result in the case of CPs where the percentage is higher, approximately 44.44% and 60.00% for PANGA and NMT solutions, respectively. Thus, these time series of the vertical component experience flattening at a low frequency, as discovered and analyzed by He et al. (2019) [25].
To further support this result, we draw a boxplot of the velocity uncertainties with the selected FL + WN and optimal noise model, as shown in Figures 8 and 9, respectively. In Figure 8, a clear trend is seen from stable monuments such as the DDBm. However, Figure 9 does not show the differences between the types of monuments. We think that a large number of receivers on CP experience GGM noise, which halves the trend uncer-tainty. The total effect is a reduction in the correlation between the type of monument and velocity uncertainties.  From Figures 8 and 9, we can also see that the RTC monument performs as well as the DDBm. To understand this result, we should look again at Tables A7 and A8 in the Appendix A; they display the amplitude of the white noise for the PANGA processing using the AIC criteria. We notice that the amplitude of the white noise is largest for the RTC monuments. Therefore, the RTC setup may partially mask small, transient, geophysical signals by increasing the amplitude of the white noise, rendering this setup disadvantageous for geophysical studies at the regional scale when looking at small amplitude short-duration signals. From Tables A7 and A8, we can also see that the estimated white noise amplitude is very different between the NMT and PANGA solutions. As discussed in Section 1, PANGA uses GIPSY-OASIS (GIPSYX) without any scale parameter in the Helmert transform, whereas the NMT solution results from the processing with GAMIT/GLOBK including a scale parameter. It results in different stochastic noise properties of the time series, with a smaller scatter than the PANGA ones (~30%), and hence, a smaller noise amplitude.
In light of this discussion, we can conclude that the higher white noise amplitude associated with the RTC monuments can mask the presence of a RW component (see Tables A9 and A10 in Appendix A, which show the relation between the monument type and estimated noise model) in the stochastic noise model than in other monuments (SDBm and DDBm). This can explain why the velocity uncertainties in Figures 8 and 9 are smaller.
From the above results and analysis, we can conclude that there is reduced discrimination as to which monument performs the best over the three coordinates. The spatial distribution of the monuments supersedes the type of monuments in the selection of the noise model for the global, not regionally filtered, GNSS time series. Previous studies, such as Herring et al. (2016) [34], have warned about the spatial distribution across North America. Williams et al. (2004) restricted their study to a small area to determine the influence of the various types of monuments [22]. Beavan (2005) [68] concluded that monument noise is not the dominant factor in the stochastic noise properties of the GPS time series; hence, we corroborate his results. Finally, we supplement CME reduction and the elimination of some stations in active areas in the Appendix A, to explore the effect of regional filtering and large postseismic relaxation events on monument performance. We can see that regional filtering does not change the overall results. The DDBm is still the best monument type, which supports what we have already previously discussed in this paper.

Conclusions
This work has investigated the noise properties of 568 GPS daily position time series with a length of 10 years in North America, computed by NMT and PANGA, together with their spatial distribution. These time series are given with respect to a global reference frame (ITRF2008) and no regional filter has been applied.
To model the noise in the GPS time series, the following models were used: FN + WN, RW + FN + WN, GGM + WN, and PL + WN. The selection of the optimal noise model was based on the log-likelihood value and the information criteria (i.e., AIC, BIC, BIC_tp) to avoid overselecting noise models with many parameters. This approach was criticized by authors of various studies who do not recommend the use of AIC or BIC. However, when we apply these information criteria to Monte Carlo simulations of synthetic time series, we can demonstrate that for noise RW + FN + WN, we have a very low percentage of false positives. Thus, if we detect RW + FN + WN noise, then we have high confidence this is correct. On the other hand, we still have a high number of false negatives, which means that PL + WN is still detected, while in reality, the underlying noise model is RW + FN + WN, in agreement with previous studies [23,25]. The simulations also demonstrate that we are able to reliably detect the GGM noise model. Statistically, we found that for around 10% of the stations, the power spectral density shows flattening at low frequencies, which can be better described with a Generalized Gauss Markov (GGM) noise model, further extending the conclusions of He et al. (2019).
For the Cascadia subduction region, which is tectonically active, it is necessary to include episodic tremor and slip in the trajectory model. Otherwise, the misfit between observations and model is absorbed as a large random walk component, which affects the value of the estimated linear motion and increases its uncertainty. In addition, our results show that the selection of the RW component in the RW + FN + WN model accounts for a non-negligible percentage when not modeling the ETS events, whereas when properly modeled, the GGM + WN and PL + WN models are more often selected. Note that AIC is generally more sensitive to the RW noise in the selection of the RW + FN + WN and GGM + WN model. Furthermore, the spatial distribution of selected noise models shows a pattern in our case study. Stations including the RW component are generally located near the coasts (distance <10 km from the shoreline) or in active tectonic areas (Cascadia subduction zone, San Andreas fault). In terms of differences between the two processing strategies, the time series from the NMT solutions are more often modeled with the GGM + WN than the PANGA solution, therefore exhibiting a flattening of the power-spectrum. In addition, the largest percentage of RW + FN + WN model is selected in SDB and DDB monuments mainly for the horizontal components. The PL + WN and FN + WN models are generally the dominant stochastic noise models for the monuments, but the GGM + WN is also often selected for the CP monuments, especially on the Up component, which means that the very long time series generated from monitoring the CP monuments are experiencing this flattening at low frequency. The difference in results between AIC and BIC_tp is relatively small for each monument in selecting the various stochastic noise models. However, we show that the amplitude of the white noise is largest for the RTC. Therefore, the RTC setup may partially mask small transient geophysical signals by mismodeling the amplitude of the noise. Moreover, when investigating the relationship between the various types of monuments and the estimation of the geophysical signals, the results emphasized that the smallest uncertainties associated with the tectonic rate are recorded by the DDBm monuments for the East and North components, especially selecting the FL + WN model. On the vertical component, the results difference between the DDBm and SDBm is not significant. Additional results obtained by removing CME and monuments located in tectonically active areas (such as the Cascadia subduction zone) show that the uncertainties decrease, especially with CP and DDBm monuments on the East and North components. CP and DDBm monuments are mostly located in tectonically active areas or close to the coastline. Overall, we recommend the use of DDBm in noisy or tectonically active areas. However, their stochastic noise properties are different to those in more stable areas. Therefore, we conclude that the location of the station supersedes these types of monuments in the estimation of the tectonic rate uncertainty, supporting previous studies [34,68].

Conflicts of Interest:
The authors declare no conflict of interest.  As mentioned in the previous section on 'GPS processing', the PANGA and NMT solutions did not implement any regional filtering. We extract the regional common mode error (CME) on the GPS network for both PANGA and NMT solution. For filtering on the CME, we implement the VBPCA method to extract CME (Li et al., 2020). The VBPCA method imposed a priori distributions on the model parameters and estimated the hyperparameters by maximizing the evidence of observed signals. This method can naturally handle missing data in the Bayesian framework and utilizes the variational expectation-maximization iterative algorithm to search the CME for the incomplete GNSS position time series. Moreover, it can automatically select the optimal number of principal components for data reconstruction and be more reliable against the overfitting problem than other statistical signal decomposition methods. The boxplots with CME filter of the velocity uncertainties with the selected FL + WN and optimal noise model are shown in Figure A2a. From Figure A2a together with Figures 8 and 9, we can see that regional filtering does not change the overall results. The DDBm is still the best monument type. The CME removal is only interesting looking at the CP monuments. However, we only have 45 CP in this study. In addition, some of the 568 analyzed sites were located in a region that was known to contain tectonic or volcanic transient events, which will affect the stochastic noise properties and, therefore, the monument performance. Thus, we removed these sites (reducing the total number to 538 sites). The boxplot figures after filtering the CME and with optimal model selection are shown in Figure A2b. Overall, the filter improves only the results with CP monuments. These monuments are located in the Pacific Northwest, with postseismic relaxation, and this result supports what we have already said in this study.   Figure A3. Boxplot with CME filter with optimal model for PANGA and NMT solution. (a) Boxplot with CME filter with optimal model for PANGA and NMT. (b) Boxplot with CME filter with optimal model for PANGA and NMT (remove large postseismic relaxation sites).

Appendix B
For the postseismic relaxation or Episodic Tremor and Slip events, we assume that the observations are the sum of a deterministic model and stochastic noise. This deterministic model is as follows: ( ) H t is the Heaviside step function and used to model offsets with amplitudes j b . An annual and semiannual signal are commonly included in the model and therefore have their own keywords, 'seasonal' and 'halfseasonal', in Hector [43]. In the equation, their angular velocities are represented by i ω . Other periodic signals need to be entered using the keyword 'PeriodicSignals', followed by a list of their periods in days. The time of the slow slip event and the delay of the postseismic deformation are Figure A3. Boxplot with CME filter with optimal model for PANGA and NMT solution. (a) Boxplot with CME filter with optimal model for PANGA and NMT. (b) Boxplot with CME filter with optimal model for PANGA and NMT (remove large postseismic relaxation sites).

Appendix B
For the postseismic relaxation or Episodic Tremor and Slip events, we assume that the observations are the sum of a deterministic model and stochastic noise. This deterministic model is as follows: where p i are the coefficients of a n p degree polynomial. By default, n p = 1 means a linear trend. H(t) is the Heaviside step function and used to model offsets with amplitudes b j . An annual and semiannual signal are commonly included in the model and therefore have their own keywords, 'seasonal' and 'halfseasonal', in Hector [43]. In the equation, their angular velocities are represented by ω i . Other periodic signals need to be entered using the keyword 'PeriodicSignals', followed by a list of their periods in days. The time of the slow slip event and the delay of the postseismic deformation are required as input parameters for the estimation of the ETS using Hector with a hyperbolic tangent by providing two parameters, t k and T k .