A Review of GNSS/GPS in Hydrogeodesy: Hydrologic Loading Applications and Their Implications for Water Resource Research.

Hydrogeodesy, a relatively new field within the earth sciences, is the analysis of the distribution and movement of terrestrial water at Earth's surface using measurements of Earth's shape, orientation, and gravitational field. In this paper, we review the current state of hydrogeodesy with a specific focus on Global Navigation Satellite System (GNSS)/Global Positioning System measurements of hydrologic loading. As water cycles through the hydrosphere, GNSS stations anchored to Earth's crust measure the associated movement of the land surface under the weight of changing hydrologic loads. Recent advances in GNSS-based hydrogeodesy have led to exciting applications of hydrologic loading and subsequent terrestrial water storage (TWS) estimates. We describe how GNSS position time series respond to climatic drivers, can be used to estimate TWS across temporal scales, and can improve drought characterization. We aim to facilitate hydrologists' use of GNSS-observed surface deformation as an emerging tool for investigating and quantifying water resources, propose methods to further strengthen collaborative research and exchange between geodesists and hydrologists, and offer ideas about pressing questions in hydrology that GNSS may help to answer.

Introduction

As Earth's climate warms at unprecedented rates, changes in the magnitude, frequency, timing, and type of precipitation and evapotranspiration have been widely observed across the globe (Barnett et al., 2005; Dyurgerov, 2021; Easterling et al., 2000; He et al., 2021; Milly et al., 2015, 2008; Mote et al., 2005; Oerlemans, 2005; F. Zhang et al., 2021). These changes are projected to intensify, impacting the global hydrologic cycle (Asadieh & Krakauer, 2017; Chikamoto et al., 2017; Dong et al., 2021; He et al., 2021; Shannon et al., 2019; Swain et al., 2018; Yoon et al., 2015), making it ever more important to monitor changes in water availability across human and natural systems. Hydrogeodesy is a promising new field that employs geodetic observations of Earth's shape, orientation, and gravitational field to quantify changes in the distribution and storage of fresh water on Earth's surface, in soils, and in groundwater reservoirs. Since traditional hydrologic methods do not directly measure total integrated water storage, hydrogeodesy has the potential to fill a key observational gap in hydrologic studies.

Earth's shape changes over time due to mass redistribution at the surface and near surface. Geodetic networks use GNSS (Global Navigation Satellite Systems that include the GPS) to observe those changes over time. GNSS position time series record the small but measurable load‐induced displacements of Earth's surface that result from the deformation of the earth under varying surface loading due to sea level changes, atmospheric circulation, glacier and ice cap evolution, snow and rainfall, and changes in water storage in rivers, lakes, soils, and groundwater aquifers.

Earth deforms both elastically (instantaneously and reversibly) and inelastically (irreversibly and over time) in response to surface loads with the largest deformation occurring nearest the load source. The elastic component of load‐induced Earth deformation can be characterized by load Love numbers (e.g., Farrell, 1972; Longman, 1963) and load Green's functions (Farrell, 1972) that depend on the density and elastic structure of the solid earth. When a load is increased (e.g., by adding more mass), the land surface depresses and nearby GNSS stations move downward and toward the load. When a load is decreased (e.g., by removing mass), the land surface rebounds and nearby GNSS stations move upward and away from the load as illustrated in Figure 1. Assuming a linear elastic Earth rheology (considered a good assumption for loading timescales of less than ∼1 year; e.g., Argus et al., 2014, 2017; Borsa et al., 2014; Chanard et al., 2014; Fu et al., 2015; Wahr et al., 2013), surface displacements recorded at a given GNSS station will reflect the instantaneous linear combination of deformation due to all loads from the near field (tens of km or less) to far field (hundreds to thousands of km) (Knappe et al., 2019).

Figure 1

Conceptual model of how Earth's surface and attached Global Navigation Satellite System stations move in response to loading, assuming strictly elastic Earth rheology. (a) Schematic uses the example of snowpack as the sole hydrologic load. Under the accumulating load of winter snowpack, the surface subsides and moves toward the load. When the snowpack melts, deformation reverses and the surface rebounds and moves away from the load. (b) Modeled displacement under a disk load of radius 50 km and an equivalent water depth of 1 m. Arrows indicate the direction and magnitude of predicted horizontal displacements and colors indicate predicted vertical displacements. Note that the magnitude of vertical displacement (downward movement represented by negative values) is greatest at the center of the load and decreases with distance; conversely, the magnitude of horizontal displacement is zero at the center of the load and greatest at the radius of the disk, beyond which it too decreases with distance.

Traditionally, geodetic studies using continuous GNSS position time series have focused on geologic studies characterizing tectonic plate motion, volcanic inflation and deflation, earthquake deformation, and natural disasters like landslides. At broader scales, GNSS helps constrain models of glacial isostatic rebound, which can be used to quantify absolute sea level rise and infer properties of Earth's mantle and crust (see Bock & Melgar, 2016; Segall & Davis, 1997 and references therein). More recently, geodetic data sets have also revealed the contribution of hydrologic loading to observe Earth deformation. For instance, an exploration of seasonal drivers of interseismic strain buildup on faults in Japan found that the surface loading from snow is largely responsible for the surface deformation observed by GNSS (Heki, 2001). Furthermore, Dong et al. (1997) and van Dam & Wahr (1998) demonstrated that the redistribution of water on Earth's surface shifts Earth's geocenter and causes variation in its gravitational field, respectively.

Ongoing increases in the number of GNSS ground stations and GNSS satellites have correspondingly improved the resolution, sensitivity, and availability of surface displacement observations over the globe, offering advantages over other geodetic methods. Early analysis predicted that global changes in continental water storage could produce tens of millimeters of surface deformation in agreement with and responsible for a significant fraction of GNSS‐observed surface displacement from 147 GNSS stations (van Dam et al., 2001). Studies of hydrologic loading signals have significantly expanded in number, scope, and geodetic techniques over the last few decades, popularizing hydrogeodesy as an emerging and impactful field of scientific research.

While not a focus of this review, several other space and airborne geodetic tools can be applied to study total water storage on land and its various components. For instance, satellite and airborne radar altimetry can be used to measure snow depth (Treichler & Kääb, 2017) and to estimate soil moisture content (e.g., the Soil Moisture Active Passive (SMAP) project; Chan et al., 2018). Interferometric synthetic aperture radar (InSAR) techniques can measure subsidence associated with groundwater depletion and provide a means of inferring elastic storage (e.g., Amelung et al., 1999; Bawden et al., 2001). Subsidence rates calculated with InSAR are also commonly combined with GNSS to verify subsidence rates or estimate three‐dimensional deformation from groundwater extraction (e.g., Bui et al., 2021; Hinderer et al., 2020; Lei et al., 2021; Yang et al., 2019 ). Light detection and ranging (LiDAR) data can be used to estimate and scale evapotranspiration (Sutherland et al., 2017). Space‐based gravitational methods, like the Gravity Recovery and Climate Experiment (GRACE) mission, quantify large‐scale water storage on land at approximately 300 km spatial and monthly temporal resolution (Dahle et al., 2019; Frappart & Ramillien, 2018; Ramillien et al., 2008; Save et al., 2016; Soltani et al., 2021; Watkins et al., 2015). GNSS complements and extends these other techniques, capturing changes in terrestrial water storage (TWS), defined as the total storage of water on land, with both high spatial and high temporal resolutions in many regions, over almost the same multi‐decade timescale provided by GRACE. Therefore, this review specifically focuses on the hydrologic application of GNSS to observe hydrologic loading and its implications for water resource research.

While geophysical applications of GNSS geodesy have been reviewed previously (e.g., Blewitt, 2015; Bock & Melgar, 2016; Riguzzi et al., 2020), little focus has been placed on hydrogeodesy and its societal implications for water resources. Therefore, this review presents the current state of hydrogeodesy with a specific focus on applications of GNSS that are most immediately relevant and available to the hydrologic science and water‐management communities. To facilitate hydrologists' use of GNSS‐observed surface deformation as a tool for investigating water resources, we begin with a brief explanation of GNSS and hydrologic loading. We provide an assessment of recent hydrologic applications of GNSS‐based hydrogeodesy, such as links between GNSS position time series and climatic drivers, TWS estimates from GNSS load displacements, and the impact of meteorological events, such as hurricanes and atmospheric rivers on hydrologic loading. We stress the importance of GNSS observations and methods in hydrogeodesy and discuss recent advancements and limitations of the GNSS analysis of hydrologic loads. Finally, we consider future directions in hydrogeodesy, offer ideas about pressing questions in the hydrologic sciences that GNSS may help to answer, and propose methods to further strengthen collaboration between geodesists and hydrologists.

Fundamentals of GNSS/GPS

This section provides an overview of GNSS measurements and basic GNSS processing and is aimed at hydrologists looking to understand the fundamentals of the technique. Experts in geodesy and GNSS may wish to skip this section. GNSS is the umbrella term for constellations of satellites that send positioning and timing data from their orbits around the globe, allowing GNSS receivers to calculate their own locations. Several countries have their own GNSS satellite constellations. The United States operates GPS, the oldest GNSS satellite constellation, dating back to its first satellite launch in 1978. Russia operates GLONASS, the European Union operates Galileo, China operates BeiDou (also known as Compass), Japan operates the Quazi‐Zenith Satellite System (QZSS) with limited coverage of the Asia‐Pacific region, and India operates the Indian Regional Navigation Satellite System (IRNSS, also operationally referred to as NavIC) with limited coverage extending approximately 1500 km past the borders of India. With the recent expansion of GNSS satellite constellations, multiple satellite constellations can now be used in combination, depending on the location and capabilities of GNSS receivers and the number of visible satellites. Here, we use the term GNSS in a general sense to refer to position, navigation, and timing (PNT) data from any single satellite constellation or combination of multiple constellations.

GNSS works by transmitting a microwave signal (frequencies within 1–2 GHz range), containing the precise time and orbital position of the satellites, to antennas and receivers on Earth, typically on the ground. The antennas require unobstructed views of at least four separate satellites for calculation of the antenna's latitude, longitude, altitude, and timing error from the transmitted signal. Using precise satellite orbits and the speed of light, the receiver calculates the pseudorange, an approximation of the distance between a satellite and GNSS station, from the measured travel time from each satellite to the GNSS antenna after adjustment for delays in signal propagation through Earth's atmosphere and ionosphere. The pseudorange also naturally includes clock synchronization errors that must be corrected using differencing of data. Trilateration of the pseudorange is then used to determine the station position in three dimensions with millimeter‐level positional accuracy (Blewitt, 2015).

GNSS used for scientific research differs from GNSS commonly used in everyday devices like smartphones and car navigation, which have positional accuracies on the order of a few meters. This difference is primarily the result of specialized processing algorithms used in scientific applications that leverage dual‐frequency satellite signals and precise orbit and clock estimates that are available for post‐processing to more accurately correct for signal‐propagation delays through the troposphere and ionosphere. For instance, scientific applications of the United States' GPS satellite constellation commonly use both the L1 (1575.42 MHz) and the L2 (1227.60 MHz) frequencies (new GPS receivers can now also use the modernized L1C, L2C, and L5 frequencies), while most civilian GPS only uses signals from the L1 frequency. For more information on the technical and theoretical details of GNSS, see Blewitt (1993, 2015) and Misra and Enge (2006), and see Segall and Davis (1997), Bock and Melgar (2016), and Herring et al. (2016) for the history and evolution of GNSS applications to geophysical processes and geohazards.

GNSS Networks and Data Access

Thousands of continuous GNSS stations have been deployed worldwide (Figure 2), each consisting of a GNSS antenna and optional protective radome cover, a GNSS receiver, a power supply, and often a means of transmitting data back to the station operator. Fixed GNSS antennas are anchored to bedrock or attached to stainless steel rods driven multiple meters beneath land surface (Figure 3) to ensure that they are not sensitive to shallow subsurface processes and move with the ground, therefore measuring elastic Earth deformation. Several networks of fixed continuous or semicontinuous GNSS stations operate as part of large‐scale observatories from which data are publicly available, allowing for a range of geodetic applications (Blewitt et al., 2018). For instance, several networks of GNSS stations are located across North America (Murray et al., 2019), the largest of which is the EarthScope Plate Boundary Observatory (PBO) across the contiguous United States (CONUS) and Alaska (Herring et al., 2016). PBO stations are now part of the greater Network of the Americas (NOTA) with approximately 1278 continuous GNSS stations across the CONUS, Alaska, Canada, Mexico, and the Caribbean (Murray et al., 2019; https://www.unavco.org/projects/major-projects/nota/nota.html). Other networks around the world include the Crustal Movement Observation Network of China (CMONOC), the EUREF Permanent Network (EPN) in Europe, the GNSS Earth Observation Network System (GEONET) in Japan, and numerous smaller, regional networks. The International GNSS Service (IGS), a cooperative global network, generates key precise GNSS products, such as orbits and clock information for GNSS data processing. Using IGS data products, several analysis centers process raw GNSS data to estimate time series of GNSS positions in the east, north, and up directions that are made publicly available for scientific research (e.g., Blewitt et al., 2018; Bock et al., 2021; Herring et al., 2016); however, GNSS data from the CMONOC are not publicly available in a raw Rinex form.

Figure 2

Global map showing Global Navigation Satellite System stations (red triangles) for which processed data product time series can be retrieved from the Nevada Geodetic Laboratory of the University of Nevada Reno. Note that these stations span variable time ranges with variable degrees of data completeness and data latencies.

Figure 3

Examples of continuous Global Navigation Satellite System (GNSS) stations mounted directly into a rock outcrop (a) and anchored into the ground in a clearing with braced mounts (b and c). All stations are firmly anchored in open sky areas to ensure a clear line of sight between antennas and GNSS satellites. A radome covers the antenna shown in (b). Each station is also equipped with a receiver and power supply, such as the solar panels and batteries shown in (c). Image (b) was taken from http://sopac-csrc.ucsd.edu/index.php/eps-em/.

Raw satellite data and derived data products (e.g., time series of position estimates) from continuous GNSS stations can be mined for various scientific research purposes. In this paper, we focus on applications involving hydrologic loading from surface water, snow, soil moisture, and groundwater. The wide range of GNSS‐based geodetic applications necessitates high‐quality and consistent site performance, data management, and initial processing to create GNSS data products that are precise and versatile (Blewitt, 2015; Blewitt et al., 2018; Herring et al., 2016). Satellite data are transformed into estimates of three‐dimensional (east, north, and up) site position through time relative to global or regional reference frames (Altamimi & Collilieux, 2009; Altamimi et al., 2002; Blewitt, 2003; Herring et al., 2016). Most of these processed time series of GNSS site position with applied standard corrections are publicly available for scientific research from several analysis centers (Table 1) via FTP or HTTP services (e.g., Blewitt et al., 2018; Herring et al., 2016; Zumberge et al., 1997).

Table 1

Details of Papers That Use Global Navigation Satellite System Displacement Time Series to Study Hydrologic Loading Induced Surface Deformation

Citation	Journal	Study location	Number of stations	Pre‐processing from/with	Reference frame
(Blewitt et al., 2001)	Science	Global	66	IGS Analysis Centers	CF
(Heki, 2001)	Science	NE Japan	75	Geographical Survey Institute (GSI) via Bernese	TSKB IGS station
(van Dam et al., 2001)	Geophysical Research Letters	Global	147	IGS Analysis Centers	Not found
(Dong et al., 2002)	Journal of Geophysical Research	Global	128	SOPAC via GAMIT/GLOBK, JPL via GIPSY/OASIS, & GSI via Bernese	ITRF1997
(Elósegui, 2003)	Geophysical Research Letters	Great Salt Lake	5	Not found	Not found
(Wu et al., 2003)	Geophysical Research Letters	Global	200	JPL via GIPSY/OASIS	ITRF2000
(Bevis et al., 2004)	Physics of the Earth and Planetary Interiors	Lago Laja, Chile	1	GAMIT/GLOBK	Attached to stable core of S American plate
(Davis, 2004)	Geophysical Research Letters	Amazon Basin	12	JPL via GIPSY/OASIS	Fiducial free
(Heki, 2004)	Book Series—The State of the Planet: Frontiers and Challenges in Geophysics	Japan	GEONET	GSI via Bernese	TSKB IGS station
(Bevis et al., 2005)	Geophysical Research Letters	Manaus, Brazil on Amazon River	1	GAMIT/GLOBK	Not found
(Grapenthin et al., 2006)	Geophysical Research Letters	Iceland	4	Bernese	REYK station
(van Dam et al., 2007)	Journal of Geophysical Research	Europe	51	IGS Analysis Centers	current ITRF
(Bettinelli et al., 2008)	Earth and Planetary Science Letters	Himalaya	Not found	Bernese	ITRF2000
(Nordman et al., 2009)	Journal of Geodynamics	Fennoscandia	7	JPL via GIPSY/OASIS	ITRF2005
(Tregoning et al., 2009)	Geophysical Research Letters	Global	80	GAMIT/GLOBK	ITRF2005
(Steckler et al., 2010)	Geophysical Research Letters	Ganges‐Brahmaputra Delta region	2	MIT via GAMIT/GLOBK	ITRF2000
(Tesmer et al., 2011)	Journal of Geodesy	Global	115	Bernese	IGS05
(Davis et al., 2012)	Journal of Geophysical Research	Alaska	1	SOPAC via GAMIT/GLOBK	Not found
(Fu & Freymueller, 2012)	Journal of Geophysical Research	Nepal, S side of Himalaya	32	UNAVCO via GIPSY/OASIS	ITRF2008
(Fu et al., 2012)	Geophysical Research Letters	S Alaska	64	GIPSY/OASIS	ITRF2008
(Nahmani et al., 2012)	Journal of Geophysical Research	West Africa	10	University of La Rochelle Analysis Center Consortium (ULR) via GAMIT/GLOBK	ITRF2005
(Dill & Dobslaw, 2013)	Journal of Geophysical Research: Solid Earth	Global	53	IGS Analysis Centers	CF
(Fu et al., 2013)	Geophysical Research Letters	Amazon Basin and SE Asia	Not found	JPL via GIPSY/OASIS	IGS08
(Ouellette et al., 2013)	Water Resources Research	Western US (CA/NV/UT/ID/WY/WA)	6	SOPAC via GAMIT/GLOBK	Not found
(Wahr et al., 2013)	Journal of Geophysical Research: Solid Earth	Lake Shasta (CA) and Greenland	4 (CA) & 1 (Greenland)	Plate Boundary Observatory Analysis Centers (CA) and JPL via GIPSY/OASIS (Greenland)	ITRF2008
(Argus et al., 2014)	Geophysical Research Letters	CA/NV/OR/WA	922	JPL via GIPSY/OASIS	ITRF2008
(Borsa et al., 2014)	Science	Western US	771	NMT	Not found
		109°W‐Pacific Coast
(Chanard et al., 2014)	Journal of Geophysical Research: Solid Earth	Nepal	31	GAMIT/GLOBK	ITRF2005
(Chew & Small, 2014)	Geophysical Research Letters	High Plains of US	15	UNAVCO	Not found
(Döll et al., 2014)	Surveys in Geophysics	Global	175	Not found	Not found
(Birhanu & Bendick, 2015)	Journal of Geophysical Research: Solid Earth	Ethiopia and Eritrea	16	GAMIT/GLOBK	ITRF2008
(Fu et al., 2015)	Journal of Geophysical Research: Solid Earth	Washington and Oregon	Not found	JPL via GIPSY/OASIS	IGS08
(Jin & Zhang, 2016)	Surveys in Geophysics	SW USA (28–43°N & 92–117°W)	Not found	JPL via GIPSY/OASIS	Not found
(Silverii et al., 2016)	Journal of Geophysical Research: Solid Earth	Apennines, Italy	174	JPL via GIPSY/OASIS	Eurasian TF aligned with IGS08
(Zhang et al., 2016)	Sensors	SW China (21.5–29.5°N & 97.5–105.5°E)	34	GAMIT/GLOBK	ITRF2008
(Argus et al., 2017)	Journal of Geophysical Research: Solid Earth	California	1276	JPL and NGL via GIPSY/OASIS	ITRF2008
(Han, 2017)	Journal of Geophysical Research: Solid Earth	Australia	120	NGL	Consistent with GRACE
(Han & Razeghi, 2017)	Journal of Geophysical Research: Solid Earth	Australia	114	NGL	ITRF2008
(Jiang et al., 2017)	Scientific Reports	Yunnan, China	26	GAMIT/GLOBK	IGb08
(Johnson et al., 2017)	Science	California	661	JPL via GIPSY/OASIS	ITRF2008
(Li et al., 2017)	Advances in Space Research	Eurasian plate	32	GAMIT/GLOBK	ITRF2005
(Wang et al., 2017)	Hydrology and Earth System Sciences	North China Plain	29	GIPSY/OASIS	ITRF2008
(Zhan et al., 2017)	Geophysical Journal International	Yunnan, China	27	GAMIT/GLOBK	ITRF2008
(Enzminger et al., 2018)	Water Resources Research	Western CONUS (25–53°N and 95–125°W)	1395 locations	NA	NA
(Ferreira et al., 2018)	Science of the Total Environment	Brazil	39	GIPSY/OASIS	IGS08
(Larochelle et al., 2018)	Journal of Geophysical Research: Solid Earth	Arabian Peninsula and Nepal Himalaya	14 for each site	NGL	ITRF2008
(Milliner et al., 2018)	Science Advances	South Texas and Louisiana	198	JPL via GIPSY/OASIS	ITRF2008
(Adusumilli et al., 2019)	Geophysical Research Letters	CONUS	1336	GAGE	Not found
(Enzminger et al., 2019)	Geophysical Research Letters	Sierra Nevada Mtns (31–50°N and 125–103°W)	924	GAGE	Not found
(Ferreira et al., 2019)	Remote Sensing	South America	397	NGL via GIPSY/OASIS	IGS08
(Fok & Liu, 2019)	Remote Sensing	SW China	34	GAMIT/GLOBK	ITRF2008
(Knappe et al., 2019)	Water Resources Research	ID/MT/WY	41	NGL via GIPSY/OASIS	ITRF2008
(Springer et al., 2019)	Journal of Geodesy	Europe	230	Italian Space Agency (ISA)	IGb08
(Argus et al., 2020)	Journal of Geophysical Research: Solid Earth	Great Lakes	3658	NGL	IGS14
(Hsu et al., 2020)	Earth and Planetary Science Letters	Taiwan	207	GAMIT/GLOBK	ITRF2008
(Knowles et al., 2020)	Journal of Geophysical Research: Solid Earth	Amazon Basin	23	NGL	IGS08
(Koulali & Clarke, 2020)	Journal of Geodesy	Antarctica	7	GAMIT/GLOBK	ITRF2014
(Lai et al., 2020)	Remote Sensing	Taiwan	Not found	GAMIT/GLOBK	ITRF2008
(Yin et al., 2020)	Journal of Geophysical Research: Solid Earth	Great Basin and Upper Colorado	77 (GB) & 10 (UC)	NMT	CF
(Jiang et al., 2021)	Remote Sensing of Environment	Yunnan, SW China	53	China Earthquake Data Center via GAMIT/GLOBK	CF
(Nespoli et al., 2021)	Journal of Hydrology	Nespoli, Italy	57	GAMIT/GLOBK	ITRF2014
(Ray et al., 2021)	Acta Geodaetica et Geophysica	NE India and Nepal Himalaya	36	GAMIT/GLOBK	ITRF2008
(Xue et al., 2021)	Geophysical Journal International	Great Lakes region	57	NGL	IGS14
(Yin et al., 2021)	Water Resources Research	Great Basin & Upper Colorado	77 (GB) & 10 (UC)	NMT	CF
(Young et al., 2021)	Journal of Geophysical Research: Solid Earth	Great Salt Lake	17	NGL	IGS14
(Zhan et al., 2021)	Journal of Geophysical Research: Solid Earth	Japan	∼1300	GSI via Bernese and NGL via GIPSY/OASIS	IGb14 & IGS14

Note. Table includes the journal where the study was published, the study location, the number of GNSS stations used, either the processing center(s) from which the data products were accessed or the software used to create the time series, and the GNSS reference frame. Frequently used acronyms: International GNSS Service (IGS), International Terrestrial Reference Frame (ITRF), Center of Surface Figure (CF), Geodesy Advancing Geosciences and EarthScope (GAGE), Scripps Orbit and Permanent Array Center (SOPAC), Jet Propulsion Laboratory (JPL), New Mexico Institute of Mining and Technology (NMT), Massachusetts Institute of Technology (MIT), and Nevada Geodetic Laboratory (NGL).

GNSS data processing typically removes signals that are large and well‐characterized or unrelated to solid Earth processes. Corrections for Earth's body and load tides are a standard procedure in GNSS processing, typically applied according to International Earth Rotation and Reference System Service (IERS) conventions (Petit & Luzum, 2010). Corrections for signal delays through the troposphere and ionosphere are also estimated although methods and models vary (e.g., Boehm et al., 2007, 2006; Martens et al., 2020; Steigenberger et al., 2009; Tregoning & Herring, 2006). Some derived data products also remove offsets associated with earthquakes and equipment maintenance (e.g., Bock et al., 2021; Herring et al., 2016).

Significant differences in data products can arise due to differences in processing software (Bertiger et al., 2020; Dach & Walser, 2015; Herring et al., 2010), choice of reference frame, treatment of tropospheric delays, integer phase ambiguity resolution, process‐noise settings, and the cutoff angle for satellites near the horizon (e.g., Herring et al., 2016; Martens et al., 2020; Springer et al., 2019). For instance, reference frames (i.e., the coordinate system used for measuring positions and velocities on Earth that are realized through a model of time‐dependent coordinates for a large number of reference stations across the globe) are continually improving as new versions become more precise and accurate through time. In addition, analysis centers use various software packages for transforming raw satellite data into estimates of site position. These include GAMIT/GLOBK (Herring et al., 2010); GipsyX (Bertiger et al., 2020), which supersedes GIPSY‐OASIS II (Zumberge et al., 1997), and Bernese (Dach & Walser, 2015). Each analysis center provides estimates of uncertainty with their data products; however, they usually do not account for the variance introduced by processing assumptions (Herring et al., 2016; Johnson et al., 2021).

Of relevance to this review, processing differences can influence assessments of water storage variations with estimates of seasonal vertical deformation varying between data products by 10%–40% or more (Martens et al., 2020). Care should be given when selecting data products from different analysis centers and when comparing loading studies from similar areas that may have processed GNSS data differently. At a minimum, we recommend that the processing software, reference frame, analysis center, and product version be clearly stated in all published analyses of GNSS data for hydrogeodetic applications (Table 1). Consistent and clear documentation among studies will facilitate reproducibility and comprehensive comparison.

As GNSS processing methods have advanced (e.g., Argus, 2012; Bertiger et al., 2020; Herring et al., 2016, 2010; Rülke et al., 2008; Steigenberger et al., 2006), the identification and removal of nonhydrologic loading signals have improved, GNSS antennas have proliferated on a global scale, and the utility and popularity of hydrogeodesy have increased significantly (Table 1). With that said, comparisons between reprocessed GNSS time series and other geodetic methods, such as very long baseline interferometry (VLBI; Tesmer et al., 2009) and GRACE (Horwath et al., 2010; Tesmer et al., 2011), have shown improvement over time, indicating that advances in processing have led to reductions in systematic errors and mismodeling and that GNSS now better reflects true vertical elastic Earth deformation.

Although not covered in detail here, GNSS position time series can also be used to study land subsidence associated with groundwater withdrawal and show promise to advance hydrogeologic characterization from the study of poroelastic properties during aquifer pumping tests (e.g., Moreau & Dauteuil, 2013 and references therein). However, we note that hydrologic loading analyses commonly omit GNSS sites overlying aquifers experiencing large groundwater withdrawal (e.g., Argus et al., 2014) because the poroelastic response from groundwater withdrawal and recharge can obscure Earth's elastic response to water loads. In some locations, poroelastic deformation is larger than the whole‐Earth elastic response to hydrologic loading from rainfall (Nespoli et al., 2021).

GNSS data are also used in studies of multipath reflectometry (i.e., the combined reception of direct and reflected, diffracted, or scattered GNSS signals) from the ground surface, vegetation, and snow. Multipath reflectometry can yield inferences of soil moisture content (Chew et al., 2014; Larson, 2019; Larson et al., 2008a, 2008b, 2010) and snow depth (Larson et al., 2009; Larson & Nievinski, 2013; Nievinski & Larson, 2014; Ozeki & Heki, 2012; Wang et al., 2020; S. Zhang et al., 2021). Several review papers have concentrated on these methods (Brocca et al., 2017; Dorigo et al., 2021; Jin & Komjathy, 2010; Larson, 2019; Mohanty et al., 2017; Peng et al., 2017). While estimates of soil moisture and snow loads offer clear synergy to the methods discussed here, multipath signals are considered an unwanted source of noise to elastic hydrologic loading applications (e.g., Elósegui et al., 1995; Nievinski & Larson, 2014); therefore, applications arising from the indirect multipath signals are not a focus of this review.

Water Storage Estimation From GNSS Observations of Surface Displacement

The following sections outline foundational studies in hydrogeodesy and discuss techniques that have led to better refinement of hydrologic loading and thus improved estimates of TWS, defined as the total storage of water on land. We also summarize the current understanding of the influence of nonhydrologic loads on Earth deformation and techniques that have been adopted to reduce error in position time series and water storage estimates.

Foundational Studies in GNSS‐Based Hydrogeodesy

Early hydrogeodetic studies established that vertical displacement recorded by GNSS is correlated with changes in terrestrial water (i.e., gains and losses and the movement of water through the hydrosphere), creating the foundation of GNSS‐based hydrogeodesy. Specifically, Van Dam et al. (2001) estimated the surface deformation expected to result from predicted hydrologic loads (from snowpack, soil, and groundwater) at GNSS stations worldwide with the vertical displacement predicted to be as large as 30 mm in some regions. They further showed close agreement between estimated and observed hydrologic load‐induced surface displacements, which indicated that annual global water storage variations contribute significantly to surface displacements recorded by GNSS. Heki (2001) directly compared snow depths to vertical surface displacements observed at GNSS stations across Japan and suggested that snow loads produce most of the annual variations in GNSS vertical time series there. Moreover, Heki (2001) also estimated snow depth across Japan using vertical displacements at GNSS stations, accurately reproducing the spatial variability of snow depths measured by in situ sensors and demonstrating the power of inverse modeling with GNSS, in which hydrologic loads can be estimated from the observed Earth deformation.

Within the last decade, a few key studies have demonstrated the utility of GNSS‐based hydrogeodesy in the context of water resources, climate, and the broader hydrologic sciences. In particular, Borsa et al. (2014) employed interannual trends in GNSS displacement time series to estimate the impact of drought on water resources in the western United States. They showed that the Earth deformation recorded by GNSS could be used to estimate changes in water storage at a sufficient spatial resolution to capture the different drought response of high‐precipitation mountains and arid deserts. Argus et al. (2017) combined hydrologic model storage estimates with GNSS‐inferred storage to infer how water loss was partitioned across hydrologic stores, finding that drought depleted groundwater storage across California, USA. They further suggested that increases in TWS during years of heavy precipitation result from spring recharge of subsurface water stores in the Sierra Nevada mountains, showing that GNSS‐based hydrogeodesy can be useful for studying hydrologic processes.

Isolating Hydrologic Loading

All sources of surface loading contribute to the elastic deformation of the solid earth, which Mangiarotti et al. (2001) modeled at the global scale by summing the individual contributions of surface loads such as snow, soil moisture, atmospheric pressure, and oceanic loading. Because surface loads add constructively in the vertical component, investigating the effects of one load, such as continental water, requires isolation of that load by removing the other components. Load isolation is commonly performed following the “peering” approach from Dong et al. (2002) where the contribution of well‐determined, modeled sources is subtracted to examine the observed residual variations. Ocean tidal loading is generally removed during GNSS processing according to IERS conventions (Petit & Luzum, 2010); however, isolating the deformation from hydrologic loading requires additional post‐processing strategies.

Here, we outline common methods for removing other sources of surface loading and eliminating systematic errors and station artifacts. While most approaches are somewhat standardized across studies, specific steps depend on the particular hydrologic application, and post‐processing should be conducted on a case‐by‐case basis. First, GNSS stations are generally filtered to exclude stations with large data gaps; however, the length of the allowed gaps varies depending on the hydrologic application and the density of GNSS stations. Generally, investigation of seasonal trends requires more complete time series and thus smaller and less frequent data gaps than investigation of multiyear trends (e.g., Borsa et al., 2014; Enzminger et al., 2019).

GNSS stations located over large aquifers experiencing active groundwater withdrawal are also commonly excluded from hydrologic loading studies. In the case of GNSS stations in the Central Valley of California, for example, the poroelastic response to ongoing groundwater pumping produces large subsidence trends (Amos et al., 2014) rather than the modest uplift expected from Earth's elastic response to the decreasing groundwater load. Conveniently, GNSS stations impacted by poroelastic effects can be identified by vertical time series that are out of phase with expected motion from seasonal wetting and drying cycles (Argus et al., 2017, 2014; Borsa et al., 2014; Johnson et al., 2017). It is also common to exclude GNSS stations in known areas of recent volcanic and seismic activity or to correct those time series using modeled deformation (Adusumilli et al., 2019; Argus et al., 2017; Borsa et al., 2014; Knappe et al., 2019) as those geologic processes can produce substantial Earth deformation.

Long‐term vertical displacement can occur from nonhydrologic sources of surface loading, like glacial isostatic rebound and tectonic uplift, and must be removed to isolate hydrologic loading. Nonhydrologic secular displacement is commonly removed by detrending GNSS time series (i.e., fitting and removing a long‐term best‐fit linear trend). It is generally best to detrend using the entire data record (often a decade or longer from continuous sites, especially those within the CONUS) to best capture long‐term secular trends. Detrending is recommended for studies of seasonal hydrologic loading; however, detrending may also remove multiyear hydrologic trends, like those associated with long‐term drought or climate change. In studies of long‐term changes in water storage, viscoelastic (e.g., glacial isostatic rebound) and tectonic signals can be removed via modeled displacement (e.g., Argus et al., 2017; Borsa et al., 2014; Wang et al., 2017).

Reducing Noise and Enhancing Hydrologic Signals

Position time series from GNSS are also commonly filtered or smoothed to mitigate short‐period noise (i.e., scatter in the time series) arising from multipath reflections, orbital and timing errors, and mismodeled tropospheric and ionospheric delays. The smoothing timescale and technique vary across studies. Several studies have used rolling or moving averages to minimize daily noise (e.g., Adusumilli et al., 2019; Argus et al., 2017; Fu et al., 2013, 2015, 2012; Fu & Freymueller, 2012; Knappe et al., 2019; Nahmani et al., 2012; Steckler et al., 2010; Wang et al., 2017). Other studies have used low‐pass filters to eliminate intramonthly variability (e.g., Borsa et al., 2014; Enzminger et al., 2019; Ouellette et al., 2013). Prior to smoothing, it is common to identify and remove outliers, which are displacements that fall well outside the station or network average behavior (e.g., Borsa et al., 2014; Enzminger et al., 2019; Fok & Liu, 2019; Springer et al., 2019). This step helps mitigate anomalous position estimates linked to data or processing station artifacts or local motion unrelated to hydrologic loading.

Discrete offsets in the GNSS time series associated with earthquakes and equipment changes should also be removed (e.g., Herring et al., 2016). Analysis centers commonly maintain a list of the dates of known equipment changes and events, whose associated offsets can be removed from the time series using simple Heaviside step functions. In addition, errors in position estimates can result from noise created by snow and ice accumulation in and on GNSS antennas or radomes. Post‐processing is required to isolate and remove these snow and ice‐related artifacts (Koulali & Clarke, 2020; Larson, 2013). However, we note that there is no standard method for identifying and removing time‐series outliers and offsets. Post‐processing of the time series should be handled with care in the context of specific stations and scientific applications. In order to make studies comparable and reproducible, a description of offset and outlier handling should be included in publications.

Atmospheric pressure loading from short‐lived high‐ and low‐pressure weather systems can produce several millimeters of vertical GNSS displacements (Argus et al., 2017; Dill & Dobslaw, 2013; Martens et al., 2020; Petrov & Boy, 2004; Steigenberger et al., 2009; Tregoning & Herring, 2006; Williams & Penna, 2011). These displacements should be removed from GNSS time series, particularly for studies of short‐period hydrologic loading on the scale of days to weeks. Surface displacements induced by atmospheric loading are computed by convolving atmospheric surface pressures, such as from the European Centre for Medium‐Range Weather Forecasts (ECMWF) (Dee et al., 2011), with load Green's functions (Farrell, 1972; Longman, 1963). Global grids of surface displacement induced by atmospheric pressure loading can also be obtained from modeling centers, such as the Earth System Modeling Group of GeoForschungsZentrum Potsdam (ESMGFZ) (Dill & Dobslaw, 2013) and the School and Observatory of Earth Sciences (EOST) Loading Service (http://loading.u-strasbg.fr/). Martens et al. (2020) showed that GNSS processing methods can significantly affect the retention of atmospheric‐loading signals in GNSS time series (cf. Kouba, 2009; Steigenberger et al., 2009; Tregoning & Herring, 2006; Tregoning & Watson, 2009, 2011). Different GNSS data products and processing methods can also affect the estimated seasonal amplitudes of surface deformation and associated hydrologic storage (e.g., Martens et al., 2020) due to different assumptions about how to align GNSS stations to the terrestrial reference frame (Herring et al., 2016).

Some recent studies of TWS also model and subtract displacements due to surface loading from large reservoirs and lakes (Argus et al., 2017, 2020; Enzminger et al., 2019) in areas where lake volume data are available. These concentrated human‐controlled loads can induce large surface deformation at nearby GNSS stations, which can bias TWS estimates inverted from GNSS displacements. Similar to calculations of atmospheric loading, predictions of reservoir loading displacements typically require convolving changes in reservoir volume with load Green's functions. Convolution integrals for elastic surface loading can be computed using open‐source software packages, including LoadDef (Martens et al., 2019), SPOTL (Agnew, 2012), and REAR (Melini et al., 2015). After the removal of all of the previously described relatively well‐known loading signals, the remaining time series should isolate and better reflect the target hydrologic loading signals with errors and uncertainties minimized to the best extent possible.

Common mode error (CME), which refers to network‐coherent noise associated primarily with orbital, reference‐frame, and atmospheric‐modeling errors, is one of the largest sources of error in GNSS time series. Initially, CME was defined as a time series of daily rigid body translations (Wdowinski et al., 1997) but more recently is described as any temporally incoherent and spatially correlated error (Dong et al., 2006). CME is often considered to be spatially uniform across GNSS networks with apertures up to about 500 km, beyond which the assumption of spatial uniformity begins to break down (Dong et al., 2006). For networks that exceed that range, CME can be estimated for smaller subsets of the overall study network, or the amplitude of the CME can be allowed to vary spatially (e.g., using principal component analysis). Draconitic error, which is related to the small difference between the orbital alignments of GNSS satellites and Earth with the sun, is also thought to contribute to CME (Fu et al., 2015). Despite typical amplitudes that are significantly smaller than those from loading, the draconitic error should be considered in seasonal studies due to its nearly annual period of 351 days.

The CME can be calculated using multiple spatial filtering methods, including stacking (Wdowinski et al., 1997), principal component analysis (PCA) (Dong et al., 2006), independent component analysis (ICA) (Gualandi et al., 2016; Liu et al., 2015), correlation‐weighted spatial filtering (CWSF) (Tian & Shen, 2016), and common‐mode component (CMC) Imaging (Kreemer & Blewitt, 2021). Milliner et al. (2018) showed that stacking, PCA, and ICA methods can produce consistent estimates of CME, and elimination of CME with ICA allowed them to track daily storm water movement across South Texas and Louisiana. Subtraction of the CME from filtered positional time series has been shown to enhance signal‐to‐noise ratios and sharpen the GNSS signals from solid Earth processes (Fok & Liu, 2019; Jiang et al., 2018; Williams, 2004). For a more in depth comparison of spatial filtering methods and CME, the reader is referred to He et al. (2017). We also note that PCA and ICA methods have been shown not only to reveal CME, but also to isolate spatiotemporal surface loading components in recent hydrologic loading studies (Hsu et al., 2020; Jiang et al., 2021; Larochelle et al., 2018; Milliner et al., 2018; Nespoli et al., 2021). Following isolation of hydrologic surface loading through removal of noise and unwanted signals, fitting a linear model to the residual time series with a constant offset, linear trend, and annual and semiannual sinusoid may be used to estimate the seasonal surface displacement due to hydrologic loading (e.g., Argus et al., 2017).

Comparing Vertical and Horizontal Displacements

Hydrologic loading induces crustal deformation most noticeably in the vertical direction on the order of several millimeters to several centimeters, depending on the total mass of hydrologic loads and their proximity to GNSS stations. For instance, GNSS stations record seasonal vertical deformation up to approximately 75 mm (from peak to trough) along the Amazon River in Brazil (Bevis et al., 2005) and up to 60 mm in Bangladesh where the Ganges, Brahmaputra, and Meghna Rivers converge (Steckler et al., 2010). Hydrologic loading also induces smaller horizontal surface displacements, which can be used together with vertical displacements to better constrain the location of changing hydrologic loads (Wahr et al., 2013, Figure 1).

Horizontal displacements are generally less than half, and often much smaller than half, of the magnitude of vertical displacements, especially in cases where the load is spread over a large region (Pinel et al., 2007; Wahr et al., 2013). For example, time series of vertical and horizontal displacement from five GNSS stations across different climatic regimes (Figure 4) show that although vertical displacement varies significantly between stations, it exceeds horizontal deformation by a factor of three or more in all cases. Tregoning et al. (2009) explain that horizontal deformation is smaller because horizontal movement may experience destructive interference from multiple loads surrounding a station as the station can be simultaneously pulled in different directions.

Figure 4

Time series of vertical and horizontal displacement (where positive displacement corresponds to up, north, and east movement, respectively) across a decade at five globally distributed Global Navigation Satellite System (GNSS) stations exhibit a wide range of Earth deformation. The station name and location are shown above each panel. Seasonal trends and interannual variability in surface deformation from hydrologic loading are apparent at each site. Offsets can be seen in the north component at NAUS at the time of logged antenna substitutions and at ELDC from an earthquake in 2018. Peak‐to‐peak vertical displacement amplitude is largest at NAUS on the Amazon River in Brazil (60 mm) and smallest at WARA in the Australian desert (20 mm). Note that the vertical scale is three times larger than the horizontal scale in all plots. GNSS positions in the IGS14 reference frame were retrieved from and pre‐processed by the Nevada Geodetic Laboratory using GIPSY‐X software. Atmospheric and oceanic loading were removed in post‐processing.

Horizontal movement from hydrologic loading is generally an order of magnitude less than long‐term tectonic movement; therefore, horizontal displacements recorded by GNSS are most commonly used to estimate secular crustal velocities to constrain tectonic‐plate motion as well as coseismic and post‐seismic deformation due to earthquakes (Ray et al., 2021). Furthermore, Borsa et al. (2014) note the use of vertical‐to‐horizontal displacement ratios to identify GNSS stations experiencing tectonic or volcanic motion where horizontal motion exceeds vertical motion.

Most of the studies that have examined hydrologic loading solely using horizontal deformation were located in areas known to have large hydrologic loads (and thus large GNSS‐recorded vertical deformation), namely the Amazon Basin and Southeast Asia (Bettinelli et al., 2008; Chanard et al., 2014; Fu et al., 2013; Ray et al., 2021). However, a few studies have demonstrated the utility of analyzing both horizontal and vertical surface deformation. For example, Silverii et al. (2016) were able to discern between the poroelastic response of karst aquifers and elastic deformation due to hydrologic loading in Apennines, Italy, using the comparison of horizontal and vertical deformation with spring discharge and GRACE‐derived TWS estimates, respectively. Furthermore, Wahr et al. (2013) found the relationship between vertical and horizontal elastic response to surface loading to be relatively insensitive to the assumed structural properties of the solid earth, which is important when inverting surface displacement to derive TWS. Milliner et al. (2018) also showed that including vertical and horizontal deformation can improve the fidelity of daily inversions of water mass during flooding from a hurricane.

Inverting Surface Displacements to Water Storage Estimates

Earth surface displacements from tens to thousands of GNSS stations are commonly inverted to derive a best‐fit model for surface load distribution through time, resulting from changes in TWS over the inversion domain (e.g., Argus et al., 2017, 2014; Borsa et al., 2014; Enzminger et al., 2018; Fu et al., 2015; Ouellette et al., 2013). Load‐deformation modeling packages, such as LoadDef (Martens et al., 2019) and SPOTL (Agnew, 2012), which forward‐model Earth deformation corresponding to a given surface load on a sphere are typically used to generate a matrix of station displacement responses to individual loads over a gridded domain (i.e., a matrix of Green's functions), which is itself used to invert for a surface‐load distribution that best explains observed GNSS displacements (e.g., Argus et al., 2014; Borsa et al., 2014). Because the number of observations is always smaller than the number of grid cells in the desired surface loading model, the inversion is typically stabilized with a model size and/or model smoothness constraint (e.g., Argus et al., 2017; Borsa et al., 2014; Enzminger et al., 2018).

For surface loads greater than approximately tens of km in area, the choice of Earth rheological model (e.g., PREM, Gutenberg‐Bullen A, and ak135) produces relatively minor differences in forward modeled surface displacements (Argus et al., 2017; Enzminger et al., 2018; Na & Baek, 2011). However, Dill et al. (2015) showed that heterogeneities in solid Earth structure can significantly influence GNSS measurements of surface loading under strong localized hydrologic signals, like floods and heavy rain.

Synthetic tests that compare known input to inverted output, such as checkerboard tests with juxtaposed grids of positive and negative displacements (e.g., Borsa et al., 2014; Jiang et al., 2021) or synthetic tests that incrementally extend the inversion area like that of Fu et al. (2015), are needed to determine the optimal inversion domain and parameterization. Estimates of spatially distributed water mass or equivalent water depth, over the inversion domain and through time, are used to quantify changes in TWS. The GNSS‐derived estimate of mass gain or loss further enhances the interpretation of hydrologic loads in the context of water resources and flood volumes.

The solutions for TWS are only as good as the GNSS data used to constrain them; thus, it is important to understand the sensitivity of GNSS‐inferred surface loading and validate the estimates of TWS. Hydrogeodetic studies sometimes combine different geodetic techniques, such as GNSS, GRACE, InSAR, and VLBI, and land surface hydrologic models for comparison, validation, and correction of each method. In particular, hydrogeodetic studies using GNSS are commonly paired with GRACE. Gravity Recovery and Climate Experiment (Tapley, 2004) and the GRACE follow on mission (GRACE‐FO; Kornfeld et al., 2019) consist of two co‐orbital satellites launched in March 2002 and May 2018, respectively. As the distance between the two satellites changes in space and time, GRACE estimates the spatial and temporal variations of Earth's gravitational field. These variations are converted to changes in water mass (expressed as equivalent water height) over a spatial resolution of about 300 km, sampled monthly and extending almost two decades (from 2002 to present) (Frappart & Ramillien, 2018). Both GNSS and GRACE infer the integrated total water storage within both surface and subsurface reservoirs. These two estimates of water load, although differing in spatial and temporal resolution, tend to be highly correlated (Adusumilli et al., 2019; Argus et al., 2014, 2017; Borsa et al., 2014; Chew & Small, 2014; Davis, 2004; Fu et al., 2012, 2013; Fu & Freymueller, 2012; King et al., 2006; Nahmani et al., 2012; Ouellette et al., 2013; Ray et al., 2021; Steckler et al., 2010; Tesmer et al., 2011; Tregoning et al., 2009).

Pairing GRACE TWS and GNSS observations has been shown to provide effective additional constraints on TWS that reduce differences and move toward better agreement between the two geodetic methods. For example, in a global analysis of TWS, Tesmer et al. (2011) showed that subtracting GRACE observations from GNSS observations (referred to as GPS‐GRACE) led to a considerable reduction in GNSS time series variance. In addition, Fok and Liu (2019) significantly reduced the root‐mean‐square (RMS) error of their solutions with respect to independently modeled TWS by using GRACE‐derived vertical surface displacements as virtual GNSS stations in southwest China where greater GNSS station density was needed to constrain inversion solutions. Similarly, Adusumilli et al. (2019) noted that jointly inverting the high temporal and spatial resolutions of GNSS displacements with the lower resolution but global coverage of GRACE TWS ensures greater consistency of TWS estimates across a range of temporal scales and enhances the quality of inversions in areas with low GNSS density.

Applications

Applying GNSS to hydrologic studies provides exciting opportunities for new insight into the redistribution of surface and near‐surface water, including the hydrologic processes controlling water storage and movement. Furthermore, the emerging discipline of hydrogeodesy continues to develop from the ongoing expansion and densification of GNSS networks and improvements in the accuracy of GNSS positioning due to technical advances in processing and data analysis. In the following section, we highlight key findings and advancements from hydrogeodetic studies that leverage GNSS data sets, focusing on their importance to hydrologic science and water resource research.

Hydrologic Loading and Storage Estimation Across Temporal Scales

Hydrologic loading is directly related to changes in water storage, providing a critical constraint for hydrologic budgets. Hydrologic loading has been primarily studied at annual scales associated with seasonal fluctuations in water mass (e.g., Argus et al., 2014; Bevis et al., 2005; Blewitt et al., 2001; Enzminger et al., 2019; Fu et al., 2015; Knappe et al., 2019). Recent studies take advantage of long GNSS position time series to study the interannual timescales associated with longer duration climatic drivers, such as drought and sustained years of high precipitation (e.g., Argus et al., 2017; Borsa et al., 2014; Chew & Small, 2014; V. G. Ferreira et al., 2018; Jin & Zhang, 2016). Additionally, many studies investigate both seasonal and interannual scale changes (e.g., Adusumilli et al., 2019; Argus et al., 2020; Elósegui, 2003; Fu et al., 2015; Hsu et al., 2020; Jiang et al., 2017; Wahr et al., 2013; Wang et al., 2017; Yin et al., 2020, 2021). However, only recently has attention turned to hydrologic loading at shorter, more dynamic temporal scales associated with individual storm events and subseasonal deformation (Adusumilli et al., 2019; Ferreira et al., 2019; Han & Razeghi, 2017; Jiang et al., 2021; Knowles et al., 2020; Milliner et al., 2018; Springer et al., 2019; Yin et al., 2021; Zhan et al., 2021). In this section, we review key applications of and advancements from GNSS‐inferred hydrologic loading and associated TWS change estimates across temporal scales ranging from several years to days.

Seasonal Loading and Storage Estimation in Snow‐ and Rain‐Dominated Regions

Seasonality in precipitation and evapotranspiration rates can cause large annual fluctuations in the regional TWS. These seasonal changes are manifested as cyclic patterns with an annual period in vertical and horizontal GNSS time series (Figure 4). Seasonality in GNSS surface displacement and associated estimates of continental water storage have been shown to result from global‐scale redistribution of continental water mass (Blewitt et al., 2001) and have been validated through correlation with hydrologic land surface models (Dong et al., 2002; van Dam et al., 2001; Wu et al., 2003) and GRACE (Dill & Dobslaw, 2013; Döll et al., 2014; Tesmer et al., 2011; Tregoning et al., 2009). Early global‐scale hydrogeodetic studies using GNSS also observed regions of more pronounced loading corresponding to seasonal wetting (Blewitt et al., 2001; Dong et al., 2002; van Dam et al., 2001; Wu et al., 2003), which paved the way for studies of more focused seasonal loading at smaller spatial scales.

Seasonal trends in displacement and TWS are most pronounced in locations with large spatial and temporal gradients in precipitation. For instance, hydrologic loading has been shown to anticorrelate with changing climatic drivers, like snow (using snow depths or measurements of snow water equivalent (SWE)) (Fu et al., 2015; Heki, 2001; Knappe et al., 2019; Ouellette et al., 2013; Yin et al., 2020) and rainfall (Argus et al., 2017; Birhanu & Bendick, 2015; Borsa et al., 2014; Fu et al., 2013; Nespoli et al., 2021; Zhan et al., 2017) across terrain with a range of aridity and topography. Using a network of 16 GNSS stations across Ethiopia and Eritrea, Birhanu and Bendick (2015) discerned different seasonal loading patterns from monsoon rain as it varied spatially across a region with variable aridity. Additionally, after examining the relationship between SWE and GNSS displacement, Ouellette et al. (2013) used time series of surface displacement to predict SWE values with reasonable accuracy at GNSS stations across the western United States. While a more thorough error analysis and investigation of scaling are needed, this hydrogeodetic application serves as an example of how GNSS can enhance valuable hydrologic information that is an important predictor of spring runoff and groundwater recharge.

Hydrogeodetic studies also commonly compare measures of storage from land surface hydrologic models to hydrologic loading in order to explore seasonal trends in storage (e.g., Argus et al., 2014; Enzminger et al., 2019; Fu et al., 2015, 2013; Han, 2017). For instance, Argus et al. (2014) compared seasonal GNSS‐inferred estimates of TWS to estimates from three different hydrologic land surface models: NLDAS‐Noah (Mitchell, 2004), which estimates total water storage as the sum of soil moisture and SWE over the continental United States; GLDAS‐Noah (Rodell et al., 2004), which is essentially a coarser resolution, global version of NLDAS‐Noah and is used in the calculation of GRACE gain factors; and a composite model that includes surface water loading in artificial reservoirs across California in addition to soil moisture and SWE. While they were unable to discern relative contributions of each load to the GNSS‐recorded displacement, Argus et al. (2014) found the best agreement between seasonal GNSS‐inferred TWS change and the composite model and concluded that seasonal TWS in California was primarily controlled by snow, soil moisture, and surface water reservoirs.

In addition, recent studies have begun to examine how TWS is partitioned into separate components of the hydrologic system through temporal comparisons of storage from land surface models, precipitation, and GNSS‐inferred hydrologic loading (e.g., Argus et al., 2017; Enzminger et al., 2019, 2018; Hsu et al., 2020; Yin et al., 2020). For example, Enzminger et al. (2019) analyzed the partitioning of TWS changes in the Sierra Nevada mountains, by comparing the timing of storage changes inferred from GNSS, remote sensing measurements of SWE and rainfall, and two land surface models. They found that annual measures of accumulated snow account for only 30%–50% of wet period storage, whereas the remaining annual TWS is stored in the subsurface after peak SWE.

A few hydrogeodetic studies have used a simple water budget calculation to investigate the timing of changing TWS across seasons (Han, 2017; Jiang et al., 2017; Yin et al., 2020). As is common in hydrologic studies, the analysis of the time lags between inputs and outputs of a hydrologic system can inform the hydrologic processes controlling streamflow generation and storage dynamics. GNSS‐based hydrogeodesy offers a new perspective on this timing analysis by providing improved constraints on changing water storage; however, this also highlights limitations of hydrogeodesy. For instance, explanation of loading and storage trends and investigating partitioning within the hydrologic system often relies on modeled and scaled estimates of storage or precipitation with their own large uncertainties. Furthermore, subsurface storage is difficult to compare across studies because its definition is not always consistent due to the lack of vertical discretization between subsurface reservoirs from GNSS, and because the partitioning of storage often relies on estimates of modeled soil moisture with large uncertainties due to soil depth, type, and soil moisture spatial heterogeneity (Walker & Houser, 2004; Wood et al., 2011).

Subseasonal and Storm Event‐Scale Storage

Extreme weather events are projected to contribute a growing percentage of precipitation within the 21st century (Wuebbles et al., 2014); therefore, short‐lived storm‐scale applications of GNSS will be more important than ever for understanding water resource changes and for potentially integrating GNSS observations into hydrologic forecasting. For instance, Springer et al. (2019) used a new land surface hydrologic model with daily time steps to show that even in Europe where extreme precipitation events and concurrent daily hydrologic loading are rare, GNSS‐recorded displacement agrees better with loading from models at daily time steps compared to the traditional monthly time steps. This suggests that sub‐monthly hydrologic loading is an important target for GNSS studies of surface displacement across the globe, especially in areas where extreme precipitation events are common and/or are becoming more frequent. Furthermore, Adusumilli et al. (2019) calculated TWS anomalies from weekly inversions of GNSS vertical deformation across the CONUS and found that atmospheric river events, an important contributor to flooding in the western United States, comprised about 40% of California's total TWS during the study period from 2007 to 2017. In addition, their orographic analysis of storage along the California coast and inland mountain ranges showed that the impact of atmospheric river events is larger and longer lasting in the higher‐elevation inland ranges where precipitation falls as both rain and snow.

A few recent hydrogeodetic studies have inverted GNSS displacements to estimate high‐frequency TWS changes at timescales as short as 1 day (Han & Razeghi, 2017; Jiang et al., 2021; Milliner et al., 2018; Zhan et al., 2021). For instance, Milliner et al. (2018) quantified flood volumes as daily changes in water storage during and following Hurricane Harvey and tracked the movement of stormwaters using 198 GNSS stations in the Southern Texas and Louisiana region, USA. They showed that incorporating horizontal GNSS displacements with more commonly used vertical displacements improved high‐frequency estimates of TWS and allowed for better tracking of the movement of focused hydrologic loads that could aid water managers during the event of catastrophic floods like those from hurricanes. Similarly, typhoons in the Northwest Pacific cause extremely heavy rainfall over just a few days that leads to ground deformation, which can be recorded by GNSS (e.g., Heki, 2020; Zhan et al., 2021). For instance, Zhan et al. (2021) tracked the progression of flooding following landfall of Typhoon Hagibis in Japan. They were able to show, with comparisons of precipitable water vapor (PWV) and ground‐truthed measurements of precipitation, that GNSS accurately tracked the path of the storm. Furthermore, they were the first to study both water vapor and GNSS‐recorded ground deformation during a typhoon and were able to discern the large and fast progression of subsidence along the path of the storm (approximately 1–2 cm deformation on the day of typhoon landfall and approximately 71 Gt of surface water loading) that would have been missed with traditional longer‐scale inversions. However, their study also highlighted a topographic limitation of GNSS‐inferred hydrologic loading, related to the fact that GNSS stations in mountainous regions tend to be located at lower elevations with easier access. In these regions with significant topographic relief, GNSS methods can overestimate surface loads, as seen by Zhan et al., 2021, due to the oversampling of valleys and basins where water accumulates. Finally, Jiang et al. (2021) used daily inversions of two winter storm events in Yunnan, China, to show that while the noise in daily inversions is high, GNSS methods are still able to resolve loading TWS changes from weather events less extreme than hurricanes and typhoons. They showed that daily changes in TWS (rapidly increasing storage following winter rain that persisted for approximately 10 days) could be captured with a somewhat sparse (average interstation distance of approximately 157 km) network of GNSS stations. Therefore, with extreme weather events forecasted to happen more frequently, GNSS shows good promise for capturing the resulting TWS anomalies in near‐real time at sub‐monthly intervals, which can be hugely beneficial to water managers.

While not a focus of this review, GNSS is commonly employed in meteorological studies of integrated PWV, especially at the temporal scale of severe storm events, and is assimilated into numerical weather prediction models (e.g., Moore et al., 2015; Shoji et al., 2015 ). We posit that the integration of “GNSS Meteorology” into improved operational weather forecasting as explained by Bonafoni et al. (2019) can serve as an example of how GNSS studies of hydrologic loading may similarly benefit operational hydrology and management of water resources.

Interannual Changes in Terrestrial Water Storage

Surface water volumes in artificial reservoirs and natural lakes vary across multiple timescales due to changing temperature, precipitation, and consumptive use demands. Several studies have shown that GNSS records these changes in local hydrologic loads as surface deformation correlated with the magnitude of changes in reservoir storage (Argus et al., 2020; Bevis et al., 2004; Dill & Dobslaw, 2013; Elósegui, 2003; Wahr et al., 2013; Xue et al., 2021). For example, Elósegui (2003) showed that surface deformation from seasonal trends and multiyear decreases in water levels in the Great Salt Lake, Utah, were captured as relatively small (maximum 6 mm) vertical displacements of two nearby GNSS stations. Dill and Dobslaw (2013) also observed seasonal changes in surface deformation resulting from lake level changes in Manton reservoir in Australia. Furthermore, Wahr et al. (2013) showed that multiyear trends in lake water levels in Lake Shasta, California, resulted in vertical and horizontal elastic Earth deformation recorded by nearby GNSS stations with decreasing sensitivity to the lake load as the distance from the lake increased.

More recently, Argus et al. (2020) and Xue et al. (2021) examined loading from increases in lake waters in the Great Lakes, United States, that resulted in up to 20 mm of land subsidence. Argus et al. (2020) were able to distinguish between causes of long‐term land subsidence (lake loading, soil and groundwater loading, and viscous collapse of the Laurentide ice sheet) in the surrounding Great Lakes region. Xue et al. (2021) further expanded on the seasonal component of hydrologic loading in the Great Lakes region to show that lake loading accounted for 10%–25% of the total loading signal immediately surrounding the Great Lakes, while soil moisture constituted the majority of hydrologic loading in the region. Moreover, the importance of lake storage to the GNSS loading analysis was demonstrated by Argus et al. (2014) through a comparison of a composite hydrologic model, reservoir volumes, and GNSS‐inferred seasonal storage. The authors show that the effect of surface water reservoirs needs to be removed from GNSS displacements to avoid biasing GNSS estimates of TWS changes, especially in places where reservoirs are relied on heavily for consumptive use and irrigation (see Section 3.2.1 for more details).

Several recent GNSS loading studies have demonstrated their utility for understanding how long‐term trends in decreased (increased) precipitation and the associated crustal uplift (subsidence) offer a new perspective on multiyear changes in TWS. For instance, Borsa et al. (2014) showed that GNSS stations across the western United States captured multiyear trends in Earth deformation during wetter than average and drier than average years, consistent with spatial and temporal trends in precipitation and streamflow. Furthermore, during the severe drought in 2014, the authors found that TWS was most impacted in the mountains of California and along the Pacific coast with losses as great as the equivalent water depth of 0.5 m. In a study of periods of drought and subsequent years of heavy precipitation, Argus et al. (2017) found that the water loss from mountain ranges in California during the severe drought from 2011 to 2015 was 1.5 times greater than the water gained during the period of increased precipitation from 2010 to 2011, showing that long‐term studies using GNSS can provide important insights for water managers with greater spatial and temporal resolutions than are available from GRACE alone. Moreover, interannual TWS gains and losses from GNSS and GRACE indicate that hydrologic land surface models consistently underestimate TWS, again underscoring the value of hydrogeodetic methods like GNSS.

GNSS loading studies have also recorded the impact of changing El Niño Southern Oscillation (ENSO) patterns on TWS. ENSO is a 2–7 years pattern in warming and cooling in the tropical Pacific Ocean that is the best‐known indicator of global interannual climate variability, a strong indicator of temperature and precipitation patterns, and typically associated with extreme La Niña and El Niño events (Adusumilli et al., 2019; Ropelewski & Halpert, 1987). Using 120 GNSS stations across Australia, Han (2017) captured interannual trends in surface deformation associated with record‐breaking rainfall during the 2010–2011 La Niña event and subsequent drying. The authors showed that the La Niña event increased TWS by 9–10 cm of equivalent water depth across the continent, which took 3–4 years to return to its long‐term average. Adusumilli et al. (2019) showed that GNSS‐derived TWS in the contiguous United States is strongly correlated with the Oceanic Niño Index (ONI), which tracks sea surface temperature anomalies indicative of El Niño and La Niña conditions. In addition, Adusumilli et al. (2019) examined the ratio of interannual to seasonal TWS anomalies across the CONUS and found that high ratios indicate areas where long‐term changes in TWS largely influence water availability. Such areas, like California and the Texas Gulf, were coincident with regions where droughts are common and municipal water storage and delivery are heavily managed, again showcasing the utility of long‐term hydrogeodetic studies using GNSS to provide pertinent information for water resource management.

Although interannual changes in TWS associated with drought and heavier than average precipitation can be very useful for understanding the relationship between climate trends and extreme events with long‐lasting impacts, hydrogeodetic estimates of storage changes are accompanied by considerable uncertainty. For instance, Argus et al. (2017) estimated uncertainty in change in equivalent water thickness to be 0.08–0.15 m and found that uncertainty increased with the time period over which water change is estimated. Uncertainties in water storage calculations can vary temporally (seasonal vs. multiyear signals) and spatially (across GNSS networks, depending on station density); however, synthetic tests (e.g., Borsa et al., 2014; Fu et al., 2015; Jiang et al., 2021) now accompany nearly all GNSS studies of TWS in order to quantify this uncertainty. Moreover, understanding and reducing such uncertainty are currently a focus of geodetic research.

Drought Indices Using GNSS

Droughts are becoming more common and widespread across the globe, threatening water resources and impacting the hydrologic cycle (Asadieh & Krakauer, 2017; Chikamoto et al., 2017; Dong et al., 2021; He et al., 2021; Huntington, 2006; Swain et al., 2018; Yoon et al., 2015). Droughts are typically characterized in three main ways (i.e., hydrological, meteorological, or agricultural) and numerous drought indices exist to classify and compare drought severity and duration (Heim, 2002; Zargar et al., 2011). As previously discussed in Section 4.1.3, TWS estimates from GNSS are able to quantify mass loss during drought in near real time, providing new insight about drought's impacts and enhancing our ability to prepare for those impacts. Moreover, a few recent hydrologic loading studies have used GNSS position time series and associated TWS estimates to create new drought indices, showing that GNSS can be used to monitor drought intensity (Chew & Small, 2014; Ferreira et al., 2018; Jiang et al., 2021).

Evaluating meteorological droughts (classified by decreased precipitation) in the context of hydrological droughts (characterized by GNSS‐derived water storage loss) can be useful in predicting drought impact and recovery. For instance, the storage deficit accumulated during drought can be estimated, which can help forecast drought in the upcoming dry season due to cumulative loss of TWS across the water year or incomplete recovery from previous hydrological droughts. Furthermore, comparison of TWS anomalies after removal of surface water reservoirs (e.g., Argus et al., 2017) may help to quantify the direct impact of drought on subsurface water resources. For example, Chew and Small (2014) calculated a drought index for each GNSS station (referred to as DIGPS) using GNSS time series of vertical surface displacement as indicators of changes in TWS. Vertical position time series were converted to anomalies by subtracting the mean, and the seasonal cycle was removed by modeling the seasonal component with least squares to allow for comparison of drought intensity across hydrologically variable times of year. DIGPS indicated the likelihood of each seasonally adjusted anomaly, given the variability of the study's full time series at that GNSS station, where lower DIGPS percentiles indicate drier periods and higher percentiles indicate wetter periods within the time period of the study. It is worth noting that the value of DIGPS depends on the drought episodes within the particular time period studied and the intensity of a particular drought may be different if the study period isolated that one drought period or included multiple periods of drought. Comparison of DIGPS with two independent indices of meteorological drought across the High Plains, USA, suggested that drought characterized by depleted TWS persisted longer than indices of meteorological drought, which improved rapidly following increased precipitation because, unlike meteorological‐based indices, DIGPS also accounted for groundwater and soil moisture lost during drought. Further comparison of drought onset and termination using drought indices that encompass different components of the hydrologic system may help to understand the propagation of drought and subsequent storage timing within the hydrologic system.

Jin and Zhang (2016) compared TWS estimates from GNSS and GRACE, TWS calculated from the water balance equation, and the drought index from the US Drought Monitor across the American Southwest. They found the best agreement between GNSS‐inferred TWS and the US Drought Monitor index, while TWS from the water balance equation and GRACE lagged 1–2 months behind. Jiang et al. (2021) showed that similar to GRACE, GNSS accurately observes hydrological droughts, and they used equivalent water height calculated from GNSS displacements to classify drought severity, onset, duration, and termination in Yunnan, China, with a better spatial resolution than GRACE could provide. Furthermore, Ferreira et al. (2018) compared drought indices derived from GRACE (DITWS) and GNSS‐recorded vertical crustal deformations (DIVCD) to classify drought across Brazil using relatively short (7 years) time series. While GRACE TWS and GNSS displacement agreed well at greater than 90% of the stations, DITWS and DIVCD were correlated at only 82% of the stations. However, their work showed that vertical deformations themselves, which are widely and rapidly available across the globe, can be as useful for monitoring drought severity as inverted estimates of TWS from GNSS, which require greater station density and are subject to increased uncertainty.

Integrating GNSS‐Inferred Loading With Traditional Hydrologic Observations

Incorporating traditional hydrologic observations into the analysis of GNSS‐inferred hydrologic loading can advance our understanding of the dynamics of the hydrologic system. For example, despite the complexity of separating the poroelastic response of aquifers to subsurface water storage changes from the accompanying solid Earth response from loading (both of which are recorded by GNSS), a few recent studies have explored the relationship between the changing water table depth and GNSS time series of surface displacement. Milliner et al. (2018) found strong agreement between GNSS‐inferred hydrologic loading and groundwater depth measurements in the Southern Texas and Louisiana area, USA, suggesting that following major flood events like that of Hurricane Harvey, GNSS time series dominantly reflect hydrologic loading rather than poroelastic aquifer properties and can be used to partition surface and groundwater impacts when paired with traditional hydrologic measurements, like potentiometric maps. Moreover, using PCA separation of groundwater level, rainfall, and GNSS‐recorded vertical displacement time series, Nespoli et al. (2021) detected seasonal and multiannual signals in piezometric data and surface displacements in response to rainfall. They were able to spatially discern between seasonal elastic and poroelastic displacement in GNSS time series, wherein elastic deformation dominated in the mountainous region of the Apennines and poroelastic response dominated in the sedimentary basin of the Po Plain, Italy. Similarly, Silverii et al. (2016, 2020) used streamflow and precipitation data to explain seasonal horizontal GNSS displacements in terms of the poroelastic response of mountain aquifers to seasonal groundwater infiltration in the Apennines and Sierra Nevadas.

Pairing traditional hydrologic measurements with hydrogeodetic studies using GNSS can also inform the partitioning of water within the hydrosphere and potentially give some insight into residence times and water transport. For instance, Chew and Small (2014) paired their GNSS‐based drought index with in situ soil moisture measurements and groundwater depths to demonstrate that TWS anomalies during drought resulted, at least in part, from depleted soil moisture and groundwater storage. The temporal changes of these traditional hydrologic measurements helped explain the timing of GNSS‐based drought intensity observed in the High Plains, USA. In a detailed hydrogeodetic study of partitioning, Hsu et al. (2020) combined GNSS‐inferred TWS, TWS from GRACE, precipitation, GLDAS soil moisture, and groundwater levels to show that groundwater is a substantial component (up to 66% of TWS) of seasonal hydrologic loading in Taiwan and suggested that land surface models and GRACE underestimate TWS. In addition, further spatiotemporal analysis using ICA revealed different response times between storage estimate methods and precipitation and the authors suggested that the observed phase lags resulted from spatially variable infiltration rates. It has also been suggested that the detailed analysis of water storage changes inferred from dense GNSS networks (tens of km station separation) could provide estimates of water transit times at the watershed scale (Knappe et al., 2019); however, no such study has been completed yet.

Advancements in Spatial Scales

The number of continuous GNSS stations worldwide has grown rapidly over the past two decades, creating denser GNSS networks over wider areas. For instance, Japan's network of approximately 1300 GNSS stations, known as GEONET, has a typical interstation separation of only 15–30 km (Zhan et al., 2021). As GNSS station density increases, the spatial resolution of GNSS‐inferred hydrologic loading and associated TWS estimates also increases. For example, Enzminger et al. (2018) used synthetic loads and tests of inversion reproducibility to show that including additional GNSS stations as constraints on the inversion improved solution accuracy.

While most hydrogeodetic studies have explored hydrologic loading at regional scales (Table 1), a handful of studies have investigated global‐scale (Blewitt et al., 2001; Dill & Dobslaw, 2013; Döll et al., 2014; Dong et al., 2002; Tesmer et al., 2011; Tregoning et al., 2009; van Dam et al., 2001; Wu et al., 2003) and local‐scale (Argus et al., 2020; Bevis et al., 2004, 2005; Elósegui, 2003; Knappe et al., 2019; Milliner et al., 2018; Wahr et al., 2013; Xue et al., 2021) hydrologic loading. Furthermore, increasing attention is being paid to TWS differences in regions with variable topography (e.g., Fok & Liu, 2019; Zhan et al., 2021) and aridity (e.g., Birhanu & Bendick, 2015; Han, 2017), which requires a high spatial resolution and an understanding of the sensitivity of GNSS measurements to hydrologic loading at various spatial scales. Probing the spatial resolution and limitations of GNSS is a current focus of hydrogeodetic studies.

Advancements in GNSS precision and post‐processing are also broadening the spatial resolution of hydrologic loading studies. Bevis et al. (2005) found that the vertical elastic response observed by GNSS is dominated by loads on the scale of approximately 100 km, while Argus et al. (2014) suggested that the spatial resolution of GNSS is approximately 50 km. In a comparison of synthetic inputs and modeled displacement inverted to estimate TWS at 1/8‐degree grid cells, Enzminger et al. (2018) showed that smoothing during the inversion process introduced load leakage that caused spatially variable differences between the inverted and input loads, especially in areas with the greatest snow loads. It is important to note that smoothing is almost always required for TWS calculations to obtain a stable solution as inversions are inherently nonunique, and regularization allows us to prioritize the better constrained solutions, which often creates leakage in the areas with the fewest constraints on the inversion (i.e., where station density is lowest). Enzminger et al. (2018) decreased leakage using their synthetic test results with grid‐cell linear scaling using empirically derived gain factors, producing near‐perfect agreement between input and inverted loads and accurate area‐averaged load estimates. Furthermore, the authors point out that sufficient GNSS station density and distribution are already in place in snow‐dominated regions in the western United States to allow for accurate estimation of annual variations in TWS from GNSS displacement at the mountain‐range scale.

Knappe et al. (2019) further reduced the scale of GNSS‐inferred hydrologic loading by showing that the regional signal can be represented by the common mode, a network‐coherent signal (similar to CME) calculated by stacking and averaging the smoothed signal for each GNSS station. Through the removal of the common‐mode regional signal, they separated the regional‐ and local‐scale contributions to the elastic deformation of Earth's crust to observe hydrologic loading with GNSS at the mountain watershed scale, which can be as small as 30 km, depending largely on the interstation distance. Increases in the spatial resolution of GNSS observations of hydrologic loading have significant implications for studies of water resources from high‐elevation, remote settings as GNSS observations of hydrologic loading fill the important spatial gap between meter‐scale point measurements and the 300 km resolution of GRACE TWS estimates.

Future Directions and Research Opportunities

We conclude with a few ideas of future research directions for GNSS‐based hydrogeodesy that could advance the science and enrich its societal impact. We offer these potential explorations as research questions that we believe are within the realm of the developing field of hydrogeodesy using GNSS. We also highlight exciting potential opportunities for further integrating hydrogeodesy with traditional hydrologic techniques and facilitating collaboration between hydrologists and geodesists. In addition, we point out that fewer than 20% of GNSS studies have been published in hydrology‐specific or multidisciplinary science journals and, instead, are primarily published in geophysical and remote sensing journals (Table 1). We assert that familiarizing hydrologists with GNSS hydrologic applications and conveying their utility will require a push toward publishing hydrogeodetic studies in multidisciplinary and water resources journals.

How Do Estimates of TWS Relate to Traditional Estimates of Hydrologic Storage?

Importantly, GNSS enables estimation of TWS changes integrated over the watershed scale between 10s and 102 km, whereas the direct measurement of hydrologic storage is not possible at this scale because of the large spatial heterogeneity of subsurface storage (Staudinger et al., 2017). Instead, hydrologists commonly rely on streamflow and turnover (or transit) time to infer watershed‐scale storage. This is a great example of the utility of GNSS‐based hydrogeodetic observations and an example of their advantages over the larger scale of GRACE and the smaller scale of point measurements. Yet, more work must be done in quantitatively relating TWS to traditional hydrologic storage. We propose that hydrogeodetic studies build on previous hydrology studies that infer watershed‐scale storage through a direct comparison with GNSS‐inferred TWS.

Furthermore, we question if immobile storage contributes to TWS measured in GNSS studies since immobile storage is thought to remain constant over time (Staudinger et al., 2017)? If that changes as the climate warms and droughts intensify, would GNSS‐based hydrogeodesy capture this change? We also ask how TWS relates to stream‐connected, dynamic storage like that estimated from the streamflow recession analysis (Kirchner, 2009). Dynamic storage is thought to be less than mobile storage as catchments can retain mobile water that may not be released as streamflow. Can GNSS provide a comparison of dynamic and mobile storage, the difference between which may be vital for evaluation of a catchment's resilience in the changing climate?

What Can We Learn About Hydrologic Partitioning From GNSS?

How can we use GNSS for a more thorough analysis of partitioning? For example, Enzminger et al. (2018) used SNODAS SWE and NLDAS soil moisture to show that the fraction of TWS at the mountain range scale (excluding groundwater, surface water, and reservoirs) that can be attributed to snow and soil moisture varies markedly over the time period of snowpack accumulation and states that, even in snow‐dominated regions like the Sierra Nevada mountains, inverted TWS cannot be attributed solely to SWE. However, the exclusion of groundwater using GRACE and the simplification of soil moisture from hydrologic land surface models, which are common across GNSS studies of hydrologic loading, preclude further analysis of partitioning.

The partitioning capabilities of GNSS hinge on the advancement of land surface models to estimate soil moisture, snow, and groundwater variability at spatial and temporal scales relevant to GNSS studies. We propose that this is an excellent opportunity for enhanced collaboration between hydrologists and geodesists. How can hydrologists improve the simplification of and reliance on land surface models that hydrogeodesy currently employs? Will the combination of GNSS with other methods (i.e., groundwater depths, snow depth from SNOTEL and altimetry, in situ and remotely sensed soil moisture, etc.) help? Better constraints on groundwater storage could inform our understanding of its potential to buffer long‐term trends in warming temperatures and increased water use and to understand the impact of short‐term effects of erratic weather patterns on groundwater resources. Moreover, we ask how improved investigation of partitioning with GNSS can inform the timing of water movement within the hydrosphere.

Can GNSS Improve Estimates of Hydrologic Fluxes?

Hydrologists commonly utilize a simple mass balance equation to calculate a hydrologic system's water budget, defining hydrologic storage (S) as the difference between system inputs (precipitation) and outputs (streamflow and evapotranspiration). In fact, some hydrogeodetic studies using GNSS have utilized the water balance equation to analyze the timing between changes in GNSS‐inferred TWS and changing hydrologic inputs and outputs in an effort to discern which hydrologic components best explain the observed TWS variation (e.g., Han, 2017; Yin et al., 2020).

Importantly, hydrogeodetic studies with GNSS present an opportunity to advance our use of the simple mass balance by recasting dS/dt (the typical unknown in water balance calculations) as an independently measured observation as has been done using TWS from GRACE (e.g., Chao et al., 2020). TWS from GNSS may further advance this work as it can be measured and incorporated into the water balance at smaller watershed scales and at shorter timescales, in which other hydrologic fluxes are more accurately measured. Therefore, we ask if GNSS enables closure of the water budget with greater certainty. We suggest that the use of TWS from GNSS in water balance analyses could also help to better constrain hydrologic fluxes with large uncertainties, such as evapotranspiration.

How Can Hydrologic Loading Applications of GNSS Improve Hydrologic Models?

Hydrologic models range from data‐driven empirical models to physics‐based distributed models and are commonly used to forecast streamflow and water resources and/or predict groundwater flow or hydrologic response. Yet, regardless of the type of model or its use, they all calibrate hydrologic response and train forecasts using streamflow. However, GNSS‐inferred TWS offers an additional integrated measure of watershed response that may be useful in calibrating and training hydrologic models. Therefore, we offer the following question to inspire future hydrogeodetic research using GNSS—how can hydrologic models incorporate observed changes in total water storage to best inform forecasting?

As we have discussed throughout this review, hydrogeodetic studies commonly explore the empirical relationships between Earth deformation under hydrologic loads and climatic drivers, which are a common input for hydrologic forecasting models. We assert that integration of our understanding of these relationships into hydrologic models is a promising direction for studies of GNSS‐inferred hydrologic loading and TWS change; very recent GNSS work has moved in this direction (Yin et al., 2021). The potential to discern between subsurface stores could have substantial benefits for hydrologic forecasting and water resource management as better resolution of soil moisture could constrain antecedent soil conditions that are important for predicting streamflow and the division of infiltration and runoff. We also suggest that GNSS studies further explore antecedent moisture conditions in an effort to improve flood and drought forecasting and capitalize on the near‐real time, watershed‐scale observations of GNSS hydrologic loading.

common mode error

contiguous United States

Global Navigation Satellite System

Global Positioning System

precipitable water vapor

snow water equivalent

terrestrial water storage