The Assessment of Hydrologic- and Flood-Induced Land Deformation in Data-Sparse Regions Using GRACE/GRACE-FO Data Assimilation

The vertical motion of the Earth’s surface is dominated by the hydrologic cycle on a seasonal scale. Accurate land deformation measurements can provide constructive insight into the regional geophysical process. Although the Global Positioning System (GPS) delivers relatively accurate measurements, GPS networks are not uniformly distributed across the globe, posing a challenge to obtaining accurate deformation information in data-sparse regions, e.g., Central South-East Asia (CSEA). Model simulations and gravity data (from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO)) have been successfully used to improve the spatial coverage. While combining model estimates and GRACE/GRACE-FO data via the GRACE/GRACE-FO data assimilation (DA) framework can potentially improve the accuracy and resolution of deformation estimates, the approach has rarely been considered or investigated thus far. This study assesses the performance of vertical displacement estimates from GRACE/GRACE-FO, the PCRaster Global Water Balance (PCR-GLOBWB) hydrology model, and the GRACE/GRACE-FO DA approach (assimilating GRACE/GRACE-FO into PCR-GLOBWB) in CSEA, where measurements from six GPS sites are available for validation. The results show that GRACE/GRACE-FO, PCR-GLOBWB, and GRACE/GRACE-FO DA accurately capture regional-scale hydrologic- and flood-induced vertical displacements, with the correlation value and RMS reduction relative to GPS measurements up to 0.89 and 53%, respectively. The analyses also confirm the GRACE/GRACE-FO DA’s effectiveness in providing vertical displacement estimates consistent with GRACE/GRACE-FO data while maintaining high-spatial details of the PCR-GLOBWB model, highlighting the benefits of GRACE/GRACE-FO DA in data-sparse regions.


Introduction
Land deformation is a response of (1) geophysical processes governed by, e.g., hydrology, atmosphere, ocean mass transports associated with the earth system, and (2) anthropogenic impact associated with the human intervention [1,2]. Such processes cause the Earth's surface to subside or uplift, severely affecting the scientific and industrial sectors [3,4]. Despite being affected by different Earth system drivers, land deformation is dominated by the hydrology cycle at annual and interannual scales [1,[5][6][7]. For instance, vertical motions' annual variations in Japan [8] and the Northwestern United States [9] are induced by snow loadings. Similarly, seasonal variations in California are explained by groundwater recharge and pumping activities [10]. In monsoon regions, significant terrestrial water storage variation (∆TWS), influenced by exceptional annual precipitation (e.g., >1000 mm), is a primary driver of the land deformation's annual variations (e.g., [11][12][13] objectives. Section 2 provides the details of the study area and data processing approaches. Section 3 describes the methods for vertical displacement computation and GRACE DA. The results are presented in Section 4, where the benefits of GRACE and GRACE DA in data-sparse regions are elaborated. Section 5 benchmarks the GRACE DA's performance, and Section 6 concludes the findings of this paper.

Study Area
This study focuses on Central South-East Asia (CSEA; 11 • N-19 • N, 97 • E-106 • E; Figure 1). The study area stretches over 800,000 km 2 , encompassing five different countries (i.e., Myanmar, Thailand, Laos, Cambodia, and Vietnam) and two major river basins (i.e., the Chao Phraya River and Tonlé Sap Basins). CSEA has a tropical wet climate with an average temperature and annual rainfall of approximately 25 • C and 1800 mm, respectively. Essential sources of water supply in most countries are lakes and reservoirs, in which the surface water contributes significantly to regional water storage dynamics, e.g., >40% over the Tonlé Sap Basin [36]. The region has experienced several episodes of extreme monsoon rainfalls and flooding [37]. A recent severe flooding was seen in Thailand in 2011, where the inundated areas extended over 18,000 km 2 [38]. Substantial hydrologic mass variations in CSEA induce a seasonal vertical motion that can likely be measured by a GPS network [33]. However, GPS measurements are very sparse in CSEA. Only six GPS sites' measurements have collected data for longer than three years, and only two of them are still operational (see Table 1). This data-sparse problem limits the comprehensive spatiotemporal deformation details to be observed. The use of satellite observations or model estimates might help improve the deformation monitoring ability in CSEA. Table 1. Characteristics of Global Positioning System (GPS) sites used in this study. The sites with an asterisk (*) are still in operation. The evaluation sites represent the model grid cell location used in the evaluation process. The measurements from GPS sites inside the same PCRaster Global Water Balance (PCR-GLOBWB) model grid cell are averaged. The averaging process results in three (independent) time series of S1-S3 sites (model grid); see, e.g., Figure 3a for their locations on the map.   Table 1.

GRACE Data
The GRACE and GRACE-FO Level 2 release 06 (RL06) products provided by the Center for Space Research (CSR) at the University of Texas at Austin (http://www2.csr.utexas.edu/grace/RL06.html and http://www.csr.utexas.edu/grace-fo) are used to derive the ΔTWSs and of vertical displacement variations (Δr) between April 2002 and December 2018. The products contain monthly (geopotential) spherical harmonic coefficients (SHC) up to a maximum degree and order of 96. The GRACE data are processed as follows. First, degree 1 coefficients are restored using values provided by Swenson et al. [39]. C20 and C30 coefficients are replaced by the value estimated from the satellite laser ranging [40]. Second, the long-term SHC mean (between April 2002 and December 2018) is computed and removed from each monthly product to obtain monthly SHC variations. Third, destriping [41] and 300-km radius Gaussian smoothing filters [42] are applied to the SHC variations to alleviate high-frequency errors. Fourth, the ΔTWS or Δr is computed from the filtered SHC variations (see also Section 3.1), and the Glacial Isostatic Adjustment correction provided by Caron et al. [43] is applied. Finally, a signal restoration approach [36,44] is performed to restore the damped signal caused by the applied filters. It is noteworthy that the signal restoration used in this study is applied to the  Table 1.
In contrast to the GPS measurements, the GRACE Level 2 data do not include shortterm atmospheric and nontidal oceanic variations. Therefore, the de-aliasing atmospheric and oceanic signals provided by the GRACE Atmosphere and Ocean De-aliasing RL06 product [45] are added to the GRACE-derived vertical displacement when comparing ∆r estimates with GPS measurements.

Hydrology Model
The global distributed hydrological model, PCR-GLOBWB Version 2 [34,46], is a gridbased hydrology model that can simulate various water storage components, e.g., soil moisture, groundwater, surface water, and snow. The model has spatial and temporal resolutions of five arcmins and one day, respectively. PCR-GLOBWB-simulated TWS comprises 27 different water storage components, i.e., eight soil moisture, four interceptions, eight snow, four inundated topwater, one surface water, and two groundwater. The complete details of PCR-GLOBWB can be found in, e.g., [34,[46][47][48].
To perform a PCR-GLOBWB simulation, we use daily precipitation and 2-m air temperature from the Integrated Multi-satellite Retrievals for Global Precipitation Measurement Mission (IMERG) [49] and the hourly land product of the European Centre for Medium-range Weather Forecasts Reanalysis Version 5 (ERA5; https://confluence.ecmwf.int/display/CKB/ ERA5-Land), respectively. Both products have a spatial resolution of 0.1 • , and they are resampled to five arcmins to overlay with the model grid space. The ERA5 s hourly data is resampled to a daily interval, consistent with the model temporal resolution. The daily reference potential evapotranspiration is calculated using the Hamon method [50].

GPS Measurements
GPS data are obtained from the Nevada Geodetic Laboratory at the University of Nevada, Reno (UNR) [51]. A final solution provides an hourly vertical motion associated with the IGS14 reference frame. Multiple GPS sites are located in CSEA (see http://geodesy. unr.edu/NGLStationPages/gpsnetmap/GPSNetMap.html) but most contain significant data gaps, i.e., only a few records available in the time series. In total, only six GPS sites have sufficiently long records (i.e., >3 years) for our analysis (Table 1). At each GPS site, the hourly records of the up (vertical) component are averaged into a monthly averaged time series, consistent with the GRACE temporal resolution. The long-term trend is removed from the time series to reduce the plate tectonic effects [7]. The trend is estimated by fitting the GPS time series using the least-squares adjustment associated with the offset, trend, annual, and semi-annual parameters. Then, the long-term mean is subtracted from the time series to obtain GPS position variations. Finally, the measurements of GPS sites inside the same model grid cell (i.e., five-arcmins grid box) are averaged before performing the validation. This averaging process results in three GPS time series associated with three independent model grid cells (S1-S3; see Table 1).

Flood Extent
The flood extent is derived from surface reflectance measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors onboard NASA's Terra satellites. The MOD09A1 collection 6 product [52] provides surface reflectance in 7 different frequency bands at 500-m spatial resolution within an 8-day window. The combinations of specific frequency bands can be used to identify open water bodies larger than 0.25 km 2 . In this study, the multifrequency surface reflectance is extracted from two MODIS tiles, h27v07 and h28v07 (covering the entire study area), between April 2002 and December 2018. The pixels with filled values or severe cloud cover are masked, using data quality control flags provided with the MODIS product. The pixel associated with surface water can be classified using the normalized difference water index (NDWI) computed using measurements from green and near-infrared (NIR) channels as follows [53]: The NDWI values are between −1 and 1, where positive values represent open water, and zero or negative values represent clouds, vegetation canopy, or soil. The flooded extent is computed by multiplying the total number of open water pixels (NDWI > 0) by 0.25 km 2 (the spatial resolution of the surface reflectance dataset).

Methods
The computation and analysis follow three steps: (1) data preparation (see Section 2), (2) vertical displacement computation, and (3) validation and accuracy assessment. Figure 2 illustrates the data processing diagram. km 2 (the spatial resolution of the surface reflectance dataset).

Methods
The computation and analysis follow three steps: (1) data preparation (see Section 2), (2) vertical displacement computation, and (3) validation and accuracy assessment. Figure  2 illustrates the data processing diagram.

Computation of Terrestrial Water Storage and Vertical Displacement Variations
The GRACE monthly geopotential spherical harmonic coefficients (SHC), together with the load Love numbers [7], are used to compute the ΔTWS and Δr as follows [54]:

Computation of Terrestrial Water Storage and Vertical Displacement Variations
The GRACE monthly geopotential spherical harmonic coefficients (SHC), together with the load Love numbers [7], are used to compute the ∆TWS and ∆r as follows [54]: where (θ, λ) are the colatitude and east longitude, (k n , h n ) are the load Love numbers of degree n, P nm is the normalized associated Legendre function of the first kind of degree n and order m, (ρ ave , ρ w ) are the average density of the Earth and freshwater, a is the radius of the Earth, and ∆C nm , ∆S nm are the SHC variation, computed by removing a long-term mean SHC (e.g., calculated between April 2002 and December 2018) from the monthly C nm , S nm . The load Love numbers computed from the Preliminary Reference Earth Model (PREM) [55] are used in this study.
Vertical displacement associated with model-simulated ∆TWS can also be computed using Equation (3). However, the model provides ∆TWS in a geographical gridded format Remote Sens. 2021, 13, 235 7 of 18 (instead of SHC), and the conversion from the gridded data to the geopotential SHC is firstly performed as follows [56]: cos mλ sin mλ sin θ dθdλ (4) In this study, the maximum harmonic degree (N max ) is approximated by N max = 180 • /ψ 1/2 [57], where ψ 1/2 is a half-wavelength of a field (in degrees), i.e., a sampling distance of model-simulated ∆TWS. The atmospheric and oceanic signals from the atmospheric and oceanic de-aliasing (AOD) product (see Section 2.2) are then added to the vertical displacement estimate to obtain signals consistent with the GPS measurements.

GRACE Data Assimilation
GRACE-derived ∆TWS is assimilated into PCR-GLOBWB using the ensemble Kalman smoother (EnKS) [30]. Comprehensive details on the DA configurations are given in Appendix A, while this section describes the core processing concepts. The EnKS comprises two-step runs, a forward run, and a repeat run. In GRACE DA, a forward run is performed to obtain the monthly average TWS states (consistent with a GRACE or GRACE-FO time window), and a repeat run is used to distribute the update evenly throughout the month. A repeat run acts like a smoother to smooth the simulation's trajectory and reduce the time series' spurious discontinuity. The DA process begins with a model propagation (F ) of each ensemble i (at time step t) as follows: where x is a vector containing the state variables, f is a vector of meteorological forcing data, and u represents the model parameters. The ensemble generation process is described in Appendix A. The number of ensembles used in this study is 100. A forward run is performed for approximately one month (T) to obtain the monthly TWS state variables. When GRACE observation is available, the model states are updated as follows: where X contains 27 monthly average TWS variables (e.g., soil moisture, groundwater, and snow; see [30] for more details); y is an observation vector; H is a measurement operator relating the state to the observation vector, and K is a Kalman gain representing the weighting factor that determines the contribution of observation in the state variable's adjustment. The subscript signs (−,+) indicate the update status: before (−) and after (+) the update, respectively. The Kalman gain is expressed as where C H is the covariance matrix of model TWS predictions (H T X i T− ), C v is the error covariance matrix of the GRACE observation, and C XH is the cross-covariance matrix computed between model TWS states and model TWS predictions. The covariance matrices can be computed empirically from the ensemble members [29]. The daily TWS update is then computed by dividing the monthly update values (from Equation (6)) by the number of days in that month. Next, a repeat run is performed by repeating a forward run (Equation (5)) but adding the daily TWS update to the initial states (associated with TWS variables) at the beginning of each day. The repeated run is skipped when GRACE or GRACE-FO observation is absent.

The Estimations of TWS and Vertical Displacement Variations
The methods described in Section 3.1 are used to compute the monthly ∆TWS and ∆r associated with the GRACE data, PCR-GLOBWB model, and GRACE DA estimate between April 2002 and December 2018. The estimated time series are least-squaresfitted and the ∆TWS annual amplitude estimates are shown in Figure 3. Note that only available GRACE/GRACE-FO data are used in the amplitude computation. The PCR-GLOBWB ∆TWS estimate (Figure 3b) exhibits much greater spatial details than the GRACE estimate ( Figure 3a). Strong TWS variations are seen in river and lake areas, implying a significant contribution of surface water to the hydrological dynamics in CSEA. Although a direct comparison between GRACE and PCR-GLOBWB estimates is difficult due to the inconsistency in their spatial resolution, some distinct differences between them can still be seen, particularly in Central and Eastern Thailand, where PCR-GLOBWB has smaller annual amplitudes by approximately a factor of two (see Figure 3a vs. Figure 3b). The model errors associated with (but not limited to), e.g., the inadequate representation of physical processes, uncalibrated model parameters, or inaccurate forcing data, likely lead to underestimated ∆TWS [27]. After assimilating GRACE into PCR-GLOBWB, the ∆TWS annual amplitude estimate is closer to the GRACE result ( Figure 3c). The GRACE DA's application also preserves the high-spatial details of the ∆TWS (from PCR-GLOBWB), and the ∆TWS features over lakes/reservoirs are still clearly seen. The GRACE Δr estimate's annual amplitude exhibits a similar spatial pattern to the GRACE ΔTWS result, showing more significant variations in the northwestern and southeastern parts of CSEA (Figure 3d). Relative to the ΔTWS result (Figure 3a), the Δr spatial pattern is slightly shifted toward Cambodia's coastal areas (see, e.g., the red oval in Figure The GRACE ∆r estimate's annual amplitude exhibits a similar spatial pattern to the GRACE ∆TWS result, showing more significant variations in the northwestern and southeastern parts of CSEA (Figure 3d). Relative to the ∆TWS result (Figure 3a), the ∆r spatial pattern is slightly shifted toward Cambodia's coastal areas (see, e.g., the red oval in Figure 3d), attributed to significant atmospheric and oceanic loading effects in these semienclosed bay regions (see later discussions). The PCR-GLOBWB ∆r estimate agrees with the GRACE result in, e.g., Northern Thailand and Laos (Figure 3e). However, the model underestimates ∆r in Central Thailand and Cambodia, consistent with underestimated ∆TWS seen in Figure 3b. Significant annual ∆r are also seen in the reservoir and lake areas (Figure 3e), confirming the large surface water seasonality in CSEA. In line with the ∆TWS estimate, applying GRACE DA leads to an improved agreement between the ∆r and GRACE estimates and maintains the relatively high spatial resolution of the PCR-GLOBWB model (Figure 3f).
Land deformation is governed by multiple mass loading sources, i.e., land, atmosphere, and ocean. While the land component dominates the ∆r in CSEA, atmospheric and oceanic loadings are found to play a crucial role in certain parts of the region, seen from their notable contributions to total vertical displacement ( Figure 4). The signal contribution associated with different components is computed as the ∆r C /∆r, where ∆r C is the annual amplitude of vertical displacement induced by an individual mass loading component (i.e., land and ocean/atmosphere), and ∆r is the amplitude of the total vertical displacement estimate (i.e., the integration of all components). Figure 4 suggests that the land water storage variation is the primary driver of ∆r (see Figure 4a vs. Figure 4b), while the contribution of atmospheric and oceanic loadings is relatively small (~20%) in the inland areas ( Figure 4c). However, large atmospheric and oceanic variations are found near Thailand and Cambodia coasts (see, e.g., a red oval in Figure 4c), influenced by a strong seasonality of ocean bottom pressure variations in the Gulf of Thailand [58], which contributes up to 60% to the ∆r annual amplitude in the coastal regions. the contribution of atmospheric and oceanic loadings is relatively small (~20%) in the inland areas (Figure 4c). However, large atmospheric and oceanic variations are found near Thailand and Cambodia coasts (see, e.g., a red oval in Figure 4c), influenced by a strong seasonality of ocean bottom pressure variations in the Gulf of Thailand [58], which contributes up to 60% to the Δr annual amplitude in the coastal regions.

The Impact of Flooding on Vertical Displacement Estimates
CSEA has experienced several flood events induced by extreme monsoon rainfalls.  Figure 5a), the surface water areas of CSEA were increased by a factor of 2.5 relative to the regionally average value (Figure 5b). The flooding severities were apparent in the Chao Phraya River Basin, Tonlé Sap Basin, and most of Southeastern Cambodia, where several major cities were inundated for more than one month [38].

The Impact of Flooding on Vertical Displacement Estimates
CSEA has experienced several flood events induced by extreme monsoon rainfalls. Figure 5 illustrates the MODIS-derived NDWI during the flood event in 2011 (Figure 5a; the most severe flood since 2002) compared with the average NDWI value computed between 2002 and 2018 (Figure 5b). The positive NDWI values reveal approximate surface water extents in the region, in which the lake, reservoir, and flooding features are clearly seen. During the flood event in 2011 (July-December 2011; Figure 5a), the surface water areas of CSEA were increased by a factor of 2.5 relative to the regionally average value (Figure 5b). The flooding severities were apparent in the Chao Phraya River Basin, Tonlé Sap Basin, and most of Southeastern Cambodia, where several major cities were inundated for more than one month [38].  Table 2). The discrepancy between PCR-GLOBWB and GRACE is seen in the annual amplitude estimates, particularly in the Tonlé Sap Basin, where the model underestimates the ΔTWS and Δr by >10% (Figure 6b; see also Table 2). After applying the GRACE DA, both the ΔTWS and Δr estimates move toward the GRACE observations, and the ρ and RMSD estimates (with respect to GRACE) are improved by up to 9% and 44%, respectively (see Table 2). Figure 6 also illustrates that the ΔTWS and Δr seasonal variations are governed mainly by the precipitation seasonality (Figure 6e,f). Exceptionally low or high ΔTWSs correspond to extreme precipitation and flooding (see, e.g., ΔTWS in 2011; Figure 6a,e,g), which results in noticeable land subsidence (negative Δr) as a response to the increased water mass. The basin-averaged ∆TWS and ∆r time series over the Chao Phraya River and Tonlé Sap Basins are analyzed together with the precipitation and flood extent ( Figure 6). The PCR-GLOBWB ∆TWS and ∆r estimates depict a significant agreement with the GRACE estimates (Figure 6a-d) in terms of seasonal variation (with ρ > 0.85; see Table 2). The discrepancy between PCR-GLOBWB and GRACE is seen in the annual amplitude estimates, particularly in the Tonlé Sap Basin, where the model underestimates the ∆TWS and ∆r by >10% (Figure 6b; see also Table 2). After applying the GRACE DA, both the ∆TWS and ∆r estimates move toward the GRACE observations, and the ρ and RMSD estimates (with respect to GRACE) are improved by up to 9% and 44%, respectively (see Table 2). Figure  6 also illustrates that the ∆TWS and ∆r seasonal variations are governed mainly by the precipitation seasonality (Figure 6e,f). Exceptionally low or high ∆TWSs correspond to extreme precipitation and flooding (see, e.g., ∆TWS in 2011; Figure 6a,e,g), which results in noticeable land subsidence (negative ∆r) as a response to the increased water mass.

ΔTWS Δr Amp (mm) ρ (-) RMSD (mm) Amp (mm) ρ (-) RMSD (mm)
The   Table 2. Annual amplitude and statistical values (correlation coefficient (ρ) and root mean square difference (RMSD)) computed from the Gravity Recovery And Climate Experiment (GRACE), PCR-GLOBWB, and GRACE DA (data assimilation)derived terrestrial water storage variation (∆TWS) and vertical displacement variation (∆r) time series in the Chao Phraya River and Tonlé Sap Basins. The ρ and RMSD values are computed with respect to the GRACE estimates. The percentage given in the GRACE DA results represents the improvement with respect to the PCR-GLOBWB estimate.

Amp (mm) ρ (-) RMSD (mm) Amp (mm) ρ (-) RMSD (mm)
The The correlation analysis shown in Figure 7 reveals that the flood extent has a stronger connection with the ∆TWS than precipitation. The precipitation-flooding relationship is affected by several factors, e.g., climate conditions, land cover (urbanization), and catchment characteristics [59], which likely lessen the correlation between the two. On the other hand, the ∆TWS reflects a saturated soil and groundwater storage condition, which is an essential factor controlling flood occurrence. When soil or groundwater storage is saturated (e.g., a high ∆TWS), the throughfall cannot infiltrate through the soil but exits the system as surface runoff or a flood. This land process mechanism also leads to a phase delay of flood occurrence (relative to, e.g., the ∆TWS), where the maximum flood extent is observed after the ∆TWS peak by about one month (see, e.g., the Chao Phraya River Basin in 2011; Figure 6g). It is noteworthy that the ∆r estimates' anticorrelations reflect the nature of vertical displacement responses to water mass loads, which results in about one cycle (π) phase difference between ∆r and ∆TWS (and, also, the flood extent) (see, e.g., Figure 6a vs. Figure 6c). The correlation analysis shown in Figure 7 reveals that the flood extent has a stronger connection with the ΔTWS than precipitation. The precipitation-flooding relationship is affected by several factors, e.g., climate conditions, land cover (urbanization), and catchment characteristics [59], which likely lessen the correlation between the two. On the other hand, the ΔTWS reflects a saturated soil and groundwater storage condition, which is an essential factor controlling flood occurrence. When soil or groundwater storage is saturated (e.g., a high ΔTWS), the throughfall cannot infiltrate through the soil but exits the system as surface runoff or a flood. This land process mechanism also leads to a phase delay of flood occurrence (relative to, e.g., the ΔTWS), where the maximum flood extent is observed after the ΔTWS peak by about one month (see, e.g., the Chao Phraya River Basin in 2011; Figure 6g). It is noteworthy that the Δr estimates' anticorrelations reflect the nature of vertical displacement responses to water mass loads, which results in about one cycle (π) phase difference between Δr and ΔTWS (and, also, the flood extent) (see, e.g., Figure 6a vs. Figure 6c).

Validation with GPS Measurements
The Δr estimates are evaluated against the GPS measurements at the S1-S3 sites (Figure 8). The long-term mean of vertical displacement estimates (GRACE, PCR-GLOBWB, and GRACE DA) with respect to the GPS data period is removed from the time series prior to the comparison. The GRACE, PCR-GLOBWB, and GRACE DA Δr estimates depict consistent seasonal variability with the GPS time series in all three evaluated locations. Discrepancies are seen during extreme events (e.g., the flood event in 2011), where GPS measurements exhibit higher variability than that observed from the GRACE and model estimates (see, e.g., the blue ovals in Figure 8b). The higher Δr variations are likely associated with local signals that coarse spatial resolution satellite observations (~300 km) or model estimates (~10 km) cannot effectively capture [6,60]. At S2 and S3 (see Figure 8c and Figure 8e), the excessive land subsidence (negative Δr) during the flood event in 2011 is clearly seen. The subsidence magnitude in July-December 2011 is increased by up to

Validation with GPS Measurements
The ∆r estimates are evaluated against the GPS measurements at the S1-S3 sites (Figure 8). The long-term mean of vertical displacement estimates (GRACE, PCR-GLOBWB, and GRACE DA) with respect to the GPS data period is removed from the time series prior to the comparison. The GRACE, PCR-GLOBWB, and GRACE DA ∆r estimates depict consistent seasonal variability with the GPS time series in all three evaluated locations. Discrepancies are seen during extreme events (e.g., the flood event in 2011), where GPS measurements exhibit higher variability than that observed from the GRACE and model estimates (see, e.g., the blue ovals in Figure 8b). The higher ∆r variations are likely associated with local signals that coarse spatial resolution satellite observations (~300 km) or model estimates (~10 km) cannot effectively capture [6,60]. At S2 and S3 (see Figure 8c,e), the excessive land subsidence (negative ∆r) during the flood event in 2011 is clearly seen. The subsidence magnitude in July-December 2011 is increased by up to 90% compared with the average subsidence computed from all July-December ∆r estimates (see, e.g., Figure 8d,f). The S1 location is located beyond the flood zone (see the S1 location in, e.g., Figure 5b), and no noticeable land subsidence feature is detected in 2011 (Figure 8a). 90% compared with the average subsidence computed from all July-December Δr estimates (see, e.g., Figure 8d,f). The S1 location is located beyond the flood zone (see the S1 location in, e.g., Figure 5b), and no noticeable land subsidence feature is detected in 2011 ( Figure 8a).  The notable agreement between ∆r estimates and GPS measurements is also seen in the annual amplitude and phase estimates ( Table 3). The annual amplitude estimates from GRACE, PCR-GLOBWB, and GRACE-DA agree with the GPS measurements withiñ 1 mm, which is reasonable, considering the 1-cm accuracy of the GPS data in the up component [16]. Similarly, the phase estimates agree within~0.3 months, consistent with the evaluated time interval (one month). The GRACE results show the highest agreement with the GPS estimates, providing a~17% closer annual amplitude than the PCR-GLOBWB results (e.g., 9.11 mm vs. 7.64 mm with respect to 8.67 mm; see Table 3). As seen in Section 4.2, the GRACE DA updates the ∆TWS estimates toward the GRACE observations, resulting in an improved agreement with the GPS data compared with the model estimates alone (e.g., the GRACE DA increases the agreement in the annual amplitude by 8%). Table 3. The annual amplitude and phase estimates from the GPS, GRACE, PCR-GLOBWB, and GRACE DA-derived ∆r at three evaluated sites (S1-S3). The average amplitude and phase values are also given in the right column (Avg).

Amplitude (mm)
Phase ( Validating against the GPS measurements, the GRACE DA ∆r estimates provide ρ and RMSD values of approximately 0.7 and 5 mm, respectively ( Table 4). The RMSD is smaller than the GRACE errors by a factor of two to three [61], denoting the effectiveness of GRACE data usage in data-sparse regions. Like previous analyses, the ∆r estimates associated with GRACE observations (GRACE and GRACE DA) outperform the model estimate, i.e., by providing higher ρ and RMSD values by~3%. The model, however, still presents one advantage in providing a continuous time series, in which ∆r can still be estimated despite the absence of GRACE, GRACE-FO, or GPS data (see, e.g., Figure 8a). Table 4. The ρ, RMSD, and RMS reduction values (with respect to GPS data) estimated from the GRACE, PCR-GLOBWB, and GRACE DA-derived ∆r at three evaluated sites (S1-S3). The statistical values associated with the inclusion (top part) and exclusion (bottom part) of atmospheric and oceanic loadings are shown. The average amplitude and phase values are also given in the right column (Avg). Our final analysis investigates the impact of atmospheric and oceanic loadings on the vertical displacement accuracy. The assessment is performed by comparing two ∆r estimate scenarios, including and excluding atmospheric and oceanic loadings in the ∆r computation. Table 4 illustrates that excluding atmospheric and oceanic signals from ∆r estimates worsens the ρ, RMSD, and RMS reduction estimates by~10%. This analysis reveals that, although land hydrology explains about 90% of the ∆r accuracy, the result is likely suboptimal when considering the land component alone. The inclusion of atmospheric and oceanic signals should be considered to achieve an accurate ∆r estimate.

Discussion
The performances of GRACE, PCR-GLOBWB, and GRACE DA-derived ∆r in datasparse regions (here is CSEA) are assessed in this study. GRACE provides the most accurate ∆r estimates, closest to the GPS measurements, with ρ and RMSD values of > 0.8 and~5 mm, respectively. A similar agreement is seen in most monsoon regions [7,11,13,33]. GRACE fails to capture local signals, due likely to the gravity data's limited spatial resolution (e.g., 300 km), causing a noticeable mismatch between the GRACE and GPS time series [60]. However, the mismatch is negligible compared with GPS uncertainties [16] and does not significantly affect the GRACE's overall performance in CSEA. A notable agreement in annual amplitude and phase estimates is also found. Previous studies (e.g., [11,33,62]) discussed the difficulty of aligning the amplitude or phase between GRACE and GPS data, which was caused by ineffective implementations of earth models. Such an issue is not observed here, likely due to the use of the PREM model, which was found to be the optimal earth model in the Southeast Asia region [33].
The PCR-GLOBWB model can also deliver an accurate ∆r estimate, despite its inferior performance to GRACE. The model underestimates the ∆r annual amplitudes but provides greater spatial details than the GRACE estimates. Underperformed ∆r is attributed in part to incorrect model ∆TWS estimates associated with, e.g., imperfect model physics, ineffective parameter calibration, or inaccurate meteorological forcing input [27]. Assimilating GRACE into PCR-GLOBWB apparently improves ∆r estimates (compared to PCR-GLOBWB results). It is noticed that GRACE DA does not outperform GRACE observations in terms of achieving a high ∆r accuracy. This behavior can be explained by the GRACE DA's nature that optimizes the contributions from the model and observation [29], resulting in the ∆r accuracy falling between the model and GRACE results. However, GRACE DA estimates present one clear advantage (over that of GRACE observations) in providing high spatiotemporal resolution ∆r estimates (e.g.,~10 km daily). In particular, accurate vertical deformation can always be obtained from GRACE DA, despite the GRACE observation absence. The GRACE DA's effectiveness found from this study also supports the concept of Springer et al. [26] in utilizing DA estimates to isolate high-frequency hydrologic signals from GPS measurements.
The GPS, GRACE, PCR-GLOBWB, and GRACE DA results depicted elevated land subsidence during flood events, particularly in 2011, where the subsidence intensity was increased by up to 90%. During flooding, saturated soil and groundwater conditions appear to cause large-scale land subsidence that can be detected from ground, model, and satellite data. Additionally, we found that the CSEA flood extent had a much stronger relationship to the ∆TWS than precipitation. This finding is in line with Tangdamrongsub et al. [36], who utilized such a robust connection for predicting the ∆TWS (from the MODIS-derived surface water extent) in the Tonlé Sap Basin. This study's analysis may underline the potential usage of ∆TWS estimates in large-scale flood monitoring applications.
It is important to include atmospheric and oceanic loadings in the ∆r computation. Omitting such components leads to an accuracy reduction of about 10%. The impact of neglecting the ocean loading component may be intensified in the semi-enclosed bay regions, e.g., the Gulf of Thailand, where the ocean bottom pressure effect is significant [58,63]. Note that this study obtains the atmosphere and ocean signals from the AOD product consistent with GRACE Level-2 processing while using different de-aliasing products or geophysical models, e.g., [64], which may have a different impact on the ∆r estimate. The performance comparison is beyond the scope of this study but can be considered in future analyses.

Conclusions
This study assessed large-scale land deformation in CSEA using GRACE/GRACE-FO, PCR-GLOBWB, GRACE/GRACE FO DA, and GPS data. It was found that the ∆r estimate in CSEA is primarily governed by hydrologic variations and extreme events (i.e., monsoon flooding), which can be detected by the GRACE/GRAE-FO satellites, hydrology model, and ground measurements. The flooding analysis associated with MODIS data illustrated a significant impact of extreme weather conditions on the ∆r estimates, where a degree of land subsidence can be intensified by almost 100% during flooding. Evaluating against the GPS data, the GRACE/GRACE-FO, PCR-GLOBWB, and GRACE/GRACE-FO DA effectively observed regional ∆r annual variations (e.g., spatial resolution > 10 km), while only the GPS measurements could capture local ∆r signals. The evaluation also exhibited the advantage of GRACE/GRACE-FO DA in improving ∆r estimates and maintaining great spatiotemporal details. A similar GRACE DA performance can likely be expected in other data-sparse regions, particularly ones with similar monsoon climates like CSEA.
Future experiments in different regions can be performed to confirm the GRACE DA's effectiveness in ∆r computations.