## 3 Results and Discussion  
### Acquired Results at the Country and Water District Levels  
The cross-correlation analysis of the GRACE TWSA and GPM IMERG precipitation was conducted at the pixel level, after resampling GRACE TWSA from its native \(1^{\circ}\) at 0.1\({}^{\circ}\), to match the resolution of GPM IMERG. Results of cross-correlation analysis at various lags and the country averaged correlation coefficients are shown in Figure 3. Thus, it can be seen that the highest correlation of 0.57 between the two variables was observed at a lag of 1 month, implying that it takes approximately 1 month for precipitation to be observed as changes in TWS. The second-highest correlation coefficient of 0.41 was observed at lag 2. Considering the inherent differences of those two datasets, but most importantly the difference in the spatial resolution, the correlation coefficient at lags 1 and 2 constitutes GPM IMERG, a promising auxiliary variable to downscale GRACE TWSA. Concerning the spatial patterns of correlation, it seems that the southwestern parts of the country demonstrate higher correlation values. The northeastern parts, where most agricultural activities and groundwater abstractions take place, have lower values  
Figure 3: Cross-correlation graphs between GPM precipitation and GRACE-derived TWSA. Data correspond to the GRACE-FO period from June 2018 to March 2021. The graph on the (**bottom**) demonstrates the country-averaged correlation at various lags. of correlation of GPM IMERG and GRACE data. Therefore, the pixel-based regression formula from Equation (1) was then arranged as follows:  
\[T_{1^{\circ}}^{t}=a_{1^{\circ}}+b_{1,1^{\circ}}P_{1^{\circ}}^{t-1}\,+b_{2,1^{ \circ}}P_{1^{\circ}}^{t-2} \tag{5}\]  
and it was then applied at the \(0.1^{\circ}\) spatial resolution, by using the GPM IMERG dataset at its native resolution, i.e., \(0.1^{\circ}\), and resampling the intercept a, and the precipitation coefficients \(b_{1}\) and \(b_{2}\) at the finer resolution using a bilinear interpolation:  
\[T_{0.1^{\circ}}^{t}=a_{0.1^{\circ}}+b_{1,0.1^{\circ}}P_{0.1^{\circ}}^{t-1}+b_{ 2,0.1^{\circ}}P_{0.1^{\circ}}^{t-2} \tag{6}\]  
Further, the residuals were computed from Equation (2) and were added to the TWSA from Equation (6), applying Equation (4), to obtain the downscaled and residually corrected TWSA. A comparison of the downscaling outcome before and after residual correction is demonstrated in Figure 4. The scatterplot on the left of Figure 4 shows the ensemble of GRACE TWSA data against the downscaled outcome before residual correction, indicating that most of the data points are above the y = x line (i.e., 1:1 line, black line on Figure 4) when TWSA is lower than \(-0.03\) and \(r^{2}\) is 0.59. The right-hand panel on the same figure demonstrates that after residual correction, the data points are distributed evenly above and below the 1:1 line and \(r^{2}\) is 0.87, considerably higher compared to the r\({}^{2}\) estimated before residual correction, implying improvement with the residual correction applied. Analogous improved results have been found by applying the residual correction approach while downscaling other remotely sensed data products, such as soil moisture [32]. To the best of our knowledge, it is the first time that such an approach has been tested in the downscaling of GRACE TWSA. Results of the downscaling outcome are provided in Figures 5 and 6. Figure 5 demonstrates the original GRACE-FO time series from June 2018 to March 2021, covering the study area. Figure 6 shows the downscaled outcome of TWSA in Greece. Please note that due to the time lags incorporated in Equations (5) and (6), the downscaled time series has a slightly different time coverage, i.e., from October 2018 to March 2021 (for details on GRACE-FO dates, please refer to GRACE dates in the Supplementary Table S1).