Spatiotemporal Downscaling of GRACE Total Water Storage Using Land Surface Model Outputs

: High spatiotemporal resolution of terrestrial total water storage plays a key role in assessing trends and availability of water resources. This study presents a two-step method for downscaling GRACE-derived total water storage anomaly (GRACE TWSA) from its original coarse spatiotemporal resolution (monthly, 3-degree spherical cap/~300 km) to a high resolution (daily, 5 km) through combining land surface model (LSM) simulated high spatiotemporal resolution terrestrial water storage anomaly (LSM TWSA). In the ﬁrst step, an iterative adjustment method based on the selfcalibration variance-component model (SCVCM) is used to spatially downscale the monthly GRACE TWSA to the high spatial resolution of the LSM TWSA. In the second step, the spatially downscaled monthly GRACE TWSA is further downscaled to the daily temporal resolution. By applying the method to downscale the coarse resolution GRACE TWSA from the Jet Propulsion Laboratory (JPL) mascon solution with the daily high spatial resolution (5 km) LSM TWSA from the Ecological Assimilation of Land and Climate Observations (EALCO) model, we evaluated the beneﬁt and effectiveness of the proposed method. The results show that the proposed method is capable to downscale GRACE TWSA spatiotemporally with reduced uncertainty. The downscaled GRACE TWSA are also evaluated through in-situ groundwater monitoring well observations and the results show a certain level agreement between the estimated and observed trends. priori input uncertainties of the GRACE TWSA. indicates EACLCO data Further analysis and quality control needed the quality of the model Similar pattern changes or improvement variations of the estimated uncertainties appear not spatially but also temporally. These results demonstrated well the role of the SCVCM as a quality control and analysis tool for the LSM EALCO TWS simulation and GRACE observations.


Introduction
Accurate knowledge on total or terrestrial water storage (TWS) and its spatiotemporal variability play a crucial role in the assessment of climate variation and water resource availability [1]. The accuracy of TWS simulated by land surface models (LSMs) at high spatial and temporal resolution is commonly limited by uncertainties in meteorological forcing, model parameter calibration, and land surface process representation [2][3][4][5]. To reduce the model uncertainties, researchers have been developing new LSMs that may account for interactions among the different components (e.g., soil moisture, groundwater, surface water, and snow) of TWS. For instance, some models introduced the impact of groundwater on surface processes [6][7][8][9] and others added surface runoff algorithms based on empirical relationships [6,10,11]. Including groundwater in an LSM enables a more complete simulation of the terrestrial water cycle, but it is also subject to large uncertainties since the hydrogeological conditions of aquifers are largely unknown in many circumstances. Due to differences in model physics and parameter values, estimates from various LSMs often exhibit large discrepancies even when models are run using identical forcing data [12][13][14][15]. These discrepancies impose large limitations in many hydrological applications and assessments. To improve quality of the TWS simulated by LSMs, researchers have been developing new models and algorithms that can combine or assimilate the temporal water storage changes or anomalies detected by the joint NASA/DLR Gravity Recovery and Climate Experiment (GRACE) satellite [16] and its follow-on missions. Since the GRACE system provides the only observation of the Earth's TWS which includes all water storage components, e.g., snow mass, surface water, vegetation water, soil moisture, and groundwater [16], it gives us a valuable calibration to the LSM TWS outputs. At the same time, we can also use it to improve LSM TWS accessibility through combination of both datasets.
Since the launch of the GRACE satellite in 2002, many researchers have been developing new models and algorithms for combining or assimilating GRACE observations and LSM TWS outputs. For example, Zaitchik et al. [17] developed an ensemble Kalman smoother (EnKS) to assimilate GRACE TWS into the NASA Catchment model in the Mississippi basin with promising results in 2008. Since then, different applications based on the EnKS were assessed and evaluated [4,18,19]. In recent years, a number of researchers reported their efforts to optimize GRACE data assimilation performance by investigating the spatial scales of assimilation [20,21], the representation of horizontal error covariance of models and GRACE observations [22][23][24], the choice of GRACE products and scaling factors [25,26], and the choice of data assimilation techniques [23,[26][27][28][29]. Most recently, the GRACE data assimilation was extended to multi-sensor land data and/or multi-mission satellite data assimilation [30][31][32][33]. For a more detailed review, readers may refer to [34]. Most of these studies showed that the GRACE data assimilation based on the EnKS improved model outputs by decreasing errors and increasing correlations between the simulated and observed variables.
The key challenge in developing a GRACE data combination or assimilation model with LSM TWS outputs is how to deal with their significant differences in spatial and temporal resolutions. The spatiotemporal resolution of GRACE TWS data available is about 300 km monthly while the spatiotemporal resolution of LSM TWS outputs can be much higher, for instance, typically 5 km daily. The GRACE data assimilation based on the EnKS works well with high spatial resolution LSM TWS data at the same monthly temporal resolutions. It is not necessary to downscale the GRACE TWS data resolution from its original coarse resolution to the fine resolution of LSM TWS data before the assimilation. A drawback of the EnKS is that it treats both GRACE and LSM TWS data as comparable observations by ignoring their differences such as missing groundwater and surface water storage (GSWS) in the LSM TWS data. In addition, the EnKS considers their uncertainties only by a given a priori value or representation of horizontal error covariance of models.
In recent years, some new downscaling methods using advanced statistical methods such as multivariate regression, artificial neural network, machine learning, and extreme gradient boosting techniques, etc. were developed for different applications. For example, Humphrey et al. [35] proposed an approach which statistically relates anomalies in atmospheric drivers to monthly GRACE anomalies. Yin et al. [36] tested statistical downscaling of GRACE-derived groundwater storage (GWS) using evapotranspiration (ET) data in the north China plain. Seyoum et al. [37] tried to downscale GRACE-derived total water storage anomaly (GRACE TWSA) into high-resolution groundwater level anomaly using machine learning-based models in a glacial aquifer system. Ahmed et al. [38] used an artificial neural network to forecast GRACE data over the African watersheds. Sun et al. [39] investigated reconstruction of GRACE total water storage through automated machine learning. Chen et al. [40] evaluated downscaling of GRACE-derived GWS based on the random forest model. Sahour et al. [41] developed optimum procedures to downscale GRACE Release-06 monthly mascon solutions based on statistical applications. Milewski et al. [42] trained a boosted regression tree model to statistically downscale GRACE TWSA to monthly 5 km groundwater level anomaly maps in the karstic upper Floridan aquifer (UFA) using multiple hydrologic datasets. Different from the EnKS that assimilates GRACE observations into existing process-based models dynamically, these new downscaling techniques based on statistical methods attempted to derive empirical relationships between GRACE observations and smaller scale quantities of interest. This typically involves training a non-parametric statistical model to predict local TWSA from GRACE TWSA and other environmental variables that impact groundwater storage [42][43][44]. A shortfall of these statistical downscaling methods is that few of them presented spatially downscaled GRACE data in a gridded format which resembles aquifers more closely and is perhaps the most useful for investigating spatial patterns in groundwater storage and abstraction [42].
In an attempt to downscale GRACE TWSA in a gridded format which closely resembles a high resolution TWS product simulated by a LSM, Zhong et al. [34] presented an iterative adjustment method based on the self-calibration variance-component model (SCVCM) for generating high spatial resolution total water storage anomaly (TWSA) through combination of GRACE observations and the TWSA simulated by the LSM Ecological Assimilation of Land and Climate Observations (EALCO). Using the SCVCM, the GRACE-derived total water storage anomaly (GRACE TWSA) can be spatially downscaled from its original coarse resolution (3-degree spherical cap/~300 km) to a high spatial resolution (5 km) grid [34].
After the GRACE TWSA is downscaled to a high spatial resolution grid spatially, another challenge is how to further downscale it from monthly to daily frequency temporally. For the very first time, we extended the iterative adjustment method based on the self-calibration variance-component model (SCVCM) to downscale spatiotemporally the monthly coarse resolution GRACE TWSA to daily 5 km resolution of the EALCO TWSA in this study. This new method combines the spatially downscaled GRACE TWSA with daily EALCO TWSA to generate a high spatiotemporal resolution TWSA product. In following, we first introduce the newly developed SCVCMs for the spatiotemporal downscaling of GRACE TWSA in Section 2. Then, a test case including data processing and model evaluation methods is presented in Section 3. After that, the results are evaluated and discussed in Section 4. Finally, we conclude this study with a summary in Section 5.

Problem Description
The objective of this study is to downscale GRACE TWSA from its original coarse resolution (monthly,~300 km) grid to the high spatiotemporal resolution (daily, 5 km) grid of the LSM EALCO TWSA through a data combination and adjustment model. Let G denote the GRACE TWSA observation of a mascon, i.e., the mascon average of a period of M days (about one month) and t be the time in days from 1 to M. For a grid point i (e.g., a 5 km × 5 km grid cell used for the LSM EALCO simulation within the mascon), let L i (t) be the EALCO TWSA at time t. Then, all observations available for a mascon and an observation period of about a month are where t ij denotes the day j for the grid point i, M is the number of the daily EALCO TWSA for the grid point i (i.e., the total number of days used for deriving the GRACE TWSA G, usually about a month), and N is the total grid cells within the mascon. Now, the objective is to build a model or to develop a process that can downscale the G from its original coarse resolution of mascon grid to a high spatiotemporal resolution grid of the EALCO TWSA L i (t) through data combination and adjustment of the observations G and L i (t).
For a very first attempt, we suggest a two-step method. In the first step, an iterative adjustment method based on the self-calibration variance-component model (SCVCM) is used to downscale the monthly GRACE TWSA G spatially to the high spatial resolution grid of the EALCO TWSA L i (t). In the second step, the spatially downscaled GRACE TWSA (SD TWSA) G i at each grid point i is further downscaled temporally to daily by using a second SCVCM for combination of the SD TWSA G i and the daily EALCO TWSA L ij .

The SCVCM for Spatial Downscaling of the GRACE TWSA
For the spatial downscaling of the GRACE TWSA from its coarse resolution mascon grid to the high-resolution grid of the EALCO TWSA, Zhong et al. [34] presented an iterative combination adjustment method based on the SCVCM [45]. The method takes the monthly GRACE TWSA G of a mascon and the monthly averaged EALCO TWSA at each grid point i, which includes soil water content, snow water equivalent, and plant water as input observations with unknown uncertainties, and then establishes an observation system to estimate the unknown TWSA at the high-resolution grid points. To consider the missing groundwater storage (GWS) and surface water storage (SWS), which are not included in the EALCO TWSA, or simply the systematic differences between the GRACE TWSA and the EALCO TWSA, an additional systematic parameter for each grid point is introduced to the system and estimated through an iterative adjustment process with a posteriori variance-component estimation technique. Following the process described in Zhong et al. [34], the final model outputs are spatially downscaled TWSA (SD TWSA) at the high spatial resolution grid. We take the same model and process directly as the first step of our downscaling processing in this study.

The SCVCM for Temporal Downscaling of the GRACE TWSA
After the GRACE TWSA is spatially downscaled to the high spatial resolution (5 km) grid in the first step, the next objective of this study is to downscale the SD TWSA G i at each grid point i further from its temporal resolution monthly to the daily resolution of the EALCO TWSA L ij . This can be realized by applying a similar method for the spatial downscaling of the GRACE TWSA in the first step.
Similar to the SCVCM development for the spatial downscaling, we have to take the following three particular challenges into account to build up a new SCVCM for temporal downscaling of the SD TWSA. First, the model must be able to distribute the information from a single monthly coarse-temporal observation onto the M daily finertemporal model elements properly. Second, the damped signal or leakage errors of the SD TWSA introduced through the monthly average smoother must be handled properly prior to or during the combination adjustment process. Third, the systematic differences between the EALCO TWSA and SD TWSA must be considered or compensated. Otherwise, the data combination may introduce errors to or reduce the accuracy of the model outputs.
Assume that the SD TWSA G i at a grid point i equals to the average of the M daily observations G ij which are to be determined, then our objective is to restore G ij from the monthly averaged signal G i . Similar to the gridded scale gain factor k i introduced to restore the leakage errors [46] or the damped signal from the mascon averaged GRACE TWSA G in the SCVCM modeling for the spatial downscaling [34], we define a scale gain factor k ij for each day to restore the leakage errors or the damped signal from the monthly averaged SD TWSA G i . Then we can restore G ij from G i as follows where ε ij are the misfit between the daily (i.e., unaveraged) G ij and the smoothed (i.e., averaged) monthly signal G i . From Equation (2), the gain factor k ij plays a role of distributing the single monthly coarse-temporal observation G i onto the M daily finer-temporal model elements G ij . If we think that the monthly average smoother used for averaging both the daily SD TWSA and the daily LSM EALCO TWSA for their monthly signal values are same, we can determine the scale gain factor k ij using the daily LSM EALCO TWSA L ij as follows Applying the scale gain factor k ij from Equation (3) to restore the damped signals or leakage error of the SD TWSA G i and introducing unknown parameter x ij to denote the spatiotemporally downscaled GRACE TWSA (SPTD TWSA) at time t ij , then we can have the following observation equation for each daily SPTD TWSA G ij = k ij G i at time t ij and Remote Sens. 2021, 13, 900 5 of 19 the monthly averaged observation G i = G i (as a constraint condition they should meet after the downscaling): Here,σ 2 G i is a variance estimate of the monthly averaged signal G i , which is actually the spatially downscaled SD TWSA G i resulted from the first step, i.e.,σ 2 For the daily EALCO TWSA L ij , we can establish the following observation equations by introducing an additional unknown parameter s ij to model the missing groundwater and surface water storage anomaly (GSWSA) component in L ij , L ij is a variance estimate of the daily EALCO TWSA L ij . If we treat the x ij in Equations (4) and (5) as unknown model parameters and the s ij as additional stochastic systematic parameters to be determined from the observations G ij and L ij , Equations (4) and (5) form a typical SCVCM with two different observation types (i.e., G ij and L ij ) and one group of stochastic systematic parameters (i.e., s ij ) [45].
By defining following unknown parameter vectors, we obtain the following model coefficient matrices and observation vectors as well as the functional model of the SCVCM L = Aδx + Hδs + ε (8) or equivalently as the error equations Following the same solution method used to the SCVCM for the spatial downscaling [34], we obtain the final estimates of the parameters δx, and δs and their variance and covariance estimates as follows: o is an estimate of the unit weight variancê From Equation (6),x = x o + δx is an estimate vector of the unknown spatiotemporally downscaled TWSA x, and the corresponding diagonal variances in the covariance matrix (11) are their accuracy or uncertainty estimates. Similarly,ŝ = s o + δŝ is an estimate vector of the GSWSA parameters s and the corresponding diagonal variances in the covariance matrix (11) are their accuracy or uncertainty estimates.

Study Area
To verify the SCVCM's performance, we need independent TWSA observations at the high spatiotemporal resolution grid points for a comparison. However, such TWSA observations are not available. Alternatively, assuming that the surface water storage anomaly (SWSA) is ignorable in a selected dry area, then we may think that the difference between the GRACE TWSA and the EALCO TWSA, i.e., the GSWSA, is caused mainly by the groundwater storage anomaly (GWSA). Hence, if we can determine the GWSA from in-situ groundwater monitoring well (GMW) observations, then the estimated GSWSA parameter,ŝ i , should be comparable with the observed GWSA at the same location. Therefore, for the model evaluation, we selected a study region located in the dry region of Canada's Prairie where surface water is a less significant factor in the TWS than other components. The test area covers 7 full and 2 partial mascons in the province of Alberta, Canada ( Figure 1).
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 20 Figure 1. Map for the test region. There are 32 unconfined or semi-unconfined groundwater monitoring wells (dots) in the 9 mascons (numbered from 1 to 9) over the study region located in the province of Alberta, Canada.

Data and Data Processing
Same as the GRACE TWSA data used for Zhong et al. [34], we take the monthly Figure 1. Map for the test region. There are 32 unconfined or semi-unconfined groundwater monitoring wells (dots) in the 9 mascons (numbered from 1 to 9) over the study region located in the province of Alberta, Canada.

GRACE TWSA Data
Same as the GRACE TWSA data used for Zhong et al. [34], we take the monthly GRACE TWS product (RL06 V1.0) from the JPL mascon solution [47][48][49] as one of our input observations for this study. In order to make it comparable with the EALCO TWSA data, a new baseline value, the mean from 5 April 2002 to 31 December 2016, is deducted from the raw data downloaded from the JPL TWS product website (https://grace.jpl.nasa. gov/data/get-data/jpl_global_mascons/, last access: on 28 July 2020). That means that the input GRACE TWSA observations for this study are the monthly TWSA derived from the JPL mascon solution by deducting the new baseline value. Its temporal resolution is monthly, and its spatial resolution is the JPL mascon grid (~300 km).

EALCO TWSA Data
The EALCO TWSA data used for this study are also similar to the data used in Zhong et al. [34]. Details about the LSM EALCO refer to [50][51][52][53][54][55][56][57][58][59][60][61][62]. For this study, we first calculated the daily EALCO TWS from the 30-minute simulations for soil water, snow water equivalent, and plant water for each EALCO grid cell at 5 km × 5 km for the period corresponding to GRACE observations from 5 April 2002 to 31 December 2016. The daily EALCO TWSA is obtained by deducting its mean (baseline) values from 5 April 2002 to 31 December 2016. The EALCO simulation doesn't include the groundwater and surface water. Hence, the input EALCO TWSA observations used in this study only include soil water, snow water equivalent, and plant water.

GWSA from GMW Observations
Similar to the GWSA data used in Zhong et al. [34], we selected 32 unconfined or semi-confined GMWs (see the black and red dots in Figure 1) in the study area and obtained their daily groundwater water level observations from the GMW network of the Alberta Environment and Parks, Government of Alberta (http://environment.alberta.ca/apps/ GOWN/#, last access: 14 November 2020). The 32 GMWs are selected due to their longest monitoring period and their greatest overlap with the GRACE data. The minimum number of days with available observations is 3620 of 5385 days in total (67.2%). Among them, 28 of the 32 GMWs (87.5%) have more than 3993 of 5385 days (74.2%) with observations. More information about these GMWs is available from the website of the Alberta Environment and Parks, Government of Alberta.
By matching the GMW observations for the same periods used for deriving the GRACE TWSA and then deducting its baseline value, we obtained the daily GMW level height variations. As the GMW level height observations are not directly comparable with the GRACE TWSA [63], they are multiplied by a specific yield of 5% for all 32 GMWs to convert them to GWSA. The adopted specific yield 5% is a close to the mean value derived by performing a least-squares fit between the daily GMW level height variations δh(t ij ) and the corresponding estimated parameter value s(t ij ) for the individual GMW i. That is, where ε s ij are the misfit between s(t ij ) and δh(t ij ). M i is the number of the GMW observations available at time t ij in the grid cell i.

Evaluation Methods
As discussed in Zhong et al. [34], the determination of the GWSA from in-situ GMW level observations is a technical challenge of the model evaluation because it depends on accurate specific yield information required to convert the groundwater level changes to GWSA [63]. The aquifer settings and human activities such as withdrawal, injection, mining, etc. around a GMW can largely affect the accuracy of the local GWSA determination, too. In addition, the simulation errors of the input EALCO TWSA and seasonal surface water influences, etc. may cause significant differences of the estimated GSWSA for a representation of the actual GWSA. Therefore, it is difficult and not practical to compare the estimated GSWSA against the observed GWSA from GMW observations qualitatively for the model evaluation. A more practical way is to quantitatively compare the trend changes between the time series of the estimated GSWSA and the observed GWSA at the GMW's location [34]. For convenience, we name the estimated GSWSA as EGSWSA and the GSWA derived from in-situ GMW level observations as observed GWSA (OGWSA) in following model evaluation respectively.
For evaluation of the trend changes agreement between two time series X = x 1 , x 2 , . . . , x n , Y = y 1 , y 2 , . . . , y n , simple and effective comparison indicators are the root mean square difference (RMSD) and Person's correlation coefficient (PCC) [64]: Pearson's correlation where x i and y i are the X and the Y time series respectively, n is the total number of available observations in the time series, and x and y are the mean of x i and y i , respectively. As pointed out in Zhong et al. [34], we use the RMSD to measure the agreement goodness of the EGSWSA time series against the OGWSA time series. Although its definition is same as the root mean square error (RMSE) used to measure the error of a model in predicting quantitative data, its meaning is different. We cannot see it as the RMSE because none of the EGSWSA and the OGSWA can be treated as a true value or standard for the measurement error.
Following the same model evaluation methods used in Zhong et al. [34], we can also compare the EGSWSA and the OGSWA not only in temporal domain but also in spatial domain if we see the EGSWSAs at different GMW's locations at a same observation time as a data series. In temporal domain, we compare the EGSWSA and the OGSWA time series for a GMW along the entire observation period. In spatial domain, we compare the EGSWSA and the OGSWA for different GMWs at a specific observation time. In addition, we can evaluate the model output through the estimated uncertainty by the model self too. If the estimated uncertainty of the SPTD TWSA from Equation (11) is comparable with the a priori input uncertainty of the GRACE TWSA and the EALCO TWSA, this means the model outputs are reasonable. From the estimated uncertainty information, we can also infer the quality of input observations and how well they fit with each other at an observation time within a mascon. For the time series of the SPTD TWSA, we can also compare their uncertainties at different observation times.
Since the main factors that influence the quality of the EALCO TWSA include the disparity in forcing data, modelling methods, and parameters used for the EALCO simulation, and the uncertainty of the resultant EALCO TWSA is the major error source that impacts the uncertainty of the SPTD TWSA, a posteriori uncertainty analysis of the SPTD TWSA can also help us better inform accuracy of the EALCO TWSA and identify where and when the input data needs to be improved for a better model output.

Downscaled TWSA
For a visual inspection on the downscaled GRACE TWSA, we plotted the input data, i.e., the monthly GRACE TWSA and the daily EALCO TWSA, against the output data, i.e., the spatiotemporally downscaled daily TWSA (SPTD TWSA) and the estimated daily EGSWSA in Figure 2.

Downscaled TWSA
For a visual inspection on the downscaled GRACE TWSA, we plotted th i.e., the monthly GRACE TWSA and the daily EALCO TWSA, against the ou the spatiotemporally downscaled daily TWSA (SPTD TWSA) and the estima SWSA in Figure 2.  From Figure 2, we can see that the suggested method is capable to downscale the GRACE TWSA from its original coarse-resolution (monthly~300 km) to the high spatiotemporal resolution (daily 5 km) grid of the EALCO TWSA. The spatial patterns of the SPTD TWSA are similar to the EALCO TWSA. This indicates that the SPTD TWSA is just an adjusted EALCO TWSA by the EGSWSA. For a quick validation, it is clear that the SPTD TWSA had been changing from 2002 to 2016 and this change is especially significant in the southwest region marked as mascon 3, 6, and 8. The reasons for these changes are multiplex. Some changes can be explained as the effects of extreme droughts and floods occurred in this region, for instances, 2013 Alberta Floods [65][66][67] caused the TWSA increases in Southern Alberta in June 2013. In comparison to the TWSA's changes, the EGSWSA are relatively small and stable (Figure 2d). Some visible changes happened in the region marked as mascon 4, 7, and 8 from 2013 to 2016. Those relatively small changes of the EGSWSA are consistent with the OGWSA derived from the in-situ GMW level observations (more details in Section 4.3).

Quantitively Comparison Results
To assess the model effectiveness for spatiotemporally downscaling GRACE TWSA, we compared three time series of the GSWSA derived from different model outputs against the 32 GMWs observations. The first time series is derived by differencing the EALCO TWSA from the original GRACE TWSA without any downscaling process. We named this time series as EGSWSA1. The second time series is derived by differencing the EALCO TWSA from the spatially downscaled GRACE TWSA, i.e., the SD TWSA without conducting the temporal downscaling process. We named this time series as EGSWSA2. The third time series are derived by conducting both the spatial and temporal downscaling processes described in Section 3. We named it as EGSWSA3. Figure 3 shows comparison plots of the RMSD (a) and the PCC (b) in temporal domain for the 32 GMWs. It is clear that the EGSWSA3 derived from the spatiotemporal downscaling show significantly smaller RMSD than EGSWSA1 and EGSWSA2, which are derived from the original GRACE TWSA and the spatially downscaled GRACE TWSA, i.e., SD TWSA, respectively. The comparison between the EGSWSA1 and EGSWSA2 shows that the EGSWSA2 has smaller RMSD than the EGSWSA1. However, the PCC of all three times series did not show significant differences. They all show a very weak correlation level. These results tell us that both the spatial and temporal downscaling processes are helpful in reducing the RMSD between the EGSWSA and the OGWSA. They functioned like a smoothing filter in adjusting the observation errors of the EALCO TWSA and the GRACE TWSA.

Visual Comparison Results
For a visual comparison, we first plotted scatter plots of all pair points for the three time-series of EGSWSA1, EGSWSA2, and EGSWSA3 against the OGSWA derived from the 32 GMWs in Figure 5, respectively. Figure 5 shows that the scatter plot of the EG-SWSA3 looks improved and gives a smaller RMSD than the EGSWSA1 and EGSWSA2. The PCC looks small for all three time-series and no improvement is made by both spatial and temporal downscaling processes. This indicates that the iterative adjustment used for downscaling process did improve the RMSD but not the PCC. This is consistent with the results of the quantitative evaluation above.

Visual Comparison Results
For a visual comparison, we first plotted scatter plots of all pair points for the three time-series of EGSWSA1, EGSWSA2, and EGSWSA3 against the OGSWA derived from the 32 GMWs in Figure 5, respectively. Figure 5 shows that the scatter plot of the EGSWSA3 looks improved and gives a smaller RMSD than the EGSWSA1 and EGSWSA2. The PCC looks small for all three time-series and no improvement is made by both spatial and temporal downscaling processes. This indicates that the iterative adjustment used for downscaling process did improve the RMSD but not the PCC. This is consistent with the results of the quantitative evaluation above. For visual comparisons of the trend changes of the EGSWSAs and the OGWS series, we plotted the EGSWSA1, EGSWSA2, and EGSWSA3 against the OGWSA selected 8 GMWs in Figure 6. These 8 GMWs are located in the different mascons m as 1, 2 ,4 ,5 ,6 ,7 ,8, and 9 (see the red dots in Figure 1). There is no GWMs availabl mascon 3.
The direct trend changes comparisons of the EGSWSAs and the OGWSA tim in the left Colum of Figure 6 look very noisy because of significant seasonal infl simulation, or observation errors, etc. For a better trend comparison and analysis, w a 1-year (365 days) moving average smooth filter to deseasonalize some seasona ences and smooth out some noise signals or observation errors. The smoothed tim are compared in the right Colum of Figure 6. For visual comparisons of the trend changes of the EGSWSAs and the OGWSA time series, we plotted the EGSWSA1, EGSWSA2, and EGSWSA3 against the OGWSA for the selected 8 GMWs in Figure 6. These 8 GMWs are located in the different mascons marked as 1, 2, 4, 5, 6, 7, 8, and 9 (see the red dots in Figure 1). There is no GWMs available in the mascon 3. For visual comparisons of the trend changes of the EGSWSAs and the OGWSA time series, we plotted the EGSWSA1, EGSWSA2, and EGSWSA3 against the OGWSA for the selected 8 GMWs in Figure 6. These 8 GMWs are located in the different mascons marked as 1, 2 ,4 ,5 ,6 ,7 ,8, and 9 (see the red dots in Figure 1). There is no GWMs available in the mascon 3.
The direct trend changes comparisons of the EGSWSAs and the OGWSA time series in the left Colum of Figure 6 look very noisy because of significant seasonal influences, simulation, or observation errors, etc. For a better trend comparison and analysis, we used a 1-year (365 days) moving average smooth filter to deseasonalize some seasonal influences and smooth out some noise signals or observation errors. The smoothed time series are compared in the right Colum of Figure 6.  From Figure 6, we can see that the variation ranges of three unfiltered EGSWSAs are much bigger than the OGWSA. Their amplitude of individual observation looks not comparable. The reasons for the big differences may be multifaceted. As mentioned in Section 4.3, it might be caused by inaccurate specific yield used to convert the GMW level height changes to OGWSA or possible simulation errors in the input EALCO TWSA or significant seasonal SWSA caused by something such as floods, etc. This is why it is still not practical to compare the amplitudes of the EGSWSA and the OGWSA qualitatively for the model evaluation. However, if we compare their long-term trend changes, the EGSWSA3 shows a better trend change agreement than the EGSWSA1 and EGSWSA2 for all selected GMWs except for very few GMWs, especially after deseasonalizing by the 1-year (365 days) moving average filter. This conclusion is valid for most of the 32 GMWs. Only a few of the GMWs showed some reversed results. An exception among the 8 GMWs is the well W229 in mascon 7, which shows that the EGSWSA2 is closer to the OGWSA. The well W287 and W117, which show negative correlations in Figure 3, have similar results. The reason for this result is not clear yet. Possibly, these exceptional wells have a different specific yield from others or its OGWSA may not adequately represent either the groundwater response or the EGSWSA because of unknown aquifer settings and human activities such as withdrawal, injection, mining, etc.

Uncertainty Comparison Results
To evaluate the model uncertainty, we compared the a priori input uncertainties of the monthly GRACE TWSA against the uncertainties of the downscaled daily SPTD TWSA in Figure 7.  The direct trend changes comparisons of the EGSWSAs and the OGWSA time series in the left Colum of Figure 6 look very noisy because of significant seasonal influences, simulation, or observation errors, etc. For a better trend comparison and analysis, we used a 1-year (365 days) moving average smooth filter to deseasonalize some seasonal influences and smooth out some noise signals or observation errors. The smoothed time series are compared in the right Colum of Figure 6.
From Figure 6, we can see that the variation ranges of three unfiltered EGSWSAs are much bigger than the OGWSA. Their amplitude of individual observation looks not comparable. The reasons for the big differences may be multifaceted. As mentioned in Section 4.3, it might be caused by inaccurate specific yield used to convert the GMW level height changes to OGWSA or possible simulation errors in the input EALCO TWSA or significant seasonal SWSA caused by something such as floods, etc. This is why it is still not practical to compare the amplitudes of the EGSWSA and the OGWSA qualitatively for the model evaluation. However, if we compare their long-term trend changes, the EGSWSA3 shows a better trend change agreement than the EGSWSA1 and EGSWSA2 for all selected GMWs except for very few GMWs, especially after deseasonalizing by the 1-year (365 days) moving average filter. This conclusion is valid for most of the 32 GMWs. Only a few of the GMWs showed some reversed results. An exception among the 8 GMWs is the well W229 in mascon 7, which shows that the EGSWSA2 is closer to the OGWSA. The well W287 and W117, which show negative correlations in Figure 3, have similar results. The reason for this result is not clear yet. Possibly, these exceptional wells have a different specific yield from others or its OGWSA may not adequately represent either the groundwater response or the EGSWSA because of unknown aquifer settings and human activities such as withdrawal, injection, mining, etc.

Uncertainty Comparison Results
To evaluate the model uncertainty, we compared the a priori input uncertainties of the monthly GRACE TWSA against the uncertainties of the downscaled daily SPTD TWSA in Figure 7. Figure 7 shows that the uncertainties of the SPTD TWSA estimated by the SCVCM are usually improved for the most times and the most parts in the test region in comparison to the a priori input uncertainties of the original GRACE TWSA. This result indicates that the GRACE TWSA and the EALCO TWSA fit each other well for the most times and the most parts in the test region and the proposed method is beneficial to reduce the uncertainties. From Figure 7, we can also see that the uncertainties of the SPTD TWSA are enlarged sometimes at some parts to a significant extent. For instance, in May and June of 2006, the estimated uncertainties for the region of mascon 6 are significantly larger than the a priori input uncertainties of the GRACE TWSA. This indicates very possibly that there are significant differences between the GRACE TWSA and the EACLCO TWSA input data during that period in the area in question. Further analysis and quality control are needed to improve the quality of the model outputs. Similar pattern changes or improvement variations of the estimated uncertainties appear not only spatially but also temporally. These results demonstrated well the role of the SCVCM as a quality control and analysis tool for the LSM EALCO TWS simulation and GRACE observations. of the GMWs showed some reversed results. An exception among the 8 GMWs is the w W229 in mascon 7, which shows that the EGSWSA2 is closer to the OGWSA. The w W287 and W117, which show negative correlations in Figure 3, have similar results. T reason for this result is not clear yet. Possibly, these exceptional wells have a differe specific yield from others or its OGWSA may not adequately represent either the groun water response or the EGSWSA because of unknown aquifer settings and human activit such as withdrawal, injection, mining, etc.

Uncertainty Comparison Results
To evaluate the model uncertainty, we compared the a priori input uncertainties the monthly GRACE TWSA against the uncertainties of the downscaled daily SPT TWSA in Figure 7.  From our test computations with different a priori variances (σ 2 E , σ 2 G and σ 2 S ) assigned to the observation L ij , G ij , and the apriori value s o (zero as virtual observation of the parameter s i ) in Equation (9), we found that the uncertainties of the model outputs vary to some extent. These a priori variance settings will impact the convergence speed of the iterative adjustment too. For a fast convergence, equal weighting with Q EE =σ 2 E I, Q GG =σ 2 G I and Q SS =σ 2 S I is recommended. For improvement of the model outputs, more advanced weighting options based on the improved accuracy estimation of the GRACE TWSA (e.g., [46]) and the EALCO TWSA may be explored in future studies.

Conclusions
In this study, we presented a two-step method for downscaling GRACE-derived total water storage anomaly (GRACE TWSA) from its original coarse resolution (monthly, 3degree spherical cap/~300 km) to a higher spatio-temporal resolution (daily, 5 km) using the land surface model (LSM) simulated high spatiotemporal resolution terrestrial water storage anomaly (LSM TWSA). In the first step, an iterative adjustment method based on the self-calibration variance-component model (SCVCM) is used to downscale the monthly GRACE TWSA spatially to the higher resolution grid used for the LSM simulation. In the second step, the spatially downscaled monthly GRACE TWSA is further downscaled using the daily LSM TWSA to generate finer temporal resolution TWSA. By applying the method to produce high spatiotemporal resolution TWSA through combination of the monthly coarse resolution (~300 km) GRACE TWSA from the JPL (Jet Propulsion Laboratory) mascon solution and the daily high resolution (5 km) LSM TWSA from the Ecological Assimilation of Land and Climate Observations (EALCO) model, we evaluated its benefit and effectiveness by comparing the estimated groundwater and surface water storage anomaly (EGSWSA) against the observed groundwater storage anomaly (OGWSA) by assuming that the surface water in a dry test region is ignorable. From the evaluation results, we can draw following conclusions.
First, the proposed method is capable to downscale GRACE-derived total water storage anomaly (GRACE TWSA) from its original monthly coarse resolution (monthly, 300 km) to the higher spatiotemporal resolution (daily, 5 km) by combining or integrating the high resolution EALCO TWSA. The evaluation results show that the EGSWSA by the proposed method agrees better with the OGWSA in both temporal and spatial domain than the GSWSA derived from the original GRACE TWSA and the spatially downscaled GRACE TWSA. These results indicate that the proposed method makes sense and is beneficial.
Second, in comparison to most of the developed GRACE TWSA combination or assimilation methods, the proposed method will deliver not only high spatiotemporal resolution TWSA and GSWSA estimate at each grid point but also their internal accuracy, i.e., estimated variance or uncertainty. From the uncertainty information, we can inform potential quality issues in the LSM EALCO TWSA and/or the GRACE TWSA, and how well they fit with each other in both temporal and spatial domain within a mascon. This means the proposed method can be used as a quality control and analysis tool for the LSM EALCO TWS simulation and the GRACE observations.
From this study, we recognized that it is still a challenge to evaluate the model's performance qualitatively because of the lack of accurate and comparable TWSA observations. We used the groundwater monitoring well (GMW) level observations to evaluate the model's performance quantitatively through comparison of the trend changes of the estimated GSWSA and the observed GWSA. However, this evaluation is only quantitative and not qualitative. It is also based on the assumption that the surface water in a dry test region is ignorable. To improve the quality of the final TWSA product, the key is to improve the quality of the input observations, i.e., either enhancing the quality of the LSM EALCO TWS simulation or introducing more observations such as the observed GWSA and SWSA, the simulated TWSA by different LSMs, etc. directly into the SCVCM to constraint the LSM simulation errors. The proposed SCVCM offers such a model for combining different observations. Funding: This study is a part of the research and development project "The Cumulative Effects Project and Climate Change Geoscience Program" funded by Natural Resources Canada.