Data Assimilation of Terrestrial Water Storage Observations to Estimate Precipitation Fluxes: A Synthetic Experiment

: The Gravity Recovery and Climate Experiment (GRACE) mission and its Follow-On (GRACE-FO) mission provide unprecedented observations of terrestrial water storage (TWS) dynamics at basin to continental scales. Established GRACE data assimilation techniques directly adjust the simulated water storage components to improve the estimation of groundwater, streamﬂow, and snow water equivalent. Such techniques artiﬁcially add/subtract water to/from prognostic variables, thus upsetting the simulated water balance. To overcome this limitation, we propose and test an alternative assimilation scheme in which precipitation ﬂuxes are adjusted to achieve the desired changes in simulated TWS. Using a synthetic data assimilation experiment, we show that the scheme improves performance skill in precipitation estimates in general, but that it is more robust for snowfall than for rainfall, and it fails in certain regions with strong horizontal gradients in precipitation. The results demonstrate that assimilation of TWS observations can help correct (adjust) the model’s precipitation forcing and, in turn, enhance model estimates of TWS, snow mass, soil moisture, runoff, and evaporation. A key limitation of the approach is the assumption that all errors in TWS originate from errors in precipitation. Nevertheless, the proposed approach produces more consistent improvements in simulated runoff than the established GRACE data assimilation techniques.


Introduction and Research Questions
The Gravity Recovery and Climate Experiment (GRACE) mission and its Follow-On (GRACE-FO) mission (hereinafter "GRACE" will be used to refer to both GRACE and GRACE-FO) have been providing observations of changes in the Earth's water storage since 2002 ( [1,2]). These mass change observations are temporally averaged, typically monthly, measurements of anomalies (departures from the long term mean at a given location) of the total mass of water stored on and beneath the land surface. The horizontal resolution is about 300 km at mid-latitudes ( [3]). The measured quantity is often referred to as terrestrial water storage (TWS), which encompasses snow, groundwater, soil moisture, and surface water in lakes, wetlands, and rivers. Alternatively, TWS can be seen as the net gain (or loss) of water storage from the balance of precipitation, evaporation, and runoff fluxes.
The assimilation of TWS data into a land surface model has also been used to downscale the coarse resolution of GRACE observations and to separate the TWS components (e.g., [4][5][6]). Data assimilation is a statistical technique that enables the optimal combination of model estimates and observational data and offers the potential to generate gridded information that is more accurate than that simulated by a model alone (e.g., [7,8]). For TWS, data assimilation has proven valuable for a myriad of hydrology applications such as drought monitoring (e.g., [9][10][11][12]) and forecasting ( [13]), groundwater depletion estimation (e.g., [14,15]), assessment of flood potential ( [16]), streamflow and lake storage estimation (e.g., [17,18]), identification of human driven hydrological mechanisms (e.g., [19,20]) and improving the water budget analyses (e.g., [21][22][23][24]). A key motivation for land data assimilation is to mitigate error and uncertainty in the parameterizations, parameters and forcing data used by land surface models to simulate complex water and energy cycle processes.
Precipitation drives the water balance at the land surface and is therefore essential to the accurate model representation of TWS, runoff and streamflow generation. Yet, obtaining accurate global estimates of precipitation remains a major challenge ( [25][26][27][28][29][30]). Data assimilation can help adjusting for precipitation uncertainties within models (e.g., [31,32]). As an example, the work by [33] illustrates several soil moisture data assimilation techniques to improve rainfall-runoff model estimation. In their work, the best streamflow estimates are attained when the assimilation technique adjusts precipitation fluxes [34] along with the model soil moisture states.
Observations based datasets show a good agreement between precipitation and TWS ( [35]), especially in terms of correlations in the high latitudes of Siberia, Europe, Russia, and North America ( [36]). Because of this agreement, two recent studies used GRACE TWS to estimate errors due to precipitation under-catch ( [37,38]) at high latitudes (or snowy regions), where precipitation is particularly uncertain, climatic conditions limit the influence of evaporation and runoff on TWS ( [37,39]), and the relationship between precipitation and TWS is stronger than at mid-latitudes. The main hypothesis behind these studies is that TWS uncertainties are, to first order, dominated by errors in precipitation.
Using the same hypothesis, in this paper we explore if precipitation errors can be estimated from the assimilation of TWS observations. In this first, exploratory study we assume that the main and only source of TWS uncertainty is error in the precipitation forcing data. We perform a synthetic assimilation experiment where GRACE-like TWS observations are used to correct errors in precipitation, thus leading to a better precipitation input to the land surface model. If, for example, the model's precipitation input at a specific location and time was less than the synthetic truth, the model would underestimate TWS. Our proposed assimilation scheme would then use the difference between the synthetic TWS observations and the modeled (underestimated) TWS to increase the estimated precipitation and, in a second repeat model simulation, generate improved TWS estimates. In reality, errors in simulated TWS and runoff could also be caused by errors in the model parameterization of evapotranspiration or the partitioning into infiltration and runoff ( [25,33,40]). Thus, the assimilation of actual GRACE TWS observations into the model may result in further uncertainties. Nonetheless, the analysis presented here is a first step towards improving precipitation and runoff estimates along with TWS estimates through the assimilation of TWS observations. A more in-depth investigation of other sources of errors and the design of more sophisticated assimilation approaches is left for future work.
The proposed work builds upon previous TWS assimilation methods (e.g., [4,5,9,41]). These efforts directly adjust the model's TWS prognostic state variables such as groundwater, soil moisture, and snow; model parameters and input meteorological forcings do not change. This implies that the direct adjustments (or increments) of the TWS prognostic variables disrupt the simulated internal model mass balance. Differently from previous TWS assimilation, here we aim to calculate and apply increments directly to the precipitation errors that is used in the land model to generate the predicted TWS. The main advantages of adjusting precipitation errors rather than prognostic model state estimates is that model fluxes such as runoff and evaporation are also automatically modified directly in response to the adjusted precipitation inputs and according to model physics, thus respecting the internal model mass balance.
The main goal of this study is to explore if the accuracy of precipitation information can be increased by minimizing the differences between modeled and observed TWS. This study focuses on TWS data assimilation, broadly similar to previous experiments but with some technical differences that highlight another way of leveraging GRACE observations to further improve land surface model output. The specific research questions we aim to address are: (1) can precipitation errors be inferred from the proposed assimilation scheme? (2) Can hydrologic states and fluxes also be improved?

Models and Data
All experiments conducted in this study used the Catchment Land Surface Model (CLSM), developed by [42]. CLSM has been used in several TWS assimilation studies (e.g., [4,5,20,43]) because it is one of the few global land surface models capable of representing shallow groundwater storage changes. Briefly, CLSM solves the surface water and energy balance in response to surface meteorological forcing. The model first determines the equilibrium profile soil moisture from the shallow groundwater table to the surface based on the total subsurface water storage deficit (referred to as the "catchment deficit"). The profile soil moisture content (prmc) is then determined from this deficit and from two additional excess moisture variables that describe deviations from the equilibrium profile in a 0-5 cm surface layer and a 0-100 cm root zone layer. Moreover, CLSM horizontally redistributes soil water within each hydrological catchment based on the statistics of the catchment topography. The model's estimate of TWS is the sum of soil moisture, shallow groundwater, snow and canopy interception water.
The model domain used here included all land regions north of 20 degrees north latitude. All experiments used a 36 km resolution Equal-Area Scalable Earth (EASE) version 2 grid ( [44]) and a model time step of 7.5 min. Our period of investigation was July 2003 to July 2016, which encompassed most of the period of the GRACE mission. We used the model parameters defined for Version 4 of the SMAP Level-4 Soil Moisture product ( [45]). Surface meteorological forcing data from the Modern Era Retrospective analysis for Research and Applications Version 2 (MERRA-2, [46]) were used as input to CLSM. The precipitation forcing used here was the precipitation generated by the Atmospheric Circulation Model (AGCM) within the MERRA-2 system following the assimilation of atmospheric variables (e.g., temperature, humidity, among others). That is, the MERRA-2 precipitation used here was not corrected with gauge-based data ( [47]).

Methods
We developed a synthetic data assimilation experiment to test our research questions, which allowed us to test the newly proposed technique in a controlled environment where all sources of errors and their characteristics are known. Four steps were involved in the synthetic experiment: (1) generate a set of synthetic observations using an independent model run (Section 3.1); (2) perform the model-only (or open-loop) ensemble simulation (Section 3.2); (3) assimilate the synthetically generated observations (Section 3.3); and (4) assess the performance of the proposed assimilation framework with respect to the synthetic truth (Section 4).

Synthetic Truth
Several methods are available to produce a synthetic (imaginary) truth. Following a so-called "identical twin experiment" approach, we adopted the same model and the same forcing datasets used in the open-loop and data assimilation experiments (MERRA-2) except for the precipitation, which was instead chosen from an 8-year shifted period. That is, we used MERRA-2 precipitation data from 1995 to 2008 as the imaginary synthetic precipitation (P True ) for 2003 to 2016, as illustrated in Figure 1. The rationale behind the 8-year shift was to construct a synthetic precipitation dataset that differs from the one used in the ensemble and assimilation runs (Sections 3.2 and 3.3). The resulting "truth" and simulated TWS, soil moisture, and SWE data were unbiased, exhibit consistent seasonality, and did not include patterns that were not modeled (e.g., groundwater depletion). This approach restricted all errors in the system to errors in precipitation, so that we could directly address the hypothesis tested in this study. Using this forcing setup we generated a suite of synthetic prognostic and diagnostic estimates (soil moisture, terrestrial water storage, runoff, etc.). As an example, Figure 1 shows the synthetic TWS obtained from the deterministic model run that uses synthetic precipitation input. We used the suite of synthetic prognostic and diagnostic estimates to verify if the assimilation of synthetic TWS observations could recover other aspects of the hydrology, beside constraining the precipitation fluxes.

Generation of Synthetic TWS Observations
The suite of synthetic diagnostic variables includes estimates of TWS at the model spatial resolutions (i.e., 36 km). From these, we generated synthetic TWS observations (Z) via temporal and spatial averages. Temporal-and spatial-averaging and re-gridding were necessary to reproduce GRACE-like synthetic observations. Temporally, we aggregated TWS values to 1-month 1-month ( Figure 2a) and spatially, apply a Gaussian smoother with a 300-km half-width on the 36-km modeled TWS ( Figure 2b) and finally sub-sampled the smoothed TWS values at 3 × 3 degrees spacing ( Figure 2c) to approximately match the spatial resolution of GRACE operational products (e.g., [2,3,48,49]). Note, that GRACEderived TWS retrievals are usually expressed in terms of anomalies with respect to a specific reference time. Prior to data assimilation, we converted the TWS anomalies into absolute TWS values based on a model estimate of TWS at the reference time. This common practice is equivalent to confronting the anomalies of the model estimates with the original anomaly TWS retrievals in that only TWS anomaly information is assimilated. In the proposed synthetic assimilation approach (Section 3.3), sources of errors were known and only attributed to precipitation fluxes, thus we arbitrarily used a measurement error equal to 5 mm.

Ensemble Simulation
In general, any land surface model determines TWS from the mass balance equation: This means that for a specific model time step t the uncertainties associated with the estimation of TWS depended upon errors/uncertainties in the precipitation (P), evaporation (E), and runoff (R) fluxes on the right-hand side of Equation (1). In this study, we focused on the error dependencies between TWS and precipitation fluxes and assume dthat the only source of modeling errors was inaccurate precipitation input. While precipitation errors were likely to dominate in GRACE data assimilation, this assumption was nevertheless a simplification that was not fully representative of reality. The assumption was, however, useful for testing the feasibility of the proposed assimilation method in our synthetic experiment. In our ensemble-based data assimilation system, the precipitation errors could be represented through the generation of an ensemble of equiprobable model realizations. We generated an ensemble of 24 members ( [5]) of possible precipitation realizations using a multiplicative error scheme as follows: where the subscript j = 1, 2, . . . , 24 indicates the ensemble member, P M2,m is the nominal precipitation forcing (in this case, MERRA-2 precipitation, Section 2), and b j,m is the precipitation error coefficient for a given ensemble member j and month m. These error coefficients were sampled at the beginning of each month m from a lognormal distribution with mean = 1.0 and coefficient of variation = 0.5, which implies that the perturbations corresponded to 50% of the magnitude of the nominal MERRA-2 precipitation. This approach has been extensively used in previous efforts to generate an ensemble of precipitation forcings (e.g., [4,[50][51][52][53]). Horizontal correlation of the precipitation errors (b j,m ) is isotropic with a 2-degree correlation length ( [5]). We assumed that errors in the precipitation forcing were perfectly correlated within 1 month. This assumption helped to assess the feasibility of using GRACE for improving precipitation estimates by narrowing its uncertainties and allowed systematic updating of the precipitation errors (b j,m ) using the monthly GRACE-like synthetic observations. With these settings we generated an ensemble of model realizations, hereafter referred to as the open-loop (or model only) simulation, which included monthly estimates of TWS, along with soil moisture at different soil levels (surface, rootzone), and all associated fluxes (e.g., precipitation, runoff, evaporation) for July 2003-July 2016.

Data Assimilation Approach
The proposed data assimilation approach was broadly similar to existing GRACE data assimilation schemes (e.g., [4,5,9,41]), except that we update errors in precipitation fluxes instead of the TWS prognostic states (i.e., soil moisture and snow). That is, we employed a two-step, ensemble smoothing approach in which the land model integration was performed twice over the course of the same month: first to collect monthly TWS observations minus predicted differences, and a second time to update that month's simulated TWS. Details of the assimilation methodology can be found in [5]. As in [5], CLSM computed, for each model grid cell and every 7.5 min, an ensemble of TWS estimates that were driven by the ensemble of precipitation forcings (Equations (1) and (2)). At the end of each month, these estimates were aggregated in space and time to compute monthly TWS observation predictions (M(b − j,m ))), which were the model's best guess of the synthetic observations. Both the spatial-and the temporal-averaging procedures were the same as those used for the generation of the synthetic observations (Z, Section 3.1.1). The ensemble-based analysis used these observation predictions to compute precipitation error increments: where the superscripts "+" and "−" in the monthly precipitation error coefficients (b j,m ) refer to the assimilation (or posterior) and open-loop (or prior assimilation) ensembles, respectively. In Equation (3), GRACE-like synthetic TWS observations (Z j ) were also perturbed according to [54] using the measurement error (Section 3.1.1). The array K represents the Kalman gain, which controlled the magnitude of the precipitation error increments (b + j,m − b − j,m ), and was calculated as: where R is the measurement error covariance (a diagonal matrix with non-zero elements equal to the measurement error), and C MM is the sample error covariance of the observation predictions. C MM is simply obtained as a covariance matrix from the observation predictions M(b − j,m ). Thus, the key technical difference with respect to previously developed GRACE data assimilation efforts was in the calculation of the cross-covariance matrix C bM , the error covariance term computed between the forecast observation predictions (M(b − j,m )) and the precipitation error coefficients (b − j,m ). This was different from previous studies which instead computed this cross-covariance directly from the ensemble of the TWS prognostic variable states (e.g., groundwater, soil moisture, snow, etc.) and the observation predictions while updating the prognostic states directly.
After computing the precipitation error increments (Equation (3)), we rewound the model to the beginning of the month and re-ran the same month with the updated precipitation errors. Under the assumption that TWS errors retained information of the precipitation errors, the rewound run should include a more realistic precipitation input, as the assimilation step minimized the differences between modeled and observed (synthetic) TWS. The data assimilation run produced a series of estimates of precipitation errors, along with fluxes (precipitation, evaporation, runoff), and prognostic variables (snow, soil moisture) over the same period as the open-loop, from July 2003 to July 2016. Thus, the adjusted precipitation was included among the outputs of the assimilation step.

Metrics for Performance Evaluation
We assessed the performance of the proposed assimilation scheme using the Pearson correlation coefficient (R) and the unbiased root mean-square error (ubRMSE, [55]), computed for both the open-loop and data assimilation experiments with respect to the set of synthetic hydrologic variables (Section 3.1). These two metrics are commonly used to evaluate modeled vs. observed hydrological estimates in terms of their dynamic variability (R) and random error (ubRMSE). Bias metrics were not computed as they were close to zero by the design of the synthetic experiment, and RMSE metrics were effectively equivalent to ubRMSE metrics here.

Verification of the Method
To verify that the new assimilation scheme behaves as designed, we first investigated the R and ubRMSE performance skill metrics for TWS at the spatial and temporal scales of the assimilated synthetic TWS observations. The domain average values for the R and When the synthetic observations of TWS were assimilated, the R skill increased everywhere over the domain region, with some areas receiving larger increases than others. The largest improvements (darker blues in Figure 3b) correspond to regions where the openloop TWS performs poorly (blue areas in Figure 3a). These improvements were expected because the assimilated TWS should, by design, agree better with the observed (synthetic) TWS (Figure 3c). For the sample location shown in Figure 3c

Assessment of Precipitation Errors
The open-loop (Section 3.2) precipitation errors were updated by the assimilation scheme (Section 3.3). The mechanism behind the precipitation error updates is illustrated for a sample location in Figure 4. The black solid lines in the figure report true-error values, determined as: b True,m = P True,m P M2,m where P M2,m is the MERRA-2 precipitation for a given month m and the synthetic "true" precipitation P True,m was obtained from the 8-year shifted MERRA-2 precipitation (Section 3.1). The above ratio inevitably created some instabilities for those times when the nominal MERRA-2 precipitation (P M2,m ) approached zero. To limit these instabilities, b True,m we calculated and showed only for the months when P M2,m was greater than its only minimally (because of sampling noise); its mean remained close to one (red solid line in Figure 4a) and did not follow the true precipitation errors.
For the assimilation case, the precipitation error distributions from the ensemblebased analysis (Section 3.3) followed the true errors much more closely. Nonetheless, the assimilation results could not quite retrieve some of the high values (extremes) reported in the synthetic precipitation error time-series as could be seen, for example, in September 2007 and in June 2008. In general, high values of b True,m corresponded to times when synthetic precipitation (i.e., the 8-year shifted MERRA-2 precipitation) was much higher than the nominal one from MERRA-2. The ensemble-based analysis performed best under a Gaussian linear assumption which may not be met for some of these instances given the non-linearity of the land surface model. Finally, note that the ensemble spread of the assimilation was much smaller than that of the open-loop. In fact, as it combined multiple sources of information, the assimilation allowed for increased robustness (i.e., reduced uncertainty) of the estimates around the mean value.
Finally, Figure 4b Figure 5 and Table 1 summarize the domain-average performance skill in terms of R and ubRMSE for all storage and flux variables of the water budget for the open-loop and assimilation cases. The assimilation results (blue bars in Figure 5) indicated consistently improved statistics relative to the open-loop results (red bars) for both the R and ubRMSE metrics because of the precipitation correction in the TWS assimilation scheme. These bulk statistics demonstrated that GRACE-like TWS assimilation could help correct (adjust) the precipitation forcing, which in turn enhanced TWS, snow mass, soil moisture, runoff, and evaporation estimates. Next, we investigated the spatial pattern dynamics related to fluxes and terrestrial water storage components.

Precipitation
The synthetic data assimilation experiment allowed us to understand if TWS observations (in this case, synthetically derived, but similar to GRACE observations) had the capability to inform precipitation fluxes, especially in regions where precipitation observations were scarce or otherwise uncertain. Figure 6 illustrated the difference in correlation skill (vs. the synthetic observations) between the data assimilation and open-loop estimates of total precipitation (i.e., rainfall plus snowfall), snowfall only, and rainfall only. Domain-mean differences were 0.08, 0.07, 0.05 for the total precipitation, rainfall, and snowfall respectively (Table 1). Domain-mean differences in ubRMSE were −0.13 mm/d, −0.11 mm/d, and −0.03 mm/d for the total precipitation, rainfall, and snowfall respectively (Table 1). Although small, these differences (increased R and decreased ubRMSE) confirmed an improvement thanks to data assimilation.
In general, the data assimilation scheme improved total precipitation correlations everywhere in the domain, except for a few localized regions where the skill was considerably degraded (red areas in Figure 6a). Specifically, there were large discrepancies between assimilated precipitation and the synthetic observations in some areas like around the Gulf of Alaska and in the Yenisei River Delta (Russia), mainly during the summer months. Because of this, the spatial pattern of rainfall correlation differences (Figure 6b) largely resembled that of total precipitation (Figure 6a). The Gulf of Alaska received abundant precipitation along the maritime coastal mountain range, but several kilometers inland it decreased dramatically. It is unlikely that this strong gradient in precipitation was resolved by the assimilation of coarse scale TWS observations. The assimilation yielded improvements nearly everywhere when we looked at snowfall correlation differences (Figure 6c). This might be related to the fact that, when snowfall occurred, nearly all of the snowfall resulted in a change in TWS due to the accumulation of SWE and the fact that runoff and evaporation were typically small when the temperature was below freezing. When it rained, a substantial fraction of the rainfall evaporated and/or is converted into runoff, thus the resulting change in TWS was generally smaller than the rainfall total. The correlation differences demonstrated that, in general, TWS synthetic observations could be used to constrain precipitation if the latter was unknown or uncertain.

Evaporation and Runoff
The differences in correlation between the open-loop and the assimilation results (∆R) for runoff and evaporation are reported in Figure 7a Runoff fluxes were improved everywhere in the study domain. By and large, the spatial patterns in these maps corresponded to those seen in the difference map for the precipitation fluxes ( Figure 6). Not surprisingly, improved precipitation fluxes corresponded to better runoff estimates, which would only have been marginally possible with the traditional GRACE data assimilation approaches. In the latter, runoff changes would have been limited to improvements in baseflow and would not capture surface runoff, which, for the most part, dominates the total runoff signal. Regions of degradation include a few sporadic areas where the runoff ∆R is negative (red area in Figure 7b), such as the coast of the Gulf of Alaska, where precipitation (rainfall) fluxes experienced degradation as well (Section 4.3.1).

TWS Storage Components
Here, we focused on analyzing the TWS components in terms of snow water equivalent (SWE) and prmc. The prmc could be used as a proxy for both groundwater table and unsaturated zone soil moisture dynamics (Section 2).
The In general, TWS skill improvement (Figure 3b) resembles that obtained in the prmc (Figure 7d). For the prmc, the open-loop domain-average R skill equaled 0.54 and it was 0.73 for the assimilation runs. For ubRMSE, these were 0.020 m 3 m −3 for the open-loop and 0.014 m 3 m −3 for the assimilation. The assimilation yielded improvements nearly everywhere in the study domain, and the patterns of these changes were very similar to those seen in TWS. This is not surprising, given that the prmc encompasseed shallow groundwater and unsaturated soil moisture.

Summary and Conclusions
The GRACE and GRACE-FO satellite missions have been measuring TWS globally for nearly two decades at spatial and temporal resolutions that are often considered too coarse for practical hydrological applications. State-of-the-art data assimilation methods have proven effective for spatial and temporal downscaling through direct updates to the TWS prognostic variables (i.e., soil moisture, groundwater, snow states, and surface water). This work presents an alternative approach in which the TWS data assimilation scheme adjusts errors in precipitation fluxes instead of the land model prognostic variables that make up TWS. The proposed method has the main advantage of adjusting water fluxes and states according to physically-based principles encoded in the land model as opposed to the statistically-based water increments obtained in the traditional GRACE assimilation schemes. In this way, this alternative assimilation method ensures the closure of the model internal water budget. Additionally, model fluxes such as runoff and evaporation are automatically and directly modified in response to the adjusted precipitation inputs and according to model physics.
To assess the potential of the proposed approach, we conducted a synthetic experiment where the ensemble of precipitation errors is sampled from a lognormal distribution. Under the assumption that TWS errors retain information of the precipitation errors, the assimilation scheme updates these errors using an ensemble-based analysis approach. Our main findings can be summarized as follows:

1.
On average, the assimilation can retrieve reasonable precipitation errors. However, when the true errors are high, the assimilation does not always recover high enough error estimates. This might be during times when the nominal MERRA-2 precipitation is very low and the synthetic precipitation error is high.

2.
The improvements in precipitation estimates are more robust for snowfall than for rainfall. Degradations sometimes occur over regions where there is a strong horizontal gradient in precipitation. Nonetheless, our results suggest that the proposed TWS assimilation technique can improve precipitation flux estimates. This is consistent with earlier studies ( [37,38]) that concluded that GRACE observations could be beneficial in constraining precipitation fluxes in regions with large uncertainties (e.g., due to a scarcity of in situ data).

3.
The assimilation procedure improves estimates of evaporation, runoff, SWE and profile soil moisture, thus demonstrating the feasibility of the proposed approach as an alternative to directly updating the TWS states in land surface models. This mechanism may work best where the precipitation fluxes are the main driver of errors in the water budget.
In conclusion, the results presented here demonstrate that the assimilation of synthetic TWS observations can correct (adjust) the model's precipitation forcing and in turn enhance model estimates of TWS, soil moisture, snow mass, runoff, and evaporation. Because of this, the proposed approach has a more positive impact on simulated runoff than the standard GRACE data assimilation techniques in which TWS is updated directly. Our findings further demonstrate the benefits of GRACE assimilation on land surface model simulations, thus adding to the myriad ways that GRACE TWS retrievals have been harnessed to improve our knowledge of the terrestrial hydrologic cycle. Nonetheless, this work only presents a synthetic experiment and future work will have to test these findings using real data. Finally, an inherent assumption of the proposed synthetic data assimilation scheme is that all model errors arise from errors in the precipitation forcing, which are the largest source of error. Future work will need to address other unaccounted sources of errors (e.g., in the model inputs and parameters) for estimating TWS and related water storages and fluxes. Funding: This study was supported by the NASA Terrestrial Hydrology and GRACE and GRACE Follow-On Science Team programs.