Mass Balance of the Greenland Ice Sheet from GRACE and Surface Mass Balance Modelling

The Greenland Ice Sheet (GrIS) is losing mass at a rate that represents a major contribution to global sea-level rise in recent decades. In this study, we use the Gravity Recovery and Climate Experiment (GRACE) data to retrieve the time series variations of the GrIS from April 2002 to June 2017. We also estimate the mass balance from the RACMO2.3 and ice discharge data in order to obtain a comparative analysis and cross-validation. A detailed analysis of long-term trend and seasonal and inter-annual changes in the GrIS is implemented by GRACE and surface mass balance (SMB) modeling. The results indicate a decrease of −267.77 ± 8.68 Gt/yr of the GrIS over the 16-year period. There is a rapid decline from 2002 to 2008, which accelerated from 2009 to 2012 before declining relatively slowly from 2013 to 2017. The mass change inland is significantly smaller than that detected along coastal regions, especially in the southeastern, southwestern, and northwestern regions. The mass balance estimates from GRACE and SMB minus ice discharge (SMB-D) are very consistent. The ice discharge manifests itself mostly as a long-term trend, whereas seasonal mass variations are largely attributed to surface mass processes. The GrIS mass changes are mostly attributed to mass loss during summer. Summer mass changes are highly correlated with climate changes.

To estimate the GrIS mass balance changes, three basic methods can be employed. The total mass balance (MB) has two main components that deal with the surface processes of precipitation, melt, and runoff (which is called surface mass balance (SMB), which is the net balance between the processes of accumulation and ablation on a glacier's surface) and the dynamic component related to glacier flow, including calving processes (i.e., ice discharge (D)). The first method is based on computing the diffirence of all the processes of accumulation and ablation by MB = SMB−D [1,17]. Unfortunately, the surface mass balance and ice discharge are affected by uncertainties from model simulations [18][19][20].
Unfortunately, the surface mass balance and ice discharge are affected by uncertainties from model simulations [18][19][20]. The second method involves the analysis of time-variable satellite gravity data for monitoring surface mass changes, such as GRACE. The application of this method is, however, hindered by a low spatial resolution and a leakage effect, in addition to large uncertainties in the glacial isostatic adjustment (GIA) correction [21,22]. The third method uses satellite altimetry for the monitoring of elevation changes. The measured elevation changes are then converted to equivalent mass changes by taking into consideration the ice/snow density [5,23,24]. In this case, the accuracy is affected especially by density uncertainties of snow and firn ice.
To complement the above deficiencies of each of these techniques [25,26], uncertainties could partially be mitigated by their combined usage. In this study, we used only the first two methods (described above) to obtain a comparative analysis and cross-validation while disregarding the third method that uses altimetry data. Because the elevation changes in the ice sheet measured by altimetry could be converted into corresponding mass changes based on adopting a density model of ice/snow. Unfortunately, the large uncertainties in the adopted density model introduce large errors in the estimates of ice mass changes. The modeled SMB variations need observational validation, such as GRACE. This is because the space-geodetic data do not exclusively reflect the SMB, but also ice dynamics, glacial isostatic adjustment, and other phenomena. Mass changes in the GrIS can be caused by variations in the glacial dynamic mass balance (ice discharge) and the SMB. However, the GRACE observations cannot directly separate these two physical causes, as they detect the sum of entire mass changes with a limited spatial resolution (around 330 km at the equator). The GrIS inter-annual variations are mainly attributed to SMB fluctuations that can be estimated more accurately by using a high-resolution SMB model, such as the SMB model RACMO2.3p2 [27] provided with a 1 km resolution that is used in this study.
In this paper, we explore the combination of GRACE and SMB-D modeling. Firstly, we investigate and compare estimates of the GrIS mass balance changes obtained from GRACE timevarying gravity field and the SMB model and the ice discharge from 2002 to 2017. We then provide a detailed analysis of seasonal changes, inter-annual and long-term trends from GRACE. We compare the results with annual mass balance of SMB-D (SMB from RACMO2.3p2 and ice discharge data from Mankoff et al. 2019 [28]) and GRACE to determine annual mass changes. Based on this, we conclude on the nature and statistical significance of the changes in the ice mass balance of the Greenland Ice Sheet.

Materials and Methods
The input data and models used for estimating the mass balance changes in Greenland and their further classification in terms of various climate variables are briefly described in this section.

Materials and Methods
The input data and models used for estimating the mass balance changes in Greenland and their further classification in terms of various climate variables are briefly described in this section.

GRACE Data
GRACE observables provide information about spatio-temporal variations in the gravitational field that reflect the Earth's mass transport. Since GRACE's launch on the 17th of March in 2002, the official GRACE Science Data System continuously releases monthly gravity solutions from different processing centers, namely Geoforschungs Zentrum Potsdam (GFZ), the Center for Space Research (CSR) at the University of Texas in Austin, the Jet Propulsion Laboratory (JPL), and the Institute of Geodesy in Graz (ITSG). The ITSG-GRACE2018 is the latest GRACE-only gravity model that provides unconstrained monthly and Kalman-smoothed daily solutions [29]. To reduce the noise in gravity field solutions within the available scatter of solutions, we employed the Level 2 monthly spherical harmonic coefficients from these four data centers' solutions (CSR-RL06, GFZ-RL06, JPL-RL06, and ITSG-2018). These coefficients are compiled with a spectral resolution complete to the spherical harmonic degree of 60. The monthly gravity coefficients are provided from April 2002 to June 2017.
The monthly GRACE solutions consist of the (fully normalized) spherical harmonic coefficients C nm and S nm of degree n and order m. The terrestrial water storage anomalies (TWS) over land (in terms of the equivalent water thickness) can directly be computed from these coefficients for a particular time period t (typically choosing a monthly period) according to the following expression [30]: where ∆σ(θ, λ) is the mass anomaly (in terms of the equivalent water thickness) at point (θ, λ), of which horizontal position is described by spherical co-latitude θ and longitude λ; ρ is the Earth's mean density; ρ w is the freshwater density; R is the Earth's equatorial radius; P nm are the (fully normalized) Legendre-associated functions of degree n and order m; k n are the (degree-dependent) Love numbers; and W n denotes the (degree-dependent) kernel functions of the Gaussian filter. The low-degree terms of spherical harmonic coefficients of GRACE solutions need to be pre-processed. The first-degree spherical harmonic coefficients that could not directly be detected by GRACE were determined from combining the GRACE data with numerical ocean models [31]. The second-zonal spherical harmonic coefficient C 20 (due to the Earth's flattening) was determined based on the analysis of the Satellite Laser Ranging measurements [32].
To reduce the noise in gravity field solutions, we adopted mean estimates of the GrIS mass changes based on averaging values from the CSR-RL06, GFZ-RL06, JPL-RL06, and ITSG-2018 models (more details are given in Section 3.1.1). We then estimated the terrestrial water storage changes from the monthly GRACE coefficients between April 2002 and June 2017, while interpolating missing monthly data over this period from the two adjacent months. The residual Stokes coefficients were obtained after removing the mean gravity field from January 2003 to December 2016. In order to reduce systematic and correlated errors in the GRACE measurements, we applied a Gaussian filter for the 500 km swath width as well as the P4M6 destriping [33]. Since the GRACE monthly solutions have a limited spatial resolution and are spatially filtered in order to reduce data noise, a leakage effect suppresses the actual gravitational signal attributed to mass variations. This effect causes inaccurate results, particularly in ocean-land areas [34]. It significantly attenuates the amplitudes and biases of mass loss estimates in Greenland [10,[34][35][36][37]. To reduce the ocean-land leakage effect, we applied a forward modeling technique [34] of GRACE monthly water storages changes results. By taking into account annual and semi-annual signals, we used an ordinary least-squares fit at each grid point to estimate the long-term change rate and seasonal changes.

Surface Mass Balance Model and Ice-Discharge Data
We used the SMB monthly means with a 1 km resolution from the Regional Atmospheric Climate Model (RACMO) 2.3p2 to obtain the surface mass changes in Greenland. The RACMO was prepared by the Institute for Marine and Atmospheric research Utrecht (IMAU) over the period from 1958 to 2016 [38]. The SMB is defined as the snow accumulation from which the surface ablation (erosion and sublimation, runoff, wind, and snow) is subtracted. We note that those seasonal mass variations in Greenland mainly reflect signals related to SMB. It is also worth mentioning that the SMB anomalies plus the ice discharge signal represent the mass changes detected by GRACE. The spatial resolution of the RACMO 2.3p2 model is about 11 km. In Greenland, the SMB represents about 50% of the ice loss signal [39]. Within the frame of RACMO2.3p2 products, the physical snow/ice surface model is used to calculate the time-varying surface albedo as a function of the ice cover performance. This allows a better representation of the SMB processes, such as the melt water penetration or re-freezing. The mass loss due to a solid ice discharge across the ice sheet grounding line estimates used in this paper were prepared by Mankoff et al. (2019) [28].

GIA Model
In order to estimate the mass balance trends from GRACE and interpret them as changes in the water content of hydrologic basins, ocean bottom pressure, or ice sheet mass, the GIA trend has to be taken into consideration. The two main constituents in any GIA model are the ice (deglaciation) history and the viscosity profile in the mantle. In this study, we used the GIA model of Geruo et al. (2013) [40] to correct the GRACE estimates. This model was compiled based on the ICE-5G deglaciation history model, the VM2 viscosity model, and the Preliminary Reference Earth Model (PREM)-based elastic structure [41].

Secular Mass Change
As stated above, after replacing the low-degree terms, applying spatial filters, and correcting the leakage and GIA effect to the GRACE solutions, we obtained the monthly mass changes of the GrIS. The total mass change (TWS-R) of the GrIS represents the sum of mass changes in each grid point in GrIS. We analyzed the monthly GRACE estimates from four different centers. As seen in Figure 2, these solutions do not differ significantly before 2013. Unfortunately, the GRACE observations deteriorated since 2013 due to aging of instruments on satellites. This resulted in large differences between individual solutions caused by differences in the processing strategies, including differences in the introduced models and inversion methods. To mitigate the uncertainties of individual GRACE solutions from the GFZ, CSR, JPL, and ITSG processing centers, we used a simple average of these four solutions.  As seen in Figure 2, the rates of mass change before and after 2013 differ, partly due to the aforementioned reasons. Apparent transient changes during the entire time span exist. To get a longterm trend, we first deducted the annual and semi-annual terms and 161-day cyclical fluctuations (which is due to GRACE tidal alias period for S2) in time series using Equation (1). Then, we used the R Package (segmented) software to estimate the regression models with unknown break points of the As seen in Figure 2, the rates of mass change before and after 2013 differ, partly due to the aforementioned reasons. Apparent transient changes during the entire time span exist. To get a long-term trend, we first deducted the annual and semi-annual terms and 161-day cyclical fluctuations (which is due to GRACE tidal alias period for S2) in time series using Equation (1). Then, we used the R Package (segmented) software to estimate the regression models with unknown break points of the time series of mass change [42]. It allows estimating linear and generalized linear models (and virtually any regression model) having one or more segmented relationships in a linear predictor [43]. to June 2017, the melting rate slowed down to only −148 Gt/yr due to the influence of atmospheric circulation in the northern hemisphere. This finding could be verified by published results that indicate a growing trend of glaciers in Greenland in recent years. The Jakobshavn glacier, for instance, has been Greenland's fastest-flowing and fastest-thinning glacier for the last 20 years. It is now flowing more slowly, thickening and advancing towards the ocean, instead of retreating further inland [44]. The apparent mass change rates derived from GRACE with low-degree terms correction and spatial filters, but without the leakage effect correction, are shown in Figure 4a. Figure 4b shows the recovered mass change rates derived from the monthly recovered water storage changes after applying a forward modeling technique to the monthly water storage changes. We note that a detailed explanation of the forward modeling technique can be found in [34]. To validate the forward modeling results, we estimated the GrIS mass change rates by applying the GRACE CSR RL05 mascon solution [45]. This result is shown in Figure 4c. As seen in Figure 4a,b, the leakage effect to GRACE results is significant. After applying the forward modeling technique to reduce the oceanland leakage, we observe considerable changes in the spatial distribution pattern of mass change rates. The forward modeling and mascon solutions exhibit similar spatial patterns (Figure 4b,c). Before applying the leakage correction (Figure 4a), due to the use of truncation and filtering methods, the mass loss occurs mainly in the northwest and southeast parts of Greenland, including the south inland area, with mass changes at the rate varying from −10 to about −5 cm/yr. After applying the leakage correction (Figure 4b), the mass ablation signal leaking to the adjacent ocean and inland ice sheets in the southeast and northwest coastal regions are characterized by many overflowing glaciers and fast ice flow rates. Along the southeast and northwest coastal regions, the mass ablation reached an alarming rate of −20 cm/yr, with some localized mass changes even reaching −30 cm/yr. Despite the mass change decreases along the northern coast, the trend kept at about 15 cm/yr. The mass change rate inland is significantly smaller than that detected along coastal regions. Coastal areas are also experiencing a more rapid decrease in total mass. The ocean-land leakage thus represents about  Over the whole period from April 2002 and June 2017, we estimated a melting at the rate of −267.7 ± 8.7 Gt/yr based on a least-squares fitting, while taking seasonal changes into account. This amount corresponds to a global sea level rise at the rate of 0.74 mm/yr. We also estimated a quadratic trend (i.e., acceleration), yielding an acceleration rate of −12.1 ± 2.6 Gt/yr 2 between 2002 and 2017, while taking the secular trend (of −266.8 ± 8.0 Gt/yr) and seasonal changes into account.
The apparent mass change rates derived from GRACE with low-degree terms correction and spatial filters, but without the leakage effect correction, are shown in Figure 4a. Figure 4b shows the recovered mass change rates derived from the monthly recovered water storage changes after applying a forward modeling technique to the monthly water storage changes. We note that a detailed explanation of the forward modeling technique can be found in [34]. To validate the forward modeling results, we estimated the GrIS mass change rates by applying the GRACE CSR RL05 mascon solution [45]. This result is shown in Figure 4c. As seen in Figure 4a,b, the leakage effect to GRACE results is significant. After applying the forward modeling technique to reduce the ocean-land leakage, we observe considerable changes in the spatial distribution pattern of mass change rates. The forward modeling and mascon solutions exhibit similar spatial patterns (Figure 4b,c). Before applying the leakage correction (Figure 4a), due to the use of truncation and filtering methods, the mass loss occurs mainly in the northwest and southeast parts of Greenland, including the south inland area, with mass changes at the rate varying from −10 to about −5 cm/yr. After applying the leakage correction (Figure 4b), the mass ablation signal leaking to the adjacent ocean and inland ice sheets in the southeast and northwest coastal regions are characterized by many overflowing glaciers and fast ice flow rates. Along the southeast and northwest coastal regions, the mass ablation reached an alarming rate of −20 cm/yr, with some localized mass changes even reaching −30 cm/yr. Despite the mass change decreases along the northern coast, the trend kept at about 15 cm/yr. The mass change rate inland is significantly smaller than that detected along coastal regions. Coastal areas are also experiencing a more rapid decrease in total mass. The ocean-land leakage thus represents about 40% of the GRACE total signal.

Seasonal Mass Changes
The GRACE gravity datasets were corrected for tidal and non-tidal effects. Since glaciers cover 87% of the whole Greenland landmass, the GRACE time-varying gravitational field reflects the seasonal and long-term mass variations in ice and snow. Figure 5 shows the annual and semi-annual changes in ice and snow in Greenland between April 2002 and June 2017. The annual amplitude along coastal areas is relatively large, especially along the southeast and southwest coastal regions (Figure 5a,b). There are many outbound glaciers along the south coast. The coastal mass balance is also influenced by the inland mass balance, particularly by the outflow glaciers along the south coast, yielding larger annual amplitudes in these areas. The annual amplitude along the north coast is less than 10 cm. This is due to the lower average temperatures at higher latitudes as well as the smaller seasonal temperature variations. Since the inland area of the southern Greenland is closer to the coast than in the north, the annual amplitude there is less than 2 cm while about 5 cm in the northern inland. In the interior of Greenland, the accumulation of snow and the snow/ice melting largely cancel out, maintaining a relatively stable mass balance. Similarly, the semi-annual amplitude along coastal areas is much larger than inland. changes in ice and snow in Greenland between April 2002 and June 2017. The annual amplitude along coastal areas is relatively large, especially along the southeast and southwest coastal regions ( Figure  5a,b). There are many outbound glaciers along the south coast. The coastal mass balance is also influenced by the inland mass balance, particularly by the outflow glaciers along the south coast, yielding larger annual amplitudes in these areas. The annual amplitude along the north coast is less than 10 cm. This is due to the lower average temperatures at higher latitudes as well as the smaller seasonal temperature variations. Since the inland area of the southern Greenland is closer to the coast than in the north, the annual amplitude there is less than 2 cm while about 5 cm in the northern inland. In the interior of Greenland, the accumulation of snow and the snow/ice melting largely cancel out, maintaining a relatively stable mass balance. Similarly, the semi-annual amplitude along coastal areas is much larger than inland.
We see that the semi-annual amplitude along the west coast is larger than that along the east coast. This pattern differs from a spatial distribution pattern of the annual amplitude. Figure 5c,d shows the annual and semi-annual phase changes in the mass of ice and snow. Overall, the mass accumulation reaches the maxima in May and June, followed by the mass ablation until September and October. The phase (in Figure 5c,d) indicates in which month the mass reaches the maxima. As seen in Figure 5c, the mass along the coastal areas in the southeast and west reaches the maxima between June and July, while inland ice mass reaches the maxima between October and December. This pattern is explained by the fact that inland the annual temperature is mostly below 0 °C [46,47]. Hence, the total mass budget there largely reflects a snow accumulation, while the contribution of ice discharge is minor. The total mass thus reaches the maxima (June−July) in coastal areas before reaching the maxima inland (October−December). As seen in Figure 6, the average mass budget in Greenland reaches the maxima in May, followed by the snow and ice melting until December. This surface process propagates into the overall decrease in mass budget, reaching the minima in October. The annual average mass accumulation is much (a) 8   We see that the semi-annual amplitude along the west coast is larger than that along the east coast. This pattern differs from a spatial distribution pattern of the annual amplitude. Figure 5c,d shows the annual and semi-annual phase changes in the mass of ice and snow. Overall, the mass accumulation reaches the maxima in May and June, followed by the mass ablation until September and October. The phase (in Figure 5c,d) indicates in which month the mass reaches the maxima. As seen in Figure 5c, the mass along the coastal areas in the southeast and west reaches the maxima between June and July, while inland ice mass reaches the maxima between October and December. This pattern is explained by the fact that inland the annual temperature is mostly below 0 • C [46,47]. Hence, the total mass budget there largely reflects a snow accumulation, while the contribution of ice discharge is minor. The total mass thus reaches the maxima (June−July) in coastal areas before reaching the maxima inland (October−December). As seen in Figure 6, the average mass budget in Greenland reaches the maxima in May, followed by the snow and ice melting until December. This surface process propagates into the overall decrease in mass budget, reaching the minima in October. The annual average mass accumulation is much smaller than the corresponding mass loss. The average change in ice sheet mass is positive in spring and negative during the rest of the year (Figure 6b). We consider September as the start of a new annual cycle for the GrIS. Meanwhile, the mass of ice sheet decreases from spring to autumn and then increases from autumn to winter, which is consistent with the monthly ice sheet mass changes. The GrIS largely gains snow from September or October, accumulating ice through autumn, winter and into spring. As the year warms up in late spring, the ice sheet begins losing more ice through surface melt than it gains from fresh snowfall. The mass loss normally begins at the end of May. This melting season generally continues until August or September.
Water 2020, 12, x FOR PEER REVIEW 9 of 16 smaller than the corresponding mass loss. The average change in ice sheet mass is positive in spring and negative during the rest of the year (Figure 6b). We consider September as the start of a new annual cycle for the GrIS. Meanwhile, the mass of ice sheet decreases from spring to autumn and then increases from autumn to winter, which is consistent with the monthly ice sheet mass changes. The GrIS largely gains snow from September or October, accumulating ice through autumn, winter and into spring. As the year warms up in late spring, the ice sheet begins losing more ice through surface melt than it gains from fresh snowfall. The mass loss normally begins at the end of May. This melting season generally continues until August or September.

Mass Balance from RACMO2.3 and Ice Discharger Modeling
Regional atmospheric climate models, derived from passive microwave data, have indicated a changing ratio between surface melt and precipitation, leading to an increasing imbalance between accumulation and ablation over the past two decades. We extracted the SMB information from the regional atmospheric climate model (i.e., RACMO2.3p2) for comparison with that from GRACE. The mean surface mass changes from RACMO2.3p2 are plotted in Figure 7. The resulting changes in surface mass are characterized by a significant decrease in the ablation and accumulation zones, which mirrors the change in run-off, as well as a partial increase in the interior owing to an increase in snowfall. We observe an increase in gradients in surface mass changes along the perimeter of the ice sheet. Significant ablation zones are present along the southwest and north coastal areas. The ice discharge is definitely negative. Hence, the situation in which the SMB becomes persistently negative in these ablation zones leads to a definite negative mass balance even if the ice discharge is zero. A large snow accumulation occurs along the southeast coast. A small snow accumulation occurs extensively inland.

Mass Balance from RACMO2.3 and Ice Discharger Modeling
Regional atmospheric climate models, derived from passive microwave data, have indicated a changing ratio between surface melt and precipitation, leading to an increasing imbalance between accumulation and ablation over the past two decades. We extracted the SMB information from the regional atmospheric climate model (i.e., RACMO2.3p2) for comparison with that from GRACE. The mean surface mass changes from RACMO2.3p2 are plotted in Figure 7. The resulting changes in surface mass are characterized by a significant decrease in the ablation and accumulation zones, which mirrors the change in run-off, as well as a partial increase in the interior owing to an increase in snowfall. We observe an increase in gradients in surface mass changes along the perimeter of the ice sheet. Significant ablation zones are present along the southwest and north coastal areas. The ice discharge is definitely negative. Hence, the situation in which the SMB becomes persistently negative in these Water 2020, 12, 1847 9 of 16 ablation zones leads to a definite negative mass balance even if the ice discharge is zero. A large snow accumulation occurs along the southeast coast. A small snow accumulation occurs extensively inland. The surface mass balance (SMB) was computed as the sum of precipitation and evaporation from which the runoff was subtracted. The precipitation includes rainfall and snowfall (while rainfall is minor in Greenland). The runoff generally refers to surface water runoff produced by surface ablation, the iceberg calving caused by the dynamic process of the ice sheet, and the bottom melting of the ice frame caused by ocean thermal processes. The positive evaporation represents the snow condensation, while negative evaporation represents the snow evaporation. Table 1 gives the monthly surface mass changes from RACMO2.3p2 during the investigated period of 15 years. The amount of mass accumulation or ablation per month varies due to climate variability, not to mention the amount of mass change during the period. The annual fluctuations in surface mass changes are also inconsistent for different years. The SMB in summer is negative, meaning that the amount of mass loss is greater than that of mass accumulation, with the largest mass loss in July. Throughout the investigated period, the smallest mass accumulation occurred in 2012 and the largest in 2013.

Comprarison between GRACE and SMB-D
We compared the GRACE time series of mass change with that from the RACMO2.3p2 cumulative SMB anomalies and ice discharge changes ( Table 2). The mass loss by solid ice discharge across the ice sheet grounding line (D) (column 9 in Table 1) estimates used in this paper are from Mankoff et al. (2019) [28]. As mentioned above, the ice sheet mass balance (MB) is controlled by the difference between net mass gains at the surface, expressed by the SMB, and mass loss by solid ice discharge across the ice sheet grounding line (D), so that MB = SMB−D. The annual mass loss due to ice discharge was without significant changes between 2002 and 2016, with an average of annual ice discharge of 491 Gt each year. The ice discharge thus manifests itself mostly as a long-term trend, whereas seasonal mass variations are largely attributed to surface processes. The surface mass balance (SMB) was computed as the sum of precipitation and evaporation from which the runoff was subtracted. The precipitation includes rainfall and snowfall (while rainfall is minor in Greenland). The runoff generally refers to surface water runoff produced by surface ablation, the iceberg calving caused by the dynamic process of the ice sheet, and the bottom melting of the ice frame caused by ocean thermal processes. The positive evaporation represents the snow condensation, while negative evaporation represents the snow evaporation. Table 1 gives the monthly surface mass changes from RACMO2.3p2 during the investigated period of 15 years. The amount of mass accumulation or ablation per month varies due to climate variability, not to mention the amount of mass change during the period. The annual fluctuations in surface mass changes are also inconsistent for different years. The SMB in summer is negative, meaning that the amount of mass loss is greater than that of mass accumulation, with the largest mass loss in July. Throughout the investigated period, the smallest mass accumulation occurred in 2012 and the largest in 2013.

Comprarison between GRACE and SMB-D
We compared the GRACE time series of mass change with that from the RACMO2.3p2 cumulative SMB anomalies and ice discharge changes ( Table 2). The mass loss by solid ice discharge across the ice sheet grounding line (D) (column 9 in Table 1) estimates used in this paper are from Mankoff et al. (2019) [28]. As mentioned above, the ice sheet mass balance (MB) is controlled by the difference between net mass gains at the surface, expressed by the SMB, and mass loss by solid ice discharge across the ice sheet grounding line (D), so that MB = SMB−D. The annual mass loss due to ice discharge was without significant changes between 2002 and 2016, with an average of annual ice discharge of 491 Gt each year. The ice discharge thus manifests itself mostly as a long-term trend, whereas seasonal mass variations are largely attributed to surface processes. The tendencies of mass balance estimated from GRACE and from the SMB minus ice discharge (SMB-D) are consistent. The mass balance results from SMB-D also indicate that there was a rapid and sharp mass loss over the period 2009-2012. This result agrees well with the GRACE results shown in

Discussions and Summary
We have estimated the GrIS mass changes based on combining the GRACE monthly solutions (over the whole period between April 2002 and June 2017) released by four different agencies in order to reduce the model errors. We used the high-resolution SMB product RACMO2.3p2 [28]. This newest SMB data product provides a more detailed representation of the mass changes in narrow ablation zones, glaciers, and disconnected ice caps. It complements GRACE's inability to detect small glaciers. These results, based on the combined processing of the GRACE and SMB data, revealed that the GrIS total mass within the entire investigated period decreased at an annual rate of −267.77 ± 8.68 Gt/yr, and this accelerated by −12.12 ± 2.59 Gt/yr 2 , which is comparable with previous studies (e.g., the range −252 Gt/yr to −274 Gt/yr from Groh et al. (2019) [48]; −255 Gt/yr from the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE) team (2019) [49]). During the same period, the Antarctic ice sheet's mass change rate was −162.91 ± 5.09 Gt/yr, with an acceleration of 10.12 ± 4.3 Gt/yr 2 . These figures correspond to the global sea level rise at a rate of 1.16 mm/yr (0.72 mm/yr due to the GrIS and 0.44 mm/yr due to the Antarctic ice sheet.) The uncertainty of the GRACE results include errors in the forward modeling technique, formal errors of the GRACE estimated mass rate from the least-squares fit, and the GIA model errors. The SMB-D is limited by the lack of information on outlet glacier ice thickness, as well as uncertainty in interior surface mass balance models.
We used the R Package (segmented) software to analyze the unknown break points of the time series of mass change. In this way, we were able to analyze more realistically the fluctuations and trends in mass changes in different basins in Greenland than by using only simple linear and quadratic trend estimates [1,13,50]. Our result shows that the GrIS total mass during the whole investigated period has gone through a rapid decline (2002)(2003)(2004)(2005)(2006)(2007)(2008) The IBMIE team (2019) also found that since 1992, the ice loss progressively increased and the trend has reversed after 2012, with a progressive reduction in the rate of mass loss during the subsequent period.
The mass change rate inland is significantly smaller than that detected along coastal regions, especially in the southeastern, southwestern, and northwestern regions. Coastal areas are also experiencing a more rapid decrease in total mass. Our mass balance estimates from a regional climate model and ice discharge indicate that ice discharge manifests itself mostly as a long-term trend, whereas seasonal mass variations are largely attributed to surface mass processes.
We compared the GRACE time series of mass change with that from the cumulative SMB anomalies and ice discharge changes. The tendencies of mass balance estimated from GRACE and from the SMB minus ice discharge (SMB-D) agree with each other. The annual mass loss due to ice discharge was without Our results indicate that a continuous mass loss in Greenland is mainly attributed to the melting of ice sheets in summer. We found throughout 15 years that the largest mass loss occurred in 2012, while the smallest mass loss was detected in 2013. Therefore, we further focused our analysis on the summer mass loss caused by climate change. The North Atlantic oscillation (NAO) index is an important atmospheric indicator reflecting the state of the general atmospheric circulation over Greenland. It has been shown to be connected to extreme melting events in Greenland [51][52][53]. We used the monthly tabulated NAO index from the National Oceanic and Atmospheric Administration (NOAA) (https://www.ncdc.noaa.gov/teleconnections/nao/) that is standardized by the 1981-2010 climatology. By analyzing the relationship between these investigated phenomena of mass changes from GRACE and SMB-D and NAO index, we were able to understand better the mechanism of how climate changes influence the GrIS mass. As seen in Figure 8, the negative phase of the summertime NAO (sNAO) index (June-July-August), such as in 2008 and 2012, increased the prevalence of high pressure, enhancing the surface absorption of solar radiation and decreasing precipitation. It caused the migration of warm air from southern latitudes into western Greenland. These changes promoted higher air temperatures, a longer ablation season, and enhanced melt and run-off. The positive phase of the sNAO index, such as in 2013, increased precipitation and decreased air temperatures, leading to a shorter ablation season and smaller mass loss. An increased mass loss has been related to increased occurrences of the negative phases of NAO during summer, which favors an anticyclonic circulation at upper levels and warmer conditions over Greenland at lower levels through the advection of warm air along its western coast. NAO is a large-scale climate driver, and it influences mass changes in Greenland by controlling regional temperature, precipitation, and atmospheric circulation. We also found that the highest summer temperature of Greenland was in 2012 over the whole investigated period of 16 years, and the amount of precipitation after 2013 was higher than in previous years. This explains, to a large extent, why in 2013 we detected a slowdown in the mass loss. As seen in Figure 8, the largest mass loss in Greenland occurred in 2012, which was followed by the smallest loss in 2013.