Using GRACE Data to Study the Impact of Snow and Rainfall on Terrestrial Water Storage in Northeast China

: Water resources are important for agricultural, industrial, and urban development. In this paper, we analyzed the inﬂuence of rainfall and snowfall on variations in terrestrial water storage (TWS) in Northeast China from Gravity Recovery and Climate Experiment (GRACE) gravity satellite data, GlobSnow snow water equivalent product, and ERA5-land monthly total precipitation, snowfall, and snow depth data. This study revealed the main composition and variation characteristics of TWS in Northeast China. We found that GRACE provided an e ﬀ ective method for monitoring large areas of stable seasonal snow cover and variations in TWS in Northeast China at both seasonal and interannual scales. On the seasonal scale, although summer rainfall was 10 times greater than winter snowfall, the terrestrial water storage in Northeast China peaked in winter, and summer rainfall brought about only a sub-peak, 1 month later than the maximum rainfall. On the interannual scale, TWS in Northeast China was controlled by rainfall. The correlation analysis results revealed that the annual ﬂuctuations of TWS and rainfall in Northeast China appear to be inﬂuenced by ENSO (EI Niño–Southern Oscillation) events with a lag of 2–3 years. In addition, this study proposed a reconstruction model for the interannual variation in TWS in Northeast China from 2003 to 2016 on the basis of the contemporary terrestrial water storage and rainfall data.


Introduction
Terrestrial water storage (TWS) is an important component of the global hydrological cycle system and plays a key role in the Earth's climate system [1][2][3]. Water resources are the basis of agricultural, industrial, and urban development. Therefore, the study of water reserves has a high social value and great scientific research significance. In terms of the vertical spatial distribution of water storage, TWS is the sum of surface water, snow, ice, soil water storage, and groundwater [4]. In addition, changes in TWS are an integrated reflection of activities such as precipitation, evaporation, runoff, and groundwater [5]. Traditional means of observing TWS mainly rely on ground-based hydrological stations (such as observation wells) and terrestrial hydrological models (such as Global Land Data Assimilation System (GLDAS) model [6]). Hydrological station observations are susceptible to the limitations of uneven spatial distribution of stations. TWS measurements taken at a single location can hardly be representative of the entire interest area [7]. Therefore, observation wells can provide high-resolution estimates of groundwater storage (GWS) in localized areas, but they struggle to Remote Sens. 2020, 12, 4166 3 of 21 To date, the studies of snow cover in China have focused on the following areas: the temporal and spatial distribution characteristics of snow [44][45][46]; the effects of topography, elevation, and location on snow cover variation [47]; and the combination of GRACE and hydrological model data to study groundwater storage changes and regional drought phenomena [35,48]. Northeast China (Figure 1), with winter temperatures below −11 • C, is a region with extensive and stable snow cover. The area of Northeast China is approximately 1.45 million km 2 , which is well within the nominal resolution of the GRACE data (about 0.1 million km 2 ). In winter, snowfall is frequent and evenly distributed, and snow cover lasts from November until April of the following year [44]. There have been few systematic studies of rainfall and snowfall on TWS changes in Northeast China thus far.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 20 distributed, and snow cover lasts from November until April of the following year [44]. There have been few systematic studies of rainfall and snowfall on TWS changes in Northeast China thus far.
In this paper, we combined GRACE satellite data, snow cover, snowfall, rainfall, and snow depth data to investigate the ability of GRACE data to monitor large-scale snow cover and the impact of rainfall and snowfall on terrestrial water storage in Northeast China. The main research objectives are as follows: (1) to understand the impact of rainfall and snowfall on seasonal changes in terrestrial water storage anomalies (TWSA), (2) to understand the impact of ENSO on rainfall and TWSA, (3) to simulate TWS changes from 2003 to 2016 in Northeast China on the basis of the new insights, and (4) to understand the challenges in investigating snow cover mass change by using GRACE data and uncertainties in the impact of ENSO on rainfall and TWS.

GRACE Data
We used the monthly Release 06 solutions, with spherical harmonic coefficients truncated to degree 60, provided by JPL (the Jet Propulsion Laboratory), CSR (the Center for Space Research at the University of Texas at Austin), and GFZ (the German Geosciences Research Center) (http://icgem.gfz-potsdam.de/series). The GRACE monthly Release 06 solutions provided by JPL, CSR, and GFZ were averaged January 2003 -December 2016). The spherical harmonic coefficients were treated as follows: the C20s coefficients were replaced by the satellite laser ranging solutions [49]; the degree-1 coefficients were added back according to GRACE Technical Note 13 [50], which were computed on the basis of the data of Sun et al. [51]; the GIA (glacial isostatic adjustment) effect was corrected by a three-dimensional model on the basis of the data of Geruo et al. [52]. The spherical harmonic coefficients were filtered using the DDK4 method [53]. Then, Equation (1) [12] was used to calculate the gridded TWSA. The time series of TWSA in Northeast China was calculated by using an area-based weighting method.  In this paper, we combined GRACE satellite data, snow cover, snowfall, rainfall, and snow depth data to investigate the ability of GRACE data to monitor large-scale snow cover and the impact of rainfall and snowfall on terrestrial water storage in Northeast China. The main research objectives are as follows: (1) to understand the impact of rainfall and snowfall on seasonal changes in terrestrial water storage anomalies (TWSA), (2) to understand the impact of ENSO on rainfall and TWSA, (3) to simulate TWS changes from 2003 to 2016 in Northeast China on the basis of the new insights, and (4) to understand the challenges in investigating snow cover mass change by using GRACE data and uncertainties in the impact of ENSO on rainfall and TWS.

GRACE Data
We used the monthly Release 06 solutions, with spherical harmonic coefficients truncated to degree 60, provided by JPL (the Jet Propulsion Laboratory), CSR (the Center for Space Research at the University of Texas at Austin), and GFZ (the German Geosciences Research Center) (http: //icgem.gfz-potsdam.de/series). The GRACE monthly Release 06 solutions provided by JPL, CSR, and GFZ were averaged January 2003 -December 2016). The spherical harmonic coefficients were treated as follows: the C20s coefficients were replaced by the satellite laser ranging solutions [49]; the degree-1 coefficients were added back according to GRACE Technical Note 13 [50], which were computed on the basis of the data of Sun et al. [51]; the GIA (glacial isostatic adjustment) effect was corrected by a three-dimensional model on the basis of the data of Geruo et al. [52]. The spherical harmonic coefficients were filtered using the DDK4 method [53]. Then, Equation (1) [12] was used to calculate the gridded TWSA. The time series of TWSA in Northeast China was calculated by using an area-based weighting method.
where ∆h is the variation of equivalent water height (EWH), ρ e and a are the mean density and radius of the Earth, θ and λ are geocentric colatitude and geocentric longitude, l represents the degree of the spherical harmonic coefficients, m represents the order of the spherical harmonic coefficients, k l is the Load Love number at degree l [54], P lm is the fully normalized Legendre function of degree l and order m, and ∆C lm and ∆S lm are the residual fully normalized spherical harmonic coefficients relative to the average gravity field between 2005 and 2010.

GlobSnow Snow Water Equivalent Products
The GlobSnow project, funded by the ESA (European Space Agency), has well-known precision characteristics and aims to create temporally and spatially extensive snow products [55]. The framework of GlobSnow project produces 2 snow parameter products, snow water equivalent (SWE) and snow extent (SE) [55]. The SWE information for terrestrial covers the area between 35 • N and 85 • N, excluding Greenland and glaciers [55]. The GlobSnow SWE product is projected onto the Equal Area Scalable Earth Grid (EASE-Grid) with a nominal resolution of 25 × 25 km for a single pixel [55]. It has a higher accuracy than the AMSR-2 and FY-3B products [56]. In this paper, we used monthly GlobSnow v3.0 SWE data ( January 2003 -December 2016) to validate the snow depth data (http://www.globsnow.info/swe/archive_v3.0/L3B_monthly_SWE/). The SWE data were treated as follows: (1) the projected SWE product was described in a plane coordinate system (x, y, f : x and y are the locations of SWE in plane coordinate system, f is the value of SWE), while the TWSA calculated by GRACE was described in a geodesic coordinate system, and thus the projection transformation method was needed to convert SWE product to be described in geographic coordinates (L, B, f : L and B are the locations of SWE in the geographic coordinate system, f is the value of SWE); (2) the cubic polynomial interpolation method was used to impute the SWE product to a 1 • × 1 • grid; (3) to unify the spatial resolution of SWE and TWSA, we calculated SWE data in spherical harmonics using the same cut-off degree (60) and DDK4 filter method as for GRACE, and then derived it into product expressed as equivalent water height. The time series of SWE in Northeast China was calculated by using an area-based weighting method.

Average Monthly Data of ERA5-Land
ERA5-Land is a reanalysis dataset produced by the ECMWF ERA5 climate reanalysis model. The dataset includes monthly dynamic data since 1981, which contains 50 indicators characterizing temperature, lakes, snowpack, soil moisture, radiation and heat, evaporation and runoff, wind speed, barometric pressure, precipitation, and vegetation. Compared to ERA5, it provides a consistent view of the evolution of land variables over several decades at a much higher resolution (approximately 9 × 9 km, 0.1 • × 0.1 • ) [57]. It mainly consists of 2 datasets, ERA5-Land hourly data and ERA5-Land monthly averaged data. We used monthly average reanalysis data of snowfall, snow depth, and total precipitation (TP; includes snow and rainfall) downloaded from the Climate Data Store (CDS) (https://cds.climate.copernicus.eu/cdsapp#!/home) for the period 1981-2016 [58]. All data provided by CDS are expressed as EWH. As mentioned in Section 2.2, we should unify the spatial resolution with the rest of the paper. Snowfall, snow depth, and total precipitation data were treated as follows: (1) the monthly average reanalysis snowfall, snow depth, and total precipitation data provided by CDS are the daily average for the month, and thus they should be multiplied by the number of days in the month; (2) we used the average accumulation method to impute snowfall, snow depth, and total precipitation data to a 1 • × 1 • grid; (3) we calculated these data in spherical harmonics using the same cut-off degree (60) and DDK4 filter method as for GRACE, and then derived them into products expressed as equivalent water height. The time series of snowfall, snow depth, and rainfall (TP minus snowfall) in Northeast China was calculated by using an area-based weighting method.

SOI
The Southern Oscillation Index (SOI) shows the development and intensity of El Niño and La Niña events in the Pacific Ocean. It can be calculated using the pressure differences between Tahiti and Darwin. Sustained negative SOI values below −7 typically indicate an El Niño phenomenon, while sustained positive SOI values greater than +7 indicate a typical La Niña phenomenon [59]. In this paper, the SOI data were provided by the Bureau of Meteorology of the Australian Government (ftp://ftp.bom.gov.au/anon/home/ncc/www/sco/soi/soiplaintext.html). We used it to investigate the effect of the ENSO events on interannual variations in Northeast China.

Fitting Method
Prior research has generally assumed that TWS and precipitation are the sums of linear variations, interannual variations, and semi-annual variations [32]. To obtain the phase and amplitude of TWSA, snowfall, rainfall, snow depth, and SWE, we used the least squares method to fit these time series. The annual and semi-annual cycles were also considered in the following equation [60]: where a is the constant, b is the trend, c and d are the annual and semi-annual amplitudes, x is the time, T 1 and T 2 are the annual and semi-annual cycles, ϕ 1 and ϕ 2 are the initial phases, and ε is the residual error.

Methods for Separating Interannual and Seasonal Signals
TWSA, SWE, snowfall, snow depth, and total precipitation (TP; includes snowfall and rainfall) data were all expressed in EWH. Figure 2a shows the time series of TWSA estimated using CSR, JPL, and GFZ solutions and their average. Figure 3 shows the time series of SWE, snowfall, snow depth, and rainfall. As shown in Figure 2a, we can tentatively conclude that TWSA series contained significant interannual variations with a greater amplitude than seasonal variations. In order to investigate interannual and seasonal components separately, we used a sliding average method to separate interannual and seasonal signals [4]. The original TWSA, SWE, snowfall, snow depth, and TP time series (S 0 ) were averaged over a sliding 12-month window to obtain the interannual variation time series (S inter ). For terminal observations, we used all available observations for the average. The seasonal variation (S season ) time series were obtained through S 0 − S inter , as shown in Equation (3). Figure 2b shows that TWSA contained evident interannual and seasonal signals and confirms that signal separation is necessary.
Remote Sens. 2020, 12, 4166 6 of 21 time series (Sinter). For terminal observations, we used all available observations for the average. The seasonal variation (Sseason) time series were obtained through S0 − Sinter, as shown in Equation (3). Figure 2b shows that TWSA contained evident interannual and seasonal signals and confirms that signal separation is necessary.     Figure 4a shows the time series of seasonal TWSA, rainfall, snowfall, and snow depth in Northeast China. Combined with the spectral analysis results (Figure 4b-f), it can be seen that the time series of seasonal changes of snowfall, snow depth, and rainfall had strong annual and semi-  Figure 4a shows the time series of seasonal TWSA, rainfall, snowfall, and snow depth in Northeast China. Combined with the spectral analysis results (Figure 4b-f), it can be seen that the time series of seasonal changes of snowfall, snow depth, and rainfall had strong annual and semi-annual cycles. To analyze seasonal signals in detail, we used Equation (2) to fit the TWSA, SWE, snowfall, snow depth, and rainfall seasonal variation time series. On the basis of the adjustment of the indirect observations model [61], we estimated the coefficients of Equation (2) and then estimated the annual amplitudes and annual phases of seasonal time series of TWSA, rainfall, snowfall, snow depth, and SWE with their root mean square errors (RMSEs). Table 1 shows that the largest annual amplitude of rainfall was about 56.64 mm; the annual amplitude of TWSA was similar to snowfall, about 6-7 mm; the annual amplitude of snow cover was consistent with snow depth, about 10-12 mm; and there was good consistency among annual phase of snowfall, snow cover, snow depth, and TWSA (annual amplitude peaked around January, with RMSE of annual phases below 0.5 months). However, the seasonal amplitude of rainfall was about eight times greater than that of TWSA. We speculate that this was because the increase in rainfall is largely offset by the corresponding increase in evapotranspiration and runoff on seasonal time scales. Northeast China is characterized by a temperate monsoon climate. Rainfall is concentrated in July and August, and snowfall dominates the winter when evapotranspiration and runoff are weakened; hence, TWS is susceptible to snowfall. The difference in sensitivity to rainfall and snowfall explains the different agreement in the results of phase (~1 year) and amplitude (6-12 mm) of snow depth, snowfall, and TWSA. The RMSE values for TWSA, rainfall, snowfall, snow depth and SWE annual amplitudes were 2.16, 2.17, 0.64, 0.60, and 0.65 mm, respectively. Comparing these RMSE values with their annual amplitudes, we found the relative RMSE (RMSE/amplitude) of rainfall, snowfall, snow depth, and SWE to be about 4%, 11%, 6%, and 5%, respectively, while the relative RMSE of TWSA reached 33%. It is clear that the seasonal term of TWSA is less reliable (Figure 4c).

Impact of Rainfall and Snowfall on Seasonal Changes in TWSA in Northeast China
In order to further analyze the influence of rainfall and snowfall on TWSA at the seasonal scale, we averaged the seasonal variation time series of 2003 to 2016 ( Figure 5) by month to obtain the monthly mean variation. The monthly mean variations in rainfall, SWE, snowfall, and snow depth were generally consistent with the climate characteristics of Northeast China. Snowfall was the main contributor to TWS changes in winter in Northeast China, and it dominated the variation of snow cover and snow depth. The summer rainfall only brought about a sub-peak, the peak month of which was 1 month later than the maximum rainfall. As shown in Figure 5a, SWE and snow depth were controlled by snowfall after September and reached their maximum values in February of the following year. Then, SWE and snow depth gradually decreased with the initiation of snow melt due to the increasing temperature. The TWSA series peaked in February and August. The average annual anomaly series of TWSA was generally consistent with snowfall from October to February. Due to snow melt and increase in runoff, TWSA gradually decreased from March to June and reached its minimum in June. Then it increased again until August to reach another local maximum. As shown in Figure 5b, rainfall gradually increased after May and reached its maximum value in July; thus, rainfall can be regarded as the main reason for TWS increase in summer. However, there was a month delay between rainfall and TWS, which was due to the delay in the response of transpiration and runoff to rainfall.   In order to further analyze the influence of rainfall and snowfall on TWSA at the seasonal scale, we averaged the seasonal variation time series of 2003 to 2016 ( Figure 5) by month to obtain the monthly mean variation. The monthly mean variations in rainfall, SWE, snowfall, and snow depth were generally consistent with the climate characteristics of Northeast China. Snowfall was the main contributor to TWS changes in winter in Northeast China, and it dominated the variation of snow cover and snow depth. The summer rainfall only brought about a sub-peak, the peak month of which was 1 month later than the maximum rainfall. As shown in Figure 5a, SWE and snow depth were controlled by snowfall after September and reached their maximum values in February of the  snow melt and increase in runoff, TWSA gradually decreased from March to June and reached its minimum in June. Then it increased again until August to reach another local maximum. As shown in Figure 5b, rainfall gradually increased after May and reached its maximum value in July; thus, rainfall can be regarded as the main reason for TWS increase in summer. However, there was a month delay between rainfall and TWS, which was due to the delay in the response of transpiration and runoff to rainfall.  Figure 6 shows the spatial distribution of monthly variations of TWSA, SWE, snowfall, and snow depth in Northeast China. SWE, snow depth, and TWSA results were anomalies relative to their November values. The spatial distribution of snowfall and snow depth had a good consistency with their monthly mean variation series, as shown in Figure 5a, which confirmed the results discussed above. Snowfall in the Northeast China was greater in December and March than in January and February, with a gradual increase from January to March. Under the conditions of reduced evaporation and runoff activities in winter, snowfall was preserved primarily in the form of snow cover, resulting in increased snow depth from December to February. The spatial distributions of snowfall and snow depth in March were inconsistent, which was likely due to the impact of temperature on snow melt. However, SWE and snowfall as well as snow depth showed different spatial patterns. The GlobSnow SWE record is assimilated from space-borne passive radiometer data (SMMR, SSM/I, and SSMIS) and ground-based synoptic weather stations data, and thus it probably has large uncertainties affected by vegetation cover, topography impacts, and temperature [40]; therefore, the November SWE result may have been underestimated (Figure 5a). Consequently, when subtracting the November values from the SWE results, it will lead to higher values compared to  Figure 6 shows the spatial distribution of monthly variations of TWSA, SWE, snowfall, and snow depth in Northeast China. SWE, snow depth, and TWSA results were anomalies relative to their November values. The spatial distribution of snowfall and snow depth had a good consistency with their monthly mean variation series, as shown in Figure 5a, which confirmed the results discussed above. Snowfall in the Northeast China was greater in December and March than in January and February, with a gradual increase from January to March. Under the conditions of reduced evaporation and runoff activities in winter, snowfall was preserved primarily in the form of snow cover, resulting in increased snow depth from December to February. The spatial distributions of snowfall and snow depth in March were inconsistent, which was likely due to the impact of temperature on snow melt. However, SWE and snowfall as well as snow depth showed different spatial patterns. The GlobSnow SWE record is assimilated from space-borne passive radiometer data (SMMR, SSM/I, and SSMIS) and ground-based synoptic weather stations data, and thus it probably has large uncertainties affected by vegetation cover, topography impacts, and temperature [40]; therefore, the November SWE result may have been underestimated (Figure 5a). Consequently, when subtracting the November values from the SWE results, it will lead to higher values compared to snow depth, which was the main reason for the wide difference in the spatial patterns of snowfall, snow depth, and SWE.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20 snow depth, which was the main reason for the wide difference in the spatial patterns of snowfall, snow depth, and SWE. On the basis of the spatial distribution of the average monthly change of snowfall ( Figure 6) and the average monthly series ( Figure 5), we can see that snow cover of Northeast China generally reached its maximum in February. Winter precipitation of Northeast China occurred mainly in the form of snowfall, which dominated the TWS in Northeast China in winter. To quantify the effect of snowfall on TWS and verify the reliability of the GRACE data for detecting snow cover, we used the linear regression method to analyze the series of snow depth and water storage in February, when the snow storage reached the maximum. As shown in Figure 7, snow depth and TWSA in February were closely correlated. The correlation between snow depth and GRACE water storage reached 0.81, and the correlation between SWE and GRACE water storage reached 0.72, which can be regarded as a strong correlation [62]. The standard deviations of their differences were both less than 4.5mm, with a 99% confidence level. On the basis of the spatial distribution of the average monthly change of snowfall ( Figure 6) and the average monthly series ( Figure 5), we can see that snow cover of Northeast China generally reached its maximum in February. Winter precipitation of Northeast China occurred mainly in the form of snowfall, which dominated the TWS in Northeast China in winter. To quantify the effect of snowfall on TWS and verify the reliability of the GRACE data for detecting snow cover, we used the linear regression method to analyze the series of snow depth and water storage in February, when the snow storage reached the maximum. As shown in Figure 7, snow depth and TWSA in February were closely correlated. The correlation between snow depth and GRACE water storage reached 0.81, and the correlation between SWE and GRACE water storage reached 0.72, which can be regarded as a strong correlation [62]. The standard deviations of their differences were both less than 4.5mm, with a 99% confidence level.
In summary, rainfall and snowfall controlled the seasonal TWS changes of Northeast China in different seasons. There were two peaks in the seasonal variation of TWSA, being found in February and August. The winter TWSA peak in February was caused by snowfall, and the summer TWSA peak in August was controlled by rainfall. Although the seasonal amplitude of snowfall was only 1/10 that of rainfall, it had a much greater impact on the seasonal-scale TWSA, and summer rainfall brought about only a sub-peak with a 1-month delay. Furthermore, the annual amplitude and annual phase of snowfall and snow cover were generally consistent with winter TWSA calculated by GRACE. The time series of snow cover in February and TWSA calculated by GRACE were highly consistent, verifying that GRACE is a feasible means to monitor large-scale seasonal stable snow cover in Northeast China. In summary, rainfall and snowfall controlled the seasonal TWS changes of Northeast China in different seasons. There were two peaks in the seasonal variation of TWSA, being found in February and August. The winter TWSA peak in February was caused by snowfall, and the summer TWSA peak in August was controlled by rainfall. Although the seasonal amplitude of snowfall was only 1/10 that of rainfall, it had a much greater impact on the seasonal-scale TWSA, and summer rainfall brought about only a sub-peak with a 1-month delay. Furthermore, the annual amplitude and annual phase of snowfall and snow cover were generally consistent with winter TWSA calculated by GRACE. The time series of snow cover in February and TWSA calculated by GRACE were highly consistent, verifying that GRACE is a feasible means to monitor large-scale seasonal stable snow cover in Northeast China.

Impact of ENSO on Rainfall and TWSA in Northeast China
ENSO is an important factor for short-term climate predictions [63]. In China, ENSO events have a clear modulation effect on climate anomalies [32]. Here, we studied the response of TWS and rainfall to the ENSO phenomenon on an interannual scale in Northeast China. The interannual variation time series of TWSA, rainfall, and SOI are shown in Figure 8. When the time series of SOI was shifted back 2 years, we could tentatively conclude that the shifted SOI time series had an improved consistency with rainfall and TWSA series, as shown in Figure 8a. This preliminarily

Impact of ENSO on Rainfall and TWSA in Northeast China
ENSO is an important factor for short-term climate predictions [63]. In China, ENSO events have a clear modulation effect on climate anomalies [32]. Here, we studied the response of TWS and rainfall to the ENSO phenomenon on an interannual scale in Northeast China. The interannual variation time series of TWSA, rainfall, and SOI are shown in Figure 8. When the time series of SOI was shifted back 2 years, we could tentatively conclude that the shifted SOI time series had an improved consistency with rainfall and TWSA series, as shown in Figure 8a. This preliminarily indicated that rainfall and water storage in Northeast China are influenced by the ENSO phenomenon, but that there is a lag of about 2 years. Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 20  Figure 9, the maximum correlation coefficient between TWSA and SOI in Northeast China was found to be 0.66, and the maximum correlation coefficient between rainfall and SOI was 0.64. TWSA and rainfall had a strong correlation with SOI, with lags of ≈29 months and ≈26 months behind SOI, respectively. In addition, the maximum value of the time series correlation coefficient between rainfall and interannual variation of TWSA was 0.85, with a lag of about 7 months. It has been shown that ENSO events modulate climate mainly by means of precipitation anomalies [64,65], which then modulate TWS [33]. The strong correlation among TWSA, rainfall, and SOI series indicates that the ENSO events influenced the interannual variation of rainfall and affected the interannual variation of TWS through rainfall. Lu et al. analyzed the relationship between drought  and ENSO events, and pointed out that the ENSO events may affect regional climate by influencing atmospheric circulation. They also have a significant impact on China's precipitation and a lag of about 2 years in Northeast China [66], which is consistent with the findings presented in this paper. Thus, the impact of the ENSO events on precipitation and TWS in Northeast China is a longterm effect with a lag of about 2-3 years.  Figure 9, the maximum correlation coefficient between TWSA and SOI in Northeast China was found to be 0.66, and the maximum correlation coefficient between rainfall and SOI was 0.64. TWSA and rainfall had a strong correlation with SOI, with lags of ≈29 months and ≈26 months behind SOI, respectively. In addition, the maximum value of the time series correlation coefficient between rainfall and interannual variation of TWSA was 0.85, with a lag of about 7 months. It has been shown that ENSO events modulate climate mainly by means of precipitation anomalies [64,65], which then modulate TWS [33]. The strong correlation among TWSA, rainfall, and SOI series indicates that the ENSO events influenced the interannual variation of rainfall and affected the interannual variation of TWS through rainfall. Lu et al. analyzed the relationship between drought  and ENSO events, and pointed out that the ENSO events may affect regional climate by influencing atmospheric circulation. They also have a significant impact on China's precipitation and a lag of about 2 years in Northeast China [66], which is consistent with the findings presented in this paper. Thus, the impact of the ENSO events on precipitation and TWS in Northeast China is a long-term effect with a lag of about 2-3 years.

Simulation of the Interannual Contribution of Rainfall to TWS in Northeast China
As mentioned in Section 3.2, the interannual variation of TWS is mainly influenced by rainfall in Northeast China. As shown in Figure 8a, the maximum interannual variations of TWSA and rainfall could reach 87.61 mm and 21.90 mm, respectively. However, the interannual variation of rainfall was found to only account for one-quarter of interannual variation of TWSA. This is mainly because there is a cumulative effect of rainfall on terrestrial water storage (i.e., rainfall in previous months still affects the variation of current water storage). Then, we constructed a model to simulate the contribution of rainfall to the interannual variation of TWSA (Equation (4)) [67]. In the model, n is the maximum number of months of cumulative rainfall, and a is the efficiency of the impact of rainfall. The model simulates the interannual variation time series of rainfall by continuously adjusting the values of n and a. The RMSE, mean error (ME) [61], and correlation coefficient between the simulation results and the water storage capacity were used as indicators of the evaluation accuracy. After continuous adjustment, the values of the parameters n and a were determined, and we obtained the simulation results of the effect of rain on the interannual variation of TWSA.
where a is the scale factor (0 < a < 1), t is the time (month) of simulated rainfall, Pi is the original rainfall value at the moment i, P0 is the mean rainfall value over the whole period, and f(t) is the TWSA at moment t. n is taken as 6, 7, 8, ..., 36, and a is taken as 0.1, 0.2, 0.3, ..., 1.0. By continuously adjusting the values of n and a, the sequence of RMSE, ME, and correlation coefficients could finally be obtained, as shown in Figure 10a-c. It is shown that the correlation coefficient reached its maximum value (0.82)

Simulation of the Interannual Contribution of Rainfall to TWS in Northeast China
As mentioned in Section 3.2, the interannual variation of TWS is mainly influenced by rainfall in Northeast China. As shown in Figure 8a, the maximum interannual variations of TWSA and rainfall could reach 87.61 mm and 21.90 mm, respectively. However, the interannual variation of rainfall was found to only account for one-quarter of interannual variation of TWSA. This is mainly because there is a cumulative effect of rainfall on terrestrial water storage (i.e., rainfall in previous months still affects the variation of current water storage). Then, we constructed a model to simulate the contribution of rainfall to the interannual variation of TWSA (Equation (4)) [67]. In the model, n is the maximum number of months of cumulative rainfall, and a is the efficiency of the impact of rainfall. The model simulates the interannual variation time series of rainfall by continuously adjusting the values of n and a. The RMSE, mean error (ME) [61], and correlation coefficient between the simulation results and the water storage capacity were used as indicators of the evaluation accuracy. After continuous adjustment, the values of the parameters n and a were determined, and we obtained the simulation results of the effect of rain on the interannual variation of TWSA.
where a is the scale factor (0 < a < 1), t is the time (month) of simulated rainfall, P i is the original rainfall value at the moment i, P 0 is the mean rainfall value over the whole period, and f (t) is the TWSA at moment t. n is taken as 6, 7, 8, . . . , 36, and a is taken as 0.1, 0.2, 0.3, . . . , 1.0. By continuously adjusting the values of n and a, the sequence of RMSE, ME, and correlation coefficients could finally be obtained, as shown in Figure 10a-c. It is shown that the correlation coefficient reached its maximum value (0.82) when n was 16, and the RMSE and ME reached their minimum values (14.50 mm, 12.04 mm) when n was also 16 and a was around 0.4. These results indicate that the possible ranges for n and a were 14-18 months and 0.3-0.4, respectively. To further obtain more reasonable values of n and a, we calculated the RMSE, ME, and correlation coefficients when n = 14/15/16/17/18 and a = 0.3/0.4, as shown in Table 2. We can see that RMSE reached its minimum and the correlation coefficient reached its maximum when n = 16 and a = 0.4. Thus, the maximum impact period of rainfall on TWS was about 16 months, with a scale factor a = 0.4.  Table 2. We can see that RMSE reached its minimum and the correlation coefficient reached its maximum when n = 16 and a = 0.4. Thus, the maximum impact period of rainfall on TWS was about 16 months, with a scale factor a = 0.4.  We used n = 16 and a = 0.4 to simulate the contribution of rainfall to the interannual TWSA from 2003 to 2016 (Figure 11). To be consistent with the simulated rainfall averages over the same timeframe, we implemented the GRACE results using a sliding average of 12 months. The RMSE value between GRACE and simulated result was 14.5 mm, and its RMSE difference was lower than 17% of the TWSA interannual variation, thus verifying the reliability of the model for the influence of rainfall on interannual land water storage changes in Northeast China given in this paper. In comparison, we calculated the simulation results of TWSA in Northeast China from 2003 to 2016  We used n = 16 and a = 0.4 to simulate the contribution of rainfall to the interannual TWSA from 2003 to 2016 ( Figure 11). To be consistent with the simulated rainfall averages over the same timeframe, we implemented the GRACE results using a sliding average of 12 months. The RMSE value between GRACE and simulated result was 14.5 mm, and its RMSE difference was lower than 17% of the TWSA interannual variation, thus verifying the reliability of the model for the influence of rainfall on interannual land water storage changes in Northeast China given in this paper. In comparison, we calculated the simulation results of TWSA in Northeast China from 2003 to 2016 using data from the climate-driven water storage change model provided by Humphrey and Gudmundsson (2019) [68]. The results are shown in Figure 11, and all series were calculated on the basis of a 12-month sliding average. The model reconstructed by Humphrey and Gumundsson (GRACE-REC) showed a trend apparently larger than the long-term change of GRACE observations. We used a simpler reconstruction method than that of Humphrey and Gumundsson but yielded observations that were more consistent with GRACE.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 20 using data from the climate-driven water storage change model provided by Humphrey and Gudmundsson (2019) [68]. The results are shown in Figure 11, and all series were calculated on the basis of a 12-month sliding average. The model reconstructed by Humphrey and Gumundsson (GRACE-REC) showed a trend apparently larger than the long-term change of GRACE observations. We used a simpler reconstruction method than that of Humphrey and Gumundsson but yielded observations that were more consistent with GRACE.

Challenge of Investigating Snow Cover Mass Change Using GRACE Data
In Northeast China, snowfall and snow cover account for a large proportion of TWS under conditions of reduced evaporation and runoff activity in winter. Therefore, we could use GRACE satellite data to detect seasonal large-scale variation of snow mass in Northeast China. Although GRACE satellite data has been used to assess the evolution of snow mass at high latitude [41,42], few studies have been conducted on snowfall in Northeast China. The result of this paper validated the ability of GRACE to monitor snow mass in a middle-high latitude region. As described in Section 3.1, the seasonal series of snowfall, snow depth, SWE, and TWSA presented good consistency. However, the mass signal detected by GRACE is a combination of all the mass variations, with signals from adjacent regions influencing each other [60]. As shown in Figure 6b-d, the negative signal in the northern region of Harbin may have been influenced by the spatial resolution constraints of the GRACE data and signal leakage. While there are various methods for correcting signal leakage errors, such as the spatial constraints method [30] and the scale factor method [69,70], these methods are mainly used to calculate the spatial distribution or rate of variation of TWS and groundwater storage. In contrast, this paper investigated the relationship between TWS and precipitation (snowfall and rainfall) mainly on seasonal and interannual scales. It did not calculate the trend terms, and thus the

Challenge of Investigating Snow Cover Mass Change Using GRACE Data
In Northeast China, snowfall and snow cover account for a large proportion of TWS under conditions of reduced evaporation and runoff activity in winter. Therefore, we could use GRACE satellite data to detect seasonal large-scale variation of snow mass in Northeast China. Although GRACE satellite data has been used to assess the evolution of snow mass at high latitude [41,42], few studies have been conducted on snowfall in Northeast China. The result of this paper validated the ability of GRACE to monitor snow mass in a middle-high latitude region. As described in Section 3.1, the seasonal series of snowfall, snow depth, SWE, and TWSA presented good consistency. However, the mass signal detected by GRACE is a combination of all the mass variations, with signals from adjacent regions influencing each other [60]. As shown in Figure 6b-d, the negative signal in the northern region of Harbin may have been influenced by the spatial resolution constraints of the GRACE data and signal leakage. While there are various methods for correcting signal leakage errors, such as the spatial constraints method [30] and the scale factor method [69,70], these methods are mainly used to calculate the spatial distribution or rate of variation of TWS and groundwater storage. In contrast, this paper investigated the relationship between TWS and precipitation (snowfall and rainfall) mainly on seasonal and interannual scales. It did not calculate the trend terms, and thus the influence of signal leakage errors was not considered in this paper. Another limitation of this study was that we did not utilize the subsurface parameters (such as soil moisture and groundwater) from models. A previous study by Niu et al. [42] showed that SWE can be derived from TWS in arctic regions where subsurface observations can be accurately provided by hydrological models. However, Northeast China has a highly changeable water storage system (Figure 11), together with a mostly unknown groundwater network in a populated and agricultural region, making it hard to employ hydrological models for a reliable estimate. Furthermore, as mentioned in Section 3.1, there was a good consistency among the time seasonal series of TWSA, snowfall, and rainfall, while notable deviations still exist. For instance, the seasonal variation series of snow depth, SWE, and TWSA showed notable deviations in the years 2012 and 2015, as shown in Figure 7a. This might have been a result of the increase in missing data since 2011, leading to increased errors, as well as the extreme drought/wet events [48]. Thus, we discussed the ability of GRACE satellite data to monitor large-scale snow cover in a qualitative sense. In addition, constrained by the spatial resolution of GRACE, the TWSA shown here was only the result of regional average, without detailed spatially distributed features. It is difficult to investigate snow cover changes at a high spatial resolution only using GRACE data.

Uncertainties in ENSO Impact on Rainfall and TWS in Northeast China
As mentioned in Section 3.2, the impact of ENSO events on precipitation and TWS in Northeast China is a long-term effect with a lag of about 2-3 years. If we adjusted the time lag to 1-12 months, as shown in Figure 9, there was no significant correlation among the precipitation, TWSA, and SOI series. This is consistent with the results of Chen et al. (2020), who analyzed the impact of two strong ENSO events in 2005-2017 on regional TWS changes in China [32]. Thus, the impact of ENSO events on precipitation and water TWS in Northeast China has a feeble short-term effect.
In addition, when changing the period of correlation analysis, we obtained different results. As shown in Figure 12 Oscillation (PDO)), which affect TWS in Northeast China, and was also limited by errors in the GRACE data [32,71]. Although different periods were chosen for the analysis, the conclusions were generally consistent. The TWSA and rainfall interannual time series lagged the SOI by ≈26 and ≈29 months, respectively, while the rainfall series lagged the TWSA by about 7 months. As mentioned above, there were three strong ENSO events worldwide in the period 2003-2016, and there was a strong correlation between TWSA and SOI series when strong ENSO events predominated in a given period. It can be confirmed that GRACE is only able to effectively monitor the impact of strong ENSO events on TWS in Northeast China.

Conclusions
In this study, we investigated the seasonal variation of rainfall and snowfall, the seasonal influence of rainfall on TWS, and the relationships between rainfall and TWS and ENSO phenomena in Northeast China on the basis of GRACE data, GlobSnow SWE products, ERA5-land monthly mean total precipitation, snowfall and snow depth data. We simulated the contribution of precipitation to terrestrial water storage in Northeast China from 2003 to 2016 and arrived at the following conclusions.
(1) The annual amplitude and phase of snowfall, snow cover, and TWSA in Northeast China were found to be generally consistent, reaching a maximum annual amplitude of about 6-12 mm around January. The annual amplitude of rainfall was about eight times that of snowfall and TWSA, but rainfall only had a moderate effect on the seasonal TWSA. Rainfall and snowfall controlled the seasonal variation of TWSA in Northeast China in different seasons. There were two peaks in the annual TWSA in Northeast China, with the winter peak mainly being caused by snowfall and the summer (usually August) peak being controlled by rainfall. Although the seasonal amplitude of snowfall was only 1/10 that of rainfall, it had a much greater impact on seasonal scale, while summer rainfall only brought about a sub-peak with a delay of 1 month.
(2) GRACE provided an effective means to monitor the stable seasonal snow cover over a large area. The annual amplitude and annual phase of TWSA calculated by GRACE were generally consistent with the results of snowfall and snow cover. At the 99% confidence level, the correlations between snow depth and SWE and TWSA series in February were 0.81 and 0.72, respectively. The

Conclusions
In this study, we investigated the seasonal variation of rainfall and snowfall, the seasonal influence of rainfall on TWS, and the relationships between rainfall and TWS and ENSO phenomena in Northeast China on the basis of GRACE data, GlobSnow SWE products, ERA5-land monthly mean total precipitation, snowfall and snow depth data. We simulated the contribution of precipitation to terrestrial water storage in Northeast China from 2003 to 2016 and arrived at the following conclusions.
(1) The annual amplitude and phase of snowfall, snow cover, and TWSA in Northeast China were found to be generally consistent, reaching a maximum annual amplitude of about 6-12 mm around January. The annual amplitude of rainfall was about eight times that of snowfall and TWSA, but rainfall only had a moderate effect on the seasonal TWSA. Rainfall and snowfall controlled the seasonal variation of TWSA in Northeast China in different seasons. There were two peaks in the annual TWSA in Northeast China, with the winter peak mainly being caused by snowfall and the summer (usually August) peak being controlled by rainfall. Although the seasonal amplitude of snowfall was only 1/10 that of rainfall, it had a much greater impact on seasonal scale, while summer rainfall only brought about a sub-peak with a delay of 1 month.
(2) GRACE provided an effective means to monitor the stable seasonal snow cover over a large area. The annual amplitude and annual phase of TWSA calculated by GRACE were generally consistent with the results of snowfall and snow cover. At the 99% confidence level, the correlations between snow depth and SWE and TWSA series in February were 0.81 and 0.72, respectively. The results for the average monthly changes in the spatial distribution of rainfall, snowfall, snow depth, and TWSA from December to March also showed great consistency.
(3) TWS and rainfall in Northeast China were found to be affected by the ENSO phenomenon on the interannual scale. ENSO drove changes in TWS through precipitation, and then affected the interannual changes of TWS, which is a long-term effect with a lag of about 2-3 years. The correlation coefficient between rainfall and ENSO was 0.64, with a lag of about 26 months. The correlation coefficient between TWS and ENSO was 0.66, with a lag of about 29 months. The correlation coefficient between rainfall and TWSA was 0.85, with a lag of about 7 months.
(4) The effect of rainfall on TWS was the cumulative effect of rainfall. We constructed a model to estimate the contribution of rainfall to the interannual variation of TWS from 2003 to 2016 on the basis of TWS and rainfall data from the period 2003-2016.