High Resolution Mapping of Ice Mass Loss in the Gulf of Alaska From Constrained Forward Modeling of GRACE Data

The resolution of Gravity Recovery And Climate Experiment (GRACE) Terrestrial Water Storage (TWS) change data is too low to discriminate mass variations at the scale of glaciers, small ensemble of glaciers, or icefields. In this paper, we apply an iterative constraint modeling strategy over the Gulf Of Alaska (GOA) to improve the resolution of ice loss estimates derived from GRACE. We assess the effect of the most influential parameters such as the type of GRACE solution and the degree of heterogeneity of the distribution map over which the GRACE data is focused. Three GRACE solutions from the most common processing strategies and three ice distribution maps of resolutions ranging from 55,000 to 20,000 km2 are used. First, we present results from a series of simulations with synthetic data and a mix of synthetic/modeled data to validate the focusing strategy and we point out how inaccuracies arise while increasing the spatial resolution of GRACE data. Second, we present the recovery of the total GRACE-derived mass change anomaly at the scale of the GOA. At this scale, all solutions and distribution maps agree, showing ∼40 Gt/year of mean ice mass loss over the period 2002–2017. This result is similar to studies using GRACE solutions from the latest releases and time-series of more than 8 years. The first studies using GRACE data published during the 2005–2008 era generally overestimated the long-term ice mass loss. Third, we show results of the three resolutions tested to focus the mass anomaly. Using focusing units (mascon) of ∼30,000 km2 or larger, the focusing procedure provides reliable results with errors below 15%. Below this threshold, errors of up to 56% are observed.

The resolution of Gravity Recovery And Climate Experiment (GRACE) Terrestrial Water Storage (TWS) change data is too low to discriminate mass variations at the scale of glaciers, small ensemble of glaciers, or icefields. In this paper, we apply an iterative constraint modeling strategy over the Gulf Of Alaska (GOA) to improve the resolution of ice loss estimates derived from GRACE. We assess the effect of the most influential parameters such as the type of GRACE solution and the degree of heterogeneity of the distribution map over which the GRACE data is focused. Three GRACE solutions from the most common processing strategies and three ice distribution maps of resolutions ranging from 55,000 to 20,000 km 2 are used. First, we present results from a series of simulations with synthetic data and a mix of synthetic/modeled data to validate the focusing strategy and we point out how inaccuracies arise while increasing the spatial resolution of GRACE data. Second, we present the recovery of the total GRACE-derived mass change anomaly at the scale of the GOA. At this scale, all solutions and distribution maps agree, showing ∼40 Gt/year of mean ice mass loss over the period 2002-2017. This result is similar to studies using GRACE solutions from the latest releases and timeseries of more than 8 years. The first studies using GRACE data published during the 2005-2008 era generally overestimated the long-term ice mass loss. Third, we show results of the three resolutions tested to focus the mass anomaly. Using focusing units (mascon) of ∼30,000 km 2 or larger, the focusing procedure provides reliable results with errors below 15%. Below this threshold, errors of up to 56% are observed.

INTRODUCTION
Glaciers represent 68.9% of fresh water resources worldwide. In many regions of the world, people rely on glacier meltwater for agriculture, hydropower, industries, and municipal water requirements (Chen and Ohmura, 1990;Blanchon and Boissière, 2009). However, over the last decades, the glacier mass losses have raised concerns in and beyond the research communities. Climate change leads to important reductions in glacial water storage. Glaciers have an important influence on sea level rise; hence, their melt threatens the living environment of costal dwellings. Jin and Feng (2016) estimated the contribution of glacial melt to sea level change between 2003 and 2012 at 1.94 ± 0.29 mm/year. From 120,000 glaciers available in the World Glacier Inventory, Radić and Hock (2011) estimated that the total volume loss could be as much as 21 ± 6% by 2100, leading to a total sea level rise of 124 ± 37 mm.
Numerous studies focused on estimating the ice mass loss over specific continents, regions, or Mountain ranges. For example, Larsen et al. (2007) investigated glacier changes in southeast Alaska and northwest British Columbia over the period 1948-2000 and 1982/1987-2000, respectively. By combining the results from these periods, they estimated an average ice mass loss rate of around 15.03 ± 4 Gt/year (considering an ice density of 900 kg/m 3 ). In the Canadian Rocky Mountains, Castellazzi et al. (2019) estimated a total of 43 Gt of glacial mass loss over the period [2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015]. Over the entire Gulf Of Alaska (GOA) area, Gardner et al. (2013) found 50 ± 17 Gt/year of glacial melt by using spaceborne altimetry data (e.g., ICESat) over the period . Berthier et al. (2010 obtained 37.47 ± 8 Gt/year of glacier ice loss from Digital Elevation Models (DEM) for the period 1962period -2006period . Larsen et al. (2015 used airborne altimetry to estimate glacier mass loss rate over the period 1994-2013 and found 75 ± 11 Gt/year. Gravity data provides direct information over ice mass changes, as the link between gravity and mass is direct and requires no calibration. In situ gravity measurements are laborintensive, costly to acquire, and point-based; while satellite gravity data are limited in resolution due to the sensing distance. Since 2002, the United States (NASA) and the German (DLR) space agencies have led the Gravity Recovery And Climate Experiment (GRACE) mission. It aims at monitoring the variations of the Earth's gravity field with a temporal resolution of a few days to a month, with a spatial resolution of ± 400 km (Tapley et al., 2004;Ramillien et al., 2017). The variations in the Earth's mass distribution cause changes in the gravity field (Wahr et al., 1998). Thus, by mapping the variations of the gravitational field, GRACE monitors the Earth's mass distribution. The mass redistribution obtained from GRACE data contains changes in Terrestrial Water Storage (TWS), and oceanic mass (Wahr et al., 1998;Chen et al., 2006). Yirdaw et al. (2009) noted that the GRACE mission estimates could be used to monitor the rate of change of TWS over large spatial scales. Ramillien et al. (2017) further indicated that continental hydrology is one of the main applications of GRACE data. The variation of TWS aggregates changes in surface water, soil moisture, ground water, snow, and ice. Jin and Zou (2015) used GRACE data to estimate a high-precision glacier mass dynamics in Greenland. In the GOA region, GRACE data have already been used to estimate glacier mass loss (Tamisiea et al., 2005;Chen et al., 2006;Arendt et al., 2008Arendt et al., , 2009Arendt et al., , 2013Luthcke et al., 2008;Baur et al., 2013;Beamer et al., 2016;Wahr et al., 2016;Jin et al., 2017).
Two main types of GRACE Level-3 solutions are used in hydrological applications. There are unconstrained solutions, relying on de-stripping and Spherical Harmonics (SH) truncation, and constrained solutions often relying on regularization or stabilization. Among the later, mascon solutions, obtained after inverting the GRACE signal into a mass change for each spatial unit ("mascon") of a predefined grid, have become particularly popular over the last years (e.g., Save et al., 2016). Meanwhile, the unconstrained GRACE Spherical Harmonics (SH) solutions present errors at high spatial frequencies (e.g., N > 60 or 300 km) and North-South stripes mainly due to gravitational model corrections, instrument errors, and gaps in data coverage. Destriping, truncation, and filtering are usually applied on these Level-2 solutions to reduce high frequency errors and to eliminate stripes. The challenge of using GRACE data to investigate hydrological fluxes, such as mass variations at the glacier scale, lies in the low spatial resolution (300/400 km) and the inability to discriminate close masses. Leakage effects are inherent to the GRACE sensing strategy, and are accentuated by the truncation and filtering used to "clean" the monthly gridded mass changes. This makes GRACE data prone to large errors when considering areas below 200,000 km 2 (e.g., Longuevergne et al., 2010). To overcome this problem, several authors combined GRACE data with other sources of information (e.g., Castellazzi et al., 2018Castellazzi et al., , 2019. Different strategies such as scaling approach, additive approach, multiplicative approach, and unconstrained or constrained forward modeling approaches permit to partially restore the signal loss (Long et al., 2016). A constrained inversion method can also be applied to improve the spatial resolution of GRACE data (Farinotti et al., 2015;Long et al., 2016;Castellazzi et al., 2018). Chen et al. (2015) used forward modeling to restore the GRACE signal amplitude of Antarctic ice and glacier loss due to the noise reduction. Long et al. (2016) showed that a constrained forward modeling can recover the distribution details of the GRACE signal. Farinotti et al. (2015) estimated the glacier mass loss in the Tien Shan (China) by subtracting the non-glacier contributions from the total mass change. They used an inversion method with a priori information about the glacier spatial distribution in area subdivisions (i.e., mascon and sub-mascon). By improving GRACE spatial resolution, their results are comparable to those derived from altimetry data and glacial melt modeling (e.g., the Cold Regions Hydrological Model- Pomeroy et al., 2007-and the Distributed Enhanced Temperature Index Model-Hock, 1999. Studying groundwater depletion in Central Mexico, Castellazzi et al. (2016Castellazzi et al. ( , 2018 improved the GRACE spatial resolution using InSAR-derived ground-displacement maps to constrain the Forward Model (FM). They subtracted groundwater contributions from other components of TWS, built a distribution map of groundwater depletion using InSAR, and focused GRACE data over it. The potential of separating masses depends on the size, separation, and amplitude of the masses, as well as on the available constraining data and their efficiency in explaining the low resolution GRACE signal.
Although several authors have injected ancillary data into GRACE post-processing to improve the resolution, there is still a need to assess how constrained modeling improves the interpretation of GRACE data versus simplistic approaches such as spatial averaging (e.g., usually over watersheds or regions). While few studies already tested how constrained FM can help interpreting GRACE data, case studies are still lacking and more importantly, the limits of the procedure are still relatively unknown. In this study, we apply this procedure to build a high resolution ice mass loss map from GRACE data over a large range of melting glaciers in the GOA. In that perspective, we use several GRACE solutions, apply a spatial constraint to focus the signal, and provide a high spatial resolution map of ice mass loss over the study area.

STUDY AREA
Our study area is delimited by zone 11 (Figure 1). The area was selected to cover the glaciers of the GOA and to include a stretch of ∼300 km of the surroundings to account for spatial leakages inherent to the GRACE signal. During the last decade, different studies used GRACE data to estimate glacier mass across the same region (Figure 1).
There are discrepancies in mass losses estimated over the GOA. Most authors obtain values ranging from −47 to −110 Gt/year of glacial mass change (Figure 1 and Table 1). These studies consider different spatial extents and data timeperiods. For example, the total surface coverage considered by Tamisiea et al. (2005) is 701,000 km 2 , while Beamer et al. (2016) consider 420.300 km 2 . The first three studies show higher values of mass loss (Tamisiea et al., 2005;Chen et al., 2006;Luthcke et al., 2008) than those more recent (Baur et al., 2013;Beamer et al., 2016;Wahr et al., 2016;Jin et al., 2017). It is potentially due to the relatively short GRACE time-series available at that time and the higher uncertainties in the earlier releases of GRACE data. In addition, variations in the GRACE data processing strategy, including levels of filtering and Spherical Harmonic (SH) truncation, might also contribute to discrepancies between studies. For example, some studies used unconstrained solutions (Tamisiea et al., 2005;Chen et al., 2006;Baur et al., 2013;Jin et al., 2017). Other authors used mascon solutions at different spatial and temporal resolutions Luthcke et al., 2008Luthcke et al., , 2013Jacob et al., 2012;Beamer et al., 2016;Wahr et al., 2016). Tamisiea et al. (2005) used unconstrained solutions completed up to SH degree and order 70. Other authors used SH degree and order up to 60 (Chen et al., 2006;Baur et al., 2013;Jin et al., 2017). Jin et al. (2017) applied a Gaussian filter with a radius of 300 km while the aforementioned studies used a radius of 500 km. Using a scaling factor approach, Tamisiea et al. (2005) reduced the GRACE signal attenuation and leakage. A forward modeling approach was also used by several authors to restore signal leakages from coastal glaciers to the ocean (Baur et al., 2013;Jin et al., 2017).
FIGURE 1 | Glaciers of the GOA and footprints of the study area considered in studies using GRACE data to assess glacier melt. The study area considered here is identified as Zone 11. Some authors also used auxiliary data to complement GRACE observations. Remote sensing of ice surface elevation (e.g., ICESat NASA, airborne laser altimetry) was combined with GRACE data Jin et al., 2017). Meteorological models also contributed to GRACE science Wahr et al., 2016). Some authors combined remote sensing and meteorological data to enrich and validate glacier mass loss estimates Luthcke et al., 2013). Beamer et al. (2016) developed hydrological models and compared their results with the airborne altimetry and GRACE data.

Glacier Data and Mascon Delineation
Glaciers are delimited using the GLIMS Glacier Database released on 27/10/2017 and available online 1 (GLIMS andNSIDC, 2005, updated in 2013). It is a continuation of the World Glacier Inventory (Raup et al., 2007) which compiles different sources: satellite imagery data, historical information from maps, and aerial photography. Recently, the GLIMS glacier Database was merged with the Randolph Glacier Inventory (RGI). According to Pfeffer et al. (2014), this combination was performed to support the fifth report of the Intergovernmental Panel on Climate Change. While the GLIMS database may contain errors, it is, at least for our study area, in good agreement with the glacier inventory of Alaska and northwest Canada proposed by Kienholz et al. (2015). According to GLIMS, the total glacier cover in the GOA area was around 90,000 km 2 in 2013 (Figure 2).
To ease data handling and make this dataset compatible with other GRACE grids used in this study, the glacier distribution was resampled to a 0.25 • grid (∼28 km by ∼14 km pixels; Figure 2). The resampling is performed using the following steps: (1) creation of a global 0.25 • raster; (2) overlay of the glacier footprint (taken as polygons) from the GLIMS data on the raster; and (3) creation of the distribution map. In the latter step, we 1 http://nsidc.org/glims/ FIGURE 2 | Glacier distribution map used to focus GRACE trend maps. Three arbitrary mascon delineations are presented: (A) 5 mascons with an average area of ∼55,000 km 2 , (B) 9 mascons with an average area of ∼30,000 km 2 , and (C) 14 mascons with an average area of ∼20,000 km 2 .
consider that a pixel from step (1) contains glacial masses if its center is within a glacier polygon from GLIMS.
Separating the contribution from a set of small masses to a low resolution signal (such as GRACE) can be challenging and the success relies principally on the initial resolution, the final resolution, the size of the masses, the distance separating them, and the amplitude of the changes for each mass. Their spatial separation is important, close masses being harder to discriminate than those spatially distant. To test the separation of the signal into a set of contributing masses, we considered three different spatial distributions of mass at different resolutions ( Table 2). Each area delimiting a mass is referred to as a mascon (mass concentration area).
The GOA glacial mass losses are spatially heterogeneous due to the difference in size of glaciers, local hypsometry, and climate (i.e., temperature and precipitation; Frans et al., 2018;Huss and Hock, 2018). For simplification, we assumed a uniform distribution of ice mass loss within each mascon and test this assumption through a synthetic test (section "Simulation With Synthetic Data"). This assumption has been used in other areas prone to ice mass loss by Chen et al. (2009) and in Farinotti et al. (2015).

GRACE TWS Data
Three versions of Level-3 GRACE data in Water Thickness Equivalent (WTE) are considered. The first solution is the stabilized solution from the Space Geodesy Research Group (GRGS RL04 2 ) available up to SH degree and order 90. The second is an unconstrained solution from the Center for Space Research (CSR RL06) complete up to SH degree and order 96 and de-striped/filtered in a single step by applying the DDK8 filter (Kusche, 2007;Kusche et al., 2009). The third is the CSR RL05 regularized Mascon solution from Save et al. (2016). All solutions are truncated up to SH degree and order 90 to better perceive the differences in the data processing strategies without being influenced by the differences in resolution. The solutions are referred to as GRGS, T96DDK8, and CSR-MASC hereafter. These solutions follow different processing protocols and assumptions; hence, we consider that their discrepancy represents an estimation of the impact of the choice of the processing strategy. They comprise 150, 156, and 159 nearmonthly measurements, respectively, extending from 2002 to 2017. For consistency, the land mask applied to produce the CSR-MASC solution (Save et al., 2016) is applied to all trend maps. All processing is performed over the same 0.25 • grid.
2D trend maps are derived from the three near-monthly solutions. The trend is computed by subtracting an iteratively fitted sine curve of 1-year period. A 13-month moving average is then applied to smooth out the residual. These steps contribute 2 http://grgs.obs-mip.fr/grace to remove the seasonal variations and prevent contamination of trend estimates by seasonal signals. Finally, for each pixel, the slope of a fitted linear curve is interpreted as the mass change rate. Glacial Isostatic Adjustment (GIA) influences over the TWS trend need to be removed. CSR-MASC, as available online, is already corrected for GIA effects using estimates from Geruo et al. (2013). We used the same model (ICE-5G) to correct the two other solutions (GRGS and T96DDK8). Apart from glacier, Chen et al. (2006) also considered other water storage components such as snow and soil moisture. In our case, we assume that the influence of the seasonal snow component over the total mass loss trend is removed by fitting and subtracting a stationary curve to the GRACE TWS data. We also consider that, by using more than 15 years of TWS data, the interannual variability of the snow cover does not affect our trend estimates significantly.

Inversion Procedure
An inversion procedure is applied to focus the low resolution mass trend map produced using the GRACE data. First, it uniformly allocates a synthetic mass to all glacier pixels within each mascon; resulting in a synthetic high resolution ice mass loss map. Second, a FM is built by truncating this map similarly to GRACE data, through the spherical harmonic domain. The FM is assessed by analyzing the difference with the actual GRACE TWS trend map. The mass allocation is iteratively modified until the FM converges to the closest possible value of the actual GRACE TWS signal. An Absolute Error (AE) approach applied over the spatial domain (the trend map in WTE) is used to quantify the degree of similarity between the FM and the actual GRACE trend map.
A Pattern search solver is used to iteratively change the mass allocated within each mascon until finding the optimal set of masses which, after truncation, is the closest to the actual GRACE-derived trend map. The Pattern search algorithm is a non-gradient-based optimization algorithm. It makes adjustments to each parameter value independently in order to iteratively converge toward a stable solution (Hooke and Jeeves, 1961;Kolda et al., 2003;Hingray et al., 2014). In order to prevent the solver to stop over local minimums, two sets of initial values are tested for each calculation. Both sets of initial values are uniform, the first is largely under the expected mass distribution, and the second largely above. This allows testing two extreme starting points and verifying that the results are not affected by the arbitrary assumption of the initial values.
We tested the three GRACE solutions with each mass distribution map. In GRACE processing, it is often assumed that the highest source of uncertainty arising in TWS estimation is related to the arbitrary choice of the processing strategy. Indeed, there are numerous methods to process GRACE data ("solutions"), and there is no obvious way to decide which one to use for any given application (see, e.g., Castellazzi et al., 2018). The shape of the GRACE trend anomaly is influenced by the highest spatial frequencies of the GRACE signal, which are the most sensitive to the processing strategy and parameterization. Hence, we assume that if different solutions are similar after focusing, the results are not affected by the processing strategy and, consequently, are close to the true TWS signal. Conversely, when the focusing procedure is applied and the results are different for each GRACE solution, we consider that the trend maps are influenced by processing residues affecting the high resolution results.

Simulation With Synthetic Data
Synthetic tests are performed to assess the performance of the focusing procedure used to retrieve local or regional masses from low resolution GRACE data. In addition, as the GRACE signal is considered dominated by the ice mass change over the GOA, synthetic simulations are also used to verify if this assumption is true when considering masses at local or regional scales.
A synthetic map is created by allocating synthetic random values, close to the ice mass change rates expected across the GOA, to the distribution map (Figure 2). The mass change values are spread uniformly on glacier pixels within each mascon to reproduce a realistic high resolution glacier mass loss map. To account for spatially distributed signals, a diffuse mass trend map is derived by computing trends from monthly soil-moisture estimates (Noah model, GLDAS v2.1; Rodell et al., 2004) over the same time-period than the GRACE observations. The final map is a realistic model of cumulated signals from ice and distributed mass changes in realistic proportions. This map is truncated similarly to the three GRACE solutions, building a FM for GRACE trend map from known mass change rates.
First, the focusing procedure is applied using solely the mass distribution map presented in Figure 2. Finally, the high resolution mass loss values recovered through the procedure are compared with the values originally injected within each mascon. Second, the inversion procedure is performed after adding a background mascon covering the whole study area, with all pixel accounted for, except those corresponding to glaciers. The first test assesses the validity of the focusing procedure and the influence of diffuse masses over the local-scale mass retrieval. The second test verifies whether the local signal from glaciers is better recovered while accounting for surrounding low-frequency signals or not.

Simulation With Synthetic Data
A synthetic map is created by allocating synthetic masses to the glacier distribution map (Figure 3A) over which a diffuse mass signal is added (Figure 3B). This map, after truncation, is a low resolution GRACE-like mass trend map (Figure 3C) built from a known high resolution map. The masses injected onto the distribution map are realistic; they are of a similar order of magnitude to the masses expected across the GOA (considering both individual mascon and the entire area. Supplementary  Table S1 shows the synthetic masses allocated to each mascon, and the mass recovered after applying the inversion procedure, considering or not the effect of a diffuse mass. The main observations from this synthetic test are compiled in Table 3. Such test, consisting in retrieving synthetic masses at different resolutions, allows to observe the limits of the focusing procedure (D) Example of a high resolution glacier mass loss map retrieved by applying the inversion procedure over C and using the 5-mascon delineation shown in Figure 2A. This synthetic test allows to assess the inversion procedure and its ability to retrieve glacier mass losses at high resolution. Results are presented in Supplementary Table S1. and to better interpret the results obtained with real data presented in the subsequent sections.
We observe that the total mass is well recovered regardless of the number of mascons and of the consideration of a diffuse mass during the inversion ( Table 3). The Mean Absolute Relative Error (MARE) of the mass recovery at the mascon resolution ranges from ∼5 to ∼15%, and increases with the resolution of the mascon map. The maximum error goes up to ∼68% while recovering masses over small mascons (14-mascon delineation). We conclude that an acceptable accuracy is obtained while using mascon size up to ∼30,000 km 2 .
The differences observed in the total mass recovery, considering or not a diffuse mass, are less than 3% for the three mascon delineations. The maximum additional error related to the inclusion of a diffuse mass is ∼2% for the 5 and 9-mascon delineations and reaches ∼12% for the 14-mascon delineation. This suggests that the inclusion of a supplementary mascon expected to account for diffuse masses (low-frequency signals) did not help to recover the glacial masses. This result, along with observations from the literature (e.g., Jin et al., 2017), and the apparently strong TWS trend (Figure 3), supports the common assumption that GRACE TWS signal is largely dominated by glacier ice mass changes in the GOA. It also suggests that the Each test is performed with and without taking into account a diffuse mass, e.g., soil moisture storage changes. inversion performed with the actual GRACE TWS data does not need to consider a diffuse mass.

GOA-Scale Ice Mass Loss
GRACE TWS signal over the GOA shows a negative anomaly in all three solutions (Figure 4). They are all showing a very similar mass depletion across the GOA, in both amplitude and spatial extent (Figures 4A-C). The GRGS solution shows a slightly different spatial pattern at the eastern tip of the anomaly. Given that all solutions were corrected for GIA using the same model (Geruo et al., 2013) and that all solutions were at the same resolution (SH coefficients up to degree/order 90), the most likely reason for the perceptible differences is related to the de-striping strategy, which is different for all the three solutions considered.
The total mass loss rates were recovered by distributing mass loss values in each mascon until the FM reaches a minimum residual when subtracted from the GRACE TWS of each solution. Results are generally consistent regardless of the distribution considered ( Table 4).
The test with synthetic data shows that the error in recovering the total mass is ∼5% regardless of the mascon resolution.
On the one hand, using different solutions over the same distribution leads to a Coefficient Of Variation (COV) of ∼3.5% (Table 4). On the other hand, using different distributions with the same solution leads to a COV of ∼2.1%. The similarity between the solutions (Figure 4D and Table 4) proves a clear improvement in GRACE processing with the latest solutions from the 2016-2018 era such as CSR RL06. When considering an ample TWS anomaly like the GOA, solutions are usually leading to similar observations. However, it might not be the case when considering low amplitude signals which are closer to the noise level, and more sensitive to the processing strategy.
These results are compared with results from previous studies performed over the same area in Table 1. Jacob et al. (2012); Baur et al. (2013), and Wahr et al. (2016) found ice loss rates within the range 47-52 Gt/year, which is in good agreement with our results (Figure 5). Tamisiea et al. (2005); Chen et al. (2006), and Luthcke et al. (2008) found significantly higher ice loss rates. This can be related to the shorter TWS time-series and the larger uncertainties associated with the early releases of GRACE TWS data. Our results are close to the estimates from studies using   different methods (Berthier et al., 2010;Gardner et al., 2013;Larsen et al., 2015, values are listed in section "Introduction"). Differences in GRACE Level-3 data interpretation strategy can lead to discrepancies with estimates from other authors. The use of a constrained forward modeling strategy can lead to lower values, as it only recovers the signal spatially correlated to glaciers (see, e.g., Long et al., 2016, for an example related to groundwater depletion). The injection of a distribution constraint is known to reduce contributions from unwanted signals with contrasting spatial patterns. These signals cannot be focused onto the constraining distribution (Farinotti et al., 2015;Castellazzi et al., 2018). However, the distribution map is a simplification of the reality as it does not consider the intra-mascon heterogeneity of the ice loss signal. This might lead to uncertainties, depending on the level of unaccounted heterogeneity and its location on the distribution map.

Mascon-Scale Ice Mass Loss
In this section, we discuss the reliability of considering mass losses focused over mascon at different scales. For each inversion, the result corresponds to the best fit to reproduce the true GRACE TWS over the GOA. Total per-mascon masses are presented in Table 5, and the spatial patterns of the mean and standard deviation of glacier mass changes are shown in Figure 6. As discussed previously, we assume that the degree of similarity between the mass losses retrieved using the three solutions represents a good indicator of accuracy. The maximum value of variation within each mascon corresponds to ∼13, ∼25, and  For the 5-mascon delineation, 55,000 km 2 scale, in %   13  12  1  1  11  ---------For the 9-mascon delineation, 30,000 km 2 scale, in % 10 13 8 3 25 1 10 9 4 ----- For the 14-mascon delineation, 20,000 km 2 scale, in %   19  16  21  32  3  40  4  3  3  10  32  9 14 13 It corresponds to the standard deviation of the results obtained from each GRACE solution expressed as a percentage of the mean.
∼40% of the mean mass allocated, for the 5-, 9-, and 14-mascon delineations, respectively. It is observed that this variation tends to increase over the East/West extremes of the study area ( Table 6 and Figure 6). This corresponds to differences observed by comparing the low resolution GRACE TWS trend maps (Figure 4) which propagates into the inversion results. These variations due to differences in processing strategies increase with the spatial resolution ( Table 6); hence the estimates for the 14-mascon delineation have the highest level of uncertainty. We evaluate the accuracy of our mass loss estimates by comparison with those from other studies. This comparison is based on the area with the largest mass losses, which are the Saint Elias Mountains and Glacier Bay, according to Luthcke et al. (2013). This was confirmed by Arendt et al. (2013), who showed that the Saint Elias Mountains, Glacier Bay, and Juneau icefield regions present the largest amount of glacier mass losses. Jin et al. (2017) indicated that the Malaspina and Bering glaciers, belonging to the Saint Elias Mountains, have the largest rates of ice mass losses. Glaciers of the Saint Elias Mountains correspond to mascons 6, 7, and 10 in Luthcke et al. (2008; see Figure 7A) and mascon 6 in Arendt et al. (2013; see Figure 7B). Our estimates agree with the aforementioned studies which indicate that glaciers of the Saint Elias Mountains have the largest mass loss rates in the region (see Supplementary Tables S2, S3). Arendt et al. (2008) used the same delineation and GRACE TWS solution as Luthcke et al. (2008; see Figure 7A). By down-sampling the solution to the Saint Elias Mountains (cumulating rates from mascons 6, 7, and 10; Figure 7A in Luthcke et al., 2008), using the ratio of ice area in each mascon over the period 2003-2007, they obtained 20.6 ± 3.0 Gt/year of ice loss. Luthcke et al. (2008) calculated 36 ± 2.0 and 30 ± 2.0 Gt/year of ice mass loss by using GRACE data for 2003-2006 and 2003-2007, respectively. First, we compare the results from these two studies with our estimates over the period 2002-2017 (Table 7), and we create a temporal subset from our estimates to make our comparison insensitive to discordances in the time-period considered.
Comparing our results with two other studies Luthcke et al., 2008) after applying similar spatial FIGURE 7 | Overlay of mascon delineations from other studies over the glacier distribution map used in this study: (A) delineation used by Luthcke et al. (2008) and (B) delineation used by Arendt et al. (2013), in which Glacier Bay is included in mascon 6, corresponding to the Saint Elias Mountains (Supplementary Table S2).
Frontiers in Earth Science | www.frontiersin.org TABLE 7 | Mean mass losses in Gt/year considering three different time-periods over a spatial subset corresponding to mascons 6, 7, and 10 in Luthcke et al. (2008; see Figure 7A), which correspond to the Saint Elias Mountains. constraints, we observe that the results from Arendt et al. (2008) are similar to our estimates regardless of the mascon delineation considered. Results from Luthcke et al. (2008) are generally larger than our estimates, as do most other studies using GRACE data from the 2003-2006 time-period. However, we note similarities for mascon numbers 5, 6, and 7, at the center of the GOA, where the mass change is the strongest (Supplementary Table S3). This observation coincides with our visual observation of the spatial patterns of mass losses ( Table 6). It points out that the weaker signal at the extreme East-West side is sensitive to the processing strategy and its corresponding residual noise pattern.
Considering the ice loss rates over the same time-period allows us to understand the source of the differences observed between our results and these two studies. To do so, we calculated the ratios between the trend at the center of the Saint Elias Mountains derived from the full length of the GRACE time-series and the trend computed after selecting a temporal subset corresponding to the time scale of their study (Supplementary Figure S2). All the trends have been determined by using the method described in Supplementary Figure S1. We then reported the ratios obtained over the trend rates from our focused trend maps. The results from our study over these short periods (i.e., 2003-2006 and 2003-2007) are ∼20% larger than those over the period 2002-2017 (Table 7). We conclude that the use of a shorter time-series plays a significant role in the differences observed with results from other studies, and that these differences may not only be due to the choice of GRACE processing strategy. In other words, estimations built from a short GRACE timesseries (3-4 years from 2003) have overestimated the long-term glacial mass loss by ∼20%. The inter-annual variability of the snow cover might strongly influence the TWS while considering short time-series.
We also compared our results over the Saint Elias Mountain with the mascon solution from NASA GSFC (v2.4; Luthcke et al., 2013; see Supplementary Figure S3), which considers GRACE data from 2003 to 2016. The mass loss rate obtained by using this dataset over the Saint Elias Mountains is 24.85 ± 2.1 Gt/year. Over the same area, we estimate the glacial mass loss at 22.01 ± 0.89, 21.57 ± 1.18, and 21.41 ± 1.64 Gt/year for our three FIGURE 8 | Residual maps of the focusing procedure in WTE. These maps represent the signal remaining after subtracting the best forward model (FM) to the actual GRACE TWS signal. The labels showed in (A) correspond to all of the residuals (A-I). According to these labels, the negative values of Mass loss, in cm/year, indicate that the estimates from GRACE solutions are higher than those from the forward model. The positive values show that the Mass loss from the forward model is higher. mascon resolutions. Our results are very close to those from the NASA GSFC mascon data (v02.4), and the slight difference might be attributed to the difference in time-series length (2003-2016 versus 2002-2017).

Inversion Residuals
Residual maps are produced by subtracting the FM of the high resolution mass distribution maps from the actual GRACE TWS signal. The amplitude of the residual decreases with the mascon extent. This confirms the ability to retrieve the total mass when a large number of mascons are used. However, as previously discussed, it is at the cost of a lower accuracy at the mascon scale.
We note that the inversion residuals from the GRGS and T96DDK8 solutions have similar spatial patterns. It is negligible inland (values below ± 2 cm/year) and more important near the ocean (close to −5 cm/year). The high absolute residual value near the ocean could be due to coastal glaciers covered by the land mask, implying that a part of the signal from these glaciers is not recovered through the focusing procedure.
Following Wahr et al. (2006) and Ramillien et al. (2017), the GRACE TWS data contain noise patterns from measurements and processing errors. We compare the amplitude of the noise with the residuals of the inversion procedure to verify the performance of the inversion (Figure 8). We consider that the noise level can be approximated by observing the maximum amplitude of mass trends in the Pacific Ocean at similar latitudes than our study area. We found values of ∼1.4, 1.6, and 1.5 cm/year for the CSR-MASC, GRGS, and T96DDK8 solutions, respectively. These values are close to the residuals (values between [−2: +2] cm/year), which suggests that the remaining signal left after mass concentration might not be attributed to the glacial mass loss.

CONCLUSION
This study shows a downscaling approach of GRACE TWS data which can be used to retrieve a high resolution ice mass loss. We tested three GRACE solutions and three uniform focusing delineations, and applied an inversion method which relies on fitting iteratively a spatially constrained FM. We used the three solutions at the same truncation level (resolution) in order to isolate the effect of the processing strategy applied to the GRACE data from Level 1 to 3. From our study, we can draw three key findings.
First, synthetic simulations indicate that the forward modeling approach is efficient, with a mean error of ∼2.5% when considering the total mass loss over the GOA, and below 10% when considering masses at the mascon-scale and with mascons up to 30,000 km 2 . It also shows that our inversion procedure is relatively insensitive to non-ice masses such as soil moisture.
Second, at the scale of the GOA, the resolution of the mass concentration units (referred to as mascon) and the choice of the GRACE solution strategy only account for an uncertainty of ∼2-4% in the total mass estimation. The three solutions provide approximatively the same total mass loss (∼40 Gt/year) over the era 2002-2017. Third, at the scale of a mascon of 30,000 km 2 or larger, the focusing procedure recovers well the regional mass loss signal. At this scale, the mascon contain from 1 to 9 Gt/year of ice mass loss, the highest mass loss rates being at the center of the area, over the Saint Elias Mountains. Variations between solutions reach 40% of the mean signal when considering high resolution mascons of 20,000 km 2 . Residuals are of the same order of magnitude as the noise level of GRACE TWS solutions. Thus, the glacier signals are well retrieved when using an inversion with a uniform constrained FM, and the distribution map is efficient at concentrating the mass loss observed by GRACE.
We compared our results with estimates reported in the literature and found a good agreement with studies using GRACE solutions from recent releases and time-series of more than 8 years. The first studies using GRACE data published during the 2005-2008 era generally overestimated the long-term ice mass loss at the GOA scale due to the data time-span, the interpretation strategy, and possibly the use of early release of GRACE data. Comparing our ice loss map at the highest resolution with those of Arendt et al. (2008) and the mascon solution from NASA GSFC over the Saint Elias Mountains, we obtained relatively similar ice loss rates. This area of comparison is at the center of the GOA, where our results are the most reliable, with mass losses and signal-noise ratio at the highest.
To our knowledge, this is the first study presenting a focused and spatially constrained ice loss map over the GOA. Comparison of our findings with results from other data sources (i.e., altimetry) could help further assessing the accuracy of our estimates. Constraining the focusing procedure using a heterogeneous mass loss map from other data sources would allow to cope with the heterogeneity of the losses, which is partially unaccounted for in this study. The findings of this study could help integrating an ice mass loss module into large-scale hydrological models; providing a framework to better understand the effect of climate change on the hydrological cycle of the main river basins of the GOA area.

DATA AVAILABILITY STATEMENT
We provided in the article/Supplementary Material the links and the references to access all datasets used for this study.

AUTHOR CONTRIBUTIONS
CD performed the data analysis and wrote the manuscript. PC trained the first author to process GRACE data, verified the results, and contributed to the writing of the manuscript. AR designed and led the project, supervised CD, and contributed to the writing of the manuscript. MA helped to build a comprehensive literature review on the topic.

FUNDING
The authors wish to gratefully acknowledge the financial support of the Natural Sciences and Engineering