Daily Post-Processing Solutions in GAMIT
Derived from the post-processing in GAMIT, the daily positions of the 46 stations (33 of the study area and 13 IGS stations defining the reference frame) were evaluated based on the quality criteria described by Herring et al. (2018). The phase ambiguities resolution is considered successful for this study with 95.4% for the wide-lane (WL) and 92.1% for the narrow-lane (NL), as a resolution in phase ambiguities of more than 90% for WL and 80% for NL is recommended. Another quality criteria is the value of the Normalized Root Mean Square (NRMS) of the position fit (post fit) defined as the square root of Chi-square (χ2) divided by the degrees of freedom. Herring et al. (2018) recommend a maximum value of 0.22. During the entire GPS processing, there are no values higher than 0.22. The NRMS ranges from 0.169 to 0.193 and the average is equal to 0.179. According to Herring et al. (2018) the best sites will have values of 3–5 mm Root Mean Square (RMS) of the phase residuals, the worst sites 7–9 mm RMS. Values from 10 to 15 mm indicate a high but acceptable noise level, and values greater than 15 mm suggest poor signal reception, a high multipath environment, severe climates, convergence problems caused by few data or a priori poor coordinate. Figure 2 shows the daily frequency histograms of the global RMS present in each year. Most of them show a left-sticky behavior, where on the abscissa axis are located the RMS values of the phase residuals and can be considered as sites with non-significant noise levels. However, in most of the years there were RMS values higher than 10 mm but not higher than 12 mm. This fact shows stations with an acceptable level of noise. Nevertheless, the number of days that this occurred was less than 5 days per year, which results in the behavior of the histograms tending to decrease in frequency towards the right, with a behavior similar to a lognormal distribution. In this sense, the processing in GAMIT ensures the accuracy and consistency of coordinate solutions with daily repeatability and non-significant noise levels.
Velocity and uncertainty analysis at 33 CORS stations in northwestern Mexico
Prior to any velocity estimation, and to avoid biases, the Globk tsfit module was implemented during GPS data processing, and then the time series were cleaned by visual inspection with TSVIEW (Herring 2003), wich allows implementing the tsfit models graphically, eliminating outliers when residuals are greater than 3σ and identifying offsets in the time series derived from seismic activity and antenna changes. It is also necessary to consider seasonal signals, i.e., annual and semi-annual effects, which according to Blewitt and Lavallée (2002) generate a bias in the velocity estimation. Thus a continuous operation time of at least 2.5 years is recommended. In this sense, a comparison has been made between estimated velocities considering and not considering the seasonal signal effects. Figure 3 shows on the left vertical axis the velocity difference between considering and not considering the seasonal signal effects expressed in mm/year. As well as the bar graph corresponding to this difference in the 3D velocity at each station. On the right vertical axis, the operating time of each station has been plotted (blue line with diamonds), ordering the stations according to the time of data acquisition from the shortest to the longest. It is evident that those stations with an operating time of less than the 2.5 years recommended by Blewitt and Lavallée (2002) are the stations with the greatest difference in velocity (between considering and not considering the seasonal signal effects), as is the case of the LPAZ and IAGX stations with variations between 0.40 and 0.80 mm/year. Stations with operating times of no more than 4.5 years (Blewitt and Lavallée 2002) such as IHER, NVLX and CORX show velocity variations between 0.14 and 0.36 mm/year.
As operating time increases, the velocity differences of each station between considering and not considering the seasonal signal effects decrease and become minimal. The differences of velocity by seismic event is shown by an experimental case with CORS station HER2, which intentionally did not undergo linear adjustment due to the co-seismic displacement detected by the station, evidencing a strong difference in velocity due to not considering the correction to this displacement of -1.72 mm North, 5.98 mm East and − 2.37 mm in the vertical component. So despite having an operating time of 8.4 years, a bad adjustment in the corrections applied in the time series will bias the velocity estimation (Fig. 3). When the linear trend adjustment is applied before and after co-seismic displacement, whose displacement behavior in the time series is generally logarithmic or exponential, the velocity bias is minimized, as shown in Fig. 3. A similar behavior is shown at the station of MEXI, which identified the Cucapah earthquake of magnitude 7.2 Mw with a co-seismic displacement of -183.88 mm in the North component, 182.09 mm in the East component and 16.87 mm in the vertical component. Therefore, after adjusting the displacement the station shows a minimum variation in velocity (0.02 mm/year) and an operating time of 11.2 years
Once the processing had been completed, as well as the corresponding cleaning of the time series, identification of hardware and software updates, in addition to the inclusion of white noise and temporally correlated noise using the real-sigma algorithm, the differences in estimated velocity uncertainties were evaluated, because in the state of the art it has been identified that considering only white noise in the estimation in velocity will give a underestimation of velocity uncertainties (Williams 2004; Shen et al. 2011). Therefore, a comparison was made in the uncertainties estimated with white noise and the combination of white noise plus temporally correlated noise (random walk noise). This comparison is shown in Table 2, where it can be seen that the uncertainties (σ) of the white noise (WN) effect are lower and are considered underestimated because this noise coming from the hardware and is independent of frequency and time, while considering temporal correlation through the inclusion of random walk noise (RW) caused by instability of the monument, as well as to the correlation of the positions from one epoch to another, favoring the realistic estimation of the uncertainty of the velocity.
Table 2
Realistic velocities and uncertainties comparing White Noise (WN) and White Noise plus Random Walk Noise (WN + RW). Vn, Ve and Vu stand for north, east and up velocity components, respectively. σVn, σVe and σVu stand for the standard deviation of Vn, Ve and Vu, respectively
Station
|
RW + WN
|
WN
|
Vn
|
σVn
|
Ve
|
σVe
|
Vu
|
σVu
|
Vn
|
σVn
|
Ve
|
σVe
|
Vu
|
σVu
|
LPAZ
|
24.40
|
0.39
|
-51.03
|
0.28
|
-2.77
|
1.23
|
23.96
|
0.17
|
-50.31
|
0.14
|
-2.87
|
0.33
|
IAGX
|
-5.31
|
0.49
|
-17.79
|
0.10
|
-0.90
|
0.56
|
-6.12
|
0.12
|
-17.42
|
0.07
|
-1.39
|
0.25
|
IHER
|
-6.99
|
0.17
|
-11.63
|
0.29
|
1.10
|
0.71
|
-6.85
|
0.09
|
-11.81
|
0.08
|
1.32
|
0.24
|
NVLX
|
33.66
|
0.12
|
-27.82
|
0.09
|
-76.77
|
0.39
|
33.66
|
0.03
|
-27.81
|
0.02
|
-76.53
|
0.08
|
CORX
|
20.10
|
0.34
|
-40.94
|
0.10
|
0.01
|
0.26
|
19.84
|
0.05
|
-40.91
|
0.03
|
-0.35
|
0.09
|
PB2Y
|
16.71
|
0.11
|
-38.93
|
0.09
|
0.49
|
0.14
|
16.70
|
0.03
|
-38.92
|
0.02
|
0.50
|
0.05
|
PLCX
|
20.44
|
0.13
|
-42.61
|
0.05
|
0.42
|
0.43
|
20.35
|
0.03
|
-42.62
|
0.02
|
0.44
|
0.07
|
GUAX
|
23.80
|
0.14
|
-47.42
|
0.06
|
0.13
|
0.14
|
23.77
|
0.04
|
-47.43
|
0.02
|
0.10
|
0.07
|
TSFX
|
17.80
|
0.16
|
-41.36
|
0.07
|
0.86
|
0.23
|
17.78
|
0.03
|
-41.35
|
0.02
|
0.83
|
0.06
|
TNBA
|
20.02
|
0.22
|
-42.79
|
0.20
|
0.33
|
0.14
|
20.02
|
0.03
|
-42.77
|
0.02
|
0.30
|
0.07
|
TNPP
|
-8.46
|
0.14
|
-13.11
|
0.04
|
0.95
|
0.16
|
-8.49
|
0.02
|
-13.08
|
0.01
|
0.93
|
0.04
|
NAYX
|
14.31
|
0.20
|
-25.32
|
0.34
|
-1.21
|
0.38
|
14.36
|
0.02
|
-25.32
|
0.03
|
-1.04
|
0.05
|
TNCU
|
-6.22
|
0.15
|
-11.51
|
0.05
|
0.50
|
0.36
|
-6.20
|
0.02
|
-11.50
|
0.02
|
0.51
|
0.04
|
USMX
|
-7.10
|
0.17
|
-11.58
|
0.05
|
0.90
|
0.36
|
-7.10
|
0.02
|
-11.58
|
0.01
|
0.91
|
0.05
|
YUMX
|
16.60
|
0.24
|
-36.23
|
0.15
|
4.12
|
0.19
|
16.61
|
0.02
|
-36.21
|
0.01
|
4.03
|
0.04
|
HER2
|
-7.15
|
0.11
|
-10.55
|
0.16
|
-0.76
|
0.37
|
-7.38
|
0.02
|
-9.75
|
0.02
|
-1.09
|
0.03
|
PHJX
|
14.88
|
0.09
|
-31.03
|
0.12
|
-0.20
|
0.13
|
14.89
|
0.01
|
-31.03
|
0.01
|
-0.15
|
0.03
|
YESX
|
-6.73
|
0.19
|
-11.09
|
0.03
|
0.08
|
0.44
|
-6.72
|
0.02
|
-11.08
|
0.01
|
0.11
|
0.03
|
QUEX
|
-7.99
|
0.16
|
-16.72
|
0.11
|
-4.55
|
0.17
|
-8.03
|
0.02
|
-16.70
|
0.01
|
-4.69
|
0.03
|
IMIE
|
19.69
|
0.10
|
-41.41
|
0.13
|
-2.15
|
0.30
|
19.68
|
0.01
|
-41.42
|
0.01
|
-2.16
|
0.03
|
PSTX
|
20.82
|
0.09
|
-42.62
|
0.16
|
-0.90
|
0.13
|
20.82
|
0.01
|
-42.62
|
0.01
|
-0.91
|
0.02
|
IPAZ
|
21.92
|
0.20
|
-49.37
|
0.23
|
0.37
|
0.35
|
21.93
|
0.02
|
-49.37
|
0.02
|
0.36
|
0.04
|
DAEX
|
18.08
|
0.15
|
-37.63
|
0.48
|
2.50
|
0.09
|
18.05
|
0.01
|
-37.63
|
0.02
|
2.47
|
0.02
|
PTEX
|
18.29
|
0.08
|
-40.93
|
0.19
|
-1.94
|
0.18
|
18.29
|
0.01
|
-40.93
|
0.01
|
-1.94
|
0.02
|
PJZX
|
17.65
|
0.06
|
-40.80
|
0.25
|
-0.18
|
0.16
|
17.64
|
0.01
|
-40.79
|
0.01
|
-0.20
|
0.02
|
PLTX
|
19.46
|
0.14
|
-39.47
|
0.20
|
1.51
|
0.11
|
19.46
|
0.01
|
-39.46
|
0.01
|
1.51
|
0.02
|
PLPX
|
19.90
|
0.09
|
-37.85
|
0.28
|
0.28
|
0.19
|
19.90
|
0.01
|
-37.85
|
0.01
|
0.30
|
0.02
|
IMIP
|
-4.93
|
0.07
|
-12.57
|
0.04
|
-3.21
|
0.18
|
-4.91
|
0.01
|
-12.59
|
0.01
|
-3.20
|
0.02
|
PALX
|
20.28
|
0.08
|
-42.00
|
0.15
|
-0.96
|
0.21
|
20.29
|
0.01
|
-42.00
|
0.01
|
-0.95
|
0.02
|
IDGO
|
-5.34
|
0.07
|
-9.81
|
0.06
|
-3.10
|
0.44
|
-5.32
|
0.01
|
-9.81
|
0.01
|
-3.13
|
0.04
|
MSD1
|
20.74
|
0.56
|
-49.15
|
0.25
|
1.12
|
0.93
|
20.73
|
0.03
|
-49.15
|
0.02
|
1.15
|
0.07
|
MEXI
|
8.81
|
0.55
|
-29.74
|
0.78
|
-1.86
|
1.00
|
8.80
|
0.03
|
-29.74
|
0.03
|
-1.89
|
0.04
|
Since the treatment and analysis of the time series has been manual, there is the possibility that the velocity may have some variation in relation to human manipulation. However, to corroborate that there are no significant variations in the calculation of the velocity and therefore a biased interpretation of the displacements, the MIDAS algorithm (Blewitt et al. 2016) has been considered as well. MIDAS is an automatic velocity estimator that is insensitive to problems such as offsets, discontinuities, outliers, seasonality and skewness among others, selecting pairs of data separated by one year estimated by the median of all possible slopes of the data set and based on the distribution of the sampled slopes whose algorithm is described in detail in Blewitt et al. (2016). In this sense, the residuals of the difference of the velocities obtained by TSVIEW and MIDAS are estimated. These residuals are shown in Fig. 4, where the frequency distributions with respect to the ITRF14 show a highest accumulated frequency within the values of 0 to 0.5 mm/year for the horizontal components; the most values of the vertical component are between 0 and 1 mm/year. The residuals fixed to the North American Plate are between 1 and 2 mm/year for the horizontal components and between 0 and 2 mm/year in the vertical. Finally, the highest accumulated frequency fixed to the Pacific Plate is between 0 and 1 mm/year for the three components (E,N,U). Based on this velocity comparison (TSVIEW - MIDAS) it is certain that the velocity variations could be considerated negligible, therefore in the following sections it is considered to work with the velocity solutions generated with TSVIEW (Herring 2003).
Velocity Fields and Deformation Interpretation in Northwest Mexico
Once the velocities and their corresponding uncertainties had been estimated considering white noise plus temporally correlated noise, as well as the rigorous processing models to achieve the best accuracy, the velocity field corresponding to the solution adjusted to the ITRF14 was estimated from the 2010 to 2021 GPS data set. The velocity field for Northwest Mexico is shown in Fig. 5a, and it is evident that the behavior of the stations located on the Pacific Plate show a displacement mostly in a northwest (NW) direction, while the stations located on the North American Plate show a displacement in a southwest (SW) direction. The analysis derives from a global tectonic circuit which refers to an absolute movement, therefore both plates move together at different displacement rates, i.e., the stations located on the North American Plate have a lower displacement rate compared to those on the Pacific Plate. On average, the stations on the Pacific Plate have a velocity of 47.20 mm/year while those on the North American Plate have an average of 12.88 mm/year, these average velocities being adjusted to the ITRF14 (Fig. 5a). Figure 5b shows the velocity field fixed to the North American Plate in order to demonstrate the displacement of the stations on the Pacific Plate, which show a displacement behavior similar to that of a block in the NW (Northwest) direction, since the velocities are similar to each other both in magnitude and direction. The NVLX station shows influence due to local deformations given its proximity to the transform fault zone, moving at a lower velocity in the East component in relation to the other stations. This velocity is consistent with the velocity field to northern Baja California reported by Gónzalez-Ortega et al. (2018). Similarly, the velocity field fixed to the Pacific Plate has been estimated as shown in Fig. 5c, where the displacements of the stations on the North American Plate describe a Southeast (SE) direction, including the IAGX station located on Angel Island in the Gulf of California. The average displacement obtained is 44.72 ± 0.29 mm/year in the NW direction for the Pacific Plate fixed to the North American Plate and 45.34 ± 0.18 mm/year in the SE direction for the North American Plate fixed to the Pacific Plate.
In comparison with other studies, the velocities estimated in this study are similar specifically for the North American Plate. DeMets et al. (2010) computed through geological methods a velocity of 48.6 mm/yr for the North American Plate and 51 mm/yr for the Pacific Plate. The NUVEL1-A model (DeMets et al. 1994) estimates a velocity of 47.4 mm/yr for the North American Plate and 51.6 mm/yr for the Pacific Plate. The MORVEL model estimates a velocity of 49.4 mm/yr and 50.2 mm/yr for the North American and Pacific Plates, respectively. A computation by Platner et al. (2007) reports a velocity of 46.8 mm/yr for the North American Plate and 51.1 mm/yr for the Pacific Plate. The box and whisker diagram on Fig. 5d shows the different between the studies described and this. It is evident that the velocities for both tectonic plates are in the third quartile of the diagram and that the estimated velocities for the North American Plate show a smaller difference compared to the mean of this study (horizontal red line Fig. 5d) of between 1 and 4 mm/year, while the velocities in the Pacific Plate have differences of between 5 and 6 mm/year. However, these differences are in the range of minimum and maximum velocities for each tectonic plate, in addition to the fact that the NUVEL1-A and MORVEL models are generated through geological techniques such as magnetic anomalies in global tectonic circuits.
Crustal deformation model for the Northwest (NW) of Mexico
Derived from the calculation and estimation of velocities at CORS stations, the interpolation algorithm implemented by Sandwell and Wessel (2016) in GMT6 was applied to the set of horizontal velocities (horizontal vector data). These were manipulated using Green's functions (Equations 5 and 6), adjusting the elastic coupling under Poisson's coefficient (ν) with a value equal to 0.5, i.e. typical elasticity based on Bian et al. (2020) and Wang et al. (2022). Figure 6a shows the interpolated 2D crustal deformation model, where the red vectors show the velocities corresponding to the CORS stations calculated and analyzed in this research, with which the interpolation was implemented using Green's functions on a 15’ x 15’ grid in order to observe the deformation characteristics in areas where there are no CORS stations. The blue vectors show the interpolated cortical velocity model based on the known input data (velocities from GPS processing). This velocity model favors the continuous deformation analysis, giving evidence of the global behavior of Northwestern Mexico. Subsequently, this model will be referred to as NW-MEXVEL. Furthermore, this displacement can be evidenced in some specific locations such as Guadalupe Island in the Pacific Plate and Angel and Tiburon Islands in the Gulf of California, showing correspondence and coherence to the closest displacement vectors (Fig. 6a). The NW-MEXVEL cortical velocity model is estimated considering the crust of Northwest Mexico as an elastic surface, modified under the geodynamic behavior and interaction of the Pacific and North American plates evaluated by the horizontal displacements in the GPS stations of the study area. Moreover, the plate boundary near the San Andres fault located in the Mexicali Valley, shows a change of direction in the interpolated vectors when crossing the fault (yellow line in Fig. 6a), which demonstrates the close behavior at the transform boundary.
Figure 6b shows the interpolated velocity model based on the displacement vectors calculated for the same coordinates of the 33 CORS stations obtained from the GEODVEL 2010 model (Argus et al. 2010), taking into account that this model has been estimated under geodetic space observations (GPS, VLBI - Very Long Baseline Interferometry, SLR - Satellite Laser Ranging and DORIS - Doppler Orbitography and Radiopositioning Integrated by Satellite) for the 11 major tectonic plates. In this sense, there are areas close to plate boundaries where the GEODVEL model may have a greater uncertainty to the local tectonic behavior, one example being the velocity estimated at Angel Island in the Gulf of California (Fig. 6b). Compared to the velocity estimated with discrete data from the IAGX station, the GEODVEL model's velocity is more similar to the velocities in the Pacific Plate. This is due to the island's proximity to the transform boundary (Fig. 6b). This model (Fig. 6b) differs slightly from the model based on the velocities of each GPS station (NW-MEXVEL-Fig. 6a). The former is smoother and shows a continuous change in direction of the interpolated displacement vectors in northern Baja California near the San Andres fault.
The NW-MEXVEL model fixed to the Pacific Plate (Fig. 6c) has been computed. The movement of the North American Plate is more evident with a continuous southeasterly displacement, better defining the plate boundary as it is in the area of northern Baja California near the San Andres fault, where the velocity model vectors tend to be parallel to the transform faults along the Gulf of California and the vectors over the Baja California Peninsula (Pacific Plate) decrease in velocity.
Similarly, the crustal deformation model fixed to the North American Plate has been computed with the aim of evidencing the displacement of the Pacific Plate and particularly of the Baja California Peninsula as well as Guadalupe Island, which are the main areas of interest on the Pacific Plate in Northwestern Mexico. The model is shown in Fig. 6d, where the displacement vectors show a northwestward direction. Also, the vectors near the plate boundary are parallel to transform faulting along the entire Baja California Peninsula. The vectors near to the plate boundary on the North American Plate show a displacement in a northwest direction also parallel to the transform faulting, however being considered as a fixed plate the velocity vectors tend to zero values and even with the velocity inversion show a change of direction, which is why towards the northeast (NE) the displacement vectors have a change of direction, and given that the model considers an elastic surface subject to the forces applied by the velocities in the GPS stations, it is that a supposed formation of two blocks with different velocity and direction is observed. However, it is not advisable to make any interpretation on the North American Plate, since it is the same plate considered fixed for the interpretation of the displacements in the Pacific Plate (Fig. 6d).
Model accuracy verification and error analysis
In order to validate the values obtained in the NW-MEXVEL interpolation model a direct point-to-point validation was performed, for which the model published by Argus et al. (2010) (GEODVEL 2010) has been considered as the basis for comparison, estimating the residuals as R = RNWMEXVEL model - RGEODVEL model, as well as calculating the RMSE of the proposed model and the corresponding correlation coefficients in the North and East velocity components between the NW-MEXVEL model (this study) and the GEODVEL model.
Figure 7a and 7b show the histograms of the interpolated velocities for the NW-MEXVEL model and the GEODVEL model to evidence the similarities by component and thus calculate the correlation coefficients which are shown graphically in Fig. 7c and 7d, with a strong correspondence by component between both models. The correlation coefficient for the North component is 0.989 to the linear best fit line describing the behavior of the data under equation y = 1.2926x + 0.3121 (Fig. 7c) while the correlation coefficient for the East component describes a strong positive correlation with a value of 0.994 to the best fit line described by y = 1.1743x + 1.7351. The calculation of the RMSE of the NW-MEXVEL model for each horizontal component of the 15’ x 15’ grid is 0.052 mm/year for the East velocity component and 0.116 mm/year for the North velocity component, which allows us to infer that the error that can be obtained for the East component in the absolute range of the interpolated velocities is lower than that of the North component, showing a better correlation in the East component between the interpolated NW-MEXVEL model of this study and the GEODVEL 2010 model. In addition, the residuals of the NW-MEXVEL models have a behavior similar to a normal distribution in both components, concentrating the highest frequency accumulation between values of -2.5 to 2.5 mm/year. Based on the above, the interpolation model for Northwest Mexico (NW-MEXVEL) shows reliability in the estimated uncertainties, as well as in the correlation coefficients between both components. The direct point-to-point validation is the method used here for the error analysis, which in the interpolation has been estimated by means of the RMSE. Here the best performance is in the East component, slightly higher than the North component, so the reliability of the NW-MEXVEL model will estimate with less uncertainty in the East component than in the North component. However, the model is reliable for displacements greater than 0.127 mm/year, this being the estimated global error.
Similarly, the values obtained in the NW-MEXVEL models fixed to the Pacific Plate have been validated. The residuals also fit a normal distribution, where the highest cumulative frequency is in values from − 1 to 1 mm/year, with a correlation coefficient of 0.9968 for the best linear fit corresponding to the equation y = 0.9195x + 0.4880 in the East component and a correlation coefficient of 0.9961 for the best linear fit corresponding to the equation y = 0.9233x − 1.0509 in the North component. Moreover, the RMSE has been calculated obtaining a value of 0.046 mm/year for the East component and 0.051 mm/year for the North component. Thus, the estimation of the velocities relative to the Pacific Plate is similar in both components, being reliable for displacements relative to the Pacific Plate greater than 0.068 mm/year. Finally, for the NW-MEXVEL model fixed to the North American Plate a correlation coefficient of 0.996 was obtained for the North component and the line of best fit described by y = 0.9084x − 2.3531. The East component has a correlation coefficient of 0.994 described by the line of best fit by y = 0.9121x + 1.9161. The result of the RMSE in the North component is 0.049 mm/year and for the East component an RMSE of 0.046 mm/year. Based on the above, the NW-MEXVEL model fixed to the North American Plate shows reliability in relation to the uncertainty of the error, as well as in the corresponding correlation coefficients. This model has been shown to be reliable for displacements on the Pacific Plate, and specifically in the Baja California Peninsula and Guadalupe Island greater than 0.067 mm/yr, this being the estimated global error.