Space Geodesy SIO 237 GNSS Theory and Analysis

Space Geodesy – SIO 237 GNSS Theory and Analysis Yehuda Bock January 23, 2020 Readings: Bock, Y. and D. Melgar (2016), Physical Applications of GPS Geodesy: A Review, Rep. Prog. Phys. 79, 10, doi: 10. 1088/0034‐ 4885/79/10/106801. Teunissen, P. and Montenbruck, O. eds. , 2017. Springer handbook of global navigation satellite systems. Springer.

GNSS research at IGPP • • • • Crustal deformation/earthquake cycle Transients: fault creep, episodic tremor and slip (ETS), postseismic Transition from coseismic to postseismic deformation Static and kinematic slip modeling (e. g. , Ridgecrest earthquakes) Earthquake magnitude scaling Real‐time high‐rate GNSS seismology and seismogeodetic data analysis (earthquake and tsunami early warning) GNSS meteorology (extreme weather – monsoons, atmospheric rivers) Ocean surface positioning for bathymetry and seafloor geodesy Ice sheet movements Subsurface water/hydrologic storage/drought Sea level rise (e. g. , Venice) Stratospheric observations of temperature for climate models GNSS reflectometry (soil moisture, tide gauges, snow)

Triangulation & GPS Prawirodirdjo, L. , Y, Bock, J. F. Genrich, S. S. O. Puntodewo, J. Rais, C. Subarya, and S. Sutisna (2000), One century of tectonic deformation along the Sumatran fault from triangulation and Global Positioning System surveys, J. Geophys. Res. , 105, 28, 343‐ 28, 361.

“San Andreas Discrepancy” 20 GPS surveys 1986‐ 1992; VLBI observations 1984‐ 1991 Feigl, K. L, D. C. Agnew, Y. Bock, D. Dong, A. Donnellan, B. H. Hager, T. A. Herring, D. D. Jackson, T. H. Jordan, R. W. King, S. Larsen, K. M. Larson, M. H. Murray, Z. Shen, and F. Webb (1993), Space geodetic measurement of crustal deformation in central and southern California, 19841992, J. Geophys. Res. , 98, 21, 677 -21, 712.

Survey-mode GPS Surveys Oblique subduction of the Pacific Plate beneath the Australian Plate at the Hikurangi Margin, North Island New Zealand. GPS velocities (black vectors) from s. GPS campaigns between 1991 and 2003 are shown relative to the Pacific Plate. Red vectors show estimated long‐term convergence rates (mm/yr) at the Hikurangi trough from a combination of GPS velocities and geologic slip rates. The GPS velocities in the eastern North Island shows that 50– 60 mm/yr of convergence occurs offshore of the northeastern North Island with the rates decreasing to ~20 mm/yr in the southern North Island. The southward decrease in offshore convergence rates is accompanied by an increase in upper plate shortening and produces rapid clockwise tectonic rotation of the eastern North Island relative to the bounding Pacific and Australian Plates. The margin‐parallel component of Pacific/Australia relative plate motion is accommodated by a combination of strike‐slip faulting and clockwise rotation of the eastern North Island. The inset shows the location of the 2016 Mw 7. 8 Kaikoura earthquake. Wallace, Laura M, Martin Reyners, Ursula Cochran, Stephen Bannister, Philip M Barnes, Kelvin Berryman, Gaye Downes, Donna Eberhart‐Phillips, Ake Fagereng, and Susan Ellis. 2009. 'Characterizing the seismogenic zone of a major plate boundary subduction thrust: Hikurangi Margin, New Zealand', Geochemistry, Geophysics, Geosystems, 10.

Continuous GNSS Typical continuous GPS station. Deeply‐anchored braced Southern California Integrated Network (SCIGN) monument and antenna (under the radome) of a typical c. GPS station (SIO 5 in La Jolla, for monitoring tectonic plate boundary deformation, earthquake early warning and GPS meteorology. The small white box on the monument’s vertical leg contains a MEMS accelerometer used for seismogeodesy. In the background are equipment enclosures, solar panels, a radio antenna for real‐time transmission of data and meteorological instruments. Photo courtesy of D. Glen Offield. Bock Y. , et al. (1997), Southern California Permanent GPS Geodetic Array: Continuous measurements of crustal deformation between the 1992 Landers and 1994 Northridge earthquakes, J. Geophys. Res. , 102, 18, 013‐ 18, 033.

Continuous GPS and In. SAR Coseismic motion detected by c. GPS and In. SAR for the 2010 Mw 7. 3 Landers earthquake, southern California. Solid arrows indicate surface horizontal displacements observed by 4 stations of the Permanent GPS Geodetic Array (PGGA) (Bock et al. 1997). Open arrows show the corresponding displacements from a dislocation model (Savage and Burford 1970) consisting of 7 linear fault segments to describe the rupture geometry. The contour lines show the modeled coseismic displacement field (in millimeters). The heavy line denotes the surface trace of the fault rupture and the dashed line is the Mw 6. 5 Big Bear earthquake's subsurface trace. The earthquake occurred 3 hours after the Landers event. The superimposed ERS-1 In. SAR image shows at least 20 phase fringes representing about 560 mm in displacement in the line of sight to the satellite (Massonnet et al. 1993). Massonnet, Didier, Marc Rossi, César Carmona, Frédéric Adragna, Gilles Peltzer, Kurt Feigl, and Thierry Rabaute. 1993. 'The displacement field of the Landers earthquake mapped by radar interferometry', Nature, 364: 138‐ 42.

Evolution of Continuous GPS 1999 -2018 Klein, K. , Y. Bock, X. Xu, D. Sandwell, D. Golriz, P. Fang, L. Su (2019), Transient deformation in California from two decades of GPS displacements: Implications for a three‐dimensional kinematic reference frame, J. Geophys. Res. , DOI: 10. 1029/2018 JB 017201.

GPS Signals Original GPS: Two frequencies, L 1 and L 2 C/A‐code modulation of L 1 (L 1 C) @1. 023 MHz – 300 m wavelength P‐code modulation of L 1 and L 2 signals @10. 23 MHz – 30 wavelength L 1 = 1, 575. 42 MHz (154*10. 23) (~19 cm wavelength) L 2 = 1, 227. 6 MHz (120*10. 23) (~24 cm wavelength)

GPS Observations The basic measurement is the time required for the GNSS signal to propagate from a satellite to the receiver. This can be obtained by tracking the PRN code modulation of the signal. Within the receiver, a local copy of the PRN sequence is generated, which is continuously compared and aligned with the signal received from the satellite. This tracking loop provides continuous measurements of the instantaneous code phase and hence the transmission time corresponding to the currently received signal. By comparing this time with the local receiver time, the signal propagation time, and – upon multiplication by the speed of light – the distance or range from receiver to satellite are obtained. Pseudorange: A measure of the difference between the receiver clock at signal reception and the satellite clock at signal transmission (scaled by the speed of light). The pseudorange measures the satellite–receiver distance, the precision of which is in the dm‐range. Carrier phase: A measure of the instantaneous beat phase and the accumulated number of zero‐crossings obtained after mixing with a reference signal of the nominal frequency. Changes in carrier phase over time reflect the change in pseudorange but are ~ 2 orders of magnitude more precise. In case of interrupted tracking the accumulated cycle count is lost and the carrier‐phase measurements exhibit a cycle slip.

Phase and Group Velocity

Ionosphere-free Observables

GPS Time and Atomic Time Atomic time is the basis of a uniform time scale on the Earth and is realized by International Atomic Time (TAI); TAI is a continuous time scale, maintained by the Bureau International des Poids et Mesures (BIPM), using data from about three hundred atomic clocks in over 50 national laboratories. The fundamental interval unit of TAI is one SI second defined at the 13 th general conference of the International Committee of Weights and Measures in 1967 as the “duration of 9, 192, 631, 700 periods of the radiation corresponding to the transition between two hyperfine levels of the ground state of the cesium 133 atom. ” The SI day is defined as 86, 400 seconds Because TAI is a continuous scale it does not maintain synchronization with the solar day (universal time) since the Earth’s rotation is slowing. This problem is solved by defining Universal Coordinated Time (UTC), which runs at the same rate as TAI but is periodically incremented by leap seconds. Today (January 2020), UTC + 37 seconds = TAI. The time signals broadcast by the GPS satellites are synchronized with atomic clocks at the GPS Master Control Station in Colorado. Global Positioning System Time (GPST) is a continuous time reference that was set to 0 h UTC on 6 January 1980 and it is not incremented by UTC leap seconds. GPST+19 seconds = TAI. There have been a total of 18 leap seconds, the last one on December 31, 2016. Today (January, 2020), GPST is ahead of UTC by 18 seconds (GPST=UTC+18 seconds). GPS phase and pseudorange measurements are given in GPST.

Relativistic Effects on GNSS clocks Special Relativity: The satellite atomic clocks fall behind clocks on the ground by about 7 mas/d due to the time dilation effect of their relative motion (velocity difference). General Relativity: The curvature of spacetime make the satellite clocks appear to run faster by about 45 mas per day. Net effect: 38575. 008 ns/d; position error of about 10 km a day if not taken into account. In GPS analysis a nominal frequency offset and a periodic relativistic correction derived as a dot product of the satellite position and velocity vectors compensate for this effect. [Important for satellite clock estimation and timing: Provision must be made for departures in the GNSS satellite orbits (in particular in its semi‐major axis and its distance from the Earth’s mass), which can reach up to 10 ns/d in apparent clock rate, and for the Earth’s gravity field oblateness. Ignoring these effects can result in errors in the relativistic transformation of 0. 2 ns/d and periodic errors of 0. 1 and 0. 2 ns, with periods of about 6 h and 14 d, respectively. ]

Inertial and Earth-Fixed Frames

Satellite Motion: Inertial Reference Frame

Terrestrial Reference Frame By convention, the X and Y axes are in the Earth’s equatorial plane with X in the direction of the point of zero longitude and the Z axis in the direction of the Earth’s pole of rotation. Plate Tectonic Frame

Transformations Geodetic to Cartesian: with semimajor axis a, inverse flattening (1/f), and e the ellipsoidal eccentricity. The U. S. Department of Defense WGS 84 system used for GPS is consistent within ± 1 m with the ITRF, which is internally consistent at the sub‐cm level. The WGS 84 ellipsoidal parameters are semi‐major axis a=6378137 and 1/f=298. 257 223 563. Conversion to local system:

Precise Point Positioning (PPP) Between station differencing Network Positioning (NP) Double differencing

GPS Observation Quartet

Phase Ambiguity Resolution

Deviations from the ideal quartet • Non-dispersive neutral atmosphere (the “troposphere”) of direct signal, “dry” (modeled) and “wet” (estimated) components. Any common-mode errors eliminated by differencing between stations, total effect in PPP • Multipath off of objects near the GPS antenna • Station clock biases (eliminated by differencing between satellites) • Satellite clock biases (eliminated by differencing between stations, not possible in PPP) • Satellite orbital errors (estimated as part of global network adjustment – need for PPP) • Differential code (pseudorange) biases (DCBs), to account for group delay differences between the signals tracked by the receiver and those of the clock reference signal (for PPP analysis) • Measurement errors • Centering, leveling and orientation of GPS antennas • Receiver and antenna phase center offsets from • center of mass of satellite, variations in phase centers (PCVs) (use IGS calibrations – not perfect) • Monument stability

Modeled as part of GPS analysis • Tidal effects, gravitational potential of the Sun and Moon causes loading displacements o Solid Earth tides (up to 1 m) o Ocean tides (~100 m) o Atmospheric tides (several mm) o Pole tide (~5 mm) • Non-tidal effect, mass loadings (heating of atmosphere from Sun’s radiation) o Atmospheric loading – atmospheric pressure on solid Earth, several mm o Ocean loading – ocean‐bottom pressure causes coastal displacements, < 1 mm Dong, D. , Fang, P. , Bock, Y. , Cheng, M. K. and Miyazaki, S. I. , 2002. Anatomy of apparent seasonal variations from GPS‐derived site position time series. Journal of Geophysical Research: Solid Earth, 107(B 4), pp. ETG‐ 9.

Phase Observation Equations

PPP Observation Equations – Ionosphere-Free These are called fractional cycle biases

Modernized GPS Original: Two frequencies, L 1 and L 2 C/A‐code modulation of L 1 (L 1 C) @1. 023 MHz – 300 m P‐code modulation of L 1 and L 2 signals @10. 23 MHz – 30 m Broadcast ephemeris modulated on L 1 Triple Frequency: L 1 = 1, 575. 42 MHz (154*10. 23) L 2 = 1, 227. 6 MHz (120*10. 23) L 5 = 1, 176. 45 MHz (115*10. 23) L 1 C New L 1 C/A signal (1575. 42 MHz) L 2 C New L 2 C/A signal (1575. 42 MHz) C/A‐code modulation of L 2 for redundancy and ionospheric corrections for users with C/A code receivers. M‐code for more secure military access includes two new encrypted P(Y) signals transmitted at both L 1 and L 2.

GLONASS Observations GLONASS Complicated by the frequency division multiple access (FDMA) modulation. Makes use of slightly different signal frequencies on about 15 distinct channels. The individual channels are separated by 562: 5 k. Hz and 437: 5 k. Hz for L 1 and L 2. May result in interfrequency channel biases (IFCB) for both code and phase observations. These biases affect the generation of precise orbit and clock products as well as the use of GLONASS observations for precise point positioning and ambiguity resolution. Useful for static and real-time kinematic positioning for surveying applications especially in urban canyons and other obstructed areas

GNSS Constellations ~160 satellites ~400 signals Teunissen, P. and Montenbruck, O. eds. , 2017. Springer handbook of global navigation satellite systems. Springer.

GNSS International Standards & Infrastructure RINEX – Receiver Independent Exchange Format for phase, pseudorange & metadata RTCM – Radio Technical Commission for Maritime for real‐time transmission of data IGS – International GNSS Service for precise satellite orbits, EOP and reference frames SINEX – Solution Independent Exchange Format for estimated coordinates and their covariances Software Hatanaka compression – for compressing RINEX files for archive teqc – quality control and manipulation of RINEX files IGS Network http: //www. igs. org/network

Time Series Analysis

2003 Mw 6. 5 San Simeon earthquake Raw Time Series: Station MIDA in the Parkfield region 2004 Mw 6. 0 Parkfield earthquake

Modeled Time Series: Station MIDA in the Parkfield region 2004 Mw 6. 0 Parkfield earthquake – postseismic motion

Residual Time Series: Station MIDA in the Parkfield region

Modeled Time Series: Station ALBH in British Columbia ETS Signal

Continuous GNSS Displacement Signals Horizontal Transient Horizontal Motions. Accumulated residual horizontal displacements from the start of 2010 until mid-October, 2019 for California and Nevada (Klein et al. , 2019) showing transients and secular differences with the secular model of Zeng and Shen (2017)*. We see significant secular differences with the model, for example, in Cascadia, the Northern Basin and Range, and the Santa Maria basin, that can then be used for further secular fault slip models. Transients include postseismic motion for 4 earthquakes in this 10 -year period, and the horizontal expression of magmatic inflation at Long Valley Caldera and Lassen Peak. *Zeng, Y. , and Shen, Z. ‐K. (2017), A fault‐based model for crustal deformation in the western United States based on a combined inversion of GPS and geologic inputs. BSSA, 107(6), 2597– 2612. https: //doi. org/10. 1785/0120150362

Continuous GNSS Displacement Signals – Vertical Area Klein, K. , Y. Bock, X. Xu, D. Sandwell, D. Golriz, P. Fang, L. Su (2019), Transient deformation in California from two decades of GPS displacements: Implications for a three‐dimensional kinematic reference frame, J. Geophys. Res. , DOI: 10. 1029/2018 JB 017201. Source a. Cascadia Subduction Tectonic b. Mount Lassen Magmatic c. Sacramento Valley Tectonic d. Coastal Ranges Tectonic e. Central Valley Groundwater pumping + drought f. Sierra Nevada Ranges Tectonic g. Long Valley Caldera Magmatic h. Transverse Ranges Tectonic i. Eastern California Shear Zone 1999 Mw 7. 1 Hector Mine postseismic j. Southern terminus SAF system 2010 Mw 7. 2 El Mayor-Cucapah postseismic k. Santa Maria Basin Groundwater pumping l. Ventura Basin Tectonic/ Groundwater pumping m. Los Angeles/Santa Ana basin Groundwater pumping n. Southern California Tectonic o. South of Salton Sea Geothermal Field p. North San Francisco Bay Area Tectonic/ Groundwater pumping

Seismogeodesy Melgar, D. , B. W. Crowell, Y. Bock, and J. S. Haase (2013), Rapid modeling of the 2011 Mw 9. 0 Tohoku‐oki earthquake with seismogeodesy, Geophys. Res. Lett. , 40, 1‐ 6. doi: 10. 1002/grl. 50590

GNSS Meteorology Atmospheric River Estimating Precipitable Water from ground GPS networks 2 km 7° Δttotal = Δtgeom + Δtiono + Δttrop + … GNSS station 32 km Zenith hydrostatic delay = f(surface pressure) Zenith wet delay TD(θ) = ZHD mh(θ) + ZWD mw(θ) Total trop delay Mapping functions In solving for this we estimate this PW = κ × ZWD 1/κ = 10‐ 6 × ρRv[(k 3/Tm) + k 2’] ≈ 6. 5 Mean atmospheric temperature With surface pressure and a surface temperature, we derive PWV from zenith troposphere delay

GAMIT Software Station. info : Metadata Session. info : Model and run‐time settings Sittbl: coordinate constraints igs_08. atx : Antenna phase center mapping ut 1. usno : ERA Pole. usno : Polar Motion Solar ephemeris Lunar ephemeris Nutation tables Leap seconds http: //sopac‐csrc. ucsd. edu/index. php/gambit‐globk/