GNSS2TWS_Slepian: A software to recover daily GNSS-inverted terrestrial water storage changes based on Slepian basis functions

GNSS2TWS_Slepian: A software to recover daily GNSS-inverted terrestrial water storage changes based on Slepian basis functions | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article GNSS2TWS_Slepian: A software to recover daily GNSS-inverted terrestrial water storage changes based on Slepian basis functions Zhongshan Jiang, Miao Tang, Haiping Wen, Linguo Yuan This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4678987/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Changes in terrestrial water storage (TWS) can deform the Earth’s solid surface in the form of geodetically measurable vertical motions. Here, a new open-source Matlab software, named GNSS2TWS_Slepian, is developed to achieve the recovery of daily TWS changes from Global Navigation Satellite System (GNSS) crustal vertical positions. Differing from the widely-used spatial-domain inversion strategy based on Green's function method, our inversion modeling is implemented in the spectral domain based on Slepian basis functions, which aims to infer daily large-scale TWS changes using non-uniformly distributed GNSS vertical data. GNSS2TWS_Slepian is designed with different structured modules and the logic of the program workflow can be easily followed. To obtain daily estimates of TWS changes, the principal component analysis is integrated into our time-varying inversion model. To demonstrate the main functionalities, equivalent water height changes are investigated in the Western United States. This study aims to provide a scientific mathematical tool for resolving large-scale water mass loads, which is instrumental in broadening the applications of GNSS in hydrology. GNSS Slepian basis functions Terrestrial water storage Vertical motions Hydrological extremes Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 1 Introduction Terrestrial water storage (TWS) constantly oscillates attributed to climate changes and human activities and concomitantly stimulates various geophysical phenomena such as gravity changes and crustal deformation. Modern geodetic techniques can accurately sense these geophysical feedbacks and facilitate some advanced applications in hydrology (Chen et al. 2022 ; White et al. 2022 ). It contributes to the birth of an emerging discipline–hydrogeodesy, which aims to study spatiotemporal distribution characteristics of terrestrial water storage/fluxes using precise geodetic infrastructure. Through accurately measuring water-induced gravity changes or elastic deformation, geodetic techniques (e.g., Gravity Recovery and Climate Experiment (GRACE) and Global Navigation Satellite System (GNSS)) are used to quantify the variations in total water storage at different spatial (local, regional, or global) and temporal (monthly, annual, or multi-year) scales. The GRACE twin satellites are widely applied to the studies of polar ice sheet loss, continental water storage fluctuations, and extreme hydrological droughts by mapping detailed measurements of the changes in global gravitational fields (Chen et al. 2022 ; Harig and Simons 2012 ). However, the GRACE measurements have some inherent and unavoidable defects, such as the coarse temporal (monthly sampling) and spatial resolution (approximately 350 km) that result in the insensitivity to high-frequency and fine-scale hydrological mass variations and the data gaps (both long-term and random) which break the continuous analysis on some hydrological applications (Jiang et al. 2022b ; Yang et al. 2021 ). GNSS-measured Earth's vertical deformation in an elastic response to hydrological loads is used to infer surface mass changes at multiple spatial (global down to watershed) and temporal (multi-decadal down to daily) scales (Argus et al. 2017 ; Argus et al. 2022 ; Blewitt et al. 2001 ; Fok and Liu 2019 ; Fu et al. 2015 ; Hsu et al. 2020 ; Jiang et al. 2021a ; Knappe et al. 2019 ; Milliner et al. 2018 ; Zhong et al. 2020 ; Zhu et al. 2023 ). This inversion modeling depends on the mass loading theory that describes the physics law between surface loads and their caused crustal deformation, which can be expressed both in the spatial and spectral domains (Farrell 1972 ). Inversion of GNSS displacement data for large-scale surface mass loads is generally implemented using global spherical harmonic functions (Blewitt et al., 2001 ; Kusche and Schrama, 2005 ), which generally involves a physically motivated regularization due to the continent-rich, ocean-poor distribution of geodetic data (Kusche and Schrama 2005 ). However, this inversion of global GNSS data using the global spherical harmonic functions is generally applicable to the determination of low degree/order components of surface mass changes (Han and Razeghi 2017 ). In the medium scale (regional and continental), the Slepian basis functions (or localized spherical harmonic functions) are an alternative to traditional global spherical harmonics functions to realize the inversion of sparsely instrumented GNSS network for relatively large-scale surface mass changes (Han and Razeghi 2017 ; Jiang et al. 2021a ; Li et al. 2023 ). On the small scale (local and regional), the dense GNSS network is emerging as an additional tool to estimate fine-scale TWS variations based on the spatial-domain Green's function method (Argus et al. 2014 ; Fu et al. 2015 ; Jiang et al. 2021b ; Jiang et al. 2021c ; Jin and Zhang 2016 ). The GNSS-recorded vertical displacements are extremely sensitive to local- to regional-scale water mass changes, which is promising to realize a high-resolution spatiotemporal characterization of total TWS changes for the hydrological community. To facilitate more potential applications of GNSS in the hydrogeodetic community, we release another open-source Matlab package, named GNSS2TWS_Slepian. This numerical tool aims to invert unevenly distributed GNSS vertical data for recovering large-scale TWS changes. Compared to our previous open-source tool, GNSS2TWS (Jiang et al. 2022a ), which is appropriate to the spatial-domain inversion modeling based on Green's functions, the new inversion tool is implemented in the spectral domain based on Slepian basis functions. To demonstrate the GNSS2TWS_Slepian’s capability and functions, we take the Western United States (WUSA) as a study area and equivalent water height (EWH) estimates are characterized in combination with other hydrometeorological datasets. Our research would further broaden the hydrologic application of GNSS and is instrumental in water sustainability plans. 2 GNSS2TWS_Slepian software description GNSS2TWS_Slepian, an open-source tool developed in Matlab, is capable of inverting discrete GNSS vertical positions for large-scale TWS changes. We assume that potential readers have some background knowledge in geodesy and are familiar with the mass loading theory (Farrell 1972 ). Therefore, we adopt a scripting language that concatenates different structured modules so that we can conveniently modify the code to make it available for interested studies. Here, we introduce the design mentality and software workflow of our inversion scheme. 2.1 Inversion methodology We briefly summarize the inversion methodology of inferring large-scale TWS changes using GNSS vertical observations based on Slepian basis functions and more details are shown in Texts S1 and S2, as well as described in previous studies (Han and Razeghi 2017 ; Harig and Simons 2012 ; Simons et al. 2006 ). Taking vertical crustal displacement (VCD) as an example, the discrete VCD data on the Earth's surface can be expressed with localized spherical harmonic functions as follows: $$\begin{array}{c}{\mu }_{L}\left(\theta ,\lambda ,t\right)=a\sum _{\alpha =1}^{{\left(L+1\right)}^{2}}{s}_{\alpha }^{vcd}\left(t\right){g}_{\alpha }\left(\theta ,\lambda \right)\#\left(1\right)\end{array}$$ where \({\mu }_{L}\left(\theta ,\lambda ,t\right)\) is the time-varying vertical positions truncated with the degree of \(L\) , \({g}_{\alpha }\left(\theta ,\lambda \right)\) and \({s}_{\alpha }^{vcd}\left(t\right)\) denote the \(\alpha \text{t}\text{h}\) Slepian basis functions and corresponding VCD-related Slepian coefficient. In general, we need to truncate the Slepian basis functions with a small number of degrees, so that only a small number of Slepian coefficients are estimated using the least squares regularization. After obtaining the VCD-related Slepian coefficients, we convert them into EWH-related Slepian coefficients based on the mass loading theory (Farrell 1972 ), which can be written as follows (Han and Razeghi 2017 ): $$\begin{array}{c}{s}_{\beta }^{ewh}\left(t\right)=\sum _{\alpha =1}^{J}{s}_{\alpha }^{vcd}\left(t\right)\sum _{l=0}^{L}\sum _{m=-l}^{l}\left(\frac{{\rho }_{e}}{3{\rho }_{w}}\right)\frac{2l+1}{{h}_{l}^{{\prime }}}{g}_{\alpha ,lm}{g}_{\beta ,lm}\#\left(2\right)\end{array}$$ where \({s}_{\beta }^{ewh}\left(t\right)\) is the \(\beta\) th time-varying EWH-related Slepian coefficients, \({g}_{\alpha ,lm}\) and \({g}_{\beta ,lm}\) are the localization coefficients related to VCD and EWH, respectively. \({h}_{l}^{{\prime }}\) denotes \(l\) th-degree load Love numbers, \({\rho }_{w}\) and \({\rho }_{e}\) represent the density of water and Earth, \(J\) is the truncated degree of Slepian basis functions, which is generally determined according to the magnitude of eigenvalue \({\gamma }_{\alpha }\) . After obtaining the EWH-related Slepian coefficients, EWH estimates at any location can be calculated using the following formula: $$\begin{array}{c}{\sigma }_{L}\left(\theta ,\lambda ,t\right)=a\sum _{\alpha =1}^{J}{\gamma }_{\alpha }{s}_{\alpha }^{ewh}\left(t\right){g}_{\alpha }\left(\theta ,\lambda \right)\#\left(3\right)\end{array}$$ where \({\sigma }_{L}\left(\theta ,\lambda ,t\right)\) is the time-varying EWH estimates at the colatitude of \(\theta\) and longitude of \(\lambda\) . To make the synthesis more robust, an eigenvalue ( \({\gamma }_{\alpha }\) ) weighted strategy is generally applied according to Han and Razeghi ( 2017 ), which aims to damp the contributions from less accurately solved coefficients at higher \(\alpha\) . To improve the inversion efficiency, we design a time-varying inversion scheme by integrating the principal component analysis (PCA) method. The time-varying modeling only performs inversion several times depending on the number of selected PCs (only 6 times in our example), in contrast to the traditional time-consuming epoch-by-epoch model (e.g., 6209 times for the study period from Jan 1, 2006, to Dec 31, 2022). The PCA method embedded with the alternating least squares (ALS) algorithm can address that GNSS data inevitably have stochastic missing values and long-term data gaps due to random sensor failure, removal of outliers, or asynchronous installation dates. Our time-varying modeling strategy is briefly summarized in Fig. 1 . First, the GNSS vertical position time series at all stations are decomposed into several so-called PCs using the ALS-based PCA algorithm. Each PC consists of spatial functions and temporal functions, respectively. Afterward, the Slepian basis functions are generated using the open-source package SLEPIAN Alpha (Harig et al. 2015 ). To each PC, the spatial functions are converted into spectral-domain VCD-related Slepian coefficients and are then transformed into EWH-related Slepian coefficients. Eventually, we calculate grided EWH values for each PC and perform a reconstitution by summing the product of grided EWH values and the corresponding temporal function overall all PCs. 2.2 Software workflow The GNSS2TWS_Slepian software follows the workflow of GNSS2TWS (Jiang et al. 2022a ) by using different structured modules and readers can easily follow the logic of program execution. As demonstrated in Fig. 2 , the procedure workflow is composed of (1) loading scenario, (2) loading data, (3) PCA decomposition, (4) calculating eigenmatrix, (5) inversion modeling, and (6) displaying results. 1) Loading scenario Loading scenario is the first step, which mainly includes two aspects: parameter configuration and directory generation. Parameter configuration includes the assignment of relevant parameters and path settings of necessary files. The initialization also produces some directories that store the intermediate variables and final results. 2) Loading data All preprocessed GNSS vertical position time series stored in each file are collected to build an observation matrix with NaN for missing values. Besides, this step also collects some pertinent information regarding station name and location (i.e., longitude and latitude). Note that this software begins with the postprocessing GNSS data primarily associated with water cycles (see Section 4.1 for extracting hydrological load displacement). 3) PCA decomposition The built-in ALS-based PCA function is called to decompose the previously generated observation matrix and two groups of matrixes (i.e., spatial and temporal functions) are produced in this step. The determination of the PC’s number generally depends on that gradually increasing the PC’s number does not noticeably improve the fit to the raw data. 4) Calculating eigenmatrix This step calls SLEPIAN Alpha (Harig et al. 2015 ) to calculate the eigenvalues and eigenvectors as described in Text S1, which only rely on the specific study area and truncated degrees. The eigenvectors are used for the production of Slepian basis functions. 5) Inversion modeling In our inversion, the spatial functions of each PC are converted into spectral-domain VCD-related Slepian coefficients and are then transformed into EWH-related Slepian coefficients using load Love numbers (Wang et al. 2012 ). After that, grided EWH values for each PC are calculated and then multiplied by the corresponding temporal function, and the total EWH changes at all grids are synthesized by summing the foregoing products of all PCs. 6) Displaying results To quickly view inversion results, some figures are plotted for readers. This step displays figures of spatial and temporal functions, maps of annual EWH amplitudes, maps of selected Slepian basis functions, and figures of eigenvalues as the degree changes. 3 Application example: GNSS2TWS_Slepian applied to the WUSA region To demonstrate GNSS recovery of daily water height changes using the open-source GNSS2TWS_Slepian software, we take the WUSA region as our study area (Fig. 3 a). The WUSA region encompasses the westernmost U.S. states bounded by the Mississippi River and the Pacific Ocean. First, we extract hydrological loading displacements contaminated by other non-hydrological signals and these post-processed vertical position time series are inverted for the estimation of daily EWH values. At last, we evaluate our inversion results in comparison with other independent EWH measurements. 3.1 Extracting hydrological load displacement The GNSS position time series (Blewitt et al. 2018 ) processed with the GipsyX-1.0 software are downloaded from the Nevada Geodetic Laboratory (NGL, http://geodesy.unr.edu/ ), University of Nevada, Reno, United States. In total, 332 GNSS stations are well-chosen to infer large-scale TWS changes in the WUSA region (See more details in Text S3). To extract the water-induced vertical surface movements, the well-modeled predictions including non-tidal oceanic loading and atmospheric loading effects are subtracted from the previously downloaded GNSS time series. We then fit each GNSS time series with a constant, a linear, two sinusoids, and some Heavisideand and logarithmic functions when needed (Jiang et al. 2022b ). After the removal of constant, linear, offset, and postseismic terms, the residual time series containing sinusoid signals are dominated by seasonal and inter-annual hydrology-induced deformations (Figs. 3 b– 3 f), which feature upward motions as water runs off and subsidence when water loads. The comparison between GNSS-observed and PCA-reconstructed time series also indicates the good performance of ALS-based PCA in filling GNSS missing values and data gaps (Figs. 3 b– 3 f). 3.2 Inversion modeling for TWS changes in the WUSA region To characterize the water cycles in the WUSA region, the postprocessing GNSS vertical positions (i.e., residual time series generated in Section 3.1) are then inverted for daily EWH estimates. In the ALS-based PCA decomposition, we consider six PCs when no noticeable enhancement of the fit to the raw data is observed as the increase of selected PCs. These six PCs significantly contribute to the total data variance of the raw data, accounting for 82.5%. Each of these six PCs can explain 73.0%, 11.4%, 6.9%, 3.5%, 2.9%, and 2.4% of the filtered data (Figure S1 ). The first PC makes a substantial contribution to the variance compared to the other five PCs, and it reveals conspicuous annual and interannual hydrological loading variations and positive spatial responses at all GNSS stations (Figure S1 a). With the increase of the number of PCs, it manifests gradually weakened and more detailed spatiotemporal features of hydrological loading signals (Figures S1 b–1f). After obtaining the extended boundary data, region-specified Slepian basis functions are generated using SLEPIAN Alpha (Harig et al. 2015 ). In the inversion model, the bandwidth of the spherical harmonic degrees is chosen as \(L=80\) , corresponding to 6561 orthogonal Slepian basis functions and we only consider Slepian basis functions with eigenvalues of \({\gamma }_{\alpha }\ge 0.1\) (Figure S2). Figure S3 shows maps of the selected 67 Slepian basis functions. In our study area, the GNSS vertical data at 332 stations are converted into 67 VCD-related Slepian coefficients and then transformed into 67 EWH-related Slepian coefficients using PREM-based load Love numbers generated by Wang et al. ( 2012 ). To avoid the inversion instability of high-order coefficients, a 150-km Gaussian filter is applied in the recovery of water height changes. In case of having sufficient GNSS data, readers can try out higher/lower truncation degrees for more refined/coarser TWS estimated results. Figure S4 indicates relatively fine-scale EWH results when considering a large bandwidth value ( \(L=100\) ) in contrast to low-resolution water estimates using a lower truncation degree of \(L=60\) . In other words, fine-scale loading signals can be retrieved in a spectral-domain inversion model with high truncation orders and good GNSS spatial coverage, but it is very time-consuming and memory-consuming when calculating Slepian basis functions and corresponding coefficients. 3.3 Water height changes derived from various datasets We spatially compare annual amplitudes of GNSS-based water estimates with other independent products from GRACE Mascon solutions, North American Land Data Assimilation System (NLDAS), and Snow Data Assimilation System (SNODAS) models (Fig. 4 ). The results indicate that conspicuous seasonal water oscillations with annual amplitudes of > 150 mm are predominantly located along major large-scale mountain systems, including the Cascade Range, the Rocky, and Sierra Nevada Mountains. The NLDAS-EWH (Fig. 4 c) features higher resolution than water estimates from GNSS (Fig. 4 a) and GRACE (Fig. 4 b). The seasonal SNODAS-EWH oscillations have strong spatial consistency (a correlation coefficient of 0.94) with that in NLDAS-EWH results and both datasets reveal that the peak regions of annual water amplitudes are clustered along these major mountain systems. It demonstrates that seasonal snowfall is a key driving factor in the regional hydrologic system. Both GNSS and GRACE results underestimate the seasonal water oscillations along these main mountain systems attributed to low spatial resolution. Although we artificially lower the spatial coverage of GNSS stations, the GNSS-inferred results show more detailed seasonal changes than those derived from GRACE Mascon solutions (Figs. 4 a and 4 b). Figures 5 a– 5 e show time-varying basin-averaged water estimates from multiple datasets (i.e., GNSS, GRACE, NLDAS, and SNODAS). The GNSS-inverted water estimates feature temporal consistency with the GRACE-, NLDAS-EWH products in the Upper Colorado, Great Basin, Pacific Northeast, and California watersheds. In the lower Colorado River basin, extremely low correlations are found between GNSS, GRACE, and NLDAS results, due to the weak signals in seasonal water changes and noises in respective measurements (Fig. 5 b). The Pacific Northwest hydrologic unit is characterized by the most conspicuous seasonality in water storage variations and the annual amplitudes can reach up to 109, 117, and 138 mm for GNSS-, GRACE- and NLDAS-EWH results (Fig. 5 d). The California watershed (Fig. 5 e) has weaker seasonal oscillations (73, 81, and 93 mm for GNSS, GRACE, and NLDAS) than those in the Pacific Northeast basin. Both Upper Colorado and Great Basin watersheds present similar amplitudes of approximately 50 mm (Figs. 5 a and 5 c). All water estimates indicate interannual variation characteristics that are caused by extreme hydrometeorological events. To quantify the intensity of these extremes, we convert the GNSS-inferred water estimates into drought severity index (DSI) products (see more details in Text S4). Figures 5 f– 5 j present a comparison between GNSS-DSI and other DSI datasets derived from GRACE and NLDAS in these five river basins. In general, GNSS-DSI can temporally correspond with GRACE-DSI and NLDAS-DSI. In the California, Great, and Pacific Northeast Basins, the GNSS results are simultaneously consistent with GRACE and NLDAS estimates, with correlation coefficients between 0.51–0.85. Good agreement (0.60) is found between GNSS and GRACE results in the Upper Colorado Basin, but a bad correlation (0.23) exists between GNSS and NLDAS (Fig. 5 f). In the Lower Colorado Basin, a low correlation is seen in the GNSS-GRACE and GNSS-NLDAS results (Fig. 5 g). 4 Conclusion and recommendations This study releases a numerical calculation tool, named GNSS2TWS_Slepian, for inferring daily large-scale water storage changes using GNSS vertical position time series. This inversion modeling of GNSS data is implemented in the spectral domain based on Slepian basis functions, in contrast to the widely-used spatial-domain Green's function method. To improve the inversion efficiency, a PCA-based inversion method is designed to recover time-varying water height changes. To demonstrate the main functionalities and applications, we investigate the daily TWS changes and characterize hydrological extremes in the Western United States. Through contrastive analysis, the GNSS-inverted water estimates have good spatiotemporal correspondence with other water measurements derived from GRACE and NLDAS. The GNSS characterization result of hydrological extremes illuminates the availability of GNSS for drought monitoring and assessment. With the optimization and enhancement of continuously operating GNSS networks worldwide, inversion modeling of crustal deformation measurements is very likely to become a powerful means in the recovery of water height changes for the hydrological community. It is very promising and necessary to explore more application prospects using GNSS for more extensive hydrogeodetic studies, which would contribute to our better understanding of water cycles, hydrological extremes, hydrological dynamics, and hydroclimatic interaction. Therefore, the potential users of GNSS2TWS_Slepian are encouraged to utilize our tools to develop more hydrogeodetic applications and to achieve deeper interdisciplinary collaborations with hydrological communities. Statements and Declarations Acknowledgments Many thanks are sent for the valuable suggestions of the editors and reviewers. The open-source Slepian Alpha software is available from GitHub (https://github.com/csdms-contrib/slepian_alpha). Some figures are created using the Generic Mapping Tools (version 6.4.0) software (Wessel et al. 2019). Ethics approval and consent to participate Not applicable. Funding This work was supported by the National Natural Science Foundation of China (Grant No. 41904015). Data Availability The GNSS vertical position time series with extra columns are available from the Nevada Geodetic Laboratory, University of Nevada, Reno (http://geodesy.unr.edu/gps_timeseries/tenv3_loadpredictions). The GRACE RL06 Mascon products are downloaded from the CSR’s webpage (https://www2.csr.utexas.edu/grace/). The NLDAS model is available from the NASA GES DISC website (https://disc.gsfc.nasa.gov/datasets/NLDAS_NOAH0125_M_002/summary?keywords=NLDAS). The SNODAS products are from the NOAA National Weather Service’s National Operational Hydrologic Remote Sensing Center (https://nsidc.org/data/g02158/versions/1). Our inversion code for large-scale water storage changes based on Slepian basis functions is released on GitHub (https://github.com/jzshhh/gnss2tws_slepian). Authors' contributions Zhongshan Jiang: Conceptualization, Methodology, Software, Funding, and Writing- Original draft preparation. Miao Tang: Writing- Original draft preparation. Haiping Wen: Writing- Original draft preparation. Linguo Yuan: Supervision and Writing - Reviewing and Editing. Author Biographies Zhongshan Jiang is an associate professor at the School of Geospatial Engineering and Science, Sun Yat-Sen University, China. His research interests include seismogeodesy and hydrogeodesy. Miao Tang is a postgraduate student at the Faculty of Geosciences and Engineering, Southwest Jiaotong University, China. His research interest is hydrogeodesy. Haiping Wen is a postgraduate student at the Faculty of Geosciences and Engineering, Southwest Jiaotong University, China. His research interest is hydrogeodesy. Linguo Yuan is a professor at the Faculty of Geosciences and Engineering, Southwest Jiaotong University, China. His research interests include GPS data processing and application and solid Earth tides. References Argus, D.F., Fu, Y., & Landerer, F.W. (2014). Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion. Geophysical Research Letters, 41 , 1971-1980. doi: 10.1002/2014gl059570 Argus, D.F., Landerer, F.W., Wiese, D.N., Martens, H.R., Fu, Y., Famiglietti, J.S., Thomas, B.F., Farr, T.G., Moore, A.W., & Watkins, M.M. (2017). Sustained Water Loss in California's Mountain Ranges During Severe Drought From 2012 to 2015 Inferred From GPS. Journal of Geophysical Research: Solid Earth, 122 . doi: 10.1002/2017jb014424 Argus, D.F., Martens, H.R., Borsa, A.A., Knappe, E., Wiese, D.N., Alam, S., Anderson, M., Khatiwada, A., Lau, N., Peidou, A., Swarr, M., White, A.M., Bos, M.S., Ellmer, M., Landerer, F.W., & Gardiner, W.P. (2022). Subsurface Water Flux in California's Central Valley and Its Source Watershed From Space Geodesy. Geophysical Research Letters, 49 . doi: 10.1029/2022gl099583 Blewitt, G., Hammond, W., & Kreemer, C. (2018). Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos, 99 . doi: 10.1029/2018eo104623 Blewitt, G., Lavallee, D., Clarke, P., & Nurutdinov, K. (2001). A new global mode of Earth deformation: seasonal cycle detected. Science, 294 , 2342-2345. doi: 10.1126/science.1065328 Chen, J., Cazenave, A., Dahle, C., Llovel, W., Panet, I., Pfeffer, J., & Moreira, L. (2022). Applications and Challenges of GRACE and GRACE Follow-On Satellite Gravimetry. Surv Geophys, 43 , 305-345. doi: 10.1007/s10712-021-09685-x Farrell, W.E. (1972). Deformation of the Earth by surface loads. Reviews of Geophysics, 10 . doi: 10.1029/RG010i003p00761 Fok, H.S., & Liu, Y. (2019). An Improved GPS-Inferred Seasonal Terrestrial Water Storage Using Terrain-Corrected Vertical Crustal Displacements Constrained by GRACE. Remote Sensing, 11 . doi: 10.3390/rs11121433 Fu, Y., Argus, D.F., & Landerer, F.W. (2015). GPS as an independent measurement to estimate terrestrial water storage variations in Washington and Oregon. Journal of Geophysical Research: Solid Earth, 120 , 552-566. doi: 10.1002/2014jb011415 Han, S.C., & Razeghi, S.M. (2017). GPS Recovery of Daily Hydrologic and Atmospheric Mass Variation: A Methodology and Results From the Australian Continent. Journal of Geophysical Research: Solid Earth, 122 , 9328-9343. doi: 10.1002/2017jb014603 Harig, C., Lewis, K.W., Plattner, A., & Simons, F.J. (2015). A Suite of Software Analyzes Data on the Sphere. Eos, 96 . doi: 10.1029/2015eo025851 Harig, C., & Simons, F.J. (2012). Mapping Greenland's mass loss in space and time. Proc Natl Acad Sci U S A, 109 , 19934-19937. doi: 10.1073/pnas.1206785109 Hsu, Y.-J., Fu, Y., Bürgmann, R., Hsu, S.-Y., Lin, C.-C., Tang, C.-H., & Wu, Y.-M. (2020). Assessing seasonal and interannual water storage variations in Taiwan using geodetic and hydrological data. Earth and Planetary Science Letters, 550 . doi: 10.1016/j.epsl.2020.116532 Jiang, Z., Hsu, Y.-J., Yuan, L., Cheng, S., Li, Q., & Li, M. (2021a). Estimation of daily hydrological mass changes using continuous GNSS measurements in mainland China. Journal of Hydrology, 598 . doi: 10.1016/j.jhydrol.2021.126349 Jiang, Z., Hsu, Y.-J., Yuan, L., Feng, W., Yang, X., & Tang, M. (2022a). GNSS2TWS: an open-source MATLAB-based tool for inferring daily terrestrial water storage changes using GNSS vertical data. GPS Solutions, 26 . doi: 10.1007/s10291-022-01301-8 Jiang, Z., Hsu, Y.-J., Yuan, L., & Huang, D. (2021b). Monitoring time-varying terrestrial water storage changes using daily GNSS measurements in Yunnan, southwest China. Remote Sensing of Environment, 254 . doi: 10.1016/j.rse.2020.112249 Jiang, Z., Hsu, Y.J., Yuan, L., Tang, M., Yang, X., & Yang, X. (2022b). Hydrological drought characterization based on GNSS imaging of vertical crustal deformation across the contiguous United States. Sci Total Environ, 823 , 153663. doi: 10.1016/j.scitotenv.2022.153663 Jiang, Z., Hsu, Y.J., Yuan, L., Yang, X., Ding, Y., Tang, M., & Chen, C. (2021c). Characterizing Spatiotemporal Patterns of Terrestrial Water Storage Variations Using GNSS Vertical Data in Sichuan, China. Journal of Geophysical Research: Solid Earth, 126 . doi: 10.1029/2021jb022398 Jin, S., & Zhang, T. (2016). Terrestrial Water Storage Anomalies Associated with Drought in Southwestern USA from GPS Observations. Surveys in Geophysics, 37 , 1139-1156. doi: 10.1007/s10712-016-9385-z Knappe, E., Bendick, R., Martens, H.R., Argus, D.F., & Gardner, W.P. (2019). Downscaling Vertical GPS Observations to Derive Watershed‐Scale Hydrologic Loading in the Northern Rockies. Water Resources Research, 55 , 391-401. doi: 10.1029/2018wr023289 Kusche, J., & Schrama, E.J.O. (2005). Surface mass redistribution inversion from global GPS deformation and Gravity Recovery and Climate Experiment (GRACE) gravity data. Journal of Geophysical Research: Solid Earth, 110 . doi: 10.1029/2004jb003556 Li, X., Zhong, B., Li, J., & Liu, R. (2023). Inversion of GNSS Vertical Displacements for Terrestrial Water Storage Changes Using Slepian Basis Functions. Earth and Space Science, 10 . doi: 10.1029/2022ea002608 Milliner, C., Materna, K., Burgmann, R., Fu, Y., Moore, A.W., Bekaert, D., Adhikari, S., & Argus, D.F. (2018). Tracking the weight of Hurricane Harvey's stormwater using GPS data. Sci Adv, 4 , eaau2477. doi: 10.1126/sciadv.aau2477 Simons, F.J., Dahlen, F.A., & Wieczorek, M.A. (2006). Spatiospectral Concentration on a Sphere. SIAM Review, 48 , 504-536. doi: 10.1137/s0036144504445765 Wang, H., Xiang, L., Jia, L., Jiang, L., Wang, Z., Hu, B., & Gao, P. (2012). Load Love numbers and Green's functions for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0. Computers & Geosciences, 49 , 190-199. doi: 10.1016/j.cageo.2012.06.022 Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F., & Tian, D. (2019). The Generic Mapping Tools Version 6. Geochemistry, Geophysics, Geosystems, 20 , 5556-5564. doi: 10.1029/2019gc008515 White, A.M., Gardner, W.P., Borsa, A.A., Argus, D.F., & Martens, H.R. (2022). A Review of GNSS/GPS in Hydrogeodesy: Hydrologic Loading Applications and Their Implications for Water Resource Research. Water Resour Res, 58 , e2022WR032078. doi: 10.1029/2022WR032078 Yang, X., Tian, S., You, W., & Jiang, Z. (2021). Reconstruction of continuous GRACE/GRACE-FO terrestrial water storage anomalies based on time series decomposition. Journal of Hydrology, 603 . doi: 10.1016/j.jhydrol.2021.127018 Zhong, B., Li, X., Chen, J., Li, Q., & Liu, T. (2020). Surface Mass Variations from GPS and GRACE/GFO: A Case Study in Southwest China. Remote Sensing, 12 . doi: 10.3390/rs12111835 Zhu, H., Chen, K., Hu, S., Liu, J., Shi, H., Wei, G., Chai, H., Li, J., & Wang, T. (2023). Using the Global Navigation Satellite System and Precipitation Data to Establish the Propagation Characteristics of Meteorological and Hydrological Drought in Yunnan, China. Water Resources Research, 59 . doi: 10.1029/2022wr033126 Additional Declarations No competing interests reported. Supplementary Files Supplementary.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4678987","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":337405323,"identity":"b29155ed-dfcd-4238-93f1-95dc29436d71","order_by":0,"name":"Zhongshan Jiang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABHUlEQVRIiWNgGAWjYBACxmYwdYDBgL2BgSGhACbORowWngNALQZEaIECoBaJBCBNjBbmduZnD7/8uWNvLvk68cMDAwZ7g/NnDBg+lB1m4J/dgMNhbObGMjzPEnfOzt0sAXRY4oYbOQaMM84dZpC4cwCXX8ykJSQOJxjczt0A0pJgcIPHgJm37TDUqdi0sH+TljA4bG9w8+zmHwlQhzH/xauFx0zyQ8Jhxg03eLeBbGHccCDHgJkRv5YyaYYDhxM3nMndZpFgIJE480ZawcGec+k8EjewazHsP75N8scfoMOOn91880eFjT3f+cMbH/wos5bjn4FDSwMwoHkQfAkGhQOgaGJg4MGqHgjkQY77gSLSgEvtKBgFo2AUjFQAAOidYiI9frkSAAAAAElFTkSuQmCC","orcid":"","institution":"Sun Yat-sen University","correspondingAuthor":true,"prefix":"","firstName":"Zhongshan","middleName":"","lastName":"Jiang","suffix":""},{"id":337405324,"identity":"ea66d7b0-b360-4205-9c68-8bb73fad011e","order_by":1,"name":"Miao Tang","email":"","orcid":"","institution":"Southwest Jiaotong University","correspondingAuthor":false,"prefix":"","firstName":"Miao","middleName":"","lastName":"Tang","suffix":""},{"id":337405325,"identity":"955a3fa6-f6ce-4195-9c64-13e0a004cd6b","order_by":2,"name":"Haiping Wen","email":"","orcid":"","institution":"Southwest Jiaotong University","correspondingAuthor":false,"prefix":"","firstName":"Haiping","middleName":"","lastName":"Wen","suffix":""},{"id":337405326,"identity":"a55283d3-2017-4964-b0b3-226b97045182","order_by":3,"name":"Linguo Yuan","email":"","orcid":"","institution":"Southwest Jiaotong University","correspondingAuthor":false,"prefix":"","firstName":"Linguo","middleName":"","lastName":"Yuan","suffix":""}],"badges":[],"createdAt":"2024-07-03 08:53:32","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4678987/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4678987/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":62081513,"identity":"80868623-24de-42bf-a7a3-1a8da102a5da","added_by":"auto","created_at":"2024-08-09 05:42:44","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":128853,"visible":true,"origin":"","legend":"\u003cp\u003eInversion framework of the GNSS2TWS_Slepian software. The PCA is integrated to achieve a time-varying inversion model. This flowchart takes the WUSAas an example.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/b38ee52d9cc2b260be86befd.png"},{"id":62082049,"identity":"f2f41341-c20d-4e5d-9fe1-570b0af3cc4f","added_by":"auto","created_at":"2024-08-09 05:50:44","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":116249,"visible":true,"origin":"","legend":"\u003cp\u003eProgram workflow and main modules of GNSS2TWS_Slepian, which is modified from Jiang et al. (2022a).\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/552e363dd824eb32c768d1bd.png"},{"id":62081515,"identity":"26c851f4-2c36-4eb1-9836-da1cfa11ec78","added_by":"auto","created_at":"2024-08-09 05:42:45","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":367497,"visible":true,"origin":"","legend":"\u003cp\u003e(a) Geographic settings and GNSS observations used in our study. Red circles show the locations of GNSS stations. Red lines show the boundary of five major watersheds, i.e., the Upper Colorado (A), Lower Colorado (B), Great Basin (C), Pacific Northeast (D), and California (E) watersheds. (b)–(f) show examples of GNSS time series marked with green labels in (a). Red and blue lines show GNSS observed and PCA-reconstructed time series, respectively.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/1de5ca38300468aa0000b503.png"},{"id":62082050,"identity":"c52146d8-ab05-4421-bd93-0a39f73fd15b","added_by":"auto","created_at":"2024-08-09 05:50:44","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":458284,"visible":true,"origin":"","legend":"\u003cp\u003eAnnual EWH amplitudes derived from GNSS inversion (a), GRACE (b), NLDAS (c), and SNODAS (d). Note that the snow water equivalent products in the NLDAS model are replaced with those from the SNODAS model.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/94ab57608332a767d1756f45.png"},{"id":62081511,"identity":"933996cb-f558-4afe-b705-97ef54505c44","added_by":"auto","created_at":"2024-08-09 05:42:44","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":316691,"visible":true,"origin":"","legend":"\u003cp\u003eTime-varying water estimates (a–e) and various DSI datasets (f–j) derived from GNSS (red line), GRACE (black line), NLDAS (blue line) in the Upper Colorado, Lower Colorado, Great Basin, Pacific Northeast, and California watersheds. Green bars indicate SNODAS-based snow water equivalent results. R1 (GNSS-GRACE), R2 (GNSS-NLDAS), and R3 (GRACE-NLDAS) show the correlation coefficients between various DSI datastes.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/cf1ac1aaf1aad08c0c87ac1a.png"},{"id":62082382,"identity":"a2cb98d2-8d33-41bd-b96d-67a574bff73c","added_by":"auto","created_at":"2024-08-09 05:58:50","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1851673,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/c83f8f7d-4daf-476d-bab0-a78cb2e03b0f.pdf"},{"id":62081516,"identity":"25aefff2-25af-4d41-81ea-14d10c0e9751","added_by":"auto","created_at":"2024-08-09 05:42:45","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":2854209,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementary.docx","url":"https://assets-eu.researchsquare.com/files/rs-4678987/v1/127ec3aa805815b91d5c10c8.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"GNSS2TWS_Slepian: A software to recover daily GNSS-inverted terrestrial water storage changes based on Slepian basis functions","fulltext":[{"header":"1 Introduction","content":"\u003cp\u003eTerrestrial water storage (TWS) constantly oscillates attributed to climate changes and human activities and concomitantly stimulates various geophysical phenomena such as gravity changes and crustal deformation. Modern geodetic techniques can accurately sense these geophysical feedbacks and facilitate some advanced applications in hydrology (Chen et al. \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; White et al. \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). It contributes to the birth of an emerging discipline\u0026ndash;hydrogeodesy, which aims to study spatiotemporal distribution characteristics of terrestrial water storage/fluxes using precise geodetic infrastructure. Through accurately measuring water-induced gravity changes or elastic deformation, geodetic techniques (e.g., Gravity Recovery and Climate Experiment (GRACE) and Global Navigation Satellite System (GNSS)) are used to quantify the variations in total water storage at different spatial (local, regional, or global) and temporal (monthly, annual, or multi-year) scales. The GRACE twin satellites are widely applied to the studies of polar ice sheet loss, continental water storage fluctuations, and extreme hydrological droughts by mapping detailed measurements of the changes in global gravitational fields (Chen et al. \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Harig and Simons \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). However, the GRACE measurements have some inherent and unavoidable defects, such as the coarse temporal (monthly sampling) and spatial resolution (approximately 350 km) that result in the insensitivity to high-frequency and fine-scale hydrological mass variations and the data gaps (both long-term and random) which break the continuous analysis on some hydrological applications (Jiang et al. \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2022b\u003c/span\u003e; Yang et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2021\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eGNSS-measured Earth's vertical deformation in an elastic response to hydrological loads is used to infer surface mass changes at multiple spatial (global down to watershed) and temporal (multi-decadal down to daily) scales (Argus et al. \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Argus et al. \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Blewitt et al. \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Fok and Liu \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Fu et al. \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Hsu et al. \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Jiang et al. \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2021a\u003c/span\u003e; Knappe et al. \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Milliner et al. \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Zhong et al. \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Zhu et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). This inversion modeling depends on the mass loading theory that describes the physics law between surface loads and their caused crustal deformation, which can be expressed both in the spatial and spectral domains (Farrell \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e1972\u003c/span\u003e). Inversion of GNSS displacement data for large-scale surface mass loads is generally implemented using global spherical harmonic functions (Blewitt et al., \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Kusche and Schrama, \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2005\u003c/span\u003e), which generally involves a physically motivated regularization due to the continent-rich, ocean-poor distribution of geodetic data (Kusche and Schrama \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2005\u003c/span\u003e). However, this inversion of global GNSS data using the global spherical harmonic functions is generally applicable to the determination of low degree/order components of surface mass changes (Han and Razeghi \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). In the medium scale (regional and continental), the Slepian basis functions (or localized spherical harmonic functions) are an alternative to traditional global spherical harmonics functions to realize the inversion of sparsely instrumented GNSS network for relatively large-scale surface mass changes (Han and Razeghi \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Jiang et al. \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2021a\u003c/span\u003e; Li et al. \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). On the small scale (local and regional), the dense GNSS network is emerging as an additional tool to estimate fine-scale TWS variations based on the spatial-domain Green's function method (Argus et al. \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2014\u003c/span\u003e; Fu et al. \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Jiang et al. \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2021b\u003c/span\u003e; Jiang et al. \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2021c\u003c/span\u003e; Jin and Zhang \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). The GNSS-recorded vertical displacements are extremely sensitive to local- to regional-scale water mass changes, which is promising to realize a high-resolution spatiotemporal characterization of total TWS changes for the hydrological community.\u003c/p\u003e \u003cp\u003eTo facilitate more potential applications of GNSS in the hydrogeodetic community, we release another open-source Matlab package, named GNSS2TWS_Slepian. This numerical tool aims to invert unevenly distributed GNSS vertical data for recovering large-scale TWS changes. Compared to our previous open-source tool, GNSS2TWS (Jiang et al. \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022a\u003c/span\u003e), which is appropriate to the spatial-domain inversion modeling based on Green's functions, the new inversion tool is implemented in the spectral domain based on Slepian basis functions. To demonstrate the GNSS2TWS_Slepian\u0026rsquo;s capability and functions, we take the Western United States (WUSA) as a study area and equivalent water height (EWH) estimates are characterized in combination with other hydrometeorological datasets. Our research would further broaden the hydrologic application of GNSS and is instrumental in water sustainability plans.\u003c/p\u003e"},{"header":"2 GNSS2TWS_Slepian software description","content":"\u003cp\u003eGNSS2TWS_Slepian, an open-source tool developed in Matlab, is capable of inverting discrete GNSS vertical positions for large-scale TWS changes. We assume that potential readers have some background knowledge in geodesy and are familiar with the mass loading theory (Farrell \u003cspan class=\"CitationRef\"\u003e1972\u003c/span\u003e). Therefore, we adopt a scripting language that concatenates different structured modules so that we can conveniently modify the code to make it available for interested studies. Here, we introduce the design mentality and software workflow of our inversion scheme.\u003c/p\u003e\n\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\n \u003ch2\u003e2.1 Inversion methodology\u003c/h2\u003e\n \u003cp\u003eWe briefly summarize the inversion methodology of inferring large-scale TWS changes using GNSS vertical observations based on Slepian basis functions and more details are shown in Texts S1 and S2, as well as described in previous studies (Han and Razeghi \u003cspan class=\"CitationRef\"\u003e2017\u003c/span\u003e; Harig and Simons \u003cspan class=\"CitationRef\"\u003e2012\u003c/span\u003e; Simons et al. \u003cspan class=\"CitationRef\"\u003e2006\u003c/span\u003e). Taking vertical crustal displacement (VCD) as an example, the discrete VCD data on the Earth's surface can be expressed with localized spherical harmonic functions as follows:\u003c/p\u003e\n \u003cdiv id=\"Equa\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e$$\\begin{array}{c}{\\mu }_{L}\\left(\\theta ,\\lambda ,t\\right)=a\\sum _{\\alpha =1}^{{\\left(L+1\\right)}^{2}}{s}_{\\alpha }^{vcd}\\left(t\\right){g}_{\\alpha }\\left(\\theta ,\\lambda \\right)\\#\\left(1\\right)\\end{array}$$\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003e\u003cbr\u003e\u003c/p\u003e\n \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\mu }_{L}\\left(\\theta ,\\lambda ,t\\right)\\)\u003c/span\u003e\u003c/span\u003e is the time-varying vertical positions truncated with the degree of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(L\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({g}_{\\alpha }\\left(\\theta ,\\lambda \\right)\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({s}_{\\alpha }^{vcd}\\left(t\\right)\\)\u003c/span\u003e\u003c/span\u003e denote the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\alpha \\text{t}\\text{h}\\)\u003c/span\u003e\u003c/span\u003e Slepian basis functions and corresponding VCD-related Slepian coefficient.\u003c/p\u003e\n \u003cp\u003eIn general, we need to truncate the Slepian basis functions with a small number of degrees, so that only a small number of Slepian coefficients are estimated using the least squares regularization. After obtaining the VCD-related Slepian coefficients, we convert them into EWH-related Slepian coefficients based on the mass loading theory (Farrell \u003cspan class=\"CitationRef\"\u003e1972\u003c/span\u003e), which can be written as follows (Han and Razeghi \u003cspan class=\"CitationRef\"\u003e2017\u003c/span\u003e):\u003c/p\u003e\n \u003cdiv id=\"Equb\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e$$\\begin{array}{c}{s}_{\\beta }^{ewh}\\left(t\\right)=\\sum _{\\alpha =1}^{J}{s}_{\\alpha }^{vcd}\\left(t\\right)\\sum _{l=0}^{L}\\sum _{m=-l}^{l}\\left(\\frac{{\\rho }_{e}}{3{\\rho }_{w}}\\right)\\frac{2l+1}{{h}_{l}^{{\\prime }}}{g}_{\\alpha ,lm}{g}_{\\beta ,lm}\\#\\left(2\\right)\\end{array}$$\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({s}_{\\beta }^{ewh}\\left(t\\right)\\)\u003c/span\u003e\u003c/span\u003e is the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\beta\\)\u003c/span\u003e\u003c/span\u003eth time-varying EWH-related Slepian coefficients, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({g}_{\\alpha ,lm}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({g}_{\\beta ,lm}\\)\u003c/span\u003e\u003c/span\u003e are the localization coefficients related to VCD and EWH, respectively. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({h}_{l}^{{\\prime }}\\)\u003c/span\u003e\u003c/span\u003e denotes \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(l\\)\u003c/span\u003e\u003c/span\u003eth-degree load Love numbers, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\rho }_{w}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\rho }_{e}\\)\u003c/span\u003e\u003c/span\u003e represent the density of water and Earth, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(J\\)\u003c/span\u003e\u003c/span\u003e is the truncated degree of Slepian basis functions, which is generally determined according to the magnitude of eigenvalue \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\gamma }_{\\alpha }\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e\n \u003cp\u003eAfter obtaining the EWH-related Slepian coefficients, EWH estimates at any location can be calculated using the following formula:\u003c/p\u003e\n \u003cdiv id=\"Equc\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e$$\\begin{array}{c}{\\sigma }_{L}\\left(\\theta ,\\lambda ,t\\right)=a\\sum _{\\alpha =1}^{J}{\\gamma }_{\\alpha }{s}_{\\alpha }^{ewh}\\left(t\\right){g}_{\\alpha }\\left(\\theta ,\\lambda \\right)\\#\\left(3\\right)\\end{array}$$\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\sigma }_{L}\\left(\\theta ,\\lambda ,t\\right)\\)\u003c/span\u003e\u003c/span\u003e is the time-varying EWH estimates at the colatitude of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\theta\\)\u003c/span\u003e\u003c/span\u003e and longitude of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\lambda\\)\u003c/span\u003e\u003c/span\u003e. To make the synthesis more robust, an eigenvalue (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\gamma }_{\\alpha }\\)\u003c/span\u003e\u003c/span\u003e) weighted strategy is generally applied according to Han and Razeghi (\u003cspan class=\"CitationRef\"\u003e2017\u003c/span\u003e), which aims to damp the contributions from less accurately solved coefficients at higher \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\alpha\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e\n \u003cp\u003eTo improve the inversion efficiency, we design a time-varying inversion scheme by integrating the principal component analysis (PCA) method. The time-varying modeling only performs inversion several times depending on the number of selected PCs (only 6 times in our example), in contrast to the traditional time-consuming epoch-by-epoch model (e.g., 6209 times for the study period from Jan 1, 2006, to Dec 31, 2022). The PCA method embedded with the alternating least squares (ALS) algorithm can address that GNSS data inevitably have stochastic missing values and long-term data gaps due to random sensor failure, removal of outliers, or asynchronous installation dates.\u003c/p\u003e\n \u003cp\u003eOur time-varying modeling strategy is briefly summarized in Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e. First, the GNSS vertical position time series at all stations are decomposed into several so-called PCs using the ALS-based PCA algorithm. Each PC consists of spatial functions and temporal functions, respectively. Afterward, the Slepian basis functions are generated using the open-source package SLEPIAN Alpha (Harig et al. \u003cspan class=\"CitationRef\"\u003e2015\u003c/span\u003e). To each PC, the spatial functions are converted into spectral-domain VCD-related Slepian coefficients and are then transformed into EWH-related Slepian coefficients. Eventually, we calculate grided EWH values for each PC and perform a reconstitution by summing the product of grided EWH values and the corresponding temporal function overall all PCs.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\n \u003ch2\u003e2.2 Software workflow\u003c/h2\u003e\n \u003cp\u003eThe GNSS2TWS_Slepian software follows the workflow of GNSS2TWS (Jiang et al. \u003cspan class=\"CitationRef\"\u003e2022a\u003c/span\u003e) by using different structured modules and readers can easily follow the logic of program execution. As demonstrated in Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e, the procedure workflow is composed of (1) loading scenario, (2) loading data, (3) PCA decomposition, (4) calculating eigenmatrix, (5) inversion modeling, and (6) displaying results.\u003c/p\u003e\n\u003c/div\u003e\n\u003ch3\u003e1) Loading scenario\u003c/h3\u003e\n\u003cp\u003eLoading scenario is the first step, which mainly includes two aspects: parameter configuration and directory generation. Parameter configuration includes the assignment of relevant parameters and path settings of necessary files. The initialization also produces some directories that store the intermediate variables and final results.\u003c/p\u003e\n\u003ch3\u003e2) Loading data\u003c/h3\u003e\n\u003cp\u003eAll preprocessed GNSS vertical position time series stored in each file are collected to build an observation matrix with NaN for missing values. Besides, this step also collects some pertinent information regarding station name and location (i.e., longitude and latitude). Note that this software begins with the postprocessing GNSS data primarily associated with water cycles (see Section 4.1 for extracting hydrological load displacement).\u003c/p\u003e\n\u003ch3\u003e3) PCA decomposition\u003c/h3\u003e\n\u003cp\u003eThe built-in ALS-based PCA function is called to decompose the previously generated observation matrix and two groups of matrixes (i.e., spatial and temporal functions) are produced in this step. The determination of the PC’s number generally depends on that gradually increasing the PC’s number does not noticeably improve the fit to the raw data.\u003c/p\u003e\n\u003ch3\u003e4) Calculating eigenmatrix\u003c/h3\u003e\n\u003cp\u003eThis step calls SLEPIAN Alpha (Harig et al. \u003cspan class=\"CitationRef\"\u003e2015\u003c/span\u003e) to calculate the eigenvalues and eigenvectors as described in Text S1, which only rely on the specific study area and truncated degrees. The eigenvectors are used for the production of Slepian basis functions.\u003c/p\u003e\n\u003ch3\u003e5) Inversion modeling\u003c/h3\u003e\n\u003cp\u003eIn our inversion, the spatial functions of each PC are converted into spectral-domain VCD-related Slepian coefficients and are then transformed into EWH-related Slepian coefficients using load Love numbers (Wang et al. \u003cspan class=\"CitationRef\"\u003e2012\u003c/span\u003e). After that, grided EWH values for each PC are calculated and then multiplied by the corresponding temporal function, and the total EWH changes at all grids are synthesized by summing the foregoing products of all PCs.\u003c/p\u003e\n\u003ch3\u003e6) Displaying results\u003c/h3\u003e\n\u003cp\u003eTo quickly view inversion results, some figures are plotted for readers. This step displays figures of spatial and temporal functions, maps of annual EWH amplitudes, maps of selected Slepian basis functions, and figures of eigenvalues as the degree changes.\u003c/p\u003e\n\n\n\n\n"},{"header":"3 Application example: GNSS2TWS_Slepian applied to the WUSA region","content":"\u003cp\u003eTo demonstrate GNSS recovery of daily water height changes using the open-source GNSS2TWS_Slepian software, we take the WUSA region as our study area (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea). The WUSA region encompasses the westernmost U.S. states bounded by the Mississippi River and the Pacific Ocean. First, we extract hydrological loading displacements contaminated by other non-hydrological signals and these post-processed vertical position time series are inverted for the estimation of daily EWH values. At last, we evaluate our inversion results in comparison with other independent EWH measurements.\u003c/p\u003e\u003ch2\u003e3.1 Extracting hydrological load displacement\u003c/h2\u003e\u003cp\u003eThe GNSS position time series (Blewitt et al. \u003cspan class=\"CitationRef\"\u003e2018\u003c/span\u003e) processed with the GipsyX-1.0 software are downloaded from the Nevada Geodetic Laboratory (NGL, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://geodesy.unr.edu/\u003c/span\u003e\u003c/span\u003e), University of Nevada, Reno, United States. In total, 332 GNSS stations are well-chosen to infer large-scale TWS changes in the WUSA region (See more details in Text S3). To extract the water-induced vertical surface movements, the well-modeled predictions including non-tidal oceanic loading and atmospheric loading effects are subtracted from the previously downloaded GNSS time series. We then fit each GNSS time series with a constant, a linear, two sinusoids, and some Heavisideand and logarithmic functions when needed (Jiang et al. \u003cspan class=\"CitationRef\"\u003e2022b\u003c/span\u003e). After the removal of constant, linear, offset, and postseismic terms, the residual time series containing sinusoid signals are dominated by seasonal and inter-annual hydrology-induced deformations (Figs. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb–\u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ef), which feature upward motions as water runs off and subsidence when water loads. The comparison between GNSS-observed and PCA-reconstructed time series also indicates the good performance of ALS-based PCA in filling GNSS missing values and data gaps (Figs. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb–\u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ef).\u003c/p\u003e\u003ch2\u003e3.2 Inversion modeling for TWS changes in the WUSA region\u003c/h2\u003e\u003cp\u003eTo characterize the water cycles in the WUSA region, the postprocessing GNSS vertical positions (i.e., residual time series generated in Section 3.1) are then inverted for daily EWH estimates. In the ALS-based PCA decomposition, we consider six PCs when no noticeable enhancement of the fit to the raw data is observed as the increase of selected PCs. These six PCs significantly contribute to the total data variance of the raw data, accounting for 82.5%. Each of these six PCs can explain 73.0%, 11.4%, 6.9%, 3.5%, 2.9%, and 2.4% of the filtered data (Figure \u003cspan class=\"InternalRef\"\u003eS1\u003c/span\u003e). The first PC makes a substantial contribution to the variance compared to the other five PCs, and it reveals conspicuous annual and interannual hydrological loading variations and positive spatial responses at all GNSS stations (Figure \u003cspan class=\"InternalRef\"\u003eS1\u003c/span\u003ea). With the increase of the number of PCs, it manifests gradually weakened and more detailed spatiotemporal features of hydrological loading signals (Figures \u003cspan class=\"InternalRef\"\u003eS1\u003c/span\u003eb–1f).\u003c/p\u003e\u003cp\u003eAfter obtaining the extended boundary data, region-specified Slepian basis functions are generated using SLEPIAN Alpha (Harig et al. \u003cspan class=\"CitationRef\"\u003e2015\u003c/span\u003e). In the inversion model, the bandwidth of the spherical harmonic degrees is chosen as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(L=80\\)\u003c/span\u003e\u003c/span\u003e, corresponding to 6561 orthogonal Slepian basis functions and we only consider Slepian basis functions with eigenvalues of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\gamma }_{\\alpha }\\ge 0.1\\)\u003c/span\u003e\u003c/span\u003e (Figure S2). Figure S3 shows maps of the selected 67 Slepian basis functions. In our study area, the GNSS vertical data at 332 stations are converted into 67 VCD-related Slepian coefficients and then transformed into 67 EWH-related Slepian coefficients using PREM-based load Love numbers generated by Wang et al. (\u003cspan class=\"CitationRef\"\u003e2012\u003c/span\u003e). To avoid the inversion instability of high-order coefficients, a 150-km Gaussian filter is applied in the recovery of water height changes. In case of having sufficient GNSS data, readers can try out higher/lower truncation degrees for more refined/coarser TWS estimated results. Figure S4 indicates relatively fine-scale EWH results when considering a large bandwidth value (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(L=100\\)\u003c/span\u003e\u003c/span\u003e) in contrast to low-resolution water estimates using a lower truncation degree of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(L=60\\)\u003c/span\u003e\u003c/span\u003e. In other words, fine-scale loading signals can be retrieved in a spectral-domain inversion model with high truncation orders and good GNSS spatial coverage, but it is very time-consuming and memory-consuming when calculating Slepian basis functions and corresponding coefficients.\u003c/p\u003e\u003ch2\u003e3.3 Water height changes derived from various datasets\u003c/h2\u003e\u003cp\u003eWe spatially compare annual amplitudes of GNSS-based water estimates with other independent products from GRACE Mascon solutions, North American Land Data Assimilation System (NLDAS), and Snow Data Assimilation System (SNODAS) models (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e). The results indicate that conspicuous seasonal water oscillations with annual amplitudes of \u0026gt; 150 mm are predominantly located along major large-scale mountain systems, including the Cascade Range, the Rocky, and Sierra Nevada Mountains. The NLDAS-EWH (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003ec) features higher resolution than water estimates from GNSS (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003ea) and GRACE (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eb). The seasonal SNODAS-EWH oscillations have strong spatial consistency (a correlation coefficient of 0.94) with that in NLDAS-EWH results and both datasets reveal that the peak regions of annual water amplitudes are clustered along these major mountain systems. It demonstrates that seasonal snowfall is a key driving factor in the regional hydrologic system. Both GNSS and GRACE results underestimate the seasonal water oscillations along these main mountain systems attributed to low spatial resolution. Although we artificially lower the spatial coverage of GNSS stations, the GNSS-inferred results show more detailed seasonal changes than those derived from GRACE Mascon solutions (Figs. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003ea and \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eb).\u003c/p\u003e\u003cp\u003eFigures \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ea–\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ee show time-varying basin-averaged water estimates from multiple datasets (i.e., GNSS, GRACE, NLDAS, and SNODAS). The GNSS-inverted water estimates feature temporal consistency with the GRACE-, NLDAS-EWH products in the Upper Colorado, Great Basin, Pacific Northeast, and California watersheds. In the lower Colorado River basin, extremely low correlations are found between GNSS, GRACE, and NLDAS results, due to the weak signals in seasonal water changes and noises in respective measurements (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eb). The Pacific Northwest hydrologic unit is characterized by the most conspicuous seasonality in water storage variations and the annual amplitudes can reach up to 109, 117, and 138 mm for GNSS-, GRACE- and NLDAS-EWH results (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ed). The California watershed (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ee) has weaker seasonal oscillations (73, 81, and 93 mm for GNSS, GRACE, and NLDAS) than those in the Pacific Northeast basin. Both Upper Colorado and Great Basin watersheds present similar amplitudes of approximately 50 mm (Figs. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ea and \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ec). All water estimates indicate interannual variation characteristics that are caused by extreme hydrometeorological events. To quantify the intensity of these extremes, we convert the GNSS-inferred water estimates into drought severity index (DSI) products (see more details in Text S4). Figures \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ef–\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ej present a comparison between GNSS-DSI and other DSI datasets derived from GRACE and NLDAS in these five river basins. In general, GNSS-DSI can temporally correspond with GRACE-DSI and NLDAS-DSI. In the California, Great, and Pacific Northeast Basins, the GNSS results are simultaneously consistent with GRACE and NLDAS estimates, with correlation coefficients between 0.51–0.85. Good agreement (0.60) is found between GNSS and GRACE results in the Upper Colorado Basin, but a bad correlation (0.23) exists between GNSS and NLDAS (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ef). In the Lower Colorado Basin, a low correlation is seen in the GNSS-GRACE and GNSS-NLDAS results (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eg).\u003c/p\u003e"},{"header":"4 Conclusion and recommendations","content":"\u003cp\u003eThis study releases a numerical calculation tool, named GNSS2TWS_Slepian, for inferring daily large-scale water storage changes using GNSS vertical position time series. This inversion modeling of GNSS data is implemented in the spectral domain based on Slepian basis functions, in contrast to the widely-used spatial-domain Green's function method. To improve the inversion efficiency, a PCA-based inversion method is designed to recover time-varying water height changes. To demonstrate the main functionalities and applications, we investigate the daily TWS changes and characterize hydrological extremes in the Western United States. Through contrastive analysis, the GNSS-inverted water estimates have good spatiotemporal correspondence with other water measurements derived from GRACE and NLDAS. The GNSS characterization result of hydrological extremes illuminates the availability of GNSS for drought monitoring and assessment.\u003c/p\u003e \u003cp\u003eWith the optimization and enhancement of continuously operating GNSS networks worldwide, inversion modeling of crustal deformation measurements is very likely to become a powerful means in the recovery of water height changes for the hydrological community. It is very promising and necessary to explore more application prospects using GNSS for more extensive hydrogeodetic studies, which would contribute to our better understanding of water cycles, hydrological extremes, hydrological dynamics, and hydroclimatic interaction. Therefore, the potential users of GNSS2TWS_Slepian are encouraged to utilize our tools to develop more hydrogeodetic applications and to achieve deeper interdisciplinary collaborations with hydrological communities.\u003c/p\u003e"},{"header":"Statements and Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eMany thanks are sent for the valuable suggestions of the editors and reviewers. The open-source Slepian Alpha software is available from GitHub (https://github.com/csdms-contrib/slepian_alpha). Some figures are created using the Generic Mapping Tools (version 6.4.0) software (Wessel et al. 2019).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by\u0026nbsp;the\u0026nbsp;National Natural Science Foundation of China (Grant No. 41904015).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe\u0026nbsp;GNSS vertical position time series with extra columns are\u0026nbsp;available from the Nevada Geodetic Laboratory, University of Nevada, Reno (http://geodesy.unr.edu/gps_timeseries/tenv3_loadpredictions).\u0026nbsp;The GRACE\u0026nbsp;RL06 Mascon\u0026nbsp;products are downloaded from the CSR\u0026rsquo;s webpage (https://www2.csr.utexas.edu/grace/).\u0026nbsp;The NLDAS model is\u0026nbsp;available from the NASA GES DISC website (https://disc.gsfc.nasa.gov/datasets/NLDAS_NOAH0125_M_002/summary?keywords=NLDAS).\u0026nbsp;The SNODAS products are from the NOAA National Weather Service\u0026rsquo;s National Operational Hydrologic Remote Sensing Center (https://nsidc.org/data/g02158/versions/1).\u0026nbsp;Our inversion code for large-scale water storage changes based on Slepian basis functions is released on GitHub (https://github.com/jzshhh/gnss2tws_slepian).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026apos; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eZhongshan Jiang:\u0026nbsp;\u003c/strong\u003eConceptualization, Methodology, Software, Funding, and Writing- Original draft preparation.\u0026nbsp;\u003cstrong\u003eMiao Tang:\u003c/strong\u003e Writing- Original draft preparation.\u0026nbsp;\u003cstrong\u003eHaiping Wen:\u003c/strong\u003e Writing- Original draft preparation.\u003cstrong\u003e\u0026nbsp;Linguo Yuan:\u0026nbsp;\u003c/strong\u003eSupervision and Writing - Reviewing and Editing.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor Biographies\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eZhongshan Jiang\u0026nbsp;\u003c/strong\u003eis an associate professor at the School of Geospatial Engineering and Science, Sun Yat-Sen University, China. His research interests\u0026nbsp;include seismogeodesy and\u0026nbsp;hydrogeodesy.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eMiao Tang\u003c/strong\u003e is a postgraduate student at the\u0026nbsp;Faculty of Geosciences and Engineering,\u0026nbsp;Southwest Jiaotong University, China. His research interest is hydrogeodesy.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eHaiping Wen\u0026nbsp;\u003c/strong\u003eis a postgraduate student at the\u0026nbsp;Faculty of Geosciences and Engineering,\u0026nbsp;Southwest Jiaotong University, China. His research interest is hydrogeodesy.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eLinguo Yuan\u0026nbsp;\u003c/strong\u003eis a professor at the Faculty of Geosciences and Engineering, Southwest Jiaotong University, China. His research interests include GPS data processing and application and solid Earth tides.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eArgus, D.F., Fu, Y., \u0026amp; Landerer, F.W. (2014). Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion. \u003cem\u003eGeophysical Research Letters, 41\u003c/em\u003e, 1971-1980. doi: 10.1002/2014gl059570\u003c/li\u003e\n\u003cli\u003eArgus, D.F., Landerer, F.W., Wiese, D.N., Martens, H.R., Fu, Y., Famiglietti, J.S., Thomas, B.F., Farr, T.G., Moore, A.W., \u0026amp; Watkins, M.M. (2017). Sustained Water Loss in California\u0026apos;s Mountain Ranges During Severe Drought From 2012 to 2015 Inferred From GPS. \u003cem\u003eJournal of Geophysical Research: Solid Earth, 122\u003c/em\u003e. doi: 10.1002/2017jb014424\u003c/li\u003e\n\u003cli\u003eArgus, D.F., Martens, H.R., Borsa, A.A., Knappe, E., Wiese, D.N., Alam, S., Anderson, M., Khatiwada, A., Lau, N., Peidou, A., Swarr, M., White, A.M., Bos, M.S., Ellmer, M., Landerer, F.W., \u0026amp; Gardiner, W.P. (2022). Subsurface Water Flux in California\u0026apos;s Central Valley and Its Source Watershed From Space Geodesy. \u003cem\u003eGeophysical Research Letters, 49\u003c/em\u003e. doi: 10.1029/2022gl099583\u003c/li\u003e\n\u003cli\u003eBlewitt, G., Hammond, W., \u0026amp; Kreemer, C. (2018). Harnessing the GPS Data Explosion for Interdisciplinary Science. \u003cem\u003eEos, 99\u003c/em\u003e. doi: 10.1029/2018eo104623\u003c/li\u003e\n\u003cli\u003eBlewitt, G., Lavallee, D., Clarke, P., \u0026amp; Nurutdinov, K. (2001). A new global mode of Earth deformation: seasonal cycle detected. \u003cem\u003eScience, 294\u003c/em\u003e, 2342-2345. doi: 10.1126/science.1065328\u003c/li\u003e\n\u003cli\u003eChen, J., Cazenave, A., Dahle, C., Llovel, W., Panet, I., Pfeffer, J., \u0026amp; Moreira, L. (2022). Applications and Challenges of GRACE and GRACE Follow-On Satellite Gravimetry. \u003cem\u003eSurv Geophys, 43\u003c/em\u003e, 305-345. doi: 10.1007/s10712-021-09685-x\u003c/li\u003e\n\u003cli\u003eFarrell, W.E. (1972). Deformation of the Earth by surface loads. \u003cem\u003eReviews of Geophysics, 10\u003c/em\u003e. doi: 10.1029/RG010i003p00761\u003c/li\u003e\n\u003cli\u003eFok, H.S., \u0026amp; Liu, Y. (2019). An Improved GPS-Inferred Seasonal Terrestrial Water Storage Using Terrain-Corrected Vertical Crustal Displacements Constrained by GRACE. \u003cem\u003eRemote Sensing, 11\u003c/em\u003e. doi: 10.3390/rs11121433\u003c/li\u003e\n\u003cli\u003eFu, Y., Argus, D.F., \u0026amp; Landerer, F.W. (2015). GPS as an independent measurement to estimate terrestrial water storage variations in Washington and Oregon. \u003cem\u003eJournal of Geophysical Research: Solid Earth, 120\u003c/em\u003e, 552-566. doi: 10.1002/2014jb011415\u003c/li\u003e\n\u003cli\u003eHan, S.C., \u0026amp; Razeghi, S.M. (2017). GPS Recovery of Daily Hydrologic and Atmospheric Mass Variation: A Methodology and Results From the Australian Continent. \u003cem\u003eJournal of Geophysical Research: Solid Earth, 122\u003c/em\u003e, 9328-9343. doi: 10.1002/2017jb014603\u003c/li\u003e\n\u003cli\u003eHarig, C., Lewis, K.W., Plattner, A., \u0026amp; Simons, F.J. (2015). A Suite of Software Analyzes Data on the Sphere. \u003cem\u003eEos, 96\u003c/em\u003e. doi: 10.1029/2015eo025851\u003c/li\u003e\n\u003cli\u003eHarig, C., \u0026amp; Simons, F.J. (2012). Mapping Greenland\u0026apos;s mass loss in space and time. \u003cem\u003eProc Natl Acad Sci U S A, 109\u003c/em\u003e, 19934-19937. doi: 10.1073/pnas.1206785109\u003c/li\u003e\n\u003cli\u003eHsu, Y.-J., Fu, Y., B\u0026uuml;rgmann, R., Hsu, S.-Y., Lin, C.-C., Tang, C.-H., \u0026amp; Wu, Y.-M. (2020). Assessing seasonal and interannual water storage variations in Taiwan using geodetic and hydrological data. \u003cem\u003eEarth and Planetary Science Letters, 550\u003c/em\u003e. doi: 10.1016/j.epsl.2020.116532\u003c/li\u003e\n\u003cli\u003eJiang, Z., Hsu, Y.-J., Yuan, L., Cheng, S., Li, Q., \u0026amp; Li, M. (2021a). Estimation of daily hydrological mass changes using continuous GNSS measurements in mainland China. \u003cem\u003eJournal of Hydrology, 598\u003c/em\u003e. doi: 10.1016/j.jhydrol.2021.126349\u003c/li\u003e\n\u003cli\u003eJiang, Z., Hsu, Y.-J., Yuan, L., Feng, W., Yang, X., \u0026amp; Tang, M. (2022a). GNSS2TWS: an open-source MATLAB-based tool for inferring daily terrestrial water storage changes using GNSS vertical data. \u003cem\u003eGPS Solutions, 26\u003c/em\u003e. doi: 10.1007/s10291-022-01301-8\u003c/li\u003e\n\u003cli\u003eJiang, Z., Hsu, Y.-J., Yuan, L., \u0026amp; Huang, D. (2021b). Monitoring time-varying terrestrial water storage changes using daily GNSS measurements in Yunnan, southwest China. \u003cem\u003eRemote Sensing of Environment, 254\u003c/em\u003e. doi: 10.1016/j.rse.2020.112249\u003c/li\u003e\n\u003cli\u003eJiang, Z., Hsu, Y.J., Yuan, L., Tang, M., Yang, X., \u0026amp; Yang, X. (2022b). Hydrological drought characterization based on GNSS imaging of vertical crustal deformation across the contiguous United States. \u003cem\u003eSci Total Environ, 823\u003c/em\u003e, 153663. doi: 10.1016/j.scitotenv.2022.153663\u003c/li\u003e\n\u003cli\u003eJiang, Z., Hsu, Y.J., Yuan, L., Yang, X., Ding, Y., Tang, M., \u0026amp; Chen, C. (2021c). Characterizing Spatiotemporal Patterns of Terrestrial Water Storage Variations Using GNSS Vertical Data in Sichuan, China. \u003cem\u003eJournal of Geophysical Research: Solid Earth, 126\u003c/em\u003e. doi: 10.1029/2021jb022398\u003c/li\u003e\n\u003cli\u003eJin, S., \u0026amp; Zhang, T. (2016). Terrestrial Water Storage Anomalies Associated with Drought in Southwestern USA from GPS Observations. \u003cem\u003eSurveys in Geophysics, 37\u003c/em\u003e, 1139-1156. doi: 10.1007/s10712-016-9385-z\u003c/li\u003e\n\u003cli\u003eKnappe, E., Bendick, R., Martens, H.R., Argus, D.F., \u0026amp; Gardner, W.P. (2019). Downscaling Vertical GPS Observations to Derive Watershed‐Scale Hydrologic Loading in the Northern Rockies. \u003cem\u003eWater Resources Research, 55\u003c/em\u003e, 391-401. doi: 10.1029/2018wr023289\u003c/li\u003e\n\u003cli\u003eKusche, J., \u0026amp; Schrama, E.J.O. (2005). Surface mass redistribution inversion from global GPS deformation and Gravity Recovery and Climate Experiment (GRACE) gravity data. \u003cem\u003eJournal of Geophysical Research: Solid Earth, 110\u003c/em\u003e. doi: 10.1029/2004jb003556\u003c/li\u003e\n\u003cli\u003eLi, X., Zhong, B., Li, J., \u0026amp; Liu, R. (2023). Inversion of GNSS Vertical Displacements for Terrestrial Water Storage Changes Using Slepian Basis Functions. \u003cem\u003eEarth and Space Science, 10\u003c/em\u003e. doi: 10.1029/2022ea002608\u003c/li\u003e\n\u003cli\u003eMilliner, C., Materna, K., Burgmann, R., Fu, Y., Moore, A.W., Bekaert, D., Adhikari, S., \u0026amp; Argus, D.F. (2018). Tracking the weight of Hurricane Harvey\u0026apos;s stormwater using GPS data. \u003cem\u003eSci Adv, 4\u003c/em\u003e, eaau2477. doi: 10.1126/sciadv.aau2477\u003c/li\u003e\n\u003cli\u003eSimons, F.J., Dahlen, F.A., \u0026amp; Wieczorek, M.A. (2006). Spatiospectral Concentration on a Sphere. \u003cem\u003eSIAM Review, 48\u003c/em\u003e, 504-536. doi: 10.1137/s0036144504445765\u003c/li\u003e\n\u003cli\u003eWang, H., Xiang, L., Jia, L., Jiang, L., Wang, Z., Hu, B., \u0026amp; Gao, P. (2012). Load Love numbers and Green\u0026apos;s functions for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0. \u003cem\u003eComputers \u0026amp; Geosciences, 49\u003c/em\u003e, 190-199. doi: 10.1016/j.cageo.2012.06.022\u003c/li\u003e\n\u003cli\u003eWessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F., \u0026amp; Tian, D. (2019). The Generic Mapping Tools Version 6. \u003cem\u003eGeochemistry, Geophysics, Geosystems, 20\u003c/em\u003e, 5556-5564. doi: 10.1029/2019gc008515\u003c/li\u003e\n\u003cli\u003eWhite, A.M., Gardner, W.P., Borsa, A.A., Argus, D.F., \u0026amp; Martens, H.R. (2022). A Review of GNSS/GPS in Hydrogeodesy: Hydrologic Loading Applications and Their Implications for Water Resource Research. \u003cem\u003eWater Resour Res, 58\u003c/em\u003e, e2022WR032078. doi: 10.1029/2022WR032078\u003c/li\u003e\n\u003cli\u003eYang, X., Tian, S., You, W., \u0026amp; Jiang, Z. (2021). Reconstruction of continuous GRACE/GRACE-FO terrestrial water storage anomalies based on time series decomposition. \u003cem\u003eJournal of Hydrology, 603\u003c/em\u003e. doi: 10.1016/j.jhydrol.2021.127018\u003c/li\u003e\n\u003cli\u003eZhong, B., Li, X., Chen, J., Li, Q., \u0026amp; Liu, T. (2020). Surface Mass Variations from GPS and GRACE/GFO: A Case Study in Southwest China. \u003cem\u003eRemote Sensing, 12\u003c/em\u003e. doi: 10.3390/rs12111835\u003c/li\u003e\n\u003cli\u003eZhu, H., Chen, K., Hu, S., Liu, J., Shi, H., Wei, G., Chai, H., Li, J., \u0026amp; Wang, T. (2023). Using the Global Navigation Satellite System and Precipitation Data to Establish the Propagation Characteristics of Meteorological and Hydrological Drought in Yunnan, China. \u003cem\u003eWater Resources Research, 59\u003c/em\u003e. doi: 10.1029/2022wr033126\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"GNSS, Slepian basis functions, Terrestrial water storage, Vertical motions, Hydrological extremes","lastPublishedDoi":"10.21203/rs.3.rs-4678987/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4678987/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eChanges in terrestrial water storage (TWS) can deform the Earth\u0026rsquo;s solid surface in the form of geodetically measurable vertical motions. Here, a new open-source Matlab software, named GNSS2TWS_Slepian, is developed to achieve the recovery of daily TWS changes from Global Navigation Satellite System (GNSS) crustal vertical positions. Differing from the widely-used spatial-domain inversion strategy based on Green's function method, our inversion modeling is implemented in the spectral domain based on Slepian basis functions, which aims to infer daily large-scale TWS changes using non-uniformly distributed GNSS vertical data. GNSS2TWS_Slepian is designed with different structured modules and the logic of the program workflow can be easily followed. To obtain daily estimates of TWS changes, the principal component analysis is integrated into our time-varying inversion model. To demonstrate the main functionalities, equivalent water height changes are investigated in the Western United States. This study aims to provide a scientific mathematical tool for resolving large-scale water mass loads, which is instrumental in broadening the applications of GNSS in hydrology.\u003c/p\u003e","manuscriptTitle":"GNSS2TWS_Slepian: A software to recover daily GNSS-inverted terrestrial water storage changes based on Slepian basis functions","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-08-09 05:42:40","doi":"10.21203/rs.3.rs-4678987/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"0f3896b4-6453-436c-83fa-f5eee45026e4","owner":[],"postedDate":"August 9th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-08-09T05:42:40+00:00","versionOfRecord":[],"versionCreatedAt":"2024-08-09 05:42:40","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4678987","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4678987","identity":"rs-4678987","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}