Comment on tc-2021-278 Anonymous Referee # 1 Referee comment on " Basal Water Storage Variations beneath Antarctic Ice Sheet Inferred from Multi-source Satellite Data

I still acknowledge the innovation of the paper in facing the challenge of estimating the BWSV. The authors apply a stepwise and iterative procedure to estimate this effect in a data-driven approach. For this purpose, temporal gravity changes of GRACE are used. These gravity changes are caused by different sources. The authors try to determine all causes, except the basal mass balance (BMB), by making use of additional information (altimetry (ICESat), firn modelling, GIA modelling, and density assumptions). Those are subtracted from the total gravity change, leaving gravity changes due to BMB and and error. In the final step, BWSV is determined by using modelled basal melting rates. Nevertheless, methods for combining satellite data sets over ice sheets have been investigated in several studies in recent years. It should be noted that only the ICESat period (2003 to 2009) is considered here, although much more recent data sets such as CryoSat-2, ICESat-2, GRACE-FO, etc. are available. BWSV are certainly also relevant over the ICESat period, but it limits the novelty as well as the impact of the study. At the end the authors stated that data sets with a higher spatial and temporal resulotuion would come to a "more elaborate result". However, these data sets have been available for several years.


Introduction
Antarctic basal water is widely generated by geothermal heating, basal pressure melting and frictional heating. The basal liquid water converge on subglacial lakes or spread over ice-bed interface, and connected to each other (Wingham et al., 2006;Fricker et al., 2016), which forms basal drainage systems with complex basal water migrations (Pattyn, 2008;Carter et al., 2015).
Temporal mass variations in Antarctic basal liquid water storage are called basal water storage variations (BWSV), it 25 influences basal effective pressure and trigger changing ice sheet velocities (Bell and Robin, 2008;Fricker et al., 2007;Alley, 1992). BWSV are mainly controlled by basal melting rate and basal water migrations. In ice-bed interface, basal ice melt in high basal melting rate regions to increase liquid water storage, the liquid water refreezes when flowing through supercooling regions to https://doi.org/10.5194/tc-2021-278 Preprint. Discussion started: 5 November 2021 c Author(s) 2021. CC BY 4.0 License.
Mass redistributions lead to temporal variations in gravity field defined in geoid, and can be captured by Gravity Recovery 60 and Climate Experiment (GRACE) observations. These variations contain surface gravity variations caused by snow accumulation, sublimation/runoff and ice flow, gravity variations of BMB in ice-bed interface, and gravity variations of solid Earth that arise from the combined effect of GIA and Earth's load deformation. Besides, BMB and basal melting/refreezing also influence ice sheet's vertical movement by controlling basal mass changes and ice-water phase change process, thereby contributing to gravity variations. Accordingly, we can express the Antarctic integrated time-variable gravity field as a 65 superposition of layered gravity variations: where ∆ are Antarctic integrated time-variable gravity variations (including load effect) observed by GRACE (section 3.1); ∆ are gravity variations caused by Antarctic surface mass redistributions process, including firn layer's processes and ice flow; ∆ are BMB-derived gravity variations; ∆ are gravity variations caused by Earth's isostatic 70 deformation and mass changes in the interior of the Earth, and can be obtained from GIA model (section 3.4); ∆ are additional gravity variation caused by ice sheet's vertical movement.

Figure 1
Mass redistributions in different layers on AIS. Italics are components that contribute to gravity variations, and red fonts refer to components resulting in surface height variations. Ice sheet's vertical movement is controlled by combined effect of BMB, basal water refreezing and ice melting. Crustal vertical deformation refers to Earth's isostatic and load deformation.
Mass redistributions also lead to surface height variations of AIS: where ℎ are AIS surface height variations that can be observed by ICESat (section 3.2); ℎ are height variations 75 caused by process of firn layer including surface snow accumulation, sublimation/runoff, melting and firn compaction. ℎ https://doi.org/10.5194/tc-2021-278 Preprint. Discussion started: 5 November 2021 c Author(s) 2021. CC BY 4.0 License. are height variations caused by ice flow; ℎ are the ice sheet's vertical movement; ℎ are crustal vertical deformations that be calculated by fitting GIA predicted basal uplift rate to sparse GPS observation rate (section 3.4).

Estimation of basal mass balance and basal water storage variations
Separation of ∆ from integrated time-variable gravity field (Equation 3) is enabled by different sensitivity of satellite 80 observation with respect to gravity and height variations in different layers. According to Equation 1, ∆ can be expressed as follows: In Equation 3, ∆ and ∆ can be calculated through gravity forward modelling method and corresponding density data, providing the known height variations caused by surface process and ice sheet's vertical movement. Then mass changes 85 of BMB can be estimated from ∆ by utilizing gravity density inversion method. However, height variations caused by surface process are complicated, especially in ice flow regions where the surface height variations are dominated by ice thinning, and the combination of surface and ice density is more appropriate in gravity forward modelling process; furthermore, BMB and basal melting/refreezing cause ice sheet's vertical movement, and this movement is also expressed as height variations and result in additional gravity variation ∆ . These coupled variations are contained in gravity and surface 90 height variations that captured by satellite gravity and altimetry mission. Therefore, the key problem of separating from integrated time-variable gravity field is to distinguish different height variation components on each layer and IVM from altimetry observation data. For this purpose, we adopted an iteration method that accounts for surface density discrimination and IVM correction, to estimate the corresponding height and gravity variation (Fig. 2), details are described as follows.
Firstly, we set the initial value of ℎ to 0, and deducted ℎ and ℎ from AIS surface height variation ℎ , to 95 obtain the initial height variation from Antarctic surface mass redistributions ℎ = ℎ − ℎ − ℎ . It is noticed that ℎ only contains height variations caused by firn layer's processes ℎ and ice flow ℎ . Afterward, we utilized surface density discrimination method (Gunter et al., 2014), in combination with a semi-empirical estimated timedependent firn densification model (FDM) that describes spatio-temporal evolution of Antarctic firn layer (Ligtenberg et al., 2011) (section 3.3), to discriminate height variation component from dh . The corresponding surface density is assigned 100 as follows: where are firn density distributions (Ligtenberg et al., 2011) percolation, refreezing and runoff) and firn compaction (Ligtenberg et al., 2011). is ice density (917 kg m -3 ) and = 105 + .
As shown in Equation 4, the negative height differences between ℎ and ℎ greater than 2 are attributed to ice flow (that is, ℎ = ℎ − ℎ ) with the assignment of ice density , and surface firn density is assigned to used to estimate initial mass variations of BMB (∆ , in form of equivalent water height (EWH)), and a 300km gaussian filter was performed on ∆ to match the spatial resolution of GRACE data. Thirdly, we used basal melting rate data 115 (section 3.5) (Pattyn, 2010) to judge whether basal is in melting or freezing condition, and combined ∆ to calculate ℎ . The calculation of ℎ was performed based on the assumption of ice layer's instant vertical movement in response to BMB: in basal melting regions, ice sheet's vertical movement ℎ is replaced by ∆ ; in basal freezing regions, ℎ was replaced by ∆ * /1000. The ice layer's instant response assumption is reasonable and can be justified by strong surface height variations (up to 90m/yr) caused by subglacial lakes water changes in a short time (Fricker et al., 2016). 120 Afterward, ℎ was used to calculate the gravity correction of ice sheet's vertical movement ∆ . Calculation of ∆ can be achieved by estimating gravity difference between ice sheet before and after it moves, while this correction method is multifarious due to the complication of surface process. For simplicity, we only updated the value of ℎ in first step, and repeated the surface density discrimination method for iteration, until the ∆ is stable (with the total BMB difference between two iterations smaller than 5Gt/yr). Finally, we combined ∆ and basal melting rate data (Pattyn, 2010) to 125 evaluate basal water storage variations (∆ , in form of equivalent water height (EWH)): regions with basal melting rate greater than 0 are assumed to maintain basal water migrations, and the (∆ are the gross of ∆ and basal ice melting; regions with the melting rate of 0 are assumed to be at supercooling condition, and no BWSV occurs in these regions.

Gravity density inversion method based on forward modelling
Gravity forward modelling is used to convert the volume variations and density distribution data into gravity variations. In this 130 study, we adopted an oblique triangular prism that addresses terrain surfaces on ellipsoids without gaps. For this purpose, surface of AIS was divided into many equilateral triangles with sides approximately 50km (Fig. 3).
where are gravity variations in AIS layers that defined in geoid, N is the number of triangular prisms within the integrated radius (500km), are density of each triangular prisms, _ and _ are lower/upper surface height of each triangular prisms 140 (that is, AIS's basal and surface height that comes from BEDMAP2 (Fretwell et al., 2013)). Then a layered gravity density inversion method was developed to estimate BMB. In this method, BMB was assumed to occur on a thin layer in Antarctic ice-bed interface, with a fixed thickness and variable density. Then BMB induced gravity variations can be expressed as follow: where is the thickness of the fixed layer on Antarctic bedrock, is fixed thin layer's density distribution, is water density, is the EWH of BMB that need to be estimated.
Then, Equation 7 could be modified as follows: For M gravity points (M=10314), Equation 9 can be expressed as:

155
Then the conjugate gradient method was employed to solve the EWH of BMB (that is, X in Equation 10).

Uncertainty estimation
Estimating uncertainties of BMB/BWSV is comprehensive due to the combination of multi-source data in gravity forward/inversion including GRACE, ICESat, FDM, GPS, GIA, and basal melting rate. Among them, Uncertainties of GRACE are provided by calibrated errors for spherical harmonic coefficients, and can be estimated utilizing the method of Wahr et al. 160 (2006); Others are provided by uncertainties of volume variations or uplift rate, and can be converted to associated gravity uncertainties by using cylinder model (Heiskanen and Moritz, 1967) combined with error propagation law. The uncertainty of BMB is obtained as follows: where ∆ , ∆ , ∆ , ∆ , ∆ are uncertainty of gravity variations related to GRACE, 165 ICESat, FDM, GPS and GIA respectively.
Afterward, uncertainties of BMB ∆ are calculated from ∆ , based on the inversion of cylinder model and error propagation law, and uncertainties of BWSV ∆ are estimated as follow: where ∆ refers to EWH uncertainties of BMB, ∆ refers to standard deviations of basal melting rate (Pattyn, 170 2010; Van Liefferinge and Pattyn, 2013). for Geosciences (GFZ). Each solution is represented by fully normalized Stokes potential coefficients with degree and order to 60. The C20 coefficients of each GRACE solution were replaced by satellite laser ranging (Cheng et al., 2013), and the degree 1 coefficients are inserted by the values of Swenson et al. (2008). Striping errors were eliminated by using P4M6 filter (Chen et al., 2007) and 300km Gaussian smoothing filter. Scaling factors was multiplied to restore the amplitude dampening and leakage-out errors (Gao et al., 2015). In this study, basal mass balance signals contained in GRACE only occurs in AIS 180 region, so leakage-in errors corrections were not performed here. Linear term was abstracted through least squares adjustment of six-parameter function consisting of constant term, linear trend, annual and semi-annual periodic signals. The gravity variation trend was calculated according to the spherical-harmonic expansion of the gravity anomaly (Heiskanen and Moritz, 1967) (Equation 13), and we used the mean gravity variations trend from the three GRACE solutions to account for timevarying gravity on AIS. It is worth noting that time-varying gravity contains not only linear trend of AIS gravity variations, 185 but also load effect. Then the AIS time-varying gravity trend is expressed as follows: where is the geocentric gravitational constant, is the mean radius of Earth, is the load Love number of degree (Wahr et al., 1998), are the normalized associated Legendre functions, ∆ and ∆ are modified harmonic coefficients, and are colatitude and longitude, and ( ) is the geocentric distance related to . 190 The uncertainties were estimated through the method of Wahr et al. (2006), from the calibrated errors for spherical harmonic coefficients provided by CSR, JPL and GFZ, and coefficient errors of degree 1 and degree 2 were replaced by the corresponding standard deviations provided by Cheng and Swenson (Cheng et al., 2013;Swenson et al., 2008).

Surface height variations derived from ICESat
Height variations on surface AIS were derived from Ice, Cloud, and Land Elevation Satellite (ICESat) mission, which were provided by National Snow and Ice Data Center (NSIDC), with the mission period spanning from February 2003 to October 2009. The object of ICESat mission is to obtain the height variations with an accuracy of ≤1.5cm/yr for spatial resolution of 200 100*100km 2 on AIS (Zwally et al., 2002). Therefore, we estimated AIS height variations by using a block crossover analysis method (Brenner et al., 2007;Gunter et al., 2009) with each block size about 100*100 km 2 . In pre-processing of ICESat data, crossover points with height variations greater than 10m/yr were deleted in order to eliminate gross errors caused by the smallscale surface roughness, undetected forward scattering and the interpolation between successive footprints, then the 3σ criterion test was carried out in each block area to further eliminate the residual gross errors (Smith et al., 2005). The linear 205 surface height variations trends were abstracted through least squares adjustment of four-parameter function consisting of constant term, linear trend, and annual periodic signals.
The inter-campaign biases (ICB) have been evaluated by several research groups through a variety of approaches, while these ICB corrections vary largely due to the selection of different reference surfaces. Zwally et al. (2015) made the correction by using concurrent Envisat radar altimetry of the same surface in open water and thin ice in leads and polynyas in Antarctic sea 210 ice pack. Richter et al. (2014) and Schroeder et al. (2016) made the correction based on a near-zero surface height changes and hydrostatic equilibrium for the snow surface above Lake Vostok and its surroundings in East Antarctica (EA) that is concluded by repeat measurement of a GPS stake network and kinematic profile resurveys. Other methods for ICB estimations are performed based on the assumption of near-zero elevation changes regions in EA. To date, none of ICB corrections has been endorsed by the ICESat Science Team, NASA, or NSIDC. Therefore, uncertainties of ICB corrections also have an evident 215 influence on the estimation of BMB results. In this study, we used the mean ICB corrections value of Zwally et al. (2015), Richter et al. (2014) and Schroeder et al. (2016), to avoid the artificial near-zero elevation changes assumption.

Firn densification model
Height variations of AIS firn layer are associated with both mass-conserving (firn compaction), mass-changing (snow accumulation, sublimation/runoff) processes, and ice flow (Zwally et al., 2005;Gunter et al., 2014). Mass-conserving process only influences surface height variation, and have no effect on gravity variation. Therefore, discriminating these processes is essential for estimating volume variations from firn layer's mass-changing and ice flow, thereby contributing to converting 225 surface volume variations into gravity variations in combination of surface density distribution of Kaspers et al. (2004). For firn densification expression that is tuned to fit depth-density observation and liquid water processes (meltwater percolation, refreezing and retention). At AIS surface, the time-dependent IMAU-FDM is forced by the surface mass balance, surface 230 temperature and wind speed from the regional atmospheric climate model RACMO2/ANT, to simulate the temporal evolution of firn density and surface height (Ligtenberg et al., 2011). For simplifying calculation process, this study used a constant surface firn density model from IMAU-FDM, which would cause a mean bias of -20.4kg m -3 compared with the temporal firn density (Keenan et al., 2021), while this bias is negligible compared with the mean density of 340kg m -3 over AIS. Figure 6a-6b shows the linear surface height changes rate derived from IAMU-FDM and the associated uncertainties over the study 235 period. For the consistency of spatial resolution, a 300km Gaussian filter was also performed in IAMU-FDM-derived surface height changes rate.

Glacial isostatic adjustment and crustal vertical deformation rate derived from GPS
In our study, three GIA models were used to account for Earth's isostatic deformation and mass changes in the interior of the Earth, including ICE-6G (Peltier et al., 2015;Argus et al., 2014), IJ05_R2 (Ivins et al., 2013) and W12a (Whitehouse et al., 240 2012). ICE-6G is a global GIA modal that is constrained to fit all available geological and geodetic observation data including GPS, ice thickness, relative sea level histories and the age of marine sedimentation. IJ05_R2 and W12a are regional GIA models that are constrained by extensive geological and glaciological data. All GIA model provides geoid rate in form of stokes coefficients, as well as predicated uplift rate in grid format. Uncertainties of ICE-6G and IJ05_R2 arise from the misfit between uplift rates from GIA and GPS: for ICE-6G, the rms value of 0.89 mm/y is estimated from 42 GPS stations (Argus et 245 al., 2014); for IJ05_R2, the uniform value of 1.40 mm/y comes from 6 GPS stations with observation period over 3000 days (Ivins et al., 2013). Uncertainty of W12a is provided by the upper and lower bounds of geoid rates.
The crustal vertical deformation rate on Antarctica contains two components: Earth's isostatic deformation and load deformation of the crust due to the mass changes of ice sheet. The combined deformation rate during 1995-2013 comes from GPS observed uplift rates from a total of 118 GPS sites listed in Sasgen et al. (2016). We used 57 of them over the study period, 250 and removed sites with the errors greater than uplift rates. The sparse GPS data were interpolated by the GIA predicated spatial uplift distribution, to obtain the crustal vertical deformation throughout the whole Antarctic ice sheet.

Basal melting rate 255
Antarctic basal melting rate was used to estimate basal water storage variations (BWSV) in combination with basal mass balance (BMB). The basal melting data used in this study is determined by Pattyn (2010), through a hybrid method that combines a priori information (such as on-site measurements data), topography, accumulation, surface temperature, geothermal heat flow data with a merged ice sheet/ice stream model. The associated uncertainty is evaluated through a series of sensitivity experiments. The multi-source dataset used in estimating basal melting rate cover different period spanning from 260 1980 to 2004 (Pattyn, 2010;Van Liefferinge and Pattyn, 2013), differing from the period of this study (2003)(2004)(2005)(2006)(2007)(2008)(2009)).
Nevertheless, this period discrepancy has little influence on our estimation of BWSV, due to the stable basal conditions caused by isolation of overlying ice sheet. water mass change rate of 65 Gt/y. In this study, we used a 300km Gaussian smoothing filter on the basal melting rate to match the spatial resolution of BMB.

Basal mass balance beneath Antarctic ice sheet 270
Based on the gravity density forward/inversion iteration method in section 2.2, we calculated the Antarctic basal mass balance (BMB) that is mainly caused by basal water migrations. Figure Table 1. The drainage basins division method is employed based on the spatial similarity of simulated basal meltwater pathways and surface ice flows. As shown in Table 1, although different GIA models are used in estimating BMB, regional BMB rates related to 3 GIA models show similar values in most of the drainage basins: 280 3 drainage basins (B5, B6 and B9) exhibit obvious basal mass increases (with the absolute value of regional BMB rates greater than associated standard deviations and lager than 10Gt/yr), and 3 drainage basins (B2, B3 and B11) displays obvious basal mass decreases. Mean BMB rates (average of BMB rates related to 3 GIA models) in East Antarctic Ice Sheet (EAIS, including B1-B8, B17 and B18) and West Antarctic Ice Sheet (WAIS, including B9-B12 and B16) are 11 ± 12 Gt/yr and -31 ± 5 Gt/yr, accounting for 23% and 30% of the corresponding ice-sheet mass balance estimated by satellite gravimetry method (Shepherd 285 et al., 2018), respectively. Mean BMB rate in Antarctic Peninsula Ice Sheet (APIS, including B13-B15) is very low (-1 ± 1 Gt/yr), indicating that there are few basal mass variations. The total mean BMB rate beneath Antarctic Ice Sheet (AIS) is -21 ± 13 Gt/yr, accounting for 28% of the satellite gravimetry derived ice-sheet mass balance rate (-76Gt/yr) from Shepherd et al. (2018).  (Figure 10d-10f) also shows consistent 295 spatial distributions: the largest STD (≥15 mm/yr) are in APIS and WAIS, which mainly comes from the STD of GRACE (about 25% of total STD estimated from Equation 11), ICESat (25%, mainly from ICB), FDM (20%) and GPS (20%); the medium STD (10-15 mm/yr) are in the marginal region of EAIS, mainly coming from GRACE (40%), ICESat (30%) and FDM (20%); the low STD (<10 mm/yr) cover a large extent the interior of EAIS, mainly coming from GRACE (45%) and ICESat (35%). 300 In Figure 10a RP, SC, GVC and ASB region possess the most obvious basal mass increases, this is attributed to the low-lying basal terrain facilitating the accumulation of basal meltwater driven by basal hydraulic potential gradient (Fig 18, (Göeller, 2014)), and Wright et al. (2012) revealed the possibility of basal meltwater accumulation in ASB region by revealing extensively distributed basal pressure melting through regional airborne geophysical survey and numerical ice sheet modelling method. In SC region, the mechanism of basal meltwater accumulation is also revealed in MacAyeal and Whillans ice streams, where the 310 similar basal hydrological systems are also found (Fricker et al., 2016;Gray et al., 2005); additionally, the accumulated basal meltwater in SC region might also partially account for the dynamic thickening of Kamb ice stream. Similarly, basal mass increases in IIS regions should be also caused by low basal hydraulic potential and low-lying basal terrain that conduces to the collection of surrounding waters. Basal mass increases in RIS region arise from subglacial lake water accretion and discharging into the bedrock trench underneath the Recovery Glaciers (Bell et al., 2007), and similar basal mass increases might also occur 315 in SG region. In DML region, basal mass increases occur in basal ridge, which might be attributed to the basal supercooling https://doi.org/10.5194/tc-2021-278 Preprint. Discussion started: 5 November 2021 c Author(s) 2021. CC BY 4.0 License. condition in DML that conduce to the basal meltwater refreezing in basal ridge when flowing upward driven by basal hydraulic gradient. The possibility of basal water flowing upward to basal ridges has been verified in Gamburtsev Subglacial Mountain (GSM) through analysing a comprehensive geophysical data set (Creyts et al., 2014), and the consistent pattern is also shown in our BMB result, although the pattern is not significant enough in GSM region. 320 Among them, the most significant basal mass decreases along Amundsen Sea coast region are caused by the active ice melting in response to basal geothermal flux (Schroeder et al., 2014), this lead to a majority of basal water discharging into Amundsen Sea through basal channels (Gray et al., 2005). Basal mass decreases in the interior of EAIS are caused by the basal meltwater generated from basal pressure-melting (Augustin et al., 2007) flows to marginal regions, driven by basal hydraulic potential 325 gradient (Göeller, 2014). The decreases in EndL region lacks justification, which might also be caused by the volume scattering in satellite altimetry data that underestimate surface height variation, or caused by the error in FDM that lead to the underestimation in surface mass loss in gravity forward modelling process, thereby impeding BMB estimation in these regions, and verifying these possible conjectures need further effort.

Basal water storage variations beneath Antarctic ice sheet 330
Antarctic basal water storage variations (BWSV) were estimated based on BMB and Antarctic basal melting rate data, through the combination method in section 2.2. It reveals mass variations of liquid water storage beneath Antarctic ice sheet. Table 2 displays regional BWSV rates related to 3 GIA models and associated standard deviations. As shown in Table 2, regional BWSV rates related to different GIA models also show similar values in each drainage basin: 2 drainage basins (B5 and B6) exhibit obvious basal water increases, and 1 drainage basin (B11) displays significant basal water decreases. Mean BWSV 335 (average of BWSV rates related to 3 GIA models) rates in EAIS, WAIS are 47 ± 12 Gt/yr and -4 ± 5 Gt/yr, and no obvious BWSV occurs in APIS. The total mean BWSV rate beneath AIS is 43 ± 13 Gt/yr, accounting for 66% of the basal melting rate (65Gt/yr) (Pattyn, 2010), which indicates most of the basal melted water is stored in ice-bed interface. Figure 11a-11f displays the spatial distributions of BWSV rates on AIS related to 3 GIA models, and associated uncertainties, in form of equivalent water height (EWH). Where red colours represent basal water increases that arise from basal water migration and basal melting; blue colours represent basal water decreases that caused by basal water runoff or basal water refreezing; blue dots represent locations of active subglacial lakes that observed through surface height variations (Smith et al., 2009), grey dots are locations of definite or fuzzy subglacial lakes that detected by radio-echo sounding (RES) measurement (Andrew and Martin, 2012). The spatial distributions of BWSV shown in Figure 11a-11c is similar to that of BMB in Figure   10a-10c, and the difference situated mainly in marginal regions of EAIS, where basal water increases are more extensive.
Results also show similar spatial distributions between obvious basal water increases (with the BWSV rates greater than associated standard deviations) and active subglacial lakes (blue circles in Figure 11a), such as the active subglacial lakes in IIS, RP, GVC, ASB and SG regions. This indicates that the basal water volumes in most active subglacial lakes are increasing, 350 despite water drainage occurs frequently; among them, the most typical justification is available in RP region, where the continued surface height increases in most subglacial lakes are found (Siegfried and Fricker, 2018). Definite or fuzzy lakes (grey circles in Figure 11b) are distributed mainly in regions of EA with low BWSV rates, exceptions are Concordia Subglacial lakes in Dome C (DC) where the basal water volumes are obviously increasing, this is attributed to the enhanced basal melting in Concordia Ridge, Concordia Subglacial Lakes and Vincennes Basin regions that conduce to the replenishment of basal 355 water storage (Carter, 2007).  In most regions of AIS, obvious basal water increases appear at regions with rapid ice flows, such as Bindschadler, MacAyeal and Whillans ice streams in SC region, Ninnis and Mertz ice streams in GVC region, Totten ice stream in ASB region and Institute ice streams in IIS region. Table 3 listed the locations, names and velocities of some fast ice streams over AIS, and corresponding basal water variations estimated in this study. Among them, most flow velocities are increasing during the study 360 period. For example, flow velocities of Bindschadler and MacAyeal ice streams were accelerating during 2001-2014, which is attributed to the instability in the subglacial till and subglacial water flow (Hulbe et al., 2016); the acceleration of Ninnis and Mertz ice streams arise from the enhanced basal slip (Rignot et al., 2011); flow velocity of Totten ice stream had increased by 18% during 2000-2007(Li et al., 2016. Regions without obvious basal water variations display no obvious flow velocity increases, such as flow velocities in Rutford ice stream (Gudmundsson and Jenkins, 2009) in WA, and Fischer, Mellor and 365 Lambert ice streams along Amery ice shelf. The spatial distribution relations between obvious basal water increases and rapid/accelerated flow velocities in most regions of AIS indicates that the accumulated basal water contribute to lubricating ice-rock interface and promote the ice flow above. However, exceptional cases are situated along Amundsen Sea coast (Pine Island and Thwaites glaciers), where significant basal water decreases follow rapid/accelerated ice flow (Han et al., 2016;Kim et al., 2015). In these exceptional cases, rapid/accelerated flows arise from the combined effect of high basal geothermal flux 370 induced basal ice-rock meltwater lubrication (Gray et al., 2005;Schroeder et al., 2014), as well as the buttressing loss induced by the extension of the shear zone and progressive disintegration (Kim et al., 2015). In summary, we believe that in most regions of AIS (except the Amundsen Sea coast), basal water storage variations appear to have an important effect on ice flow, except for other factors such as the stiffness of ice shelf, and flow velocities will further accelerate if basal water continues to increase in these regions. 375 https://doi.org/10.5194/tc-2021-278 Preprint. Discussion started: 5 November 2021 c Author(s) 2021. CC BY 4.0 License. Table 3 Locations, names and flow variations of over AIS and corresponding BWSV. Columns give (1) locations of ice stream that shown in Figure 11; (2) ice stream names; (3) flow velocities, R donate rapid ice flows, S donate slow ice flows less than 500m/yr; (4) flow variations, ↑/N/↓ donate accelerated/steaty/slowed down flow; (5) basal water storage variations, ↑/N/↓ donate obvious increased/steady/decreased basal water storage; (6) (Kim et al., 2015);

Conclusions
In this study, we presented a layered gravity density forward/inversion method for estimating basal mass balance (BMB) and basal water storage variations (BWSV) rates beneath AIS during 2003-2009, by a combination of multi-source satellite observation data including satellite gravity, altimetry and GPS data, and relevant models including firn densification model (FDM), glacial isostatic adjustment and basal melting rate. The presented method used an iteration approach for separating 385 gravity components caused by surface firn process and ice flow from the integrated gravity variations, to obtain the gravity variations and mass changes caused by BMB. When estimating BMB rates, 3 monthly GRACE gravity field solutions (CSR, GFZ and JPL), 3 inter-campaign biases corrections and 3 GIA models were used to obtain the average gravity variations, height variations and Earth's isostatic adjustment, to reduce possible errors.
Plateau, Siple Coast, Institute Ice Stream regions and marginal of EAIS. Total BMB rates beneath AIS are decreasing with the mean rate of -21 ± 13 Gt/yr (EAIS: 11 ± 12 Gt/yr, WAIS+APIS: -32 ± 5 Gt/yr), accounting for 28% of the mass balance rate (-76Gt/yr) from Shepherd et al. (2018); and this indicates that the basal water migration between basal drainage systems and 395 oceans is non-negligible in estimating AIS mass changes. BWSV rates exhibit a similar spatial distribution as BMB, and the total BWSV rates are increasing with the rate of 43 ± 13 Gt/yr (EAIS: 47 ± 12 Gt/yr, WAIS+APIS: -4 ± 5 Gt/yr), which is attributed to the active basal melting rate that producing sufficient meltwater to replenish the basal meltwater loss, and the excessed meltwater is stored beneath AIS that bring about increase in basal meltwater storage. BWSV results also show similar spatial distribution between significant basal water increases and active subglacial lakes, indicating that the basal water in most 400 active subglacial lakes is increasing, despite water drainage occurring frequently. Basal water in most definite or fuzzy lakes is relatively stable, exceptions are Concordia Subglacial lakes where the basal water is obviously increasing. We also found the spatial distribution relations between obvious basal water increases and rapid/accelerated flow velocities in most of the AIS (except the Amundsen Sea coast region), which indicates that basal water storage variations appear to have an important effect on ice flow, and flow velocities will further accelerate if basal water continues to increase in these regions. 405 The main errors in estimating BMB come from the uncertainty of GRACE, FDM and inter-campaign biases of altimetry data, especially the current results of different inter-campaign biases corrections remain large discrepancy. A more elaborate result is desirable if utilizing more reliable and higher spatial/temporal resolution data, which deserves further studies. Any other specific data and code of this study is available on request to the authors.