This study quantified the contribution of each pixel to the global TWS IAV of GRACE observations and two selected predominantly data-driven models, SINDBAD and H2M, as well as its modeling errors. The results show that the global TWS IAV is mainly driven by humid tropical and semi-arid region. The hotspots of modeling errors of the global TWS IAV are mainly located in tropical regions that span across climatic regions. The study provides an improved understanding of the global TWS IAV and its modeling error. Generally, the topic is important, and the study is well written and easy to follow. My comments are as follows.

1. In the high latitudes of the northern hemisphere, glacier changes contribute to TWS, whether the SINDBAD model and the H2M model have a glacier module.

2. It needs to be further pointed out that the model is inconsistent with GRACE in typical irrigation areas, such as the western United States, northern India, etc.

3. Figure 2(a) shows that the two models are in good agreement, and they both have some differences from GRACE. Does the input of precipitation significantly affect the simulation results of the model? If other precipitation products are used as input, will the results be different?

4. The abscissa and ordinate of the scatter plot in Figure 3 have no text description

5. How much different precipitation inputs affect the modeling error of global terrestrial water storage interannual variability? Does the precipitation input or the different model structure affect the simulation error more?

The authors presented a study to diagnose the modeling errors by comparing GRACE and model TWSA based on IAV. The motivation of this study is nice, since Scanlon’s PNAS study revealed an interesting question on the discrepancy between GRACE and models. The focus on interannual is a good complementary to the focus on trend by Scanlon et al.. Generally, this study is interesting. However, I have some critical questions related to the methods used for analysis, which may largely affect the reliability of the findings.

• Generally, WGHM, PCR-GLOBWB, and maybe some other LSMs that include GW module, are more popularly used than the two models used in this study. I am not going to say the two models used here is not good enough, but I guess many researchers would be more interested on what will it like if we use WGHM, or PCR-GLOBWB, or CLSM. Besides, it is not clear how the including of GRACE in model parameter estimation and evaluation (Line 77) will impact the comparison between GRACE and the two models.
• It is not clear why using Equation (1) to derive the IAV for analysis. I cannot understand the physical meaning of subtracting long-term trend (fit ()) from monthly values. So, the question comes that what is interannual variability, and how to define it? Can we just subtracting long-term average from monthly values? I am not sure my understanding is correct or not. Please verify it.
• Since GRACE Level-3 data has been already processed by subtracting the mean of a period (2004-2009?) from monthly TWS to get TWSA. If the authors again do subtracting (2002-2017) for GRACE and models, it may lead to mismatch between GRACE and model, because different subtracting were done for GRACE (subtracting 2004-2009, and then subtracting 2002-2017) and models (subtracting 2002-2017).
• Line 128: I am not sure it is the best way to evaluate model performance by comparing the IAV derived from GRACE and models. How about compare TWSA?
• Before Figure 2, people may be interested on seeing spatial distribution map of TWSA from GRACE and models, as well as the distribution map of IAV, which both can help we better understand the difference and consistence between GRACE and models.
• Figure 3: Sorry, but I do feel difficult to understand what the exact meanings of the spatial maps are. Maybe more information can be added to the figure showing who minus who, something like that. Besides, I guess the white blank areas here are the grid cells with positive covariances, is it true?

This study models the error between the global TWS IAV observations of GRACE and two models, SINDBAD and H2M. The authors found that the global TWS IAV is mainly driven by humid tropical and semi-arid regions, and identified the hotspots of modeling errors of the global TWS IAV mainly in tropical regions that span across climatic regions. The study presents a novel way to attribute global variability to each pixel and focused on regions where hydrological cycle components in models may not be sufficiently well represented due to their complex hydrological and climatological processes.

The study in general is well-written and easy to follow. Additional to comments made by the two Anonymous Referees, which I consider important to answer, my comments are as follows:

• As the study identifies humid regions of northern South Americas as one of the main drivers of global TWS IAV,  I suggest including these references in the discussion in which global models are compared with GRACE products in a very important instrumented tropical basin.

Bolaños Chavarría, S., Werner, M., Salazar, J. F., & Betancur, T. (2022). Benchmarking global hydrological and land surface models against GRACE in a medium-sized tropical basin. Hydrology and Earth System Sciences, 26(16), 4323-4344.

Bolaños, S., Salazar, J. F., Betancur, T., & Werner, M. (2021). GRACE reveals depletion of water storage in northwestern South America between ENSO extremes. Journal of Hydrology, 596, 125687.

• I am a bit confused with Equation 1, in figure 1 I think it is clear that TWS IAV is the result of detrending and deseasonalizing TWS, but in Equation 1, I understand that only TWS is deseasonalized.
• I think is necessary to define what is the meaning of SREX Regions, I don't identify what is.
• Why the preference for the JPL mascon if there is another mascon product like the mascon CSR that has the same resolution?

Save, H., S. Bettadpur, and B.D. Tapley (2016), High resolution CSR GRACE RL05 mascons, J. Geophys. Res. Solid Earth, 121, doi:10.1002/2016JB013007.

• Figure 2 a) describes a "NSE is the Nash-Sutcliffe Efficiency", but it does not appear in the figure