Introduction

Terrestrial water storage (TWS) is closely related to human survival and the productivity of all living things. Therefore, changes in TWS have profound impacts on human social development and economic activities. Monitoring TWS changes is significant for water resource management and can also help us understand patterns of climate change. Several geodetic techniques and hydrology models have been used to explore the change in TWS. For example, the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GFO) missions track changes in TWS by detecting variations in gravity caused by mass redistribution (Tapley et al. 2004). However, the GRACE technology has limited spatiotemporal resolution, and conventional analysis based on the global time-variable gravity field using Stokes spherical harmonic coefficients contains significant uncertainties. Despite the apparent improvement in resolution provided by Mascon (Save 2019, 2020), different modeling approaches can lead to markedly different TWS changes. Global hydrological models (for example, GLDAS and WGHM) that simulate changes in TWS are mainly based on global land data assimilation analysis and have considerable uncertainties, making them difficult to use for research related to environmental change and disaster prediction (Rodell et al. 2004; Jin and Zhang 2016; Scanlon et al. 2018; Trautmann et al. 2018).

Given the technical difficulties involved in monitoring TWS for small-scale water bodies with GRACE/GFO, another geodetic measurement technology for monitoring TWS changes, the Global Navigation Satellite System (GNSS), has emerged (van Dam 2001). The principle behind GNSS monitoring of TWS is different from that of GRACE/GFO (For the sake of convenience in writing, unless otherwise specified, the term "GRACE" in the following text refers to "GRACE/GFO"). While GRACE satellite gravity observations can directly detect TWS (mass) changes, GNSS can only indirectly constrain TWS changes (Bock and Melgar 2016) through the observation of surface (geometric) deformation responses caused by crustal elastic effects (Farrell 1972). Numerous studies have shown that TWS monitoring via GNSS is broadly consistent with GRACE results (Amos et al. 2014; Borsa et al. 2014). If the GNSS network is dense enough, it can provide images of TWS changes with higher temporal and spatial resolution (Fu et al. 2015; Milliner et al. 2018). In the past decade, GNSS is now widely used to investigate TWS changes and drought and flood disasters (Johnson et al. 2017; Milliner et al. 2018; Hsu et al. 2020; Carlson et al. 2022b, a). This fully demonstrates that while GRACE can recognize the overall changes in TWS within a region, GNSS is better equipped to reflect the fine features of TWS variability within the same region, demonstrating a complementary relationship with the overall results obtained by GRACE. Therefore, some researchers jointly inverted GNSS and GRACE for TWS changes. Some of them used GRACE to generate "virtual" GNSS vertical displacement data and inverted TWS together with GNSS measurements (Liu et al. 2019; Fok and Liu. 2019; Zhu et al. 2023a), and the other part combined GNSS and GRACE into a Green's function matrix to regularize the inversion of TWS (Adusumilli et al. 2019; Razeghi et al. 2019; Carlson et al. 2022b, a; Li et al. 2023c). Yang et al. (2023) compared the joint inversion methods and considered that the former was more effective.

Droughts and floods are both related to drastic changes in TWS. Some scholars have utilized GNSS to investigate drought events in Yunnan due to the larger annual amplitude and dense distribution of stations and the significant multiyear seasonal fluctuations resulting from periodic droughts. For example, He (2017) systematically studied TWS inferred from continuous GNSS in Yunnan by using GF inversion method. Zhang et al. (2016) used Helmert variance component estimation (HVCE) and minimum mean square error (MMSE) to infer the annual equivalent water height (EWH) changes in Yunnan. Jiang et al. (2017) analyzed the monitoring ability of 26 continuous GPS stations in Yunnan Province from October 2007 to October 2015 and showed that the characteristics of regional monsoon and drought can be understood more comprehensively by combining precipitable water vapor (PWV) and vertical crustal displacements (VCD) derived from GPS with precipitation and EWH derived from GRACE. Fok and Liu (2019) speculated that the GPS-inferred TWS has larger amplitudes than those inferred from GRACE, which should be partly due to irregularly distributed GPS stations and the neglect of the terrain effect, and proposed a new method for inferring TWS by GPS. Cheng et al. (2021) inversed the TWS changes in Sichuan, Yunnan and Chongqing based on SBF method. Jiang et al. (2021a) studied the TWS changes and extreme hydrometeorology (such as heavy precipitation and drought) in Yunnan using the inversion independent component analysis method. They found that there was good consistency between the GPS-inferred TWS and water changes based on GLDAS during the winter rainstorm in 2015. Li et al. (2023a) used the SBF method to recover the TWS changes inferred from the VCD of GNSS in southwest China from 2011 to 2020 and found that the signal amplitude of SBF method was stronger than that of GF method. Li et al. (2023b) used the VCD of GNSS in China from December 2010 to February 2021 to infer TWS changes. They introduced a priori constraint matrix and estimated regularization parameters by iterative least square method. Zhu et al. (2023a) used variational Bayesian independent component analysis to reconstruct the three-dimensional crustal deformation caused by water load at 43 GNSS stations in Yunnan, and combined with the VCD inferred by GRACE, the TWS estimates of three basins in Yunnan were obtained. Zhu et al. (2023b) systematically studied the drought propagation in Yunnan by using the GNSS observation and precipitation data and found that hydrological drought was more serious than meteorological drought, and the drought propagation time was 2–7 months. Yang et al. (2023) studied the TWS changes in southwest China by integrating GNSS and GRACE/GFO observations and found that the hydrological drought events from May 2019 to December 2020 were the most serious. However, studies have revealed discrepancies between GNSS- and GRACE-derived TWS estimates in Yunnan, particularly in terms of amplitude (Li et al. 2023a). These discrepancies may arise from differences in spatial resolution and sensitivities between the two methods, inaccuracies in GNSS station time series, or the leakage problem associated with the classical processing of spherical harmonics (SH) in GRACE.

In this study, we focus on Yunnan, in southwest China, which has recently undergone several severe droughts (e.g., drought events in 2006, 2009–2010, 2011, and 2019) and floods (such as the heavy precipitation events in 2012, 2015, 2017, and 2018). In addition to the GF approach, we also apply the Slepian basis function (SBF) method to invert GNSS vertical displacement to TWS in six river basins and compare the differences in this paper. We adopt three different approaches to derive TWS from GRACE and find that the GRACE TWS derived by SBF has a similar pattern to the GNSS TWS estimates in amplitude. However, in previous studies, the amplitude of TWS, as inverted by GNSS methods, is much higher than that inverted by GRACE spherical harmonic coefficients or GRACE mascon solutions (Jiang et al. 2021a, b). This paper implements a combined analysis using GNSS, GRACE, and GLDAS to analyze the spatial and temporal characteristics of TWS and successfully detect major droughts in Yunnan.

Data and method

We choose regions with significant seasonal and interannual variations in TWS as the subjects of our study. Moreover, the study area is equipped with dense continuous GNSS stations. The data used in this research are publicly available and easily accessible, facilitating the reproducibility of the experiments and validation of the findings.

Study area

Yunnan, located in southwestern China, is bordered by the South China Sea to the southeast, the Bay of Bengal to the southwest, and the Qinghai–Tibetan Plateau to the northwest. The region experiences a typical monsoon climate, which is affected by periodic changes in atmospheric circulation and ocean currents, such as El Niño and La Niña. This results in a wet season from May to October, followed by a dry season from November to April. The precipitation distribution is uneven, with May receiving 104 mm of precipitation, almost double the dry season's average (Liu et al. 2020). Six main water networks, including the Yangtze River, Pearl River, Red River, Lancang River, Nujiang River, and Irrawaddy River, flow through Yunnan, ultimately ending up in the Pacific Ocean or Indian Ocean.

Data and preprocessing

In the study area, we gathered a substantial amount of data, including GNSS, gravity, hydrological models, and rainfall. This chapter outlines the characteristics and processing methods for each type of data. For more detailed procedures, please refer to Appendix A.

GNSS vertical positioning time series

We collect daily GNSS measurements from 46 CMONOC stations and 8 CMA stations; Fig. 1 displays the locations of all the stations. Most of the stations are observed continuously from 2010 to 2022. We obtained the CMONOC GNSS station data from the China Earthquake Data Center and processed the GNSS data using the GAMIT/GLOBK (version 10.7) software developed by Herring et al. (2018). In order to obtain precise TWS in the Yunnan region, we meticulously processed the observation data from all GNSS stations (see more details in Appendix A) and acquired time series of vertical surface displacement caused by water load changes.

Fig. 1
figure 1

Distribution of GNSS stations and the division of six major river basins in the Yunnan region. (left) The GNSS station distribution in the study area. The red rectangle in the bottom right corner indicates the location of Yunnan. Red circles and green circles indicate the Crustal Movement Observation Network of China (CMONOC) and China Meteorological Administration (CMA) GNSS stations, respectively. Blue lines represent the primary rivers in Yunnan. (right) A map of the six primary river basins: (1) Yangtze River basin, (2) Pearl River basin, (3) Red River basin, (4) Lancing River basin, (5) Nujiang River basin, and (6) Irrawaddy River basin, as defined by HydroSHEDS (https://www.hydrosheds.org/products/hydrobasins)

Full size image

GRACE data products

Changing mass loads can cause displacements, which can be expressed in terms of spherical harmonic coefficients (SHc) for the gravity field and load Love numbers (Kusche and Schrama 2005; van Dam et al. 2007). GRACE observations are mainly given in the SHc form (GRACE solutions can be downloaded from NASA PODAAC https://podaac.jpl.nasa.gov/); thus, users need to carry out a series of preprocessing steps to obtain specific physical quantities. We use an open-source GRACE MATLAB Toolbox (GRAMAT) (https://github.com/fengweiigg/GRACE_Matlab_Toolbox) to replace C20 and Degree-1, including removing the average, de-striping, obtaining the mass coefficients, and removing GIA effects with 400 km Gaussian filtering.

To confirm the variations in the EWH inversion results, we employ GRACE RL06 Mascon products provided by the Center for Space Research (CSR) (http://www2.csr.utexas.edu/grace/RL06_mascons.html), which are defined according to 0.25° rectangle grids and cover the period from January 2010 to December 2021 and include 31 months with missing data. The recently published Mascon solutions are utilized, and all well-characterized corrections, such as the Degree-1 coefficients, C20 components, oceanic and atmospheric terms, and Glacial Isostatic Adjustment (GIA) effects, are applied (Save 2019; Save 2020). Since no further empirical de-striping, smoothing, or filtering is needed for these data, we refrained from applying additional filtering or scaling factors to these CSR Mascon products.

GLDAS data product

This study utilizes a product from the Noah 2.7.1 land hydrology model from the Global Land Data Assimilation System (GLDAS). The GLDAS offers high-quality, global land surface fields and supports various water cycle investigations and water resource applications (Rodell et al. 2004). The input data are available at the Goddard Space Flight Center Distributed Active Archive Center (DISC) (http://disc.sci.gsfc.nasa.gov/hydrology/data-holdings). The integrated total water content datasets used in this investigation are obtained by summing the layers, including soil moisture at depths of 0–2 m, snow water equivalent, and canopy interception storage, and are represented as 0.25° × 0.25° rectangular grids spanning from 60°N to 90°S and 180°W to 180°E. It is important to note that groundwater variations below the soil layer are not considered in GLDAS.

Surface precipitation data

Precipitation plays a crucial meteorological role influencing water circulation and is closely linked to flood and drought events. In this study, precipitation data provided by the CMA (http://data.cma.cn/) are utilized to evaluate the coherence between precipitation and TWS changes. Daily and monthly resolved precipitation products from January 2010 to January 2022 are included, with these datasets are converted into 0.5° × 0.5° rectangular grids using records from 2,472 stations in China, 109 of which are located in Yunnan and are equally distributed. To calculate surface oscillation resulting from the thermoelastic effect, we use the skin temperature data from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis of the global climate (ERA5). These products have a horizontal spatial resolution of 0.25° × 0.25°, and they are downloaded from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/cdsapp#!/home).

Methodology

Both the GF method (GNSS-Green) (Argus et al. 2014, 2017; Fu et al. 2015) and the SBF method (GNSS-Slepian) (Simons and Dahlen 2006; Han and Razeghi 2017) can be applied for the inversion of EWH using GNSS measurements. The mathematical model for calculating EWH is elaborated in detail in the attached materials (Appendix B).

Green’s function

Green's function (GF) is a mathematical and physical equation (Farrell 1972), which represents the relationship between a specific "field" and the "source" that produces this field. For the inversion in the field of geophysics and geodesy, the primary work of linear inversion is to establish a linear function model (linear observation equations) (Wahr et al. 2013). In the study of GNSS hydrology, GF is the contribution of surface water unit load to the deformation of GNSS stations. The mathematical model of linear inversion is expressed as follows (Borsa et al. 2014).

$$\left\{\begin{array}{c}Ax=u\\ \beta Tx=0\\ {D}_{u}={\sigma }_{0}^{2}Q\end{array}\right.$$

where \(A\) is GF matrix, \(x\) is the parameter vector to be estimated (EWH); \(u\) is the observed value (GNSS deformation field); \(\beta\) is the smoothing factor of the parameter column vector; \(T\) is the smoothing constraint matrix of parameter vector (such as Laplace second-order difference smoothing operator); \({D}_{u}\) is the variance–covariance matrix of the observation column vector; \({\sigma }_{0}^{2}\) is the error in the unit weight; and \(Q\) is the cofactor matrix of the observation column vector.

According to the Tikhonov regularization principle, the TWS estimation criteria for GNSS inferred are as follows (van Dam 2001; Fu et al. 2015; He 2017; Jiang et al. 2021a, b).

$${\left(\left(Ax-u\right)/\sigma \right)}^{2}+{\beta }^{2}{\left(T\right)}^{2}\to {\text{min}}$$

Slepian basis function

The inversion results of GF are greatly influenced by the number and spatial distribution of GNSS stations. When the research area is large, the grid area is too large, which will lead to poor inversion results. Therefore, GRACE technology is mainly used in the inversion research of large-scale TWS. However, the spherical harmonic function satisfies the orthogonality in the global scope, and there will be "leakage" in the local research area (Han and Razeghi 2017). The spherical SBF can minimize the "leakage" of signals from the research area. Therefore, researchers began to study the use of SBF to invert GRACE and GNSS data to obtain large-scale EWH.

The surface water load \({\widehat{\sigma }}_{L}\left(\theta ,\lambda ,t\right)\) in the study area can be expressed as follows (Simons and Dahlen 2006; Han and Razeghi 2017; Jiang et al. 2022; Li et al. 2023a).

$${\widehat{\sigma }}_{L}\left(\theta ,\lambda ,t\right)=a\sum_{\beta =1}^{J}{\gamma }_{\beta }{s}_{\beta }^{{\text{EWH}}}\left(t\right){g}_{\beta }\left(\theta ,\lambda \right)$$

where \(a\) is radius of the earth, \(J\) is truncation order, \({\gamma }_{\beta }\) is eigenvalue (also called aggregation factor), \({s}_{\beta }^{{\text{EWH}}}\left(t\right)\) is SBF coefficient, and \({g}_{\beta }\left(\theta ,\lambda \right)\) is SBF.

Six methods for inferring EWH

We compare six different estimates of EWH time series (Fig. 2): GNSS-inferred EWH based on the SBF method (GNSS_Slepian) or the GF method (GNSS_Green), GRACE-inferred EWH based on the SH coefficient method (GRACE_SHc), or the SBF method (GRACE_Slepian), GRACE-Mascons (GRACE_Mascon), and the GLDAS NOAH hydrological model (GLDAS_NOAH). EWH can be obtained directly from GLDAS_NOAH and GRACE_Mascon products, offering users great convenience. On the other hand, GRACE_Slepian and GNSS_Slepian employ a spherical window spherical harmonic function (local orthogonal basis) to expand water storage anomalies, with the key task of determining Slepian basis function coefficients. Similarly, GRACE_SHc utilizes the load Love numbers to transform time-variable gravity spherical harmonic coefficients into water storage anomaly spherical harmonic coefficients. As to GNSS_Green, it infers water storage anomalies by leveraging the water load–deformation response relationship (Green’s function). Each of these six methods possesses distinct characteristics, with GRACE primarily focusing on identifying the overall changes in TWS within a given region (resolution ~ 400km), while the GNSS captures finer features of TWS variations within the region (resolution better than 50km) (Argus et al. 2014; Knappe et al. 2019). The GLDAS typically captures shallow groundwater while ignoring surface water (rivers, lakes, and reservoirs) and deeper groundwater sources (Scanlon et al. 2018).

Fig. 2
figure 2

Data processing procedure and data inversion strategy. The numerical values enclosed in circles in the diagram represent the following meanings: ① GNSS_Green; ② GNSS_Slepian; ③ GRACE_Slepian; ④ GRACE_SHc; ⑤ GRACE_Mascon; ⑥ GLDAS_NOAH

Full size image

Results

We conduct a statistical analysis of TWS inferred from different datasets and compare and analyze the spatial–temporal distribution patterns and variations in water storage quantities (correlation). The computed results in this chapter serve as the basis for subsequent analysis and discussion.

Analysis of spatial coherence from various water datasets

The GNSS-inferred EWH is compared with the predictions from the GRACE products and the hydrology assimilation model. The spatial variations in TWS inferred from GNSS, GRACE, and GLDAS are found to be similar (Fig. 3). The amplitude of seasonal TWS changes decreased quickly from southwest to northeast Yunnan, and the seasonal water fluctuations in the southwest are mainly attributable to the South Asian monsoon (Jiang et al. 2017; Cao et al. 2017). The largest amplitude in GNSS_Slepian is approximately 287 mm (with an average of 211 mm), which is located in southwestern and western Yunnan (Fig. 3d). GRACE_SHc exhibits an average of 84 mm and a maximum EWH of 159 mm. Although the largest GRACE_SHc-inferred water fluctuations are approximately 55.4% of the GNSS_Slepian results, the peak variations of both methods are located near the drainage basins of the Nujiang River and Lancang River (Fig. 3a and b). GRACE_Mascon and GRACE_Slepian consistently exhibit hydrological spatial distribution characteristics (Fig. 3c and f), and the amplitude is larger than GRACE_SHc, potentially due to signal leakage. The GLDAS hydrologic model (Fig. 3e) exhibits average and largest seasonal water load changes of 94 mm and 165 mm, respectively, smaller than those of GNSS and GRACE. The GNSS- and GRACE-inferred EWH variations integrate all surface and subsurface water resources, whereas the primary water component in the hydrological model is soil moisture. GNSS-EWH has a slightly higher spatial resolution than GRACE-EWH due to the density of the GNSS network in Yunnan, and the GNSS inversion produces similar spatial patterns to those inferred from other models.

Fig. 3
figure 3

Annual amplitudes of water height variations are from a GNSS_Green, b GRACE_SHc, c GRACE_Mascon, d GNSS_Slepian, e GLDAS_NOAH, and f GRACE_Slepian. Black lines are primary river basins in Yunnan

Full size image

We also discover discrepancies in the results obtained from the two methods of utilizing GNSS-inferred equivalent water heights. Specifically, the statistical analysis shows that the amplitude of the SBF inversion results is on average 25% larger than that of the GF inversion results (Li et al. 2023a). We believe that the GF inversion results are more reliable at small scales, while the SBF may suffer from insufficient spatial resolution (Appendix C). Hence, we make corrections (\(\omega\)= 0.705) to the GNSS_Slepian results, elevating the amplitude consistency between GNSS_Green and GNSS_Slepian from 0.48 to 0.88. For more details, please read the support information Appendix C and Fig. C6.

Investigation of temporal characteristics of different water estimates

We utilize multisource time-varying hydrological datasets to compare the temporal coherence of seasonal TWS across Yunnan. Figures 4 and 5 present the monthly time series of area-averaged precipitation and water storage changes, as derived from the hydrologic assimilation model (GLDAS), GRACE products, and GNSS-based solutions. Generally, precipitation attains its annual peak first, followed within 1–2 months by the seasonal maximum derived from the GLDAS-, GNSS-, and GRACE-inferred water solutions (Fig. 4 bottom panel). This time lag between precipitation and water estimates is consistent with observations made in other tropical and subtropical regions, such as Central Africa, the Amazon Basin, and Australia (Humphrey et al. 2016). The phase shift between GRACE-derived water solutions and precipitation indicates that TWS accumulates during the rainy seasons and reaches its largest value close to the end of the wet seasons (Birhanu and Bendick 2015).

Fig. 4
figure 4

Time series of EWH and precipitation. GNSS_Green (orange lines), GNSS_Slepian (green line), GRACE_SHc (red line), GRACE_Mascon (blue lines), GRACE_Slepian (black lines), GLDAS_NOAH (magenta line). The bottom panel focuses on the period of 2015−2016 from the top panel. Green bars represent the monthly precipitation

Full size image
Fig. 5
figure 5

Phase of Equivalent Water Height Time Series Inversed by Different Data and Methods. From a−f are GNSS_Green, GRACE_SHc, GRACE_Mascon, GNSS_Slepian, GLDAS_NOAH, GRACE_Slepian solutions, respectively

Full size image

We investigate the time lags and correlation coefficients between seasonal TWS solutions derived from different datasets (Table 1 and Fig. 5). It is worth noting that the phase of GNSS results lags behind rainfall by almost two months. The spatial distribution of GRACE_SHc, GRACE_Mascon, and GLDAS_NOAH aligns with the GNSS results, indicating that the northwestern to southeastern regions of Yunnan have smaller phase values compared to the central and southwestern regions. While the GRACE_SHc and the GLDAS_NOAH exhibit consistent phase values, they are smaller than those of the GNSS results, whereas the GRACE_Mascon shows larger phase values. Additionally, the GRACE_Slepian stands out as it exhibits a distinct phase maximum in northeastern Yunnan, suggesting potential outliers and notable differences compared to the other results. On the other hand, the cross-correlation coefficients between different pairs of data ranged from 0.37 to 0.98. Weak or moderate correlations (0.37–0.59) are observed between the GRACE-Slepian solution and other datasets, while significant correlations (0.65–0.98) are found between the GNSS, GLDAS, and other GRACE-based water estimates. The anomaly curve of land water reserves inferred by GNSS in six major basins in Yunnan (Figs. 1a and 6) generally correspond with the changes in the total amount of surface water resources in Yunnan, with the exception of possible abnormalities in the Pearl River Basin due to the lack of GNSS stations. Notably, Fig. 6 also reveals a 32.7% decrease in average rainfall in 2019, as reflected in the GNSS-inferred results.

Table 1 Correlation coefficients estimated between these regional monthly average solutions
Full size table
Fig. 6
figure 6

Time series of regional water storage changes around some typical GNSS stations. The orange, green, red, blue, black, and magenta lines are the monthly GNSS_Green, GNSS_Slepian, GRACE_SHc, GRACE_Mascon, GRACE_Slepian, and GLDAS_NOAH time series, respectively. Green bars illustrate the monthly precipitation. Red shaded bands indicate the drought periods. A drought is determined by a water deficit lasting for three or more continuous months

Full size image

In addition, we utilize GNSS-inferred equivalent water heights to investigate annual TWSA in six sub-basins of the Yunnan region from 2010 to 2022 (Table C1) and obtain a total annual TWSA of 72.6 ± 6.0 Gt for Yunnan province (Appendix C).

Discussion

The occurrence of drought is positively correlated with the decline in TWS. Variations in water load, whether increasing or decreasing, can induce surface deformation, whose response can be accurately captured by GNSS. Therefore, by investigating the vertical deformation observed by GNSS, it is possible to infer the magnitude of TWS and consequently reflect the severity of drought in the study area.

Characteristics of typical droughts

The distribution of precipitation in Yunnan is highly variable across seasons and regions. The rainy season lasts from May to October and accounts for 85% of the annual rainfall, while the dry season lasts from November to April of the following year and only accounts for 15% of the annual precipitation (Li et al. 2019). Although most areas receive over 1000 mm of annual precipitation, there is considerable variation in the regional precipitation distribution throughout Yunnan. The maximum annual precipitation ranges from 2200–2700 mm, while the minimum recorded precipitation is only 584 mm. Central Yunnan (Fig. 6d, e, f) is a severely arid area situated in the watersheds of the Yangtze River, Pearl River, Lancang River, and Honghe River. It receives less precipitation, and the average annual runoff accounts for only 12% of Yunnan's water resources (Liu et al. 2020; Jiang et al. 2017).

Seasonal hydrometeorological extremes, particularly during winter and spring, are frequent in Yunnan and can be observed by examining the amplitudes of seasonal water storage. For example, the winter and spring droughts in 2010 are just as severe as those in 2020. Meteorological records indicate that similar droughts have occurred periodically. In this study, we investigate the relationship between several well-documented hydrometeorological extremes and changes in TWS. Figure 6 illustrates that the GNSS, GRACE, and GLDAS estimates show lower-than-average water fluctuations during the summer of 2011, indicating significant TWS deficits resulting from insufficient precipitation. Likewise, in 2019, persistent precipitation deficiencies led to a prolonged extreme drought, which is evident from the water deficits shown by GNSS (Fig. 6). During heavy precipitation events such as the rainy summers of 2012 and 2017, the annual peak is larger in GNSS, GRACE, and GLDAS inferred EWH than during other time periods.

For further analysis, we find that drought occurred frequently in winter and spring in Yunnan in 2015 (Figs. S5 and 7). We also find that abnormally abundant rainfall occurred in January 2015, which usually falls within the dry season. The GNSS-EWH results indicate an abnormal increase during this period; however, the GRACE results cannot clearly capture this signal (Figs. S5 and 7).

Fig. 7
figure 7

Observed seasonal GNSS displacements in May (a) and September (d) in 2015, associated TWS anomalies (TWSA) from GNSS and GRACE analysis in May (b), (c) and September (e), (f)

Full size image

Drought severity index based on GNSS

The drought severity index (DSI) is widely used to assess and monitor drought events. In this study, we create a DSI by utilizing the standardized anomalies within the GNSS_Green based EWH variations (Appendix B), following the approach used in the GRACE-DSI products (Zhao et al. 2017). The GRACE-DSI products are established by calculating the standardized GRACE TWS anomalies and the standard and mean deviation of GRACE TWS for each month. Furthermore, several recent hydrologic loading studies use GNSS displacements and associated TWS estimates to develop new DSIs, indicating that GNSS is capable of monitoring drought intensity (Chew and Small 2014; Famiglietti and Rodell 2013; Jiang et al. 2021a, b).

Since 2009, severe droughts have occurred almost every year in Yunnan. The GNSS-DSI and GRACE-DSI time series (Fig. 8) exhibit obvious periodic fluctuations, corresponding to the drought cycle. From the autumn of 2009 to the spring of 2010, Yunnan suffered one of this century's worst droughts (GNSS-DSI −2.18). In 2011, Yunnan suffered a local persistent drought (GNSS-DSI −1.90), and the rainfall in 25 counties broke the records for the lowest amount of rainfall (China News Network https://www.chinanews.com/gn/2011/08-16/3261444.shtml). In 2012 and 2013, the drought in Yunnan continued to worsen (GNSS-DSI −0.88), and the water supply for tens of millions of people is affected. In 2015, Niño incident reached a strong standard, and Yunnan Province experienced high temperatures and less rain (GNSS-DSI −1.31). From 2015 to 2018, the drought is moderated, the precipitation continued to increase, and there is a sudden drop after 2019, resulting in a serious drought. In May 2019, the drought accelerated, the GNSS-DSI reached −0.55, and the proportion of severe and special drought regions increased continuously. When the intensity of the drought reaches its maximum on May 21, the area impacted by drought accounts for 94% of all Yunnan. Of this area, 66% experienced severe (or worse) drought conditions. In 2020, Yunnan suffered the worst drought (GNSS-DSI −1.77) in nearly a decade (Jiang et al. 2021a). According to the GNSS-DSI of six river basins, there were obvious spatial differences in droughts from October 2012 to May 2013, and the drought in western Yunnan was more severe than that in eastern Yunnan. In 2018, there was no obvious drought in western Yunnan, but there was a local drought in eastern Yunnan (see Fig. 8d for details).

Fig. 8
figure 8

Drought severity index for six sub-basins in the Yunnan region. a Selected 6 river basins (color-coded). b Oceanic Niño Index (red line); monthly GNSS-DSI (blue line) and scPDSI (brown line) (http://climexp.knmi.nl/); green bars show the monthly precipitation averaging anomalies over Yunnan (http://climexp.knmi.nl/selectyear.cgi); the red shaded areas indicate the existence of drought. c Annual-averaged GNSS-DSI for selected river basins. d The GNSS-DSI time series for the selected river basin

Full size image

Based on the above analysis and Fig. 8b, we find that almost all droughts in Yunnan occurred during the "trough period" within the GNSS-DSI. In addition, the self-calibrating Palmer drought severity index (scPDSI), calculated according to the precipitation and potential evapotranspiration data (Palmer 1965), is listed in Table 2 for comparison. We analyze the nine main drought events during 2010–2021 and find that all of the GNSS-DSIs are negative, which indicates that drought possibly occurred. Six scPDSIs are negative, two are positive, and one is absent. The connection between negative GNSS-DSI and drought events indicates that GNSS has the potential to calculate a reliable DSI. Both scPDSI and GNSS-DSI represent distinct types of drought indices, but scPDSI relies on the accuracy of meteorological inputs and a simple water balance model. In regions with sparse ground observations, particularly regarding precipitation and soil characteristics, the reliability of the scPDSI is reduced. The GNSS-DSI and the GRACE-DSI, on the other hand, are computed without model assimilation and account for spatiotemporal variations in local hydroclimatology. We also compare the GNSS-DSI and scPDSI with the Ocean Niño Index, which is a measure of the El Niño/Southern Oscillation and large positive values indicate El Niño and large negative values indicate La Niña; however, except for the drought in 2015 (GNSS-DSI −1.31), this does not show much correlation between the Niño index and drought events in Yunnan.

Table 2 GNSS-identified hydrological events in Yunnan
Full size table

Conclusions

We generate a TWS product by inverting daily GNSS vertical positions using the SBF method, and we utilize various geodetic and hydrological data to characterize spatiotemporal changes in water storage in Yunnan. Multiple datasets, including GNSS, GRACE, GLDAS, and precipitation, indicate a coherent spatial pattern for the annual water storage variation, which indicates decreasing seasonal variability from southwestern to northeastern Yunnan with a maximum amplitude of 100–300 mm. The maximum amplitude is approximately 287 mm (211 mm on average) in western and southwestern Yunnan. Compared to GNSS-derived EWH estimates, both the GLDAS and GRACE significantly underestimate the seasonal oscillations of TWS in Yunnan. While GRACE gravity time series appear to capture the temporal pattern of total surface mass load changes well, the small size of river basins leads to a leakage problem and severe underestimation of the amplitude of water storage cycles. The GNSS-derived time series of land water reserves provides more details regarding hydrological events in small river basins and indicates that the annual amplitude of total land water reserves in six major basins in Yunnan Province is 72.6 ± 6.0 Gt. Our GNSS-inverted EWH also captured several extreme droughts that occurred during 2010–2021, and the GNSS-derived drought severity index indicated a clear connection between drought occurrences. Our study demonstrates that GNSS can serve as an independent technique for observing water storage and provide pivotal clues for understanding hydrological dynamics. Integrating GNSS loading deformation, GRACE gravity, GLDAS, and precipitation into an integrated water monitoring network provides independent measurements of storage variations. Further investigations are needed to develop joint inversion methods that combine hydrological and geodetic observations and models.