# Two Decades of Groundwater Variability in Peru Using Satellite Gravimetry Data
## 2 Materials and Methods
### Methods
#### 2.3.1 Calculation of GWS Changes
This dual approach allows both physical and statistical attribution of changes in GWS, integrating climatic variability and anthropogenic impacts. Natural Recharge Estimation Natural groundwater recharge was estimated using a combination of satellite observations and hydrological modeling. Monthly anomalies of total water storage (TWS) derived from GRACE and GRACE-FO missions were used, which reflect changes in both surface and subsurface water masses [38]. To isolate the GWS component, a water balance approach was applied by subtracting contributions from soil moisture, snow, and surface water using global model products such as GLDAS or WGHM [39]. The residual signal represents monthly GWS variations per pixel or region. From this time series, recharge events were identified using the water table fluctuation (WTF) method, which estimates net recharge (R) as follows: R = Sy \(\times\)\(\Delta\)h where: Sy is the specific yield of the aquifer, and \(\Delta\)h is the increase in storage anomaly between a minimum (trough) and a subsequent peak. This method is particularly useful for detecting episodic recharge driven by intense precipitation events. The validity of this estimation was confirmed through comparison with ground-based methods (e.g., water balance, tracers, and well level measurements). Studies in arid zones such as Ordos (China) and Abu Dhabi have shown high agreement between GRACE and in situ observations (R\({}^{2}\) > 0.9), supporting the reliability of this approach. Finally, monthly GWS percentiles were calculated relative to the 2003-2023 baseline period, following methodologies proposed by [40]. Percentiles allow for evaluating the relative magnitude of hydrological anomalies: low values (near 0) indicate abnormally dry conditions, while high values (near 100) reflect extreme surpluses. Model Variability and Uncertainty Assessment To address interannual climate variability, the following techniques were employed: Trends and anomalies were analyzed within specific climate contexts, such as El Nino and La Nina events [41]. Seasonal effects were minimized using smoothing techniques and moving averages [42]. Hydrological models were used to estimate a range of uncertainty ([PERSON] et al., 2018) [43]. Validation was performed using in situ observations and groundwater level records from selected wells [39]. Sensitivity analyses were conducted to evaluate model responses to variations in forcing inputs [44]. This integrated approach enables a spatially continuous and robust representation of Peru's groundwater storage dynamics, providing essential inputs for water balance assessments and sustainable resource management. _Model Error Characterization and GRACE Uncertainty Propagation_GRACE products incorporate a formal covariance matrix of spherical harmonic coefficients (Stokes), capturing error contributions from instrumental noise, tide models, atmospheric corrections, and more. From this matrix, the variance for each region or pixel can be derived, representing the monthly uncertainty of the estimates. Given the complexity of propagating the full covariance matrix, spatial covariance models (homogeneous, non-stationary, or anisotropic) are used to replicate error structures and enable regional uncertainty estimations with low computational cost. _Filtering and Leakage Errors_  
GRACE uncertainty also includes components introduced by filtering (e.g., destriping) and mass concentration scaling. GRACE-Tellus products provide separate maps for measurement error and leakage error, allowing total error computation as follows:  
\[\text{Total error}=\sqrt{\text{measurement error}^{2}+\text{leakage error}^{2}} \tag{3}\]  
_Auxiliary Model Uncertainty_  
Uncertainties from auxiliary models (atmosphere, ocean, tides, glacial isostatic adjustment--GIA) are incorporated using autoregressive regression and weighted stochastic modeling, improving gravity signal estimation by adjusting for model-based error characteristics.