Using Integrated Geodetic Data to Assess Drought Characteristics Across Four Provinces and Municipalities in Southwest China

1. Introduction

Rising global temperatures and low precipitation frequency make drought a major natural disaster worldwide [1]. Due to its persistence and universality, drought seriously threatens the ecological environment and socioeconomic development. It is one of the most severe meteorological disasters that cause financial losses among many natural disasters [2]. Meteorological, hydrological, agricultural, and socioeconomic experts primarily identify drought as distinct categories. Meteorological drought mainly refers to the lack of long-term precipitation, resulting in dry climatic conditions; Hydrological drought refers to the significant decrease of water body and soil moisture due to the reduction of precipitation or insufficient groundwater supply [3].

Terrestrial water storage (TWS) mainly comprises soil moisture, groundwater, lake and river water. Although TWS accounts for only 3.5 % of the water cycle, it has an essential impact on climate and weather and is also the basis for the survival of life on Earth [4]. Therefore, monitoring TWS changes and quantifying TWS surpluses and deficits are essential for developing effective water resources management and drought response strategies.

Various geodetic techniques and hydrological models have been used to monitor changes in the TWS, including the Gravity Recovery and Climate Experiment (GRACE) and its successor GRACE Follow-On (GFO) and the Global Land Data Assimilation System (GLDAS). GRACE/GFO can monitor large-scale surface mass changes with a spatial resolution of 300-500km and a temporal resolution at monthly intervals [5,6]. However, due to the nearly one-year data gap between GRACE and GFO, GRACE/GFO is limited in monitoring small-scale and short-term TWS changes. Although the GLDAS model can simplify complex hydrological processes through mathematical modeling, it cannot simulate all water components [7].

The emergence of the Global Navigation Satellite System (GNSS) provides another option for monitoring regional TWS changes. GNSS displacement (especially GNSS vertical displacement data) allows for continuous monitoring of the Earth’s surface motion with millimeter-level accuracy to infer TWS changes [8,9,10]. Since then, GNSS vertical displacement inversion for TWS changes has been widely applied in hydrological and geodetic research [11]. Currently, researchers employ two main methods to establish the relationship between TWS changes and surface deformation [12]. One is Green’s function method based on the spatial domain [8], and the other is the Slepian basis function method based on the frequency domain [13]. Many scholars have studied the TWS changes around the world based on the above two methods [14,15,16,17,18,19], and the results show that the use of GNSS vertical displacement alone can continuously obtain the Earth’s load changes in days or months within 50-100 km.

The density of GNSS stations restricts the quality of GNSS inversion TWS change results [20]. When GNSS stations are unevenly distributed, the inversion results will be unstable. Therefore, to weaken and solve this impact, researchers combined GNSS, GRACE/GFO, and other hydrological models to obtain more accurate regional TWS changes for various studies. Fok et al. [20] and Liu et al. [21] proposed a new method of setting up virtual stations, which calculates the TWS obtained from GRACE as the vertical displacements of these virtual stations and uses them as supplementary virtual stations for inversion together with the GNSS stations to enhance the density of the GNSS stations. Li et al. [23] and Zhu et al. [24] also adopted the above method to analyze the TWS changes in the Yangtze River Basin and Yunnan Province. Adusumilli et al. [25] used GRACE-TWS as a constraint in the regularization process of GNSS inversion to obtain higher spatiotemporal resolution results. Carlson et al. [26] employed JPL Mascon to construct limitations to make them more consistent with the actual spatial resolution of GRACE. Yang et al. [27] improved the two types of joint methods (setting the virtual station method and using GRACE/GFO as the spatial constraint method), indicating that the joint methods can integrate the spatial sensitivity advantages of GNSS and GRACE/GFO.

Southwest China (including Sichuan Province, Chongqing Municipality, Yunnan Province, and Guizhou Province) has a complex topography and geomorphology structure. The precipitation here is subject to seasonal changes, resulting in uneven spatiotemporal distribution and significant differences between dry and wet seasons. With the impact of global warming, this region has experienced a substantial decrease in precipitation and an increase in temperature, leading to frequent droughts [28]. Researchers primarily focus on meteorological drought in Southwest China [29], while few studies address the interplay between comprehensive meteorological drought and hydrological drought. This study primarily employs the standardized precipitation evapotranspiration (SPEI) index as the basis for the meteorological drought index, considering multi-scale factors such as temperature change, solar radiation, wind speed, and humidity [30]. Since GNSS vertical displacement serves as a tool for monitoring water storage surpluses and deficits caused by drought and floods, it can be combined with the drought index to analyze drought events [31]. Subsequently, Jiang et al. [32] proposed the GNSS drought index (GNSS-DSI) to analyze the drought situation in Yunnan, which was calculated based on the TWS change obtained from GNSS inversion alone. Li et al. [33] adopted the results of GNSS and GRACE/GFO joint inversion to analyze the extreme flood disasters in Sichuan from 2020 to 2022. Peng et al. [34] obtained the GNSS-DSI of Poyang Lake Basin at the daily scale based on GNSS inversion of daily TWS changes and then analyzed the drought situation in the basin. Zhu et al. [35] integrated GNSS observation data with precipitation metrics to compute precipitation efficiency (PE) and drought index, facilitating a comprehensive analysis of the drought conditions in Guangdong Province. This paper uses the Joint-DSI calculated by the joint inversion of TWS changes as the hydrological drought index.

In this study, we take Southwest China (SCP, YNP, GZP, and CQM) as the research object and set up GRACE/GFO virtual stations in the whole research area. We use Akaike’s Bayesian information criterion (ABIC) to determine the relative weight factor between GNSS and GRACE/GFO, as well as the smoothing factor of prior constraints, to invert the TWS of the research area. Comparison of the results with TWS obtained using GNSS inversion alone (including Green’s function method and Slepian basis function method) and GRACE/GFO Mascon solutions reveal that the TWS derived from four methods maintain good correlation on both time and spatial scales. In addition, we also compare and validate the dTWS/dt obtained from the joint inversion with the TWS changes estimated by the water balance equation (P-ET-R). We comprehensively analyzed the drought events during the period December 2010 to June 2023 from both meteorological and hydrological perspectives. Finally, we conducted a quantitative study on drought events across four provinces and municipalities in Southwest China regarding drought frequency, severity, and classification of drought and flood.

2. Data and methodology

2.1. Study area

The study area is located in southwest China (Figure 1), encompassing Sichuan Province (SCP), Yunnan Province (YNP), Chongqing Municipality (CQM), and Guizhou Province (GZP), covering a total area of about 1.13×106 km². The topography of Southwest China is intricate, predominantly mountainous or hilly, with a terrain primarily defined by its high northwest and low southeast elevation. The prevalent climate types are the subtropical humid monsoon climate, the subtropical plateau monsoon humid climate, and the distinctive plateau climate akin to the Qinghai-Tibet Plateau [28]. Moreover, precipitation exhibits uneven spatial and temporal distribution, with more abundant in summer and autumn, scarce in winter and spring, and predominantly concentrated in the southwest region. Due to global warming, precipitation in Southwest China has significantly declined, leading to frequent drought occurrences that have severely impacted local agricultural production. Consequently, conducting drought detection and research in Southwest China is imperative.

2.2. GNSS vertical displacement time series

From the GNSS vertical displacement time series provided by the China Crustal Movement Observation Network (CMONOC) (http://www.eqdsc.com), we selected the vertical displacement time series of 79 stations in the southwest region for processing from December 2010 to June 2023. Figure 1 illustrates the distribution of stations, with an average distance of nearly 120km between 79 stations and the distribution of stations is uneven, mainly in the west (SCP and YNP) and sparsely in the east (GZP and CQM). The data provided by CMONOC performed baseline processing through GAMIT/GLOBK10.70 [36] and overall adjustment using QOCA software [37], ultimately obtaining the time series within the ITRF2008 framework.

Post-processing is needed to obtain a continuous and stable GNSS vertical displacement time series based on hydrological signals. Firstly, we employ the least square fitting method to eliminate data jumps, and then we apply the quartile range gross error detection method to remove gross errors. The Kriging Kalman filter interpolation software is used to find missing data [38]. In addition, to obtain GNSS observations dominated by hydrological signals, we use the non-tidal atmospheric load (NTAL) and non-tidal ocean load (NTOL) models under the geometric center reference framework (CF) provided by the German Research Centre for Geosciences (GFZ) to deduct the effects of atmospheric load and ocean load (http://esmdata.gfz-potsdam.de:8080/repository). The annual amplitude of NTAL is between 1.5 ~ and 5.5 mm, and the yearly amplitude of NTOL is between 0.1 ~ and 0.6 mm (Figure S1). Finally, the annual amplitude of the GNSS station is between 1 ~ and 10 mm. Figure 2 compares the vertical displacement time series of four stations in Southwest China before and after post-processing. After the above post-processing, we get a complete and clean time series. Therefore, post-processing of the original GNSS time series is significant for the subsequent use of GNSS data to obtain TWS changes.

2.3. GRACE/GFO Mascon solutions

This paper used the GRACE/GFO Mascon solutions released by the Center for Space Research (CSR) at the University of Texas at Austin to estimate the TWS changes in Southeast China [39] (https://www2.csr.utexas.edu/grace/RL06_mascons.html). The CSR Mascon solutions have applied all appropriate corrections, including C20 and C30 replacements, first-order geocentric corrections, and glacial isostatic adjustment [40]. Since the GRACE data was cut off in June 2017 and GFO was launched in May 2018, nearly one year of missing data occurred. This paper did not fill in the corresponding missing data; the remaining missing data were completed using the cubic spline interpolation method.

2.4. Hydrometeorological data

Multi-source hydrometeorological data can help us to verify the inversion results using geodetic data. In this study, monthly precipitation (P), evapotranspiration (ET), and runoff (R) data with a spatial resolution of 0.1° × 0.1° are utilized, which are provided by ERA5-Land (https://cds.climate.copernicus.eu/ ). ERA5-Land is produced by redrafting the land portion of ERA5 climate reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) [41]. GLDAS NOAH provides monthly P, ET, and R data with a spatial resolution of 0.25° × 0.25° [42] (https://disc.gsfc.nasa.gov/ ). The Global Precipitation Climate Project (GPCP) provides monthly P data with a spatial resolution of 2.5°×2.5° derived from the combinations of satellite estimation and land rain gauge measurements [43] (https://www.ncei.noaa.gov/). The Global Land Evaporation Amsterdam Model (GLEAM) provides monthly ET data with a spatial resolution of 0.1° × 0.1°, which can maximize the extraction of evaporation information from climate and environmental variables observed by different satellites [44] (https://www.gleam.eu/). The monthly P data from the three datasets are consistent between December 2010 and June 2023, with GLEAM providing significantly smaller ET data than ERA5 and NOAH, and ERA5 providing considerably smaller R data than NOAH. As shown in Figure 3, the P, ET, and R data in southwest China exhibit obvious seasonal variations, all peaking in the summer of each year. However, the magnitude of the P time series is much larger than that of ET and R, indicating that the P data significantly impacts the TWS change here more than the other two data.

In order to verify the correctness of the TWS change results, this paper quantifies TWS changes based on the water balance equation [27,45].

$$dTWS/dt=P-ET-R$$

(1)

The TWS changes obtained from GNSS, GRACE/GFO, and joint inversion are calculated as the first-order differential to time, as follows:

$$dTWS/dt=TWS(t+1)-TWS(t)$$

(2)

The TWS change calculated by P, ET, and R represents the change in the middle of each month. In contrast, the TWS first-order differential calculated using Equation (2) represents the TWS change at the beginning or end of the month. To ensure the consistency of time, this paper adopts the time series filtering method proposed by Landerer et al. [46] to filter P-ET-R:

$$\widehat{F}(t)={\displaystyle \frac{1}{4}}F(t-1)+{\displaystyle \frac{1}{2}}F(t)+{\displaystyle \frac{1}{4}}F(t+1)$$

(3)

where $F$ denotes P, ET and R, and F denotes the filtered data.

The SPEI selected in this paper is provided by the Climate Research Center (CRU) at the University of East Anglia (https://digital.csic.es/).

3. Methodology

3.1. Green’s function

According to the mass loading theory, Farrell proposed that TWS changes can cause the elastic deformation of the surface and their relationship can be expressed by a linear equation as follows [8]:

$$d_1=G_1{x}_1+e_1$$

(4)

where $x_1$ represents the mass loading vector expressed as the equivalent water height (EWH); $d_1$ represents the GNSS displacement vector; $e_1$ represents the error term; $G_1$ represents the calculated Green’s function. Argus et al. proposed the least squares model for GNSS vertical displacement inversion of TWS as follows [14]:

$${\Vert {{\sigma }_1}^{-1}(G_1{x}_1-d_1)\Vert }^2+{\beta }^2{\Vert L(x_1)\Vert }^2\to \min$$

(5)

$$E={{\sigma }_1}^2$$

(6)

where ${\sigma }_1$ denotes the standard deviation; Since the GNSS observations are less than the EWH grids, it is necessary to add the Laplace spatial smoothing constraint to the inversion model, so $L$ represents the Laplace operator; $\beta$ is the regularization factor determined by the Generalized Cross Validation (GCV) method, with a value of 0.29 (Figure S2). Through the least squares algorithm, we divide the study area into 0.5° × 0.5° grids, then the EWH value of each grid can be calculated:

$${\widehat{x}}_1={\left({G^T}_1{E}^{-1}{G}_1+\beta L^{T}L\right)}^{-1}{G^T}_1{E}^{-1}{d}_1$$

(7)

3.2. Slepian basis function

Han and Razeghi proposed to estimate the spectral information of the GNSS displacement through a regional spherical basis function [13] and apply the loading Love number to obtain the EWH of the region, which can be expressed as:

$${\sigma }_L\left(\theta ,\lambda ,t\right)=a{\displaystyle \sum _{\beta =1}^J{\gamma }_{\beta }{s}_{\beta }^{EWH}\left(t\right)g_{\beta }\left(\theta ,\lambda \right)}$$

(8)

where $a$ is the Earth’s radius; $J$ is the truncation order; ${\gamma }_{\beta }$ denotes the eigenvalue; $s_{\beta }^{EWH}$ is the Slepian basis function coefficients; $g_{\beta }$ is the Slepian basis function. Figure S3 displays the concentration ratio of each Slepian basis function. Figure S4 represents the energy spatial distribution of Slepian basis functions, with the black solid line indicating the boundary of the study region after an extension of 2°; The first 28 basis functions with concentration ratios greater than 0.1 are well concentrated in Southwest China. According to research by Li et al. [17], the average distance between GNSS stations serves as a suitable radius for Gaussian filtering. Given that the average distance between GNSS stations in this study is approximately 120 km, we set the radius of Gaussian filtering to 150 km. The minimum and maximum distances between adjacent stations are about 50 km and 323 km, while the value 65 is close to the ratio of 20000/323km (wherein 20000 km represents half of the Earth’s circumference), so the maximum order in this paper is set to 65.

3.3. Joint inversion

This paper employs a joint inversion model founded on incorporating virtual stations, as advocated by Fok et al. [21] and Liu et al. [22]. By integrating this with the modified mathematical model proposed by Yang et al. [27], we strategically establish virtual stations across the entire study region and introduce a weighting factor to ascertain the optimal balance between GNSS and GRACE/GFO datasets, thereby mitigating the risk of overfitting. TWS derived from GRACE is inverted to the vertical displacement of the GNSS virtual station and then performs inversion together with the GNSS station to augment its density. The following outlines the displacement inferred from the forward model of GRACE/GFO:

$$d_2=G_2{x}_2+e_2$$

(9)

where $d_2$ is the displacement vector calculated by GRACE/GFO; $G_2$ is the Green’s function of GRACE/GFO virtual station; $x_2$ is the corresponding GRACE/GFO solution in the study area; $e_2$ denotes the error term. In order to determine the location and number of GRACE/GFO virtual stations, this paper uses the checkerboard test to investigate the spatial resolution of TWS inversion by virtual stations [20]. The results show that we can recover the grid signal with a spatial resolution of 3.5° × 3.5° by deploying GRACE/GFO virtual stations at an interval of 2.5° in Southwest China (Figure S5), which matches the original spatial resolution of GRACE/GFO (~ 350 km). Therefore, this paper evenly distributes 56 GRACE/GFO virtual stations in the study area at 2.5° intervals. The virtual station distribution method is also applicable to the virtual station settings in Sichuan and Yunnan [27].

The displacement and Green’s function of GRACE/GFO are obtained from the forward model, composing new observation data as follows:

$$d={\left(\begin{array}{cc}{d}_1& d_2\end{array}\right)}^T$$

(10)

$$G={(\begin{array}{cc}{G}_1& G_2\end{array})}^T$$

(11)

GNSS and GRACE/GFO are two independent data types, and then introduce the weighting factor $\alpha$ to obtain the variance matrix:

$$E({\widehat{\alpha }}^2)=\left(\begin{array}{cc}{\alpha }^2{\sigma }_1^2& 0\\ 0& {{\sigma }_2}^2\end{array}\right)$$

(12)

Bring the above observations into the least squares model of Equation (5), we can obtain the EWH value of each grid calculated by the joint inversion:

$$\widehat{x}={(G^{T}E{\left({\widehat{\alpha }}^2\right)}^{-1}G+\widehat{\beta }{L}^{T}L)}^{-1}{G}^{T}E{\left({\widehat{\alpha }}^2\right)}^{-1}d$$

(13)

This paper employs the ABIC method to ascertain the weight factors within the joint inversion model. ABIC is an objective method based on the principle of entropy maximization, integrating the Akaike information criterion with the Bayesian theorem [47], and is commonly used to determine the weight factor in geodetic inversion [21,27,33,48]. By solving the minimum value of the ABIC function, we can effectively estimate the optimal weights between GNSS and GRACE/GFO data.

$$\begin{array}{c}ABIC({\alpha }^2,{\beta }^2)=n\log f(\widehat{x})+\log \Vert G^{T}E{({\widehat{\alpha }}^2)}^{-1}G+{\beta }^T{L}^{T}L\Vert -\log \Vert {\beta }^T{L}^{T}L\Vert +\\ \log \Vert E({\widehat{\alpha }}^2)\Vert +C\end{array}$$

(14)

$$f(\widehat{x})={\left(Gx-d\right)}^{T}E{({\widehat{\alpha }}^2)}^{-1}\left(Gx-d\right)+x^T{\beta }^2{L}^{T}Lx$$

(15)

where $n$ is the number of observations; $\widehat{x}$ is calculated by Equation (13); C is a constant. By seeking the minimum value of the ABIC function, the optimal weighting factors $\alpha$ and $\beta$ can be determined. Figure 4 illustrates the distribution of ABIC values in May 2023. The minimum value of ABIC occurs at $\log \left({\alpha }^2\right)=0.69$ , $\log \left({\beta }^2\right)=-2$ , and then the optimal weighting factor is obtained.

3.4. Drought index

GNSS-DSI is developed by Jiang et al. [32] based on GRACE-DSI [49]). Table S1 lists the classification of drought severity.

The calculation of GNSS-DSI is as follows:

$$GNSS-DS{I}_{i,j}={\displaystyle \frac{EW{H}_{i,j}-\overline{EW{H}_j}}{{\sigma }_j}}$$

(13)

where $EW{H}_{i,j}$ represents the water storage in the j-month of the i-year; $\overline{EW{H}_j}$ represents the average of the $j$ -month of each year; ${\sigma }_j$ represents the standard deviation of the j-month of each year.

The technical flowchart of this paper is as follows:

4. Result

4.1. Spatial distribution characteristics of TWS

Figure 6 illustrates the annual amplitudes of TWS derived from GNSS-Green, GNSS-Slepian, GRACE, and Joint methods in Southwestern China from December 2010 to June 2023. All methods consistently pinpoint the peak annual amplitudes to the southwest of our study area (i.e., southwest of YNP), with a maximum amplitude of about 284 mm for GNSS-Green, 325mm for GNSS-Slepian, and 208 mm for GRACE. GNSS-Green and GNSS-Slepian exhibit strong spatial distribution coherence, given that both methods rely solely on GNSS data for inversion (Figure 6a,b). Their spatial distribution patterns reveal that the TWS annual amplitude generally tapers off from west to east, but variations in magnitude occur, with GNSS-Slepian results 25% higher than those of GNSS-Green [17,50]. According to Argus et al. [14], the GNSS vertical displacement results from the combined effect of all hydrological loads within a radius of 50 km around the GNSS station; however, the mathematical model of Slepian basis function theoretically concentrates all hydrological influences on the GNSS station, resulting in biased inversion results. Since GNSS-Green is primarily influenced by the density and distribution of GNSS stations, GNSS-Slepian demonstrates a closer alignment with GRACE-TWS than GNSS-Green, particularly in regions with sparse GNSS stations, such as GZP.

The spatial distribution of GNSS inversion results diverges from that of GRACE-TWS. Although the peak of the TWS annual amplitude is located in the southwest of YNP, the nuanced spatial distribution patterns differ markedly. GRACE-TWS demonstrates an escalating amplitude trend from north to south, while GNSS inversion results exhibit an amplitude increase from east to west, aligning with the trend in GNSS vertical displacement’s annual amplitude variation (Figure S1c). This may be due to the sparse distribution of GNSS stations in the northeast of the study area. In contrast, GNSS stations are more densely clustered in the central part of SCP and the northern region of YNP, enabling GNSS to capture TWS signals in these localized areas more accurately. In these regions, GRACE-TWS struggles to discern TWS changes. Consequently, GNSS inversion results are better equipped to accurately capture TWS changes in areas with dense GNSS station coverage, particularly in higher-altitude plateau regions. Conversely, GRACE-TWS, constrained by its relatively low resolution [16], proves more reliable in areas with sparse GNSS station distribution. Besides, the spatial distribution of Joint-TWS closely resembles that of GRACE/GFO overall. In sparse GNSS station regions, Joint-TWS mitigates the effects of limited GNSS station availability on inversion results, thereby enhancing their reliability. At the same time, it retains the advantages GNSS provides in areas with a dense station distribution.

Figure 7 depicts the annual phases derived from GNSS-Green, GNSS-Slepian, GRACE, and Joint methods spanning from December 2010 to June 2023. The annual phases obtained by these four methods manifest commendable spatial coherence, except for GNSS-Green, which shows significant divergence in the east of the study area, such as CQM.

4.2. Temporal distribution characteristics of TWS

To further delve into the validity of the joint inversion results and the distinctions among various inversion methods, we investigate the dTWS/dt time series derived from GNSS-Green, GNSS-Slepian, GRACE, Joint methods, precipitation data, and water balance equations. Figure 8a illustrates the average time series of TWS inferred by GNSS-Green, GNSS-Slepian, GRACE, Joint methods, and precipitation within Southwest China, which manifests obvious seasonal patterns. The results reveal that the two-time series of TWS changes obtained by pure GNSS inversions exhibit remarkable consistency with a correlation coefficient of 0.98. Nevertheless, GNSS-Slepian results are quantitatively more prominent than those derived from GNSS-Green. Additionally, the time series of TWS changes extracted from pure GNSS inversion is approximately threefold that of GRACE/GFO, with correlation coefficients of 0.82 and 0.79, respectively. The value of the joint inversion results resides between the two, and the trend aligns more closely with that of GRACE/GFO, exhibiting a correlation coefficient of 0.83. Compared to the results derived from pure GNSS inversion, the trajectory of the joint inversion results approximates more faithfully that of GRACE/GFO. In contrast to the trend of precipitation data, the geodetic datasets demonstrate commendable consistency. Nevertheless, the intricate transport mechanisms of terrestrial water potentially cause a delay of 1 to 2 months in the apex value of the geodetic dataset relative to the precipitation data [51]. Since 2023, TWS fluctuations across various geodetic datasets have diminished substantially. Subsequent analysis reveals that this occurrence is predominantly a consequence of the sharp decline in precipitation during this period, which plummeted to a mere 45% of the levels recorded in corresponding periods of prior years. The impact of precipitation on the results of GNSS inversion is markedly more pronounced than that exerted by GRACE/GFO data. Figure 8b shows the time series of dTWS/dt derived from GNSS-Green, GNSS-Slepian, GRACE, and Joint methods. Table 1 presents the correlation coefficients between dTWS/dt derived from various datasets, revealing a robust correlation among the geodetic datasets. Nevertheless, we could ascribe the more pronounced variation in dTWS/dt yielded by using GNSS exclusively for TWS inversions to the heightened sensitivity of GNSS in discerning anomalous TWS changes.

The trend of dTWS/dt derived from geodetic data is similar to that of P-ET-R, but the amplitude is noticeably more significant than that of P-ET-R, suggesting that the P-ET-R calculated via the water balance equation may somewhat underestimate the changes in TWS within the region. This discrepancy may arise because P-ET-R does not account for the influence caused by deep groundwater composition or human activities [17]. Compared to the P-ET-R determined using NOAH data (0.49 ~ 0.76), the TWS differentials obtained from geodetic data exhibit a stronger correlation with those calculated from ERA5 data (0.69 ~ 0.79). This elevated correlation could be attributed to the superior spatial resolution of ERA5 data, providing more detailed and precise meteorological information, whereas the NOAH model is based on the simulation of land surface processes, constrained by the model’s underlying assumptions and parameter configurations. Compared with the results from utilizing GNSS alone to infer TWS changes, the correlation coefficient between Joint methods and other alternatives demonstrates a higher degree of association. In the subsequent investigation of drought events in Southwest China, this paper adopts the Joint-DSI results as the hydrological and meteorological drought index to facilitate comparative and analytical endeavors.

4.3. GNSS-based drought index

We delineate an event as a hydrological drought when the water deficit persists for three months or more. Our joint inversion analysis results derive the data in Table 2. Our findings reveal five distinct drought events in the southwestern region from December 2010 to June 2023. The most prolonged event transpired from January 2016 to December 2018, lasting 36 months. The most acute single drought event occurred in September 2022, exhibiting an average deficit of 86.133 km².

Given the pronounced topographic differences in Sichuan Province (SCP), Yunnan Province (YNP), Chongqing Municipality (CQM), and Guizhou Province (GZP), coupled with the substantial variations in the spatial distribution of precipitation, we have elected to undertake a meticulous analysis to evaluate the drought attributes of these areas more thoroughly. This analysis aims to provide a more precise depiction of the specific drought conditions across four provinces and municipalities in Southwest China, thus allowing us to understand the impact of topography and precipitation distribution on drought conditions. Figure 9 presents the time series of Joint-DSI, SPEI, and precipitation deficits for SCP, YNP, CQM, and GZS. While there are variations between the meteorological drought index (i.e., SPEI) and the hydrological drought index (i.e., Joint-DSI) throughout the study period, the trend remains congruent during instances of severe drought (when the drought index falls below -1.3). The changes in Joint-DSI marginally lag behind SPEI, as the meteorological drought index is typically more responsive than its hydrological counterpart. The meteorological index offers a direct reflection of alterations in precipitation, whereas the hydrological index necessitates intricate computations via processes such as the hydrological cycle to ascertain drought occurrences. As a result, SPEI demonstrates a stronger correlation with fluctuations in precipitation than Joint-DSI. Compared to the results in Table 2, Figure 9 affords more details on the temporal patterns of drought occurrence in Southwest China. Table 2 scrutinizes the Southwest region and estimates five drought incidents based on fluctuations in TWS. Nonetheless, juxtaposing this with Figure 9 unveils significant disparities across the different regions. Specifically, Yunnan remained unscathed by severe drought conditions throughout the prolonged drought period from January 2016 to December 2018. Consequently, a more nuanced analysis of the drought situation in Southwest China is imperative in conjunction with insights from Figure 9.

This investigation concentrates on the spatial evolution of drought events within Southwest China from June 2022 to June 2023. Figure S6 elucidates the comprehensive spatial developmental trajectory of Joint-DSI during this drought episode. We can discern that throughout this interval, a substantial portion of the Southwest region remains enshrouded in a state of drought. Nonetheless, the overall intensity of drought characterized by CQM remains relatively modest, an observation that resonates with Figure 9d. The locales exhibiting the most pronounced drought conditions are situated southwest of YNP, west of SCP, and at the confluence of GZP and YNP. This spatial configuration mirrors the spatial distribution of TWS annual amplitude, as shown in Figure 6. Since February 2022, the drought has escalated in severity, possibly attributable to the increasing precipitation deficiencies.

5. Discussion

5.1. Spatial resolution of joint inversion

The spatial resolution of GNSS and GRACE/GFO is different, and they are suitable for respective spatial scales. Therefore, discussions on the spatial resolution of the joint inversion results are necessary. We designed two checkerboards with different spatial resolutions (2.5°×2.5° and 3.5°×3.5°) to test the results of the joint inversion. Figure 10a,b depict that the white and brown checkerboards represent EWH signals of 0mm and 400mm, respectively. Subsequently, we calculate EWH from the GNSS station and the GRACE/GFO virtual station, as shown in Figure 10c,d.

The joint inversion effectively recovers the EWH signal in the region. We can accurately recover the 2.5° signal in areas with densely distributed GNSS stations; however, restoring the 2.5°×2.5° spatial resolution is less effective in regions with sparse GNSS stations and only GRACE/GFO virtual stations. For the checkerboard with a spatial resolution of 3.5°×3.5°, only the GRACE/GFO virtual station can recover the corresponding signal (Figure S5), and the joint inversion results can also effectively recover the signal at this spatial resolution.

The checkerboard test results indicate that the joint inversion incorporating GRACE/GFO virtual stations can compensate for the spatial resolution limitation in areas with unevenly distributed GNSS stations. This approach successfully integrates the advantages of both GNSS and GRACE/GFO satellite spatial resolutions.

5.2. Quantification of drought characteristics in Southwest China

The Joint-DSI is classified into five drought and wet categories based on Table B1 [52]. We have calculated the proportion of drought and flood areas to the total area of Southwest China during the study period, as shown in Figure 11. Due to distinct topographic conditions across four provinces and municipalities in Southwest China, there are different distributions of droughts and floods derived from the hydrological drought index. Following Peng et al. [34], we define drought frequency as the proportion of time when the drought index falls below D0. Therefore, this study defines drought frequency as the proportion of months with a drought index below -0.5 and then calculates the drought frequency distribution of Joint-DSI on the monthly scale during the periods of December 2010 to June 2023 and July 2022 to June 2023, respectively (Figure 12). This scheme enables a thorough investigation of drought characteristics in Southwest China from both temporal and spatial viewpoints.

For SCP (Figure 11a,b), exceptional drought (D4) primarily occurred from July 2017 to June 2018 and October 2022 to June 2023. Approximately 39% and 32% of the SCP experienced exceptional drought during these periods. Conversely, exceptionally wet conditions (W4) were observed mainly between September 2021 and June 2022, affecting about 14% of the area. Li et al. [33] reported that during the 2022 drought event, drought was predominantly concentrated in the Sichuan Basin (SCB) and the western mountainous area (WTM), with SCB more severely impacted than WTM. From July 2022 to June 2023, moderate drought (D1) covered the highest proportion of the SCP, approximately 72%. The drought frequency was higher in the west and lower in the east, with the central part of the region experiencing continuous exceptional drought (Figure 12b).

In YNP (Figure 11c and d), exceptional drought was most widespread from July 2017 to January 2018 and August 2022 to June 2023, impacting approximately 11% and 25% of the region in these respective periods. On the other hand, exceptionally wet conditions prevailed from January to July 2022, affecting about 13% of the area. Notably, from July 2022 to June 2023, the drought-affected area in YNP escalated to its peak, covering roughly 80% of the region. We observed that the drought frequency was lowest in the central area and highest at the peripheries, with the southwest section of the plateau experiencing persistent exceptional drought (Figure 12b).

Regarding GZP (Figure 11e,f), exceptional drought was primarily evident between July 2017 and June 2018, as well as from March to June 2023, affecting about 28% and 39% of GZP during these respective timeframes. On the contrary, exceptionally wet conditions were mainly noted between December 2010 and July 2011 and from October 2021 to June 2022, with approximately 15% and 16% of the region being exceptionally wet, respectively. As a typical karst landscape, the intricate topography of GZP plays a significant role in its drought conditions [53]. Notably, from July 2022 to June 2023, the drought-affected area in GZP peaked, comprising roughly 73% of the total area. The drought frequency was higher in the western part and lower in the eastern part, primarily manifesting as exceptional drought (Figure 12b).

For CQM (Figure 11g,h), exceptional drought was particularly noticeable between July 2017 and June 2018, as well as from October 2022 to June 2023, impacting approximately 46% and 20% of CQM during these respective timeframes. Conversely, exceptionally wet conditions prevailed between October 2021 and June 2022, affecting about 55% of the region. It is worth noting that from July 2022 to June 2023, the drought-affected area in CQM escalated to its highest proportion, roughly 57%. The drought frequency was higher in the eastern part and lower in the western part, primarily marked by exceptional drought (Figure 12b).

The duration and intensity of exceptional drought and wet conditions in Southwest China exhibit significant variations. We can attain a more profound comprehension of drought characteristics within this region via rigorous computations and analyses. Notably, the primary occurrences of unusual drought across various regions correspond to drought events 4 and 5, as outlined in Table 2. Specifically, drought event 4 is the most prolonged in Southwest China, whereas drought event 5 emerges as the most severe.

6. Conclusion

This paper employs multiple methods to monitor regional TWS changes, including using pure GNSS data based on Green’s function and Slepian basis function, pure GRACE/GFO Mascon data, and joint GNSS and GRACE/GFO datasets. For the integrated GNSS and GRACE/GFO data, we adopt the ABIC approach to ascertain the optimal weighting factor. Additionally, we delve into drought characteristics across four provinces and municipalities in Southeast China in combination with hydrological data. The spatial distribution of TWS annual amplitudes derived from these four methods demonstrates consistency in Southwest China, with a discernible decrease from southwest to northeast. It is worth highlighting that the most pronounced annual amplitude is evident in Southwest Yunnan. Our investigations reveal that the GNSS-Slepian results are approximately 1.25 times more pronounced than that of GNSS-Green, and GRACE/GFO significantly underestimates the TWS changes in the Southwest compared to GNSS. It is intriguing that Joint-TWS poses a stronger spatial resemblance to GRACE-TWS and demonstrates a more pronounced correlation with GRACE/GFO and dTWS/dt calculated based on the water balance equation. This enhanced correlation significantly bolsters the reliability of results in regions with sparse GNSS station coverage. Leveraging these joint inversion results, we have successfully tracked five drought events in Southwest China from December 2010 to June 2023. Among these events, the fifth is the most severe regarding water deficit. We perform individual computations and analyses for SCP, YNP, CQM, and GZP to investigate drought characteristics better. The results reveals that all regions exhibit remarkable drought characteristics during the fifth drought event, accompanied by widespread drought indicators. Our research indicates that integrating GNSS and GRACE/GFO data for TWS change inversion can effectively compensate for the constraints posed by uneven GNSS station distribution, as well as the inherent challenges of limited resolution and incompleteness of GRACE/GFO data, and provide an alternative for monitoring regional TWS changes and drought conditions using geodetic data.