Comment on tc-2021-181 Anonymous Referee # 1 Referee comment on " Brief communication : ICESat-2 reveals seasonal thickness change patterns of Greenland Ice Sheet outlet glaciers for the first time

The authors present a preliminary analysis of seasonal patterns in Greenland outlet glacier elevation change using 1-2 years of ICESat-2 ATL06 observations. The study focuses primarily on glaciers in the southeast, central west, and northwest. For each glacier, after detrending the data to eliminate the interannual elevation change signal, they classify the seasonal elevation change patterns. They do not find that seasonal patterns are spatially or temporally consistent, nor are they correlated with seasonality in speed. The authors suggest that future work should include critique of the summary flag used to automatically identify poor data and comparison with dense time series of terminus position change and environmental parameters. The article is easy to read overall and presents interesting new data from ICESat-2. I have several relatively minor comments below and one comment regarding reference elevations that may require considerable data reprocessing. Note that the authors may have addressed some of my comments in their response to the other reviewer comment, but I intentionally do not read comments by other reviewers so that I can provide an unbiased evaluation of the work.


1
Introduction 25 Understanding the complex nature of Earth's ice sheets is of critical importance as they have undergone dynamic changes in recent decades (Church et al., 2013;Oppenheimer et al., 2019). Greenland Ice Sheet (GrIS) marine-terminating outlet glaciers, which drive dynamic ice mass change, are projected to account for 50 ± 20% of total mass loss over the 21 st century (Choi et al., 2021). While multi-year and decadal changes of ice sheet discharge via outlet glaciers have been studied before (Mouginot et al., 2019), patterns of seasonal thickness change have not yet been studied for a representative sample of GrIS outlet glaciers. Outlet glaciers exhibit seasonal fluctuations in velocity with distinct patterns (Moon et al., 2014;Vijay et al., 2019) but the lack of seasonal thickness change measurements contributes to a lack of understanding of what processes control glacier dynamics on seasonal time scales. Seasonal thickness changes of outlet glaciers are driven by both external forcings (e.g., precipitation, evaporation, runoff, terminus melt) and internal glacier dynamics (e.g., subglacial and englacial hydrology, terminus calving) and classifying their patterns of seasonal thickness change is the first step towards a more holistic 35 understanding of the processes that control them. Here, we measure dynamic ice sheet thickness near the termini of 34 GrIS outlet glaciers at seasonal resolution for the first time using the ATL06 land ice along-track altimetry dataset from the Ice, Cloud and land Elevation Satellite-2 (ICESat-2; Markus et al, 2017;Neumann et al, 2019). Large scale observational studies such as this allow for smaller, less studied, glaciers to be observed at the same time as more well-studied glaciers and comparisons to be drawn into how these lesser-known glaciers compare with the seasonal thinning of larger glaciers, which is 40 critical for better understanding the drivers of dynamic change in a changing climate across all outlet glaciers. We use each glacier's temporal pattern of seasonal dynamic thickness changes to group glaciers into 7 distinct patterns over 2019 and 2020.
Given that we present just one to two years of data, our results are not intended to definitively characterize these glaciers but, rather, to present a method for quantifying seasonal dynamic thickness changes and to highlight the heterogeneity exhibited by these glaciers over the study time period. We discuss ways in which future work could build on our results in Section 4. 45

Data and methods
We used three data sources within this study: (1) The ATLAS/ICESat-2 L3A Land Ice Height, Version 3 (ATL06) data product, acquired by the Advanced Topographic Laser Altimeter System (ATLAS) instrument on board the ICESat-2 observatory, which provides geolocated measurements of land-ice surface heights (Smith et al., 2019); (2) Making Earth System Data Records for Use in Research Environments (MEaSUREs) glacier termini dataset of annual Greenland outlet glacier locations 50 from Synthetic Aperture Radar (SAR) mosaics and Landsat 8 OLI imagery, version 1 (Joughin et al., 2015), from which we use outlet glacier locations and identifier (ID) numbers; (3) Greenland Ice Mapping Project Digital Elevation Model (GIMP DEM; Howat et al, 2014), a digital surface elevation model of the GrIS that we used as a reference height dataset to remove along-and across-track surface slopes from the ATL06 measurements.
ATL06 provides measurements of ice sheet surface elevation at an along-track spatial resolution of 20 m, which allows 55 for ample spatial sampling of the fast-flowing, dynamic portions of GrIS outlet glaciers . We use elevation data (h_li) retrieved from all six ATLAS ground tracks to achieve the highest quantity of data available. ICESat-2 has a repeat cycle of 91 days, allowing for sufficient temporal sampling to measure seasonal changes of glaciers, although we do not receive data from every satellite pass due to cloud interference. We filter out poor quality ATL06 height data using the ATL06 quality summary flag (atl06_quality_summary), keeping only data for which the flag is set to zero. 60 The MEaSUREs glacier termini dataset contains locations for 238 glaciers across the GrIS, as well as an ID number (Joughin et al, 2015). We selected 65 glaciers from the MEaSUREs dataset due to their spatial distribution across several GrIS https://doi.org/10.5194/tc-2021-181 Preprint. Discussion started: 14 July 2021 c Author(s) 2021. CC BY 4.0 License. regions and range of average ice discharges between 68 m/yr and 8141 m/yr (Rignot and Mouginot, 2012). The 65 glaciers chosen for this study also correspond to the glaciers for which a dense record of terminus positions has been generated by the Calving Front Machine (CALFIN; Cheng, 2020). The CALFIN dataset is currently the only pan-Greenland dataset of seasonal 65 terminus positions. Although we do not use this dataset in this study, our selection of glaciers will enable comparisons of seasonal thickness change with seasonal terminus position in future studies. We define glacier seasons by three-month periods of winter (Dec-Jan-Feb), spring (Mar-Apr-May), summer (Jun-Jul-Aug), and autumn (Sep-Oct-Nov). We removed glaciers that do not contain a full year (4 seasons) of ICESat-2 data from either 2019 or 2020, reducing the number of glaciers categorized to 42 (listed in supplementary spreadsheet). 70 To collect ATL06 measurements representative of near-terminus glacier thickness change, we created a 2 km x 2 km bounding box for each glacier, centered on each glacier's location in the MEaSUREs dataset, within which we aggregated ATL06 data. We manually adjusted the MEaSUREs glacier locations slightly to ensure between one and three ICESat-2 repeat ground tracks intersect each box but we kept each bounding box within 10km of the terminus for each glacier. The 4km 2 bounding box was chosen as an arbitrary size, however it was kept to this size as a larger box may include data off the main 75 fast flowing section of the outlet glacier.
The GIMP DEM represents the mean ice sheet surface elevation between 2003 and 2009 (Howat et al, 2014). The elevation data has a 90 m spatial resolution and is used as the reference ice sheet surface elevation to account for the surface slope of the glaciers. Because the repeating passes of ICESat-2 do not exactly survey the same location on the surface of the ice sheet, ATL06 measurements from season to season are affected by both the vertical component of surface elevation change as well 80 as differences in surface elevation due to surface slope. To account for this, we sampled the GIMP DEM at each ATL06 measurement and subtracted the GIMP DEM elevation from each ATL06 surface elevation measurement. This effectively changes the datum of the ATL06 measurements to the GIMP DEM, thereby accounting for the surface slope of the ice sheet within our bounding boxes, leaving just the vertical component of surface elevation change.
We use the ATL06 data within each bounding box, a surface mass balance model, and a firn model to calculate each 85 glacier's dynamic thickness change from season to season. For each glacier, we calculate the surface elevation change (dH) between ICESat-2 observations and the GIMP DEM. We then calculated the seasonal dynamic dH as the mean of the dHs within each bounding box for each year and season, and we subtracted the surface elevation change due to changes in surface mass balance (SMB) and firn air content changes using output from the Community Firn Model (CFM; Medley et al., 2020), forced by Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) climate reanalysis (Gelaro 90 et al., 2017). Over the two-year timescale of our study, we assumed bed elevation to be constant and, thus, our surface elevation change measurements are equal to thickness change. We removed the trend from each glacier's seasonal dynamic dH, calculated over the entire duration of the available data to isolate the seasonal fluctuations from the longer-term trend. We removed 8 of the 42 glaciers with measurements of seasonal dynamic dH larger than 25m over one season, assuming that these are errors, leaving 34 glaciers for which we classified seasonal dynamic dH patterns.
To account for uncertainty in seasonal dynamic dH, we propagated error through our calculations from each data source with the assumption of random, uncorrelated error. We used the error estimates provided by ATL06 to account for error on each height data point (h_sigma). Root-mean-square differences of ±8.5m between the GIMP DEM elevations and ICESat elevations were found on ice-covered terrain (Howat et al., 2014) and we assumed this to be the uncertainty on each GIMP DEM pixel's elevation value. We assume a 20% uncertainty on the thickness change due to SMB and firn components, 100 estimated by the CFM. Assuming uncorrelated and random errors in the ATL06 and GIMP DEM surface elevation measurements, we used standard error propagation rules to calculate the error on seasonal dynamic dH, . . : where ℎ_ , represents the error on each ATL06 surface elevation measurement (h_li_sigma), 8.5 m represents the error in each GIMP DEM surface elevation, n represents the number of ATL06 elevations within the bounding box for a particular 105 season, and is the absolute value of the magnitude of surface elevation change due to changes in SMB and firn air content changes from CFM. We do not account for uncertainty in the trend that is removed from each glacier's seasonal dynamic dH because the trend is removed solely to present the thickness changes more clearly in plots. Quantifying uncertainty in the dynamic thickness change trend could be done more thoroughly in future studies, given more ICESat-2 data that will be collected over the coming years. Additionally, keeping the trend in the seasonal dynamic dH has no impact on our 110 categorization of glacier behavior for all but three glaciers, as we discuss in Section 4.
Using the time series of seasonal dynamic dH for each glacier, we manually grouped glaciers into categories based on their seasonal patterns of thickness change. Because seasonal dynamic dH had not been surveyed for a representative set of GrIS outlet glaciers, we did not prescribe categories prior to generating results. Instead, we based the categories on the timing of observed seasonal dynamic thinning and thickening for our surveyed glaciers. Each year of data is individually categorized; 115 in other words, the classification for one glacier in 2019 does not influence the classification of the same glacier in 2020.

Results
We find that, over 2019 and 2020, the 34 surveyed glaciers can be categorized into seven seasonal patterns: no statistically significant seasonal change, mid-year thinning, mid-year thickening, spring and autumn thinning with summer thickening, summer thinning with spring and autumn thickening, sharp single season thickening, and full-year thickening (Fig. 1). Glaciers 120 were classified into an additional category, "no statistical seasonal change," if seasonal dynamic dH uncertainties were larger than the amplitude of seasonal change across all seasons within a given year. Sharp single season thickening includes glaciers that undergo a lone season of significant (>3 times the change between any other seasons and >3 times the uncertainties for that glacier) thickening (either spring or summer) followed immediately by a similar sharp decline in thickness. Rink Isbrae is the only glacier that we identified with repeating sharp single season thickening across two years of results, undergoing 6-10m 125 of change during this spike (Fig. 1E). Mid-year thickening refers to glaciers exhibiting two consecutive seasons of thickening https://doi.org/10.5194/tc-2021-181 Preprint. Discussion started: 14 July 2021 c Author(s) 2021. CC BY 4.0 License.
in spring and summer before thinning in autumn. Conversely, mid-year thinning glaciers exhibit spring and summer thinning with thickening in autumn. Each glacier's detrended dynamic thickness change, alongside the seasonal trend of SMB and total dH change is plotted in the supplementary materials (Figs. S1 through S34). Although we have removed the trend to better illustrate seasonal dynamic dH for each glacier, we note that keeping the trend in the data does not alter our classifications for 130 all but three glaciers: Kangerlussuup Sermia (Fig. S9), Petermann Gletsjer (Fig. S25), and Courtauld Glacier (Fig. S28).
Without the trend removed from the dynamic dH, there is a slight thinning in Autumn 2020 for Petermann Gletsjer (Fig. S25). glaciers, their seasonal dynamic thickness changes are larger in magnitude than changes due to the 1-or 2-year trend and, thus, our classification is not sensitive to the removal of the trend.

Both Kanderlussuup Sermia and Courtauld
We find that the 34 surveyed GrIS outlet glaciers are well distributed across the seven patterns. Figure 2 shows glacier classifications for both 2019 and 2020 in the table but displays the classification from the earliest available year on the map.
With each year individually categorized, there are 47 total seasonal cycles observed between 2019 (27) and 2020 (20). Of these 140 seasonal cycles, there are 12 seasonal cycles within the summer thinning and spring and fall thickening pattern, 10 seasonal cycles with mid-year thickening, 9 seasonal cycles exhibit summer thickening with spring and fall thinning, 6 seasonal cycles experience mid-year thinning, 3 seasonal cycles with sharp single season thickening, 1 seasonal cycle exhibiting full-year thickening, and 6 seasonal cycles with no statistical seasonal change pattern. Of the 13 glaciers for which we have two years of data, we find that most glaciers exhibit seasonal thickness change patterns that differ from year to year. Two glaciers exhibit 145 repeating patterns: Rink Isbrae and Ussing Braer. However, the remaining glaciers for which ICESat-2 can so-far provide two annual cycles worth of data exhibit changing patterns between 2019 and 2020.
Although there are spatial clusters of glaciers with similar seasonal thickness change patterns, there is heterogeneity within the regions that contain multiple surveyed glaciers (Fig. 2). We use the 2019 classifications to compare glaciers per region because we have more glaciers classified in that year (27 glaciers). In the NW, 6 glaciers exhibit summer thinning with spring-150 fall thickening, 4 glaciers exhibit a mid-year thinning pattern, 2 exhibit summer thickening with spring-fall thinning, 2 exhibit mid-year thickening, and 2 exhibit no statistically significant change. In the CW, 4 glaciers exhibit summer thinning with spring-fall thickening, 3 glaciers exhibit mid-year thickening, 1 glacier exhibits sharp single season thickening, and 2 glaciers exhibit no statistically significant change. Within the SE, 3 glaciers exhibit a mid-year thickening pattern, 2 glaciers exhibit summer thickening with spring-fall thinning, and 1 glacier exhibits mid-year thinning. In the N, the single surveyed glacier, 155 Petermann Gletsjer, exhibits summer thickening with spring-fall thinning in 2019, but notably in 2020 is the only glacier to exhibit full year thickening. Small clusters of neighboring glaciers with similar patterns can be seen in the NW (glacier IDs 29,30,31,32,34,and 40), the CW (glacier IDs 5, 6, and 8), and the SE (glacier IDs 147, 148, and 153) but there is no one pattern that is representative of all glaciers within each region. We find only a weak relationship between glacier speed and seasonal dynamic thickness change patterns. The 34 surveyed 160 glaciers have a variety of speeds (Rignot and Mouginot, 2012 ; Fig. 2). The fastest glaciers in this study, with speeds above 3.5 km/yr, Kangerdlugssuaq Gletsjer (8.1 km/yr), Rink Isbrae (4.2 km/yr), and Store Gletsjer (3.71 km/yr), undergo patterns of mid-year thickening or sharp single season thickening, while medium-fast speed glaciers with speeds between 2.5 and 3.1 km/yr, such as Sermeq Kujalleq (3.1 km/yr), Kong Oscar Gletsjer (2.9 km/yr), Illullip Sermia (2.7 km/yr), and Upernavik Isstrøm S (2.5 km/yr), undergo patterns of summer or mid-year thinning. Medium-slow glaciers between 1.6 and 1.9 km/yr, 165 such as Kangiata Nunaata Sermia (1.9 km/yr), Hayes Gletsjer (1.9 km/yr), Christian IV Gletsjer (1.8 km/yr), and Kangerlussuup Sermia (1.6 km/yr), undergo mid-year thickening. Slower glaciers, with speeds below 1.6 km/yr, are more divergent in their seasonal thickness responses, for instance Cornell Gletsjer (0.5 km/yr), Sorgenfri Gletsjer (

Discussion and conclusions
Enabled by 91-day repeat measurements from ICESat-2, we have developed the first classification of GrIS outlet glacier seasonal dynamic thickness change patterns for a representative sample of glaciers from around the ice sheet. We have chosen 185 to use the ATL06 data product and to account for along-and across-track surface slopes using the GIMP DEM as a reference elevation dataset. This method allowed us to aggregate surface elevation data within customized bounding boxes, representative of each glacier's behavior. Higher-level data products, such as ATL11 and ATL15, will provide estimates of surface elevation change through time and we believe it will be worthwhile for future work to compare our results against the higher-level ICESat-2 products, both to build confidence in our results and as a check on the data products themselves. 190 Our results reveal little regional coherency in seasonal dynamic thickness change patterns, indicating that atmospheric circulation patterns are not the likely driver of differences in patterns among glaciers. While we do find small clusters of similar patterns, we do not observe similar patterns across larger ice sheet regions. If atmospheric forcing were the primary driver of seasonal dynamic thickness changes, we would expect to see coherent patterns of seasonal changes across each region.
However, we do not find this to be the case, indicating that other factors that differ from glacier to glacier within each region are causing the differences in observed patterns. This finding is consistent with seasonal glacier velocity changes, which also exhibit spatial heterogeneity (Moon et al., 2014;Vijay et al. 2019). Ocean forcing may be responsible for the differences in seasonal dynamic thickness change patterns because heat transport from the continental shelf to the termini of outlet glaciers is modulated by fjord geometry, which is heterogeneous among glaciers (Carroll et al., 2017). Each glacier's unique geometry, which has been shown to govern observed differences in terminus retreat (Catania et al., 2018) and the multi-annual upstream 200 diffusion of thinning (Felikson et al., 2021), may also be responsible for the observed heterogeneity in seasonal thickness changes.
Refining the ATL06 data quality flag (atl06_quality_summary), with the goal of accepting additional good-quality measurements that are currently flagged as poor-quality, would benefit future studies of seasonal outlet glacier change. Because ICESat-2 has a repeat cycle of 91 days, collecting good-quality data from each pass is critical to studies of the seasonal 205 thickness changes of outlet glaciers. The current set of parameters used by the ATL06 quality summary flag may remove goodquality measurements over rough topography, high surface slopes, or low-reflectivity surfaces under clouds (Smith et al., 2021). In the course of our study, we found that 12 additional glaciers, of the subset of 65 glaciers we initially selected from the MEaSUREs dataset, could be included in our results, had we ignored the quality summary flag entirely. Of course, some of the measurements that are removed by the quality summary flag are unusable and we do not advocate ignoring data quality 210 checks entirely. However, we suggest that further inspection of the parameters used for the quality summary flag to potentially reduce the strictness by which data is eliminated may prove useful and would allow additional glaciers to be considered.
As ICESat-2 continues data collection, future work should build on our two-year assessment of seasonal dynamic thickness changes by extending our record and comparing with other glacier variables and external forcings. The MEaSUREs dataset identifies 239 total outlet glaciers around the ice sheet and, by adding more outlet glaciers and extending the record 215 forward in time, future studies can examine how consistent the patterns are from year to year, identify new patterns not exhibited by the glaciers in our study, and better identify glaciers that exhibit the same or different patterns through time. With a longer and more comprehensive classification of seasonal thickness changes, future work can focus on compiling a holistic record of seasonal glacier dynamics by investigating thickness changes together with terminus position and velocity changes.
The subset of glaciers that we have selected for study are ones that have a temporally rich dataset of terminus position changes 220 from the newly developed CALFIN automated deep learning extraction method (Cheng et al., 2020) allowing our results to be directly compared with seasonal terminus positions. Finally, to advance our understanding of the processes that drive seasonal glacier behavior, future work should compare seasonal dynamic thickness changes with external forcings such as seasonal ocean temperature changes and surface meltwater runoff estimates. Our study provides the first classification of seasonal dynamic thickness changes of outlet glaciers around the GrIS to complement previous classifications of seasonal velocity 225 change (Moon et al., 2014;Vijay et al. 2019), bringing us one step closer to a holistic understanding of seasonal glacier dynamics.