Variability of glacier albedo and links to annual mass balance for the Gardens of Eden and Allah, Southern Alps, New Zealand

The Gardens of Eden and Allah (GoEA) are two of New Zealand’s largest icefields. However, their remote location and protected conservation status has limited access and complicated monitoring and research efforts. As a result, the behaviour and dynamicsTo improve our understanding of the spatial and temporal changes in mass balance of these unique glaciers are 10 poorly understood. We use annual minimum glacier-wide albedo (?̅?yr ) derived fromicefields, observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors aboard the Terra and Aqua satellites are used to monitor annual minimum glacier-wide albedo (?̅?yr ) over the period 2000-2018 to assess linkages between the spatial and temporal variability of ?̅?yr min and glacier mass balance.. The ?̅?yr min for 12 individual glaciers, identified using a new objective glacier masking approach, ranges between 0.42 and 0.70, and can occur as early as mid-January and as late as the end of April. There is 15 evidence that The evolution of the timing of ?̅?yr min has shifted?̅?yr min indicates a shift to later in the summer over the 19-year period on 10 of the 12 glaciers. However, there is only a weak relationship between the delay in timing and the magnitude of ?̅?yr , which implies that albedo is not necessarily lower if it is delayed. The largest negative departures in ?̅?yr min (lower than average albedo) are consistent with high snow line altitudes (SLA) relative to the long-term average as determined from the End-of-Summer Snowline (EOSS) survey, which has monitoredbeen the benchmark for monitoring glaciers in the Southern 20 Alps for over 40 years. While the record of ?̅?yr min for Vertebrae Col 25, an index glacier of the EOSS survey and one of the GoEA glaciers, explains less than half of the variability observed in the corresponding EOSS SLA (R = 0.43, p = 0.003), the relationship is stronger when compared to other GoEA glaciers. Angel Glacier has the strongest relationship with EOSS observations at Vertebrae Col 25, accounting for 69% of its variance (p = 1.8×10< 0.001). A key advantage in using the ?̅?yr min approach is that it enables the variability in the response of individual glaciers to be explored, revealing that topographic setting 25 plays an importanta key role in addition to the regional climate signal. The observed variability of individual glacier response at the scale of the GoEA contrasts with the high consistency of responses found by the EOSS record across the Southern Alps, and challenges the suggestion that New Zealand glaciers behave as a “unified climatic unit”. MODIS imagery acquired from the Terra and Aqua platforms also provide insights about the spatial and temporal variability of clouds. The frequency of clouds in pixels west of the Main Divide is as high as 90% during summer months, and reach a minimum of 35% in some 30 locations in winter. These complex cloud interactions deserve further attention as they are likely a contributing factor in controlling the spatial and temporal variability of glacier response observed in the GoEA. The observed variability of individual glacier response on the GoEA demonstrate that glaciers in the Southern Alps do not behave as a single climatic unit.


The manuscript provides very little in the way of new science. No comparison to actual or modelled mass balance was conducted. A proxy for mass balance was compared to another proxy for mass balance -this can be worked out from first principles.
We disagree that this manuscript does not provide new science. There is a fundamental question that needs to be answered here. It must be recognised that there is no published information on the behaviour of the glaciers in the Gardens of Eden and Allah, despite them being two of the largest ice fields in New Zealand. Thus, one needs to ask whether we want to provide the glaciological community with new and detailed insights about these glaciers, using a state-of-the-art remote sensing approach, or whether to accept the view of the reviewer that no new science is contained in the manuscript. If we accept the latter, is the glaciological community really better off not having any new insights or information about these glaciers?
We would argue that even documentation of the areal extent of these glaciers in Fig. 1, which is not currently available in any detail, is valuable new knowledge. However, the absolutely crucial contribution of this research is that it is a stepping stone to providing a new pathway to monitor glaciers remotely, and may be far more powerful than the current standard, which is the End of Summer Snowline (EOSS) approach. Importantly, it reveals the limitation of the timing of the EOSS observations and how variability in response at an individual glacier scale can be resolved in much greater detail and more widely using the albedo method. Both approaches are carefully compared to the only available direct observations of mass balance by Sirguey et al. (2016), which does not appear to have been read by the reviewer.
It is not clear to us how first principles can be meaningfully applied to understand glacier behaviour if no observations exist from a glacier of interest, whether those are sourced from actual, modelled or remote sensing data. Why is one approach better than the other? If one applies this logic, the question that begs asking is why do we observe the behaviour of glaciers at all if their changes can always be resolved from "first principles"?
Lastly, we are of the view that comparison is a cornerstone of the scientific method. As such, we believe that comparison between our method and that used by the EOSS programme is valuable to the scientific community and we disagree with the reviewer that there is no merit in doing so.

Most of the studies original shortcomings were not addressed in the revision. No climate data was added to support the idea that the glacier basins are indeed climatically different. No systematic analysis of MODIS error between Aqua and Terra was
conducted.
As noted above, we have made significant improvements to the manuscript by following in detail the feedback from the reviewers. This is articulated clearly in our responses, which include supporting figures, new analysis and a new discussion section (see Section 5.2 Implications of variability of albedo on mass balance). We would like to point out that the reviewers did not ask for climate data to be added as part of the revisions, but importantly we have refined our discussion to ensure readers understand how the hypothesis of the Southern Alps behaving as a "climate unit" has evolved from the EOSS programme, and whether this is legitimate or not. It is this point that is central to our manuscript, not the climatic assessment of individual glacier basins.

No systematic analysis of MODIS error between Aqua and Terra was conducted.
We do provide careful details about the differences between Terra and Aqua (L182-190), and it is not clear why any further systematic analysis of error between these two MODIS sensors would be necessary given that this aspect of the research is not the focus of the manuscript. It is stated in the manuscript that the use of Aqua was carefully explored to see whether it could be used to "fill in" the Terra time series but as we demonstrate, it does not add much additional benefit. This is an important result as few have explored this. However, the approach did allow us to document the spatial and temporal variability of cloud cover, which is novel and has value -this is why we report on these initial findings. Our question to the reviewer would be how would the proposed systematic analysis of MODIS error improve our current methodologywhat new insights would such an analysis bring or more importantly, how has our current approach failed readers?
Please consider ways to revise your manuscript so that the issues above are at least discussed in detail in the manuscript, if it is not possible to incorporate them in your analysis. For example,in (1) what was the reason for not implementing a simple mass balance model (using a degree-day model for melting, or a simple energy balance model), with reanalysis data as forcing, to get an estimate of glacier mass balance? What concrete evidence do you present for the links between the albedo changes and glacier mass balance (instead of changes in summer snow line or similar proxies)? In (2), what has prevented you from performing a systematic error analysis? Could more sensitivity tests (e.g. assumptions of random errors, or using  We have made substantial revisions to the manuscript and we would like our comments to the reviewers, as well as the manuscript to be assessed before making any additional changes. The former was clearly not done by the reviewer, which we believe is a necessary step before proceeding. As you are aware, there are considerable uncertainties that are introduced when modelling glacier response to climate using a simple mass balance model. These glaciers are located in a very data sparse regionthere are no in situ meteorological observations, and it is highly questionable whether reanalysis data could meaningfully resolve individual glacier response in the Gardens of Eden and Allah. That is what makes using remote sensing so attractive as it provides key insights that would otherwise not be attainable. The concrete evidence of the links between the albedo method, EOSS and mass balance is extensively documented in the literature, including glaciers in New Zealand (e.g. Sirguey et al., 2016;Rabatel et al., 2017). Importantly, the method that we use has been very carefully tested on the only glacier in the Southern Alps that has direct mass balance observations (Sirguey et al., 2016), so as we argue in the current manuscript, the next logical step is to use it on glaciers that we don't know anything about to provide new insights. We are offering the glaciological community an alternative to the EOSS approach, which clearly has some limitations. We believe there is merit in exploring a new approach to asess glacier behavioiur remotely in the Southern Alps, and to carefully document for the first time the behaviour of two of the largest ice fields in New Zealand. As clearly stated, these glaciers are located in a Wilderness Area with very limited access. We are of the view there is no better way to do this than what we offer in this manuscriptcareful remote sensing observations that provide a platform to monitor glacier behaviour more widely across the Southern Alps (L691-694). It will make much more sense to use reanalysis data and a simple mass balance model once we shift to a large-scale approach of detecting glacier response more broadly across the Southern Alps.
Having all 12 glaciers clearly defined on a map would benefit result interpretation.
Response: We modified Figure 1 to display glacier outlines and names as well as clusters.

(right). If the timing of the ablation minimum is getting later, does this mean that the minimum albedo is also decreasing due to a longer alation season? Or Not?
Response: We thank the reviewer for this constructive comment. We added a panel showing the evolution of the minimum glacier-wide minimum albedo alongside the evolution of its timing. As we revised this figure, we decided to also display yearly values to help visualise the variability better, while maintaining the 3-year rolling average as a dotted line. Furthermore, we added a scatter plot that shows the relationship between and its timing and assessed its correlation to inform the hypothesis of the reviewer. We find a weak correlation when all clusters are considered together (R 2 =0.26, p<0.001) and an even weaker correlation when all glaciers are considered separately (R 2 =0.11, p<0.001). However, within an individual cluster, the value of does not exhibit a significant correlation with its timing (cluster 1: R 2 =0.2, p=0.058; cluster 2: R 2 =0.08, p=0.248; cluster 3: R 2 =0.03, p=0.472). We have added some further insights in Section 4.4 to report on this new result, updated references to this figure in the text where appropriate, and added a glaciological interpretation and significance of these findings in a new discussion section (5.2 Implications of variability of albedo on mass balance).
As noted above, the lack of a detailed location map hinders spatial thinking, as does the organisation of Figure 6. While clearly the authors have opted to organise both graphs in Figure 6 by the scale of the x-axis values in doing so the reader is left with no clues as to how these glaciers are actually related in space, again is there a spatial influence on the data presented? It is very difficult to compare results of the left and right graph for an individual glacier as the order of the y-axis (by giving priority to the x-axis value) are different for each graph. While it is appreciated that a 'progressional' x-axis approach might be the 'neatest' presentation, something is lost in regards to the actually physical process or characteristics that might be driving the patterns being presented. For example can something more be said about W/E or N/S trends? If one colour-codes the class sizes some patterns appear, for example class 1 glaciers tend to have lower albedo and a later minimum timing, whereas class 2 (n=2, should potentially include Eve) have higher albedo and earlier minimum timing.
Response: Thanks for the constructive comment. We found the suggestion of colour coding glaciers by class in Figure 6 very helpful and well-thought out. Section 4.4 in the original manuscript did characterize the contrast in magnitude and timing of minimum albedo between glaciers of class 1 and 2 that now becomes more obvious with the colour coding. The reference to Fig. 6 has also been added to the point made in Section 4.4. As noted above, we have added a new section in the discussion that adds to the interpretation of our results by assessing how topography modulates the response and possible fate of the GoEA glaciers.
While it is appreciated that this paper is 'methods' focused there is missed opportunity to engage more fully with some of the glaciological findings.
Response: Our paper makes use of the albedo method to produce a comprehensive set of new glaciological observations in an area where very little data exist. The principles of this method are described for completeness, yet this section remains limited to echoing the body of literature where it was developed and validated. In particular, Sirguey et al. (2016) provides a very useful 'methods focused paper' that is in part focused on describing, calibrating and validating the albedo method.
In contrast, this study is very much an application of the method rather than a paper about the method. It provides a large amount of new glaciological observations unavailable to date, and characterizes glacier behaviour of largely undocumented icefields. It is true that this new application of the albedo method provides an opportunity to incrementally refine it, for example by testing the use of glacier outlines to facilitate masking and the opportunity to use MODIS Aqua data. Nevertheless, results associated with these "methodological" steps only involve three out of the now eleven sections of the results and discussion. All other results present new insights about the behaviour of glaciers in the study area. This paper also provides a discussion of the implications of the new findings regarding the variability of the glacier surface compared to results from the EOSS methodology traditionally used in New Zealand. We have also explored the validity of the "unified climatic unit" theory about the response of New Zealand glaciers. In response to comments from both reviewers, we have added a new discussion section (5.2 Implications of variability of albedo on mass balance) that we hope builds knowledge about the glaciology of this remote but important region in the Southern Alps.
In particular, the finding that the timing of the minimum albedo is occurring later in the summer, which could signal a later onset of the first winter snowfall (i.e. lengthening of the ablation season). This result also makes one wonder if there is any trend (across all the glaciers measured) of a decreasing minimum albedo over time, for if the ablation season is becoming more protracted then the snow surface would likely become more discoloured. Or alternatively, is the minimum albedo the same and the trend is associated with a later start to the ablation season (i.e. more spring snow delaying the onset of melting). It would be great to see a little more discussion of these important mass balance feedbacks.

Response:
We thank the reviewer for this input. We found the suggestion very useful as a prompt for us to explore this issue more deeply as our data challenges this hypothesis. As demonstrated above, the comparison between the magnitude of and its timing shows that a delayed minimum albedo does not necessarily lead to a lower albedo. We have added elements of our interpretation of this result in the new discussion section.
Should the authors wish to include a reference for the Rolleston Glacier mass balance programme, they could cite Purdie, H., Rack, W., Anderson, B., Kerr, T., Chinn, T.J., Owens, I. and Linton, M. 2015: The impact of extreme summer melt on net accumulation of an avalanche fed glacier, as determined by ground-penetrating radar. Geografiska Annaler, Series A: Physical Geography 97, 779-791.

Response:
We thank the reviewer for this suggestion and have added the citation to support the statement made in the introduction about the Rolleston Glacier mass balance programme.
Response: Accepted, we expressed the actual p value unless p<0.001, in which case we only used the statement of inequality.
The results section of the Abstract need to be rewritten to present all of the results.

Response:
We have refined and updated the abstract.

Furthermore, a statement of why you have done this analysis is required.
Response: This comment echoes the need for clearer objectives that was highlighted in the first general comment. We have modified the abstract and other parts of the manuscript to ensure the aim and objectives of this research are more explicit.

Response:
We are referring to the theory of the Southern Alps glaciers behaving as a "unified climatic unit" as proposed by Chinn et al. (2012). This is one of the main outcomes of the EOSS monitoring programme inferred from the high degree of intra-correlation found in the variability of snowline elevation between any index glacier and the average of the group. In the introduction where the term "unified climatic unit" is directly quoted from Chinn et al. (2012), we have clarified that this is referring to a "consistent response to climate variability inferred from the EOSS record". We have also reworded the abstract to clarify the origins and meaning of this term.

Introduction:
P2,Line 17: The minimum albedo method does not infer mass balance. It scales to ELA and AAR, which in turn scale to mass balance. Without measurements of mass balance to quantify the relationship to minimum albedo, ELA, or AAR, mass balance should not really enter into the discussion.
Response: We agree that the word "infer" was not appropriate in this sentence. We replaced it by "estimate". We however respectfully disagree with the reviewer that mass balance should not enter into the discussion. This paragraph is introducing the general concepts of remote sensing approaches for the very purpose of capturing glacier mass balance, whether its relative variability, or its absolute value when calibration data can be provided. The statement made remains factually correct in view of the supporting literature provided in the paragraph.

Response:
We have rephrased this sentence as "The annual minimum glacier-wide albedo ( ) retrieved from MODIS imagery has been found to scale to annual glacier mass balance…" P2,Line 20-23: There are more citations for this methodsee below.
Response: We thank the reviewer for these suggestions. Among those, we find that only Zhang et al. (2018) is directly relevant and added it to the Introduction section.
P2,Line 25: why is the glacier contribution "globally significant" and define (i.e., X m sea level increase).

Response:
We must stress to the reviewer that the meaning and definition of significant is not unique and solely reserved to characterise a quantitative value subject to statistical significance testing. In context, significant means "Sufficiently great or important to be worthy of attention; noteworthy" (Oxford Dictionary). As such, New Zealand is invariably represented in global studies as its own glacier region. In view of this we maintain that our statement that New Zealand Glaciers "are regarded as globally significant" is factually correct and has adequate context to suggest the meaning of significant in this sentence as "noteworthy". The alternative would be that New Zealand glaciers would be otherwise ignored in such global studies, which is clearly not the case. We respectfully disagree with the reviewer that using significant in this context necessarily calls for a specific quantified statement.

Response:
The point the reviewer is raising is not clear, nor what she/he means by "minimum ELA" in this context. We understand ELA stands for Equilibrium Line Altitude, which has meaning in surface mass balance. Braithwaite & Raper (2009) define "minimum ELA values corresponding to balance years with highly (...) positive mass balances". In view of this, we believe the reviewer's use of "minimum ELA" seems misplaced and makes a connection between ELA and glacier outlines that is irrelevant in our context. Nonetheless, we assume the reviewer is trying to say that glacier outlines should be derived from imagery timed precisely at the maximum ablation. If so, we dispute this strict requirement because glacier outlines are not a direct consequence of the surface mass balance for that year due to glacier response time. In practice, one aims for late summer imagery mostly because it corresponds to the melting of as much transient snow as possible to facilitate interpretation of glacier boundaries. At the same time, tradeoffs must be made given the availability of imagery, and mapping glacier outlines remains often heavily reliant on photo interpretation and experience. We indicated that we used the NDSI from the S2 image on 14 March as a base to derive glacier outlines. We then refined them manually (e.g., NDSI cannot map debris-covered parts and transient snow needs to be removed) using the additional data mentioned in the text. We believe this is a very common and valid approach, yet it is not within the scope of this manuscript to discuss this in more detail. For the reviewer's information, the Sentinel-2 image from 30 April 2016 shown in Figure 9 in the original manuscript would be closer to the end of ablation and the latest before winter snowfall that year. On the one hand it exhibits less residual snow, which may help. On the other hand, it can be seen in the comparison below, how topographic shading affects the image on 30 April compared to 14 March. Mapping solely from this image is not making the task easier nor more accurate. The NDSI is sensitive to shadow with high values in the shade greatly complicating the snow/ice segmentation. We maintain that the mid-March image is a better compromise despite the remaining snow patches, to derive suitable glacier outlines for this study. The consistency between snow patches with the late April image and across a few years with other Sentinel images (e.g., those shown in Fig.  9) and Pleiades 2017 images help discriminate them from glaciers at the refinement stage. Braithwaite, R. & Raper, S. (2009). Estimating equilibrium-line altitude (ELA) from glacier inventory data. Annals of Glaciology 50 (53), 127--132. (Doi: 10.3189/172756410790595930.)  NDSI (B3-B11)/(B3+B11) stretched between 0 (black) and 1 (white) NDSI (B3-B11)/(B3+B11) stretched between 0 (black) and 1 (white) P5,Line 12: What is "field knowledge"? Actually the whole section should be rewritten. What are the NDSI wavelengths used?

Response:
We have modified this section and specified the wavelength used for the NDSI. As explained, we took advantage of a two-day field trip to take close-up aerial photos of the glacier. This provided us with valuable direct evidence of this remote area that supplemented our interpretation of satellite imagery used to refine some glacier outlines. To avoid any confusion, we removed "field knowledge" from the revised manuscript.
There is almost two months difference between image capture dates. How is this reconciled?

Response:
We have reworked this section to ensure there is no confusion for readers, and we hope this addresses the reviewer's concerns about our approach to mapping glacier outlines in this study. To respond to the specific question, we acknowledge that we did not specifically mention in the manuscript that we used the 29 April 2016 image to refine our mapping, which we understand the reviewer is referring rather than the 14 March 2016 image. It should also be noted that we used 0.5m imagery from Pleiades in March 2017, as well as photos from our field work in 2018. Thus, we compiled images across two years to help interpret glacier outlines from the March 2016 Sentinel-2 imagery. For example, the Pleiades imagery and high resolution surface supported much better interpretation of debris-covered glaciers. It is also important to keep in mind that the glacier boundaries used for this study are derived to extract 250-m resolution pixels from the MODIS albedo map and that glaciers are not expected to retreat hundreds of meters per year in this region (it would need calving in terminal lakes to achieve such rates). The mapping also compounds uncertainties in the order of tens of meters due to the geolocation of Sentinel-2 imagery and the variable interpretation at pixel level. In this regard, we don't believe the two-year period between these sources of information and our approach to mapping are any cause for concern. Further evidence of the benefit of exploiting all of these data are shown below, which provides a sequence of the terminus of Lambert Response: T (resp. A) superscripts indeed stands for MODIS Terra (resp. Aqua). We believed this was self-explanatory. We are aware of the MOD/MYD convention for naming MODIS products, and we did refer to this usage in P6L9 of the original manuscript. In context, when clearly referring to either sensor rather than the products from them, we do not believe that using the acronyms of MOD or MYD are appropriate. We therefore maintain our use of the superscript to refer to either one or the other sensor, and clarified this convention in Section 2.2 to avoid confusion.
P5, Section 2.2: You should read: https://nsidc.org/data/modis/terra_aqua_differences There is some speculation that the Terra MODIS sensor is degrading. This degradation is largely correct for in the collection 6 data. There are several citations that should be considered in relation to this issue. These citations have been listed at the end of the review.

Response:
We thank the reviewer for this reading suggestion. This also echoes the general comment made earlier about "The error analysis of MODIS albedo and the potential for sensor degradation (post 2016 on the Terra bus) need to be addressed." We have provided further information about MODIS calibration issues in Section 4.8.1. However, we raise two important points about this issue in response to the reviewer's comments.
First, the main point being made on webpage given by the reviewer is about the differences in processing MODIS snow products between Terra and Aqua sensors due to the band 6 failed detectors of Aqua. There is no mention about the degradation of MODIS detectors. We are very aware of the differences in processing MOD10/MTD10 C6 products. However, these are not relevant for our study as we are not relying on MOD10 data. Instead, we complete all processing and albedo retrieval directly from L1B data, as is clearly explained in the manuscript, and in even more detail in the cited literature about our processing technique. As far as we know, the Quantitative Image Restoration (QIR) technique to restore Aqua band 6 is only used in production of MYD10 products but not in the production of the MOD03/MOD02 C6 products we are using. Please also note that we have acknowledged the failed detector of band 6 and demonstrate that it does not significantly affect our retrieval, which we believe is a valuable outcome of our study. This was expected as band 6 in the SWIR is an absorption band for snow and ice, and therefore contributes only marginally to the albedo signal.
Second, we are aware of the sensor calibration degradation in MODIS C5 data as demonstrated by Lyapustin et al. (2014). The latter concludes very clearly that "The new C6 calibration approach removes major calibrations trends in MODIS Level 1B data.". From the citations suggested by the reviewer, Casey et al. (2017) and Polashenski et al. (2015) both point out and document calibration issues with C5 data and conclude that the use of C6 data is preferable, if not necessary. Sayer et al. (2015) also confirmed the better performance of the C6 calibration. The reviewer points out that the degradation is largely corrected for in the C6 data. As we use the C6 data, which are the highest quality data available and do not have any severe calibration issues, we are of the view that further justification for using them is not necessary.
In this context, we shall stress Figure 8 and 10 to the reviewer which show albedo time series. The reviewer can assess a lack of visible trends in winter albedo which we believe should lift concerns about the existence (speculated or not) of any major calibration issue with C6 affecting our study and results. The alternative hypothesis would be that winter albedo has a trend that is precisely matched by a calibration issue. We do not believe that this hypothesis is likely. P5, Line 29: orbits are usually described as ascending and descending.
Response: Agreed, we modified the revised manuscript accordingly.
P6, Line 2: What does "exceptionally cloudy" mean? Provide a description with statistics.

Response:
We replaced "exceptionally" by "extraordinary" as described by Wardle (1986), and we reworked the sentence to make the origin of this statement unambiguous and invite the reader to consult the citation for more information. We also agreed with the reviewer and added specific statistics from Wardle (1986) that are specific to two sites surrounding our region of study. It should be noted that Cropp River, west from the GoEA, exhibits the most frequent cloud cover of all sites examined by Wardle (1986).
Why were only four days used for comparison between Terra and Aqua? Provide justification.

Response:
We clarified that in addition to clear sky acquisition on the same day, such comparison needed similar sensorzenith to minimise the effect of contrasting panoramic distortion on pixel footprint. As our assessment of cloud cover below demonstrates, there are not many opportunities to get clear-sky conditions for both Terra and Aqua. This opportunity is further reduced when adding the limitation of sensor zenith. We also stated the number of samples provided by the four images amounts to 2196 pixel pairs, which we believe is a sufficient sample size to complete this comparison.
The MODIS methods section does not indicate which albedo product was used, or how albedo was produced at 250 m. A much more clear description of data processing is required. O.k., I see the following section on albedo processing. Perhaps a sentence here to indicate that MODImLab will be detailed (and why) later in the paper.

Response:
We added a sentence in the previous paragraph that the albedo retrieval is described in the methods section.

Section 2.3 requires a much better description of what exactly was done.
Response: This comment suggests to us the reviewer is assuming we completed the mapping of the EOSS snowlines and generated the EOSS SLA records for Vertebrea Col 25 as part of this study. This is not the case. The EOSS programme is a national programme run by the National Institute of Water and Atmospheric Research or NIWA. An annual report with SLA values for each index glacier is normally published sometime after the observations are taken. Appendices about the method used to map SLA on each year are also provided on request. They indicate what method was used to derive SLA for any particular year. We have clarified that the EOSS record was obtained rather than generated by this study, summarised the methodology used by the programme to derive SLA and supported it with relevant citations in Section 2.3. We also substantially reworked the section to stress the outcome of the EOSS programme that led to the "climatic unit" theory that our study is testing.
P6, Line 21: Sentinel was used to "approximate" the timing and elevation of SLA. A great deal of effort should be spent on what approximate means in this instance.

Response:
We considerably reworked this paragraph to clarify how we used additional Sentinel-2 images to observe the evolution of the glacier surface and assess its consistency with relative variations of EOSS and albedo in those years. We also revised Section 4.6 to stress the consistency of these observations with the measured evolution of the minimum surface albedo, and related those variations to documented heatwaves that affected New Zealand glaciers.
P8, Line 30: why was 50% used? How was this value determined? Should 100% be used if you want to remove debris-covered glaciers?
Response: We agree that in theory 100% snow/ice cover should be used. This however assumes that the spectral unmixing is perfect. This is not the case and relying only on pixels truly achieving 100% snow and ice membership from unmixing excludes too many to retrieve a useful albedo signal. We use the 50% threshold because it is the one generally used for binary discrimination of snow/ice pixels (see Hall et al. 2002;Hall & Riggs, 2007;Sirguey et al., 2009;Masson et al., 2018), while being conservative enough to preserve a sufficient amount of pixels. It is true that it may include mixed pixels into the signal and it is the reason why we tested this approach with the conservative masking technique (M2). We did report and discuss in Section 4.1 that the albedo signal with the objective masking (with 50% threshold) was slightly darker on average by 1.2%, likely due to some mixing in the debris-covered pixels. The differences between both approaches however demonstrated sufficient agreement to accept this choice. We added justification for the choice of this threshold and added a citation in support. Masson Page 9, Section 3.3: The requirement of using only cloud free glaciers is very restrictive and substantially reduces the amount of data used in the analysis. Provide statistics on the amount and duration of cloud cover. When there is cloud cover will likely be as important as where there is cloud cover.
Response: Accepted, we included that the cloud-free criteria excluded about 66% of images.
Pixel is a picture element found on a display screen. Grid cell is a better usage.

Response:
We respectfully disagree with the reviewer on such a strict use of "pixel". The use of "pixel" to denote image data is widespread in the literature and remote sensing textbooks. For example, Richards (2013) defines "We talk about the recorded imagery as image data, since it is the primary data source from which we extract usable information. One of the important characteristics of the image data acquired by sensors on aircraft or spacecraft platforms is that it is readily available in digital format. Spatially it is composed of discrete picture elements, or pixels. Radiometrically-that is in brightness-it is quantised into discrete levels." Richards, J. (2013). Remote Sensing Digital Image Analysis. An Introduction. Springer-Vertag.
Page 10, L1: Artefact should be replaced by error.

Section 4.3 should be presented in the Introduction section. Typical values for the different glacier states should be provided.
This is an example of the systematic problem with this paper, in that there really aren't that much in terms of reportable results.
Response: As noted earlier we have clarified the atmospheric controls on mass balance in the Introduction, and have removed and changed some of the text in this section to provide clarity. As suggested by the reviewer, we also moved two sections of discussion (5.2.1 & 2) to the results, which now contains 9 sections of results (4.1-4.8a,b). Thus, we are of the view there is no shortage of reportable results.
Page 11, L5: The classification of glaciers into four groups is accomplished by geographical variables. Without providing climate data, it is a mystery to me how the authors determine the glaciers are not behaving as a single climate unit.

Response:
We classified glaciers into three groups, rather than four as suggested. We have responded to an earlier comment to clarify that the "unified climatic unit" refers to the theory proposed by Chinn et al. (2012) in reference to the consistent response of glaciers inferred from the EOSS record. This is based on the high degree of intra-correlation between EOSS SLA across the Southern Alps, not on the study of climate data. Since the snowline and albedo methods aim to capture a similar glacier response, our study can apply a similar approach to test this hypothesis by assessing the degree of intra-correlation between minimum albedo across the GoEA. The outcome of this test and our conclusion that it challenges the "unified climatic unit" inferred by EOSS, is unrelated to the clustering of glaciers. The classification provides additional insights that help characterise what may control the consistency or contrasts in glacier responses across this area. It does not preclude the existence of local contrasts in climate that ultimately govern glacier mass balance. Our results do raise the question about the extent of local climate variability, which itself may be complicated by the complex terrain, but can hardly be answered without reliable climate data that are not available in this area.
Page 11, L28: Either the declines are monotonic or they are not.

Response:
We respectfully disagree with the reviewer, "monotonic decreasing" has mathematical meaning that "decreasing" alone has not. A signal can be decreasing over a period albeit punctuated by increases. This is the definition of an increasing or decreasing function. A monotonic decreasing function however shall not exhibit any increase within a period of decrease. We have replaced decline by decrease to better fit the mathematical definition.

Page 13, L15 (and elsewhere) Rˆ2 values should be presented with significance values.
Response: We added the significance value.

Page 13, L29-30. What is a maximum summer snow line? When exactly was this observed within the summer period?
Response: We added the word "altitude" and the years when those observations were made. The exact date can be found in the figure caption.
P14,L23: Using only 12 glaciers, which are argued to occupy different climactic regions, can size and elevation be reasonably discounted as predictors of the min albedo and SLA?
Response: Indeed, this is what our results show. We have added the correlation analyses and significance of these as part of an earlier comment, which should help clarify and avoid any further confusion.

Response:
We reworked this sentence to be more factual with reference to supporting citations.

P15,L23: I wouldn't have characterised the Rˆ2 values found in the Results section as strong to moderate.
Response: To reiterate here, the R2 value between EOSS and for Vertebrae Col is 0.43, while for Angel glacier it is 0.69. Overall, Class 1 glaciers yield an average R2=0.51, while Class 3 glaciers yield R2=0.36. While we appreciate that there is no unique approach to qualifying the strength of the relationship, Akoglu (2018) compiled the following with all approaches compatible with a qualification of those R2 as indicative of moderate to strong relationship.
Adapted from Table 1  Section 5.2 is results and should not be presented in the discussion section.
Response: Accepted. We moved these sections to results. In relation to characterising cloud cover frequency, we added a statement in Section 2.2 that explicitly indicates that we re-used the classification of clouds from the albedo processing chain of the complete inventory of MODIS imagery for this purpose. The new discussion section (5.2) maintains key points relevant to the control of clouds on the spatial variability of glacier response.

P18,L7: "likely having as significant influence" means you don't know if it is significant or not.
Response: We replaced "significant" with "substantial".

Section 5.3. This section should be expanded. Is the correlations between minimum albedo to SLA related to the amount of cloud cover?
Response: In response to an earlier comment we added a discussion about the limitation of MODIS spatial resolution with respect to the size of glaciers. The section referred to has been moved to results and we have created a new section to address in part the potential role clouds have played on albedo and SLA.

Response:
We did indicate in Section 2.2 that this comparison used 4 pairs of clear-sky imagery for the purpose of testing the consistency of the albedo retrieval between Aqua and Terra, it is therefore unrelated to cloud cover. We must stress to the reviewer that this point is discussed in detail in Section 4.8.1 (former 5.2.1). This is related to the contrasting shadowing configuration between the morning and afternoon image acquisitions and the difficulty of retrieving an accurate observation of albedo in the shade. and protected conservation status has limited access and complicated monitoring and research efforts. As a result, the behaviour and dynamicsTo improve our understanding of the spatial and temporal changes in mass balance of these unique glaciers are 10 poorly understood. We use annual minimum glacier-wide albedo ( ̅ yr min ) derived fromicefields, observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors aboard the Terra and Aqua satellites are used to monitor annual minimum glacier-wide albedo ( ̅ yr min ) over the period 2000-2018 to assess linkages between the spatial and temporal variability of ̅ yr min and glacier mass balance.. The ̅ yr min for 12 individual glaciers, identified using a new objective glacier masking approach, ranges between 0.42 and 0.70, and can occur as early as mid-January and as late as the end of April. There is 15 evidence that The evolution of the timing of ̅ yr min has shifted ̅ yr min indicates a shift to later in the summer over the 19-year period on 10 of the 12 glaciers. However, there is only a weak relationship between the delay in timing and the magnitude of ̅ yr min , which implies that albedo is not necessarily lower if it is delayed. The largest negative departures in ̅ yr min (lower than average albedo) are consistent with high snow line altitudes (SLA) relative to the long-term average as determined from the End-of-Summer Snowline (EOSS) survey, which has monitoredbeen the benchmark for monitoring glaciers in the Southern 20 Alps for over 40 years. While the record of ̅ yr min for Vertebrae Col 25, an index glacier of the EOSS survey and one of the GoEA glaciers, explains less than half of the variability observed in the corresponding EOSS SLA (R 2 = 0.43, p = 0.003), the relationship is stronger when compared to other GoEA glaciers. Angel Glacier has the strongest relationship with EOSS observations at Vertebrae Col 25, accounting for 69% of its variance (p = 1.8×10 -5 < 0.001). A key advantage in using the ̅ yr min approach is that it enables the variability in the response of individual glaciers to be explored, revealing that topographic setting 25 plays an importanta key role in addition to the regional climate signal. The observed variability of individual glacier response at the scale of the GoEA contrasts with the high consistency of responses found by the EOSS record across the Southern Alps, and challenges the suggestion that New Zealand glaciers behave as a "unified climatic unit". MODIS imagery acquired from the Terra and Aqua platforms also provide insights about the spatial and temporal variability of clouds. The frequency of clouds in pixels west of the Main Divide is as high as 90% during summer months, and reach a minimum of 35% in some 30 locations in winter. These complex cloud interactions deserve further attention as they are likely a contributing factor in controlling the spatial and temporal variability of glacier response observed in the GoEA. The observed variability of individual glacier response on the GoEA demonstrate that glaciers in the Southern Alps do not behave as a single climatic unit.

Introduction
The retreat of mountain glaciers during periods of the 20 th and 21 st centuries has been widespread and unprecedented (Zemp et al., 2015(Zemp et al., , 2019. One approach to document this change is to monitor glaciological mass balance, which enables the response of glaciers to climate forcing to be assessed (Hock et al., 2019;Zemp et al., 2019). However, global records of mass balance 40 are sparse, as traditional in situ glaciological measurements are resource intensive (Cogley et al., 2011), with only 260 glacierwide results available from an estimated 200,000 glaciers worldwide (Pfeffer et al., 2014;Zemp et al., 2015). In addition, the ongoing cost of maintaining the repetitive measurements means that few long-term records exist, with only 30 of the world's glaciers having an uninterrupted series dating back to 1976 (Zemp et al., 2009). Mass balance records are therefore biased towards the more accessible glaciers in the Northern Hemisphere, while remote and relatively inaccessible glaciers in South 45 America and New Zealand have remained largely out of reach. In this context, there is a need to monitor the state of more mountain glaciers to capture their response to changes in regional and large-scale climate.
Recent developments in the capabilities of satellite and airborne sensors, and improved processing techniques, have provided several alternatives to the in situ glaciological approach, to derive individual or regional mass balance signals for out of reachremote glaciers (Rabatel et al., 2017). For example, mapping the end-of-summer snowline altitude (SLA) provides an 50 estimate of the glacier equilibrium-line altitude (ELA) and accumulation area ratio (AAR). These glacier properties have been used as proxies for annual mass balance and are often less resource intensive to monitor (Chinn et al, 2012;Rabatel et al., 2016;Salinger et al., 2019b). Alternatively, the 'albedo method' uses glacier surface albedo, for example retrieved from satellite imagery captured by the Moderate Resolution Imaging Spectroradiometer (MODIS), to inferestimate glacier mass balance. and/or monitor its variations (Rabatel et al., 2017). Surface albedo plays a governing role in the mass balance of 55 glaciers due to its control on absorbed short-wave radiation (Klok and Oerlemans, 2004;Oerlemans et al., 2009;Pope et al., 2016). Importantly, relationships have been found between theThe annual minimum glacier-wide albedo ( ̅ yr min ) retrieved from MODIS imagery andhas been found to scale to annual glacier mass balance in the French Alps (Dumont et al., 2012;Davaze et al., 2018), the Arctic (Greuell et al., 2007), the Himalayas (Brun et al., 2015;Zhang et al., 2018) and New Zealand's Southern Alps on Brewster and Park Pass glaciers (Sirguey et al., 2016;Rabatel et al., 2017). 60 The mountain glaciers of the Southern Alps of New Zealand comprise the largest ice mass in the Southern Hemisphere outside of Antarctica and South America, and are regarded as globally significant (Chinn et al., 2012;Zemp et al., 20162015). The Southern Alps are unique as they are surrounded by vast areas of ocean and are subject to both subtropical and polar air masses that are embedded in a prevailing westerly airflow. This unique setting contributes to the humid, maritime climate that influences the glaciers in the Southern Alps, which is reflected most starkly by the exceptionally high precipitation rates that 65 exceedexceeding 10 m a --1 in many alpine regions (Fitzharris et al., 1999). However, itMeteorological experiments on glaciers in the Southern Alps have shown the main source for melt energy is net radiation, driven primarily by net shortwave radiation (e.g. Gillett and Cullen, 2011;Cullen and Conway, 2015), but fluctuations in mass balance over time are sensitive to changes in air temperature that plays the most important roleand precipitation. Given the importance of air temperature in controlling mass balance variability in the Southern Alps both melt and the phase of precipitation, the advance and retreat of glaciers in 70 the wet climate of the Southern Alps are particularly sensitive to changes in air temperature (Oerlemans, 1997;, though the). For example, the advance of some fast-responding glaciers in the Southern Alps between 1983 and 2008, during three of the warmest decades in the instrumental record, have been attributed to regional cooling controlled by changes in large-scale atmospheric circulation in the Southern Hemisphere . The influence of synoptic scale circulation on air mass variability and changes in cloud properties remains importanthave also been shown to influence 75 mass balance in the Southern Alps Cullen et al., 2019). Changes in large-scale circulation have the potential to enhance or suppress the effects of global warming in the Southern Alps, which has been most recently demonstrated by the advance of some fast-responding glaciers between 1983 and 2008 .
ComprehensiveThe first comprehensive glaciological mass balance study in New Zealand's Southern Alps was completed on Ivory Glacier between 1969and 1975(Anderton and Chinn, 19731978). Today, comprehensive glaciological mass balance 80 programmes exist on Brewster Glacier (e.g. Cullen et al., 2017) and Rolleston Glacier, (Purdie et al., 2015), extending back to 2004 and 2010, respectively. Important glaciological measurements were also made on Ivory The size of Brewster Glacier between 1969(ca. 2 km 2 ) and 1975 (Andertonits comprehensive in-situ record enabled the albedo method with MODIS to be assessed and Chinn, 1973;1978).calibrated to estimate mass balance (Sirguey et al., 2016). With only 0.11 km 2 in 2012, Rolleston Glacier is however too small to be captured by MODIS imagery. In addition, the End-of-Summer Snowline (EOSS) 85 survey, which has been the benchmark for monitoring glaciers in the Southern Alps for over 40 years, has used aerial photography to map and retrieve the altitude of the annual snowline at the end of the ablation season (SLA on) of 50 'index' glaciers across the Southern Alps since 1977 (Chinn et al., 20122012). The SLA is used as a surrogate of the Equilibrium Line Altitude (ELA) whose annual variations inform about changes in the annual balance of the glaciers (Rabatel et al., 2017).
Despite is longevity, the EOSS survey continues to face two main challenges: (1) the temporal resolution is limited to a single 90 observation per year in late-summer; the success of which is dependent on transient conditions such as snowfall prior to a flight, and (2) limited resources constrain the total number of index glaciers that can be mapped. Furthermore, the focus on relatively small glaciers in the Southern Alps has meant that only three index glaciers exceed 2 km 2 . These challenges, combined with the difficulty of establishing a network of glaciological programmes, have meant thatresulted in the vast majority of glaciers in New Zealand remainremaining out of reach of current observation methods. 95 The EOSS programme indicates that the yearly response of the transient snowline for each of the index glaciers is consistent with the average of the group, with strong to very strong intra-correlations between SLA departures across index glaciers.
ThisThis consistent response to climate variability inferred from the EOSS record has led to the suggestion that glaciers in the Southern Alps behave as a "singlean "unified climatic unit" (Chinn et al., 2012, p. 115). Arguably, theThe linear models corresponding to this hypothesis have led to observations from single glaciers, such as Tasman or Brewster Glacier, being used 100 to estimate fluctuations in ice volume for the entire Southern Alps (Salinger et al., 2019a;2019b). To further explore the validity of this approach, more long-term continuous signals of glacier mass balance are required across a broader topographical range of glaciers, particularly those larger than 2 km 2 . Given the available approaches to monitor glaciers and the lack of in situ data, the albedo method is a desirable alternative to objectively resolve mass balance variability and trends on large, un-monitored New Zealand glaciers at near-daily temporal resolution. 105 The Gardens of Eden and Allah (GoEA), located on the Main Divide of the Southern Alps, northeast of Aoraki/Mt Cook National Park, comprise two of the largest icefields in New Zealand. (Fig. 1). They form an interesting target to test the climatic unit hypothesis, as the icefields consist of a network of smaller, interconnected glaciers, placed in a critical climatic zone straddling the Main Divide. Their position on the Main Divide, combined with their protected conservation status has greatly constrained the possible methods for data collection and extraction. As a result, these glaciers are yet to be the target of focused 110 research, and their behaviour over the past decades is largely undocumented and poorly understood.
To address this gap in knowledge, the aim of this paper is to characterise the spatial and temporal variability of glacier-wide surface albedo on the GoEA, and to investigate the linkages between this variability and glacier mass balance. By resolving and documenting spatial patterns in the evolution of glacier surface albedo across two of New Zealand's largest icefields, the single climatic unit hypothesis for glaciers in the Southern Alps can be tested, and the governing physical processes controlling 115 mass balance explored in more detail.To achieve this aim, this research has the following objectives: (1) to determine the timing and magnitude of glacier-wide surface albedo on 12 different outlet glaciers on the GoEA; (2) to compare the spatial and temporal changes in surface albedo to variability in SLA as determined from the EOSS programme; and (3) to characterise variability in cloud cover using MODIS imagery acquired by the Terra and Aqua platforms. Remote sensing with the albedo method provides an opportunity to examine the significant, yet remote and protected GoEA icefields, and provides an 120 opportunity to further assess the skill of the albedo method to remotely monitor the physical processes governing glacier mass balance of glaciers in the Southern Alps.
The next section describes the location and significance of the icefields in the New Zealand context. The methods section provides a description of the remote sensing datasets and processing techniques used to derive glacier-wide surface albedo.
Two developments to the albedo method are proposed, namely (1) a new objective way to define glacier boundaries for albedo 125 analysis, and (2) the processing and assessment of albedo from MODIS imagery acquired by the Aqua platform to increase the temporal resolution of albedo monitoring. The results section presents the records and analysis of glacier-wide surface albedo for the GoEA, including a comparison to results from the EOSS programme. Finally, the discussion compares the results to the EOSS survey, examines the contribution of MODIS Aqua to resolve albedo and its usefulness in combination with Terra to characterise cloud cover variability, and describes the limitations and future opportunities of the albedo method. lifting to generate a significantlarge precipitation gradient from west to east across the Main Divide (Henderson and Thompson, 1999). The estimated mean annual rainfall in the area over the period 1972-2016 ranges from 6,000 to 8,500 mm yr -1 (Macara, 145 2017). The mean annual temperature over the glaciated region for the period 1980-2010 is estimated to range from -1.0 to 2.5°C based on an interpolation of mean normal temperature from distant surrounding weather stations and a lapse rate determined from the data of -5.7°C km -1 . Given the important rolecontribution of snow and ice into New Zealand's water resources, changes in the volume of these icefields have the potential to impact the hydrology of these rivers under future climate change (Chinn, 2001). 150 Despite their significance, research on the GoEA is limited by the remote location of the glaciers, and direct access has been further complicated in recent years by their inclusion in the 466 km 2 Adams Wilderness Area (est. 11 February 2014).
Helicopter landings are now largely prohibited, effectively ruling out a number of in situ monitoring methods. As a result, the most significantonly comprehensive data obtained from the area to date are EOSS SLA records for the two Vertebrae glaciers at the western margin of the GoEA (Willsman et al., 2018), which are described in more detail below. (Section 2.3). The 155 Garden of Eden was also briefly included in an attempt to detect changes in the ELA of New Zealand glaciers from 15 mresolution ASTER satellite imagery . The absence of data contrasting withand the largenotable extent and volume of ice locked in a potentially sensitive climatic situationthe GoEA make the GoEAregion a compelling target for research among the more than 2500 glaciers in the Southern Alps.

Glacier outlines and surface topography 160
As part of this research, the boundaries of the icefield outlets were redefined, updating the existing New Zealand Glacier Inventory outlines derived during the 1970s (Chinn, 1991), and refining the 2017 Randolph Glacier Inventory (RGI) version 6.0 outlines (Pfeffer et al., 2014).

MODIS products
This study relies primarily on a record of glacier surface albedo retrieved from a time-series of MODIS images. MODIS is one of the key sensors aboard the Terra (EOS/AM-1) and Aqua (EOS/PM-1) satellite platforms., hereafter referred to as MODIS T 185 and MODIS A , respectively. Terra was launched on 18 December 1999 as the flagship of NASA's Earth Observing System (EOS), allowing MODIS T to capture near daily, moderate-resolution images of Earth's surface since 25 February 2000. Aqua was launched in May 2002, creating a second daily opportunity to capture MODIS A images. Terra passes north to south overdescending node crosses the equator at approximately 10:30 am, while Aqua completes the reverse journey with acquisitionascending node is crossing at approximately 1:30 pm. Historically, MODIS T imagery has been preferred for the 190 albedo method due to the longer record of imagery, as well as failed detectors with MODIS A Band 6.
Following Sirguey et al. (2009) and Dumont et al. (2012), the albedo retrieval was completed with the MODImLab software as described in Section 3.1. As yet, the use of MODISA imagery to supplement the MODIST record for snow and ice albedo retrieval has not been explored. Wardle (1986) demonstrated the extraordinary cloudy sky conditions that dominate New Zealand's Southern Alps tend to be dominated by exceptionally cloudy sky conditions (Wardle, 1986),, greatly reducing the 195 data available to space-borne remote sensors. Wardle (1986) specifically reported that the western flanks of the study area near Cropp River exhibited the highest frequency of days with some clouds (92%), as well as days with heavy clouds (72%). The frequency of cloudy days tends to decrease towards the east, with some or heavy clouds occurring on 72% and 38% of days, respectively. The consideration of MODISA data in this study provides an opportunity to compare surface albedo derived from each sensor, and to inform on key climate variables such as cloud frequency and variability over the Southern Alps. If the 200 albedo derived from MODIS A is suitable, despite the degraded Band 6, it can be used to increase the temporal resolution of the albedo method to two observations per day, or to fill gaps in the MODIST record.
WeFollowing Sirguey et al. (2009) and Dumont et al. (2012), we use MODIS Level-1B Collection 6 (C6) swath data products that contain calibrated and geolocated top-of-atmosphere (TOA) radiance counts for all 36 spectral bands. The products include bands 1 and 2 supplied at 250 m nadir resolution (MxD02QKM files, x=O/Y for Terra/Aqua, respectively). Bands 3 to 7 are 205 provided at 500 m spatial resolution (MxD02HKM files), along with an aggregated 500 m version of bands 1 and 2. All 36 spectral bands are provided at 1000 m resolution (MxD021KM files), including the aggregated bands 1-7 at 1000 m. The accompanying MxD03 file contains the geolocation information and relative geometry parameters such as solar zenith (θs), solar azimuth (φs), sensor zenith (θv) and sensor azimuth (φv), required to take into account sun-target-sensor geometry in the albedo retrieval (Dumont et al., 2012;Sirguey et al., 2016).

EOSS survey programme
The national EOSS programme is run by National Institute of Water and Atmospheric Research (NIWA). Aerial photos are captured on the first suitable weather window following 1 March. The mapping of SLA for each of the index glaciers implements three methods summarised hereafter (see Willsman, 2018). (1) When the snowline is clearly visible, it is sketched from the oblique photograph onto 1:50,000 topographic contour maps of each glacier, or digitized readily from rectified photos 225 to map the accumulation or ablation zones. The snowline elevation is derived from the ablation area and the hypsometric curve of the glacier. (2) The SLA is often determined by applying an "interpolation method" whereby all EOSS photographs of a glacier are arranged into a sequence of increasing area of snow cover (descending order of SLAs); the SLA value is interpolated form that of the adjacent years. (3) When the snowline is obscured, it is interpolated from an interpretation of the degree of snow cover surrounding the glacier compared to previous years on record. 230 For each year i, EOSS SLAi records delivered by the EOSS programme for the two Vertebrae glaciers in the GoEA provide an independent dataset to make a comparison with the retrieved albedo signal. Vertebrae Col 25 is a 0.7 km 2 mountain glacier facing southwest at the western tip of the Garden of Eden (Fig. 1). Vertebrae Col 12 is a smaller cirque glacier (0.21 km 2 ) located directly adjacent, to the north. The SLA was first recorded for the Vertebrae glaciers in 1978, despite the EOSS survey beginning in 1977, and has been recorded every year since, with the exception of 1979 (no flight), 1984 (no visit), 1987 (no 235 visit), 1990 (no flight) and 1991 (no flight). Of the 35 observations in the 40-year period since 1978, the SLA was digitised eleven times (Fig. 2). In digitised years, the snowline is sketched onto the oblique aerial photographs and then digitised to derive its elevation from 1:50,000 topographic contour maps. In the remaining years, SLAi is estimated by comparing the visible extent of the summer snowline to photographs from previous years.method 1) and interpolated (method 2 and 3) all other times (Fig. 2). 240 The SLA for each year (i.e. SLAi)As part of the EOSS programme methodology (Chinn et al, 2012;Willsman et al., 2018), the SLAi for each glacier is compared against the long-term or 'steady-state' SLA, termed 'ELA0', which as determined by the EOSS survey is 1864 m a.s.l. and 1840 m a.s.l. for the Vertebrae Col 12 and 25 glaciers, respectively. Vertebrae Col 12 is not considered in this study, as it is not of sufficient size for albedo retrieval with MODIS, however the annual SLA departures from ELA0 for the two glaciers are strongly correlated (R 2 = 0.95). The SLAi departures for Vertebrae Col 25 show eight years 245 with particularly high seasonal snowlines since 1999, indicative of a negative glacier mass balance (in particular 1999, 2008, 2011, 2012 and 2016), preceded by a majority of years with snowlines located near or below ELA0, indicative of a positive mass balance.
Both SLA departure records from Vertebrae Col 12 and 25 correlate strongly to the SLA departures averaged over the remaining index glaciers photographed in a particular year (EOSSAlps; R 2 = 0.86 and 0.92, respectively)., see Willsman et al., 250 2018). These very strong correlations between individual glacier responses to EOSSAlps suggest that the Vertebrae Col glaciers, like other index glaciers in the EOSS programme, respond uniformly to climate variability (or act as a climatic unit).. This high degree of intra-correlation in the EOSS record has led to the theory that glaciers in the Southern Alps behave as a "unified climatic unit" (Chinn et al., 2012). However, we hypothesise that changes in the climate system are likely to result in greater variability in glacier response in the Southern Alps, which is arguablymay not beingbe resolved and/or detected by the EOSS 255 programme. The ability to derive and discriminate a mass balance signal from a large number of glaciers in the GoEA using the albedo method, which has a high temporal resolution, provides us with a uniquenew opportunity to further characterize the spatial variability of glacier behaviour in the GoEA, as well as test the consistency in glacier response predicted by the EOSS programme.

Sentinel-2 data 260
As EOSS SLAi records only exist on two of the smaller glaciers in the GoEA, high-resolution Sentinel-2 data are also used to support the MODIS analysis by mappingobserving the evolution of the summer snowline across the icefields. Launched on 23 June 2015, Sentinel-2 imagery are available for the 2016, 2017 and 2018 summers, and provide a secondan independent means to assess the consistency of the albedo products and EOSS survey results. For the three summer periods (1 January -30 April), sequences of cloud-free 10 m-resolution Sentinel-2A and 2B Level-1C images are used to identify the approximate elevation 265 and timing of the SLAi overobserve the evolution of the GoEA. surface until seasonal snow ends the ablation season. There is an expectation that the discolouration of the glacier surface and elevation and timing of the SLAisnowline in the Sentinel-2 images depicting maximum ablation will closely reflect the results from compare and be consistent with relative changes determined by the EOSS survey and the summer minimum glacier-wide albedo measured across the glaciers due to the relationship between the transient snowline and glacier-wide albedo (Sirguey et al., 2016). Qualitative results of these 270 observations are presented and discussed for Lambert Glacier, with the elevation of the snowline determined from LINZ 20 m topographic contours.

Retrieving snow and ice albedo
Terra and Aqua MODIS C6 Level-1B products were processed using the MODImLab toolbox . 275 MODImLab has been widely employed to perform a series of image processing techniques on MODIS time-series data, yielding snow and ice albedo products at 250 m-resolution (Dumont et al., 2011(Dumont et al., , 2012Brun et al., 2015;Sirguey et al., 2016;Rabatel et al., 2017;Davaze et al., 2018). The operations and processes of MODImLab are described by Sirguey et al. (2009) and Dumont et al. (2012), and as such, are only covered here in brief.
First, image fusion combines the higher resolution MxD02QKM bands 1 and 2 (250 m) with the lower resolution MxD02HKM 280 bands 3 to 7 (500 m) to produce seven spectral bands at 250 m spatial resolution (Sirguey et al., 2008). As a result, MODIS bands importantused for mapping snow cover (band 4 at 555 nm and 6 at 1640 nm) have a higher spatial resolution than other common MODIS products (e.g. MxD10 -500 m). The atmospheric and topographic correction module (ATOPCOR) then corrects images containing TOA radiance counts into values of ground surface reflectance (Sirguey, 2009). Topographic corrections are calculated from the NZSoSDEM along with the MxD03 product to account for the relative illumination and 285 viewing geometry. The DEM is also used to pre-process the sky-view and terrain factor to account for diffuse and terrainreflected irradiance over rugged terrain, and the effects of both self and cast shadow.
The MODImLab algorithm has the ability to resolve the sub-pixel snow cover fraction in mountainous terrain ). Spectral unmixing estimates the relative contributions of specific land cover types to the radiometry of each individual pixel (Masson et al., 2018). Values of spectral albedo (narrowband albedo) for each pixel measured from five MODIS bands 290 are compared against a Look-Up Table (LUT) generated with DISORT (Stamnes et al., 1988). The LUT contains an array of spectral albedo values simulated for snow and ice surfaces with varying snow grain size, impurity type and content, and incident zenith angle (Dumont et al., 2011). The best match of spectral albedo can then be integrated to yield the Blue-Sky (BS) and White-Sky (WS) broadband albedo for a pixel. BS albedo is the value of broadband albedo corresponding to the specific ground irradiance. However, diffuse, or WS albedo is preferred in this research, as it allows the surface albedo to be 295 studied with reduced sensitivity to seasonal changes in illumination conditions (Dumont et al., 2012). Given the small reflectance of snow and ice targets at 1600 nm and marginal contribution to broadband albedo, the performance of albedo retrieval from MODIAMODIS Aqua data is not expected to be compromised despite the degraded Aqua band 6.

Glacier masks
Having processed the MODIS granules with MODImLab, we then create 250-m raster glacier masks using the glacier outlines 300 from Section 2.1.2 under which albedo is averaged. Glacier masks are defined from glacier outlines to avoid debris-covered areas or mixed land-cover pixels (e.g. containing a combination of snow/ice and rock) so as to preserve the integrity of the snow and ice albedo retrieval. In previous studies, this has been achieved subjectively (e.g. Dumont et al., 2012;Brun et al., 2015;Davaze et al., 2018), but takes time and introduces variability between users. As a result, we consider an alternative approach to masking that can be more objectively deployed across a large number of glaciers. 305 From all pixels within the glacier outlines, we use the sub-pixel snow classification produced by spectral unmixing in MODImLab to exclude those with snow-covered-area less than 50% snow., as this threshold is generally used to segment snow from snow-free pixels ). This excludes most non-glaciated surfaces from the overall glacier outlines (e.g.,. debris-covered) and mitigates the effect of mixed-pixels in the glacier-wide albedo. We assess the effectiveness of this masking approach (Mask 1) against a conservative glacier mask created manually (Mask 2) by comparing the glacier-wide 310 albedo ̅( ) from each mask and outlet glacier over the full sequence of MODIS T imagery. The mean surface albedo series retrieved from each mask, ̅( ) 1 and ̅( ) 2 , are compared using the linear coefficient of determination (R 2 ), root mean square difference (RMSD) and mean difference (MD).

Filtering the MODIS-albedo record
Despite the near-daily capture of MODIS images over the GoEA, cloud cover in the Southern Alps significantlygreatly reduces 315 the quantity of available data. Following Sirguey et al. (2016), only MODIS images with no cloud present in any pixels within the glacier mask are retained in the analysis. Cloudy pixels were determined using MODImLab's cloud detection algorithm, based on the original MODIS MOD35 Cloud Product described by Ackerman et al. (1998). Using these products, Brun et al. (2015) and Davaze et al. (2018) demonstrated the successful application of a cloud threshold, whereby images are retained so long as a certain proportion of the glacier surface is cloud-free in the image (>20% and >30% clear surface, respectively). 320 While there is merit in this approach for larger glaciers (e.g. Chhota Shigri Glacier: 15.7 km 2 and Mera Glacier: 5.1 km 2 ), the outlet glaciers of the GoEA identified in this analysis are much smaller (average ca. 2.8 km 2 ). On small glaciers, there is a higher probability that large parts of the accumulation or ablation areas will be obscured by cloudclouds, and therefore the surface albedo across the glacier may be misrepresented. Relying on cloud-free conditions over the glacier excluded ca. 66% of the images, resulting in an average of one clear-sky image every 3 days. The classification of clouds in the complete 19-325 year long inventory of near-daily MODIS images was then leveraged to characterise the spatial variability in cloud frequency of occurrence around the GoEA.
An increase in sensor viewing zenith angle (θv) distorts MODIS pixels due to the panoramic effect. The ground sampling distance of MODIS pixels increases from 2-to 5-fold from θv = 45º to 66º, respectively (Wolfe et al., 1998). Albedo retrieval of mountain glaciers with MODIS is most accurate with θv < 30º (Sirguey et al., 2016). However, Brun et al. (2015) and 330 Sirguey et al. (2016) accepted θv < 40º and < 45º, respectively, to increase the number of images available for retrieval. Davaze et al. (2018) confirmed that, while albedo retrieval was most accurate with view angle θv < 30º, it performed well until θv < 45º. In this study, and given the various configurations and size of glaciers in GoEA, we found that this threshold remained suitable to optimise the number of images and quality of the albedo retrieval.

Calculating glacier-wide albedo 335
Glacier-wide albedo ̅( ) was then calculated for each MODIS image (at time t) as an average of all pixel values within each outlet glacier mask. Following Dumont et al. (2011), ̅( ) was calculated using the WS albedo product, under the anisotropic reflectance model. To reduce the noise of the ̅( ) signal and smooth the time-series, a three-period rolling average was applied (Sirguey et al., 2016), and was preferred to a fixed-day average, which is heavily influenced by large gaps in the series, driven by persistent cloud cover. The ̅ yr min was then extracted as the annual minimum glacier-wide value of smoothed ̅( ), from the 340 period corresponding to the end of the ablation season (1 January -30 April). Due to complex topography and cast shadows in the GoEA, visual confirmation revealed that when ̅ yr min is reached outside of this period, it is always due to incorrect albedo retrieval (underestimated) due to inaccurate topographic correction in the shade, as identified by Davaze et al. (2018).

Characterising topographic shading
Shading at the glacier surfaces complicates the topographic correction and has the potential to create artefactserrors in the 345 albedo retrieval algorithm. This typically occurs at low sun zenith angles, outside of the late-summer period when ̅ yr min is reached. Nonetheless, it is important to understand how a changechanges in the distribution of shortwave radiation being received at the glacier surface affects processes such as the evolution of surface albedo. Net shortwave radiation is one of the key components of the glacier surface energy balance (SEB). The intensity of solar radiation received at the surface is a function of the time of year (season) and global position (latitude). At a local scale, this relationship is complicated by 350 topography, slope and aspect. Olson and Rupper (2018) show shading from topography (both cast-and self-shading) is an importanta key factor contributing to the SEB. Therefore, in addition to quantifying glacier slope and aspect (Section 2.1.2), the ATOPCOR module from MODImLab is used to calculate fractional shading at 250 m resolution across the glacier surface at the time of image capture. We use the maximum proportion of surface shading (occurring during the winter solstice, when the solar zenith angle is at its maximum) as a simple metric to compare shading between glaciers. A high maximum percentage 355 of topographic shading indicates a topographically confined glacier, while a low maximum percentage indicates an open, unconfined surface that receives radiation year round (e.g. Fig. 3).

Assessing the mask performance
Between February 2000 and April 2018, the difference between derived glacier-wide albedo ̅( ) 1 and ̅( ) 2 derived is 360 small (RMSD = 0.037; Table 1), with Mask 1 (objective masking approach) typically yielding smaller albedo compared to Mask 2 (MD = -0.012). Importantly, theThe linear agreement between the masks is highest during the critical summer period (1 October-31 March), when the ̅ yr min ̅ yr min is expected to be reached. The increased difference during winter months is likely a result of debris-covered pixels (not considered by Mask 2) that become snow-covered beyond 50%, thus included in Mask 1 and raising ̅( ). The relatively small difference between the two methods, particularly during summer, is confirmation of 365 the suitability of the objective glacier masking approach, and is therefore used to produce the following results over the GoEA.

Glacier characteristics
Despite the two icefields being an interconnected ice mass, the glaciers of the GoEA exhibit large differences in their hypsometry ( Table 2). The majority of the ice resides slightly above 1900 m a.s.l., close to the average elevation of the plateaus.
As with many other glaciers in the Southern Alps, the average surface gradient is steep, with a number of glaciers approaching 370 20º. In addition, glacier size is relatively small, ranging between 0.44 km 2 and 4.44 km 2 . Lambert Glacier is an exception (9.44 km 2 ), comprising over one quarter of the total surface area (33.89 km 2 ). The glaciers also occupy a range of aspects with variable topographic shading. As expected, glaciers with mean north-facing aspectaspects display lower values of topographic shading than south-facing glaciers (e.g. Angel Glacier, east-northeast, 5.4%; Colin Campbell, south, 82%).
It is anticipated that the large topographic differences between outlet glaciers may drive significanta large variability in the 375 temporal evolution of glacier surface albedo. To further characterise the contrasting topography of these glaciers, we perform a K-mean cluster analysis based on mean aspect, mean slope and maximum topographic shading. Three clusters are identified derived from the mapped outlines of each glacier. The glacier characteristics (Table 2) when viewed in scatter plots revealed that the 12 outlet glaciers grouped in three identifiable clusters, with the contrasting hypsometry indicated in Fig. 4.
ClassCluster membership for each glacier is provided in Table 2. Glaciers in ClassCluster 2 are characterised by southerly 380 aspects, steep slopes and incised topography (indicated by topographic shading exceeding 80%). These glaciers contrast to the north-and east-facing, unconfined glaciers in ClassCluster 1. ClassCluster 3 shares similar topographical attributes to ClassCluster 1, although the aspect is primarily west-facing. Importantly, the glaciers in Clusters 1 and 3 account for 81% of the total surface area of the GoEA , while the two south-facing glaciers in Cluster 2 occupy the remaining 19%.

Annual evolution of MODIS-derived glacier-wide albedo 385
All 12 glaciers in the GoEA exhibit a marked seasonal evolution of MODIS-derived glacier-wide albedo ̅( ) when averaged within-cluster (Fig. 5). As air temperatures rise and incident solar radiation increases following5a), characterised by a decrease in ̅( ) at the end of the austral winter (1 September),as controlled by the onset of melt. The lowering of ̅( ) decreases. This is primarily due to the melting of snow from the previous winter, exposing the underlying glacier ice and firn. To compound this process, snow and ice that remain through the summer undergoes metamorphism. This process drives an increase in grain 390 size, liquid water content and impurity concentration, which all act to decrease surface albedo. As a result, glacier-wide albedo continues to fall through the summer, until fresh snow begins to fall (accumulate in late summer and early-March). autumn.
The duration of exposure of discoloured glacier ice and firn across the glacier during summer is a critical part of the physical processes controlling surface ablation. Although each of the three clusters follows a typical pattern of seasonal evolution, the timing and magnitude of key events, such as the ̅ yr min and the rate of change of ̅( ), appears to varyvaries between glaciers. 395 On some glaciers, the seasonality of the albedo signal is complicated by a significantsubstantial decrease in ̅( ) centred around the winter solstice (21/22 June). This trend is not unique to the GoEA, and has been identified to be an artefact of widespread surface shading when the sun is at its lowest that compromises the radiometric correction by MODImLab, resulting in incorrect surface albedo (Rabatel et al., 2017;Davaze et al., 2018). This issue affects both glaciers in classcluster 2 (observable in Fig. 55a), and some in classcluster 3 (namely, Eve Icefall and Perth Glacier). These glaciers all exhibit a 400 south/southwest aspect and topographic shading greater than 60% (Table 2). Importantly, asAs this pitfall develops in winter, it has little or no impact on the summer albedo signal and the identification of ̅ yr min . We therefore disregard much of the winter signal from these glaciers, and focus on the summer ablation period.
During the ablation period, all three classes share consistent patterns over the 19-year albedo record that deviate from a monotonic declinedecrease. Short-lived increases in ̅( ) occur in mid-October and December, suggesting synoptic-scale 405 atmospheric processes result in late snowfall events, which temporarily change the albedo signal over the glaciers. A similar pattern of late spring or early summer snowfall events were also observed on Brewster Glacier between November and December (Sirguey et al., 2016), which have a strong and coherent weather-system scale signature (Cullen et al., 2019). The occurrence of these snow events likely play an importanta considerable role in the evolution of albedo through the summer period and in turn, glacier mass balance. Large events have the potential to deposit enough fresh snow at the surface to provide 410 prolonged protection of the otherwise exposed glacier ice during the height of summer. can occur as early as mid-January and as late as the end of April. In addition to a large range in timing between glaciers, the variability in the arrival of ̅ yr min on certain glaciers is also large (e.g. σ = 29.5 days for Barlow Glacier).

Annual minimum glacier-wide albedo (̅
The Spearman's rank coefficient is used to determine the topographic controls on the median value and timing of ̅ yr min across the outlet glaciers (Table 3). To avoid the circular nature of aspect data (i.e. where 0º and 360º are equal), values of mean aspect (in degrees) are converted into Cartesian coordinates and correlated independently. A strong and significant association 420 is found between the magnitude of ̅ yr min and the proportion of topographic shading (r = 0.741, p = 0.008). Similarly, strong and significant associations are also found with the timing of ̅ yr min (north/south component of glacier aspect: r = 0.805, p < 0.002; proportion of topographic shading: r = -0.782, p = 0.003). These correlation resultscorrelations, as well as Fig. 6, suggest that glaciers with northerly aspects and low topographic shading (cluster 1) exhibit a lower and delayed ̅ yr min . Conversely, glaciers with southerly aspects and heavy topographic shading (cluster 2) typically have a higher ̅ yr min , which occurs earlier 425 in the year. Interestingly, Fig. 5 shows that the timing of ̅ yr min over the period of the study appears to have evolved towards late summer for class 1 and 3, while this timing is unchanged for the two glaciers in class 2.
Interestingly, Fig. 5d shows that despite a large interannual variability, the timing of ̅ yr min appears to have evolved over the period 2010-2018 towards late summer for clusters 1 and 3, while this timing remained relatively unchanged for the two glaciers in cluster 2. However, Figure 5b reveals no significant temporal trend for ̅ yr min within each cluster (the Pearson 430 correlation coefficient between ̅ yr min and the year yielded p = 0.20, 0.70 and 0.52 for clusters 1, 2 and 3, respectively). A delayed ̅ yr min could indicate a longer ablation duration, and consequently correspond to a lower ̅ yr min . However, we only find a weak relationship between ̅ yr min and its timing when aggregating the three clusters (R 2 =0.26, p < 0.001, see Fig. 5c). This correlation proves to be weaker when considering all glaciers individually, with the Julian Day of the minimum albedo only explaining 11% (p < 0.001) of the variability of ̅ yr min . Furthermore, this relationship loses significance within clusters (cluster 435 1: R 2 =0.20, p = 0.058; cluster 2: R 2 =0.08, p = 0.248; cluster 3: R 2 =0.03, p = 0.472).

Gardens of Eden and Allah (2000 -2018)
We use the same methodology as the EOSS survey to characterise and compare the variability of ̅ yr min across the 12 outlet glaciers. The ̅ yr min across the 12 outlet glaciers are averaged to establish ̅ 0 min (56.5%), from which yearly departures are calculated. The combined 19-year record of ̅ yr min variability averaged across the 12 outlet glaciers of the GoEA is shown in 440 Fig. 7. The use of ̅ 0 min follows the theory outlined by the EOSS survey, where the long-term average elevation of the SLA is assumed to approximate a glacier in its equilibrium state. The main limitation of ̅ 0 min is that the record only consists of 19 years of data, as opposed to almost 40-years for the EOSS survey. In addition, there is no compelling proof that glaciers in the Southern Alps have been in equilibrium during either period. Therefore, it is possible that the value of ̅ 0 min proposed here actually represents the GoEA in a state that is not in equilibrium. In this instance, ̅ 0 min is arguably limited to describing relative 445 change over the observation period, as opposed to providing a measure for quantitative mass balance in absolute terms.
Negative departures from ̅ 0 min correspond to years where glacier-wide minimum albedo is lower than average, and in turn a more negative mass balance than average is expected. The opposite is expected for positive departures. Therefore, based on the relationship between ̅ yr min and annual mass balance, it can be inferred that these shifts broadly correspond to relative shifts in glacier mass balance. In general, the glaciers of the GoEA appear to have undergone four changes to ̅ yr min since 2000. The behaviour of the glacier-wide surface albedo anomaly on the GoEA over this 19-year period couples reasonably well with EOSSAlps, with about half of the variability explained (R 2 = 0.55, p<0.001). Notably, the largest negative departures in ̅ yr min are in general consistent with high snowlines observed across the index glaciers of the Southern Alps in 2000, 2008, 2011, 460 2016 and 2018. Despite the relatively limited length of the series, the agreement between the two methods in years where the snowline departure is positive (R 2 = 0.68; n = 9, p=0.006) is much stronger than in years with a negative departure (R 2 = 0.22; n = 9, p=0.201). Salinger et al. (2019a) also highlight the overwhelmingly positive departures of the transient snowline during the 2018 summer, with EOSSAlps estimated to be 386 m above the long-term average from linear regression with the EOSS of Tasman Glacier alone. Since this estimate is derived with a very different methodology than the traditional EOSS survey, it is 465 not included in the current analysis. Nonetheless, such a large positive departure is supported by widespread reports of exceptional ablation of glaciers in the Southern Alps. The effect of this extreme summer on the GoEA glaciers is now also supported by results in this study using the MODIS driven albedo method.

Assessment of ̅ yr min on Lambert Glacier (2016-2018)
The full 19-year time-series of ̅( ) retrieved from Lambert Glacier exemplifies the annual variability in ̅ yr min (Fig. 8). The 470 raw MODIS-derived values of ̅( ) are indicated in red, along with the three-period rolling average used to illustrate the seasonality of the albedo signal. In lieu of the glaciological mass balance data needed to relate the ̅ yr min to quantitative mass balance, we use Sentinel-2 data to independently support the MODIS record. The images in Fig. 9 show the end of the ablation season and maximum altitude reached by the summer snowline observed on Lambert Glacier observed in cloud-free Sentinel-2 images during the 2016, 2017, and 2018 summer periodperiods. Figure 9a shows the snowline captured on 2930 April 2016, 475 1112 days later than the ̅ yr min observed inestimated from the MODIS imagesrecord (18 April). The snowline can be distinguished by the separation between the bright white snow and discoloured ice and firn, visible above ca. 1900 m and up to 2000 m in parts.at some locations. The snowline in the 29 March 2017 image is located at ca. 1820 m, slightly above the large icefall that separates the upper and lower glacier. A much higher proportion of snow covering the glacier surface corresponds to a 6% increase in the ̅ yr min during the 2017 summer., and is found to be 10 days later using the albedo method 480 (8 Aril). This is consistent with the lower ablation reported by the EOSS programme in 2017 (Willsman et al., 2018). Despite a relatively large snowfall event brought by cyclone Gita in mid-February 2018, the image from 2829 March 2018 reveals high ablation and a similar snowline to 2016, although the proportion of exposed ice and firn appears to be slightly larger.
These events are captured in the albedo record, with a short-lived rise in ̅( ) corresponding to Gita during February, followed by a. ̅ yr min that is reached 19 days earlier (10 March) than the chosen Sentinel-2 image, with its magnitude lower than in 2018 485 compared to 2016. These images not only act to illustrate the variabilityThe Sentinel-2 observations are consistent with the relative changes of ̅ yr min , but captured by the albedo method, as well as the effects of heatwaves on New Zealand glaciers in 2016 and 2018 (Salinger et al., 2019a;2019b;.
The Sentinel-2 images shown in Fig. 9 also demonstrateillustrate the complexity of defining the summer snowline elevation on topographically complex glaciers such as Lambert. Despite the limited number of cloud-free images, the sequence of 490 Sentinel-2 images obtained during the summer and through to the arrival of seasonal snow proved to be key in interpreting the evolution of the snowline. They were found to be available within a period of less than three weeks from ̅ yr min determined using the albedo method, while the 2016 and 2017 EOSS surveys captured Vertebrae Col glacier on 11 and 9 March, or 40 and 38 days earlier than ̅ yr min , respectively (see the albedo record for Vertebrae Col glacier in Fig. 10; note that reports from the EOSS 2018 campaign have not been released to date). 495

Links between ̅ yr min and EOSS
Previous applications of the albedo method to New Zealand glaciers have developed strong relationships (R 2 = 0.89 and R 2 = 0.87, see Sirguey et al., 2016 andRabatel et al., 2017) between the MODIS-derived ̅ yr min and EOSS SLAi. Across the GoEA, SLA records only exist for both Vertebrae Col glaciers. However, with 0.7 km 2 surface area, Vertebrae Col 25 is substantially smaller than the 2 km 2 recommended to develop a reliable glacier-wide albedo average with enough MODIS pixels (Sirguey 500 et al., 2016). Although the albedo signal for Vertebrae Col 25 is obtained from only six11 250-m resolution pixels, the timeseries appears to reflect the typical seasonal cycle of albedo, albeit with a larger range in the maximum and minimum albedos than those observed from Lambert Glacier (Fig. 10).
Between 2000 and 2017, the MODIS record of ̅ yr min for Vertebrae Col 25 explains nearly half of the variability observed in the EOSS SLAi ( Fig. 11; R 2 = 0.43, p = 0.003). Alternatively, ̅ yr min of Angel Glacier exhibits the strongest relationship with 505 EOSS observations at Vertebrae Col 25, accounting for 69% of its variance (p = 1.8×10 -5 < 0.001). Other glaciers in the GoEA, namely Lambert, East Lambert, and Eve, each capture half or more of this variance as well (52%, 50%, and 55%, respectively).
Overall, it appears that the relationship between EOSS observations at Vertebrae Col 25 and ̅ yr min of each glacier is related to topographic shading (R 2 = 0.48, p-value = 0.012), slope (R 2 = 0.42, p-value = 0.022), and the north/south component of glacier aspect (R 2 = 0.34, p-value = 0.045). In contrast, glacier size (R 2 = 0.03, p = 0.609) and elevation (R 2 = 0.20, p = 0.142) exhibit 510 no significant role in determining this relationship. From the clustering analysis of the 12 glaciers in the GoEA, it is evident that changes in ̅ yr min over the unconfined glaciers in classcluster 1 are consistent to EOSS observations (mean R 2 = 0.51). More confined west-facing glaciers of classcluster 3 exhibit a weaker, yet still significant relationship (mean R 2 = 0.36). Finally, ̅ yr min of the two most confined south-facing glaciers in classcluster 2 do not have a strong relationship to the EOSS record (mean R 2 = 0.19). 515

Comparison to the EOSS survey
The EOSS survey provides a well-documented and invaluable record of the changes to glaciers in the Southern Alps over the past 40 years. Importantly, the data collected by the survey bridges the gap in glaciological mass balance records between the termination of the Ivory Glacier programme in 1975 and commencement of the Brewster Glacier programme in 2004. In 520 addition, the correlation between the snowline record and recently developed monitoring methods has allowed past mass balance trends to be reconstructed (e.g. Sirguey et al., 2016). At a larger scale, the use of EOSSAlps has proven to be an effective measure for broadly characterising the annual variability of glacier mass balance in the Southern Alps.
Published relationships developed between the MODIS-derived ̅ yr min and the EOSS SLAi departures on Brewster Glacier and Park Pass Glacier were found to be strong (Sirguey et al., 2016;Rabatel et al., 2017). Although albedo and snowline 525 methodologies differ significantly, both aim to remotely capture the point of maximum ablation at the glacier surface and use it as a proxy for mass balance. SLAiused as an estimate of the ELAprovides an effective measure of glacier mass balance (Rabatel et al., 2005), and should therefore be closely related to ̅ yr min (Dumont et al., 2012).
However, measuring SLAi remains difficult, in particular due to uncertainties associated with the method used to identify snowlines from oblique photos. Complex topography, avalanching, fresh snow, and cloud cover all present additional 530 challenges for deriving consistent results of SLAi. To compound these challenges, limited resources mean that the EOSS survey is required to assume that SLAi occurs annually in early-March, which involves careful timing of the observation flights (Willsman et al., 2018). In most cases, SLAi is estimated manually based on the overall appearance of the glacier and snow patches compared to previous years. When SLAi is digitized, photographs are compared to historical topographic maps or from orthorectified images. Furthermore, the methodology often involves an interpretation and assessment of the appearance of an 535 individual glacier in relation to other glaciers in the programme, observed during the same or from different surveys.
Despite the coarse resolution of the MODIS sensor, it has been demonstrated in this study that the albedo method is capable of retrieving robust time series of ̅ yr min that show strong to moderate relationships to EOSS signals (glaciers of Class 1 and 3, respectively). There are some clear advantages in using the albedo method over an EOSS approach, which include the approach being repeatable and not being temporally constrained. While the EOSS programme yields an archive of invaluable images, a 540 number of photos are impacted by transient snow that impact the quality of the SLAi estimates. The timing of the surveys also mean that some SLAi estimates might not fully capture late summer melt (Sirguey et al., 2016). Monitoring albedo until the onset of the accumulation season (winter) allows sustained late season ablation to be recorded. The large variability in timing of ̅ yr min , as well as the apparent trend towards a delayed occurrence for most glaciers in the GoEA (Fig. 5 and 6), demonstrates the benefit of systematically monitoring glacier surface albedo. Thus, the albedo method provides unique insights about 545 glaciers in the Southern Alps that are arguably not captured by the EOSS programme.
The EOSS programme reports significantly high intra-correlations across the 50 index glaciers that sample the Southern Alps, with pair-wise R 2 ranging 0.22 to 0.96 and averaging 0.69. Overall, our systematic application of the albedo method on the limited geographical extent of the GoEA reveals a degree of variability in glacier response that challenges the highly consistent behaviour suggested by the EOSS programme across the Southern Alps. Strong and significant intra-correlations in ̅ yr min across 550 glaciers of the GoEA exist but do not dominate, with pair-wise R 2 ranging from 0.12 to 0.88 and averaging only 0.52. Similarly, the correlation between the average GoEA albedo signal to EOSSAlps is only R 2 = 0.55 (see Fig. 7), while the average correlation between SLAi and EOSSAlps is reported at 0.81. These contrasting results suggest the EOSS approach may have led to spatial and temporal auto-correlation of the SLA records, which in turn build very strong correlations across SLA records of EOSS index glaciers and to EOSSAlps. The outcome of this is that it may have dampened or concealed the variability of mass balance 555 response of glaciers in the Southern Alps. Figure 12a illustrates the agreement between the MODImLab 250 m WS albedo retrieved from Terra and Aqua MODIS images across four days (3 March 2009, 9 March 2010, 8 March 2012, 10 March 2013. To account for the role of shadows during the 560 albedo retrieval process, each pixel within the GoEA was assigned to one of four "shading classes,", depending on the presence of shade in each image. Pixels in Class 1 were unshaded in both images, Class 2 and Class 3 pixels were only shaded in the Aqua or Terra image respectively, and Class 4 pixels were shaded in both images. Of the total 2196 matching pixels from these image pairs, 83.5% belonged to Class 1. The linear regression for these 1834 pixels is significant (p < 0.01), and strong (R 2 = 0.58), although slightly askew to a linear 1:1 relationship (Fig. 12a). Individually, each of the four days displays a similar trend 565 and distribution, with R 2 values ranging between 0.51 and 0.71, and gradients between 0.75 and 0.85. Figure 12a also shows an increase in the variability of albedo between the sensors at higher values. Overall, this agreement is consistent with the expected 10% accuracy of the albedo retrieval method (Dumont et al., 2011(Dumont et al., , 2012Sirguey et al., 2016).

MODIS albedo
The increased variability coincides with the highest density of points; where 83% of pixel albedo values fall between 0.4 and 0.8. Due to the uneven distribution of albedo across the spectrum, a second linear regression was run on a stratified random 570 sample of 150 Class 1 pixels (Fig. 12b). The linear regression of these resampled data show a much-improved fit (R 2 = 0.88, p < 0.001) that closely approximates the 1:1 relationship, and demonstrates that the degraded band 6 of MODIS A does not compromise MODImLab albedo retrieval. The slightly lower albedo found in Aqua images compared to Terra may also be explained by the timing, as snow and ice surfaces would undergo transformation compatible with a decrease in albedo from morning to afternoon. 575 To characterise the variability between the sensors in more detail, we present the spatial distribution of averaged residuals between maps of glacier surface albedo from MODIS T and MODIS A across the GoEA (Fig. 13). Although Fig. 13 confirms the good agreement between the sensors over large parts of the GoEA, it shows large departures around the fringes and near rock outcrops (+0.23/-0.38), often close to steep, complex terrain and involving mixed pixels. In particular, two large steepsided rock outcrops in the south of Lambert Glacier and Westthe west icefall of Angel Glacier correspondsdescending from 580 Mt Farrar correspond to MODIS A albedo being substantially larger than MODIS T . Close inspection of imagery before and after correction, as well as shadow maps reveal the larger extent of cast shadows produced by outcrops in the afternoon and affecting MODIS A imagery. Issues of overcorrecting spectral reflectance in cast shadows is known to challenge MODImLab (Davaze et al., 2018), and appears again to cause overestimationsoverestimates of MODIS A albedo. Figure 12 also demonstrates this effect with overestimation of albedo by MODIS A relative to MODIS T for Class 2 pixels (MODIS A pixels in 585 the shade) and conversely for Class 3 pixels (MODIS T pixels in the shade).
MODIS A involves a delay in image capture from 10:30 am to 1:30 pm. Towards the winter solstice, this decreases the proportion of shaded pixels in the steep terrain of the Southern Alps. The change in the solar zenith angle caused by the delayed timing of the image capture means pixels on south-and southwest-facing slopes have a higher chance of receiving incident radiation. This effect was seen across the GoEA, where a much lower proportion of the total pixels in the GoEA were shaded 590 in MODIS A images captured near the winter solstice, compared to MODIS T (34.5% and 49.8%, respectively). However, in summer and toward April, steep rock outcrops cast large shadows in the afternoon that challenge albedo retrieval with MODIS A . By May, low sun zenith angles cast longer shadows that are not accurately modelled due to the resolution and accuracy of the DEM. Unpredicted shadows yieldsyield severe underestimations of surface albedo, in particular affecting confined glaciers of classcluster 3. Although MODIS A images could help capture a better albedo signalssignal in such cases, 595 Figure 5 (left)5a demonstrates that, over the period of this study, imagery from MODIS T remained suitable to capture ̅ yr min before this issue becomes problematic by May. Nevertheless, despite the computational burden, the systematic processing of MODIS A imagery may gain merit as the timing of ̅ yr min appears to be more often delayed intoto late summer or early autumn more often (Fig. 5).

MODIS cloud cover 600
Finally, Lyapustin et al. (2014) documented calibration issues with MODIS Collection 5 data, and concluded that major calibration trends were removed in C6. The stability of the C6 calibration was confirmed by Sayer et al. (2015). In the context of retrieving time series of snow and ice albedo with MODIS, Casey et al. (2017) stressed how C5 data could compromise the detection and interpretation of trends, and concluded that C6 is preferable. In the case of MODIS Terra, albedo time series derived from C6 data by MODImLab shown in Fig. 8 and 10 reveal no visible trend in winter albedo. It seems unlikely, if a 605 trend in winter snow albedo existed, that it would be matched and concealed by a calibration issue of exactly the opposite magnitude. Our results therefore support the alternative hypothesis that there is no detectable trend in winter albedo nor calibration issue over the length of the record produced by this study. This is further supported by the cross-platform agreement shown in Fig. 12.

MODIS cloud cover 610
While MODIS A potentially provideprovides a means to capture more shadow-free pixels across the GoEA, its use is inhibited by daily development of cloud cover over the Southern Alps. Figure 14 shows MODIS A images display a consistently higher proportion of cloudy pixels over the GoEA than MODIS T images over the course of a year. While this trend is consistent throughout the year, it is pronounced through the summer ablation period, at the critical time when the ̅ yr min is likely to be observed. A likely driver of this pattern is the increased surfacedaytime heating of the atmosphere through the warmer summer 615 months, driving freeresulting in convection and aiding cloud development over the course of the day. Ultimately, while MODIS A shows some promise in its ability to supplement the MODIS T dataset as discussed above, the high variability of albedo from mixed-pixels, and the daily increase in cloud cover means that its application is limited in New Zealand.
Interestingly, the cloud cover results from both MODIS A and MODIS T images are still of use, as the spatial characterisation of cloud cover over the Southern Alps is very limited. (Wardle, 1986). Current research is largely limited to point based cloud 620 data from automatic weather stations (e.g. Conway et al. 2015). Cloud cover patterns and dynamics are important for glaciers, as clouds play a key role in influencing incident solar radiation at the glacier surface . Figure 15 displays the non-uniform distribution of monthly average cloud cover across the processed area. TheOver the period of study, the frequency of cloudclouds in pixels west of the Main Divide is as high as 90% during summer months, and reaches a minimum of 35% in some areas during winter, reflecting Fig. 14. During summer, cloud spill-over is limited to within a few 625 kilometerskilometres east of the Main Divide, contrasting to relatively uniform conditions through winter. It is also possible to identify the presence of 'cloud hot-spots' that persist through the 19-year record. The spatial complexity of cloud cover is likely having a significant influence on the SEB, which would suggest the physical processes controlling glacier behaviour will be different from one individual glacier to another, which is apparent by the contrasting albedo records for glaciers of the GoEA.'cloud hot-spots' that persist through the 19-year record. 630

Comparison to the EOSS survey
The EOSS survey provides a well-documented and invaluable record of the changes to glaciers in the Southern Alps over the past 40 years. The data collected by the survey bridges the gap in glaciological mass balance records between the termination 635 of the Ivory Glacier programme in 1975 and commencement of the Brewster Glacier programme in 2004. In addition, the correlation between the snowline record and recently developed monitoring methods has allowed past mass balance trends to be reconstructed (e.g. Sirguey et al., 2016). At a larger scale, the use of EOSSAlps has allowed the annual variability of glacier mass balance in the Southern Alps to be broadly characterised (Chinn et al., 2012;Willsman et al., 2018).
Published relationships developed between the MODIS-derived ̅ yr min and the EOSS SLAi departures on Brewster Glacier and 640 Park Pass Glacier were found to be strong (Sirguey et al., 2016;Rabatel et al., 2017). Despite the large differences between the albedo and snowline methodologies, both aim to remotely capture the point of maximum ablation at the glacier surface and use it as a proxy for mass balance. SLAiused as an estimate of the ELAprovides an effective measure of glacier mass balance (Rabatel et al., 2005), and should therefore be closely related to ̅ yr min (Dumont et al., 2012).
However, measuring SLAi remains difficult, in particular due to uncertainties associated with the method used to identify 645 snowlines from oblique photos. Complex topography, avalanching, fresh snow, and cloud cover all present additional challenges for deriving consistent results of SLAi. To compound these challenges, limited resources mean that the EOSS survey is required to assume that SLAi occurs annually in early-March, which involves careful timing of the observation flights (Willsman et al., 2018). In most cases, SLAi is estimated manually based on the overall appearance of the glacier and snow patches compared to previous years. When SLAi is digitized, photographs are compared to historical topographic maps or from 650 orthorectified images. Furthermore, the methodology often involves an interpretation and assessment of the appearance of an individual glacier in relation to other glaciers in the programme, observed during the same or from different surveys.
Despite the coarse resolution of the MODIS sensor, it has been demonstrated in this study that the albedo method is capable of retrieving robust time series of ̅ yr min that show strong to moderate relationships to EOSS signals (glaciers of clusters 1 and 3, respectively). There are some clear advantages in using the albedo method over an EOSS approach, which include the 655 approach being repeatable and not being temporally constrained. While the EOSS programme yields an archive of invaluable images, a number of photos are impacted by transient snow that impact the quality of the SLAi estimates. The timing of the surveys also mean that some SLAi estimates might not fully capture late summer melt (Sirguey et al., 2016). Monitoring albedo until the onset of the accumulation season (winter) allows sustained late season ablation to be recorded. The large variability in timing of ̅ yr min , as well as the apparent trend towards a delayed occurrence for most glaciers in the GoEA ( Fig. 5 and 6), 660 demonstrates the benefit of systematically monitoring glacier surface albedo. Thus, the albedo method provides new insights about glaciers in the Southern Alps that are not captured by the EOSS programme.
The EOSS programme reports considerably high intra-correlations across the 50 index glaciers that sample the Southern Alps, with pair-wise R 2 ranging 0.22 to 0.96 and averaging 0.69. Overall, our systematic application of the albedo method on the limited geographical extent of the GoEA reveals a degree of variability in glacier response that challenges the highly consistent 665 behaviour suggested by the EOSS programme across the Southern Alps. Strong and significant intra-correlations in ̅ yr min across glaciers of the GoEA exist but do not dominate, with pair-wise R 2 ranging from 0.12 to 0.88 and averaging only 0.52. Similarly, the correlation between the average GoEA albedo signal to EOSSAlps is only R 2 = 0.55 (see Fig. 7), while the average correlation between SLAi and EOSSAlps is reported at 0.81. These contrasting results suggest the EOSS approach may have led to spatial and temporal auto-correlation of the SLA records, which in turn build very strong correlations across SLA records of EOSS 670 index glaciers and to EOSSAlps. The outcome of this is that it may have dampened or concealed the variability of mass balance response of glaciers in the Southern Alps.

Implications of variability of albedo on mass balance
Given the scarcity of surface mass balance measurements in the Southern Alps, and the spatial and temporal limitations of 675 obtaining imagery from observational flights, satellite remote sensing provides a powerful tool to increase the number of glaciers being currently monitored in New Zealand. To relate the observed variability of albedo to changes in annual mass balance, the relationship between the magnitude and timing of ̅ yr min in each of the three glacier clusters identified on the GoEA (Table 2) needs to be further examined. There is evidence that the occurrence of ̅ yr min has been delayed in time (later in the year) on all glaciers other than those located on steep, south-facing slopes (Cluster 2) (Fig. 5d). Intuitively, a delay in the 680 timing of ̅ yr min on 10 of the 12 glaciers (Clusters 1 and 3) might be expected to result in a more negative summer mass balance by virtue of the ablation season being extended. However, there is no visible trend in the magnitude of ̅ yr min in any of the clusters during the observation period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Also, there is only a weak relationship between a delay in timing and the magnitude of ̅ yr min (Fig. 5c), suggesting that glacier-wide albedo is not necessarily lower if it is delayed. We interpret this dichotomy as evidence that the atmospheric controls on mass balance are complex, and that the magnitude of ablation is in 685 part sensitive to discrete weather events (Cullen et al., 2019), which either bring snowfall (increase albedo) or enhance ablation (decrease albedo) in summer. There is also compelling evidence that clouds play an important role in modulating the energy available for ablation over glaciers in the Southern Alps , and are likely to influence both the timing and magnitude of ̅ yr min in the GoEA, thus confounding the relationship. In particular, clouds are very common to the west of the Main Divide in the GoEA region during summer (Fig. 15), and their spatial and temporal variability elsewhere may be an 690 underlying factor in controlling the observed variability in ̅ yr min . We believe there is a risk that some of this complexity is being missed using the single snapshot approach of the traditional EOSS programme. Arguably, the albedo method provides more certainty of the atmospheric controls on mass balance, as it provides a platform of continuous observations of albedo and cloud cover over a large number of glaciers throughout the ablation season.
Despite the observed complexity between the timing and magnitude of ̅ yr min over the observation period, there is evidence of 695 a shift towards more negative annual departures of the mean ̅ yr min since 2008 (Fig. 7), which is consistent with the EOSS record, implying cumulative losses in mass balance for the glaciers of the GoEA. This mass loss is consistent with glaciological observations of mass balance from Brewster Glacier (Cullen et al., 2017) and the modelling of glacier response to changes in climate in the Southern Alps Salinger et al., 2019a, b). This shift towards mass loss appears to be widespread, with the 2018 summer heatwave reportedly the most damaging for glaciers in the Southern Alps over the last half 700 century (Salinger et al., 2019a). We anticipate that the GoEA are particularly vulnerable to a warming climate, as the average elevation of the icefields is close to the regional snowline, which is estimated to be 1950 m a.s.l. (Chinn, 2001). Cullen and Conway (2015) showed that the majority of precipitation at the altitude of Brewster Glacier falls at an air temperature that is very close to the threshold between snow and rain. Thus, small increases in air temperature controlled by warming is likely to impact the amount of precipitation falling as snow on the GoEA, which in turn will impact surface albedo and mass balance, 705 potentially driving a rapid decline under a sustained warming trend. While the observations we have obtained using the albedo method suggest that most glaciers are vulnerable, the contrasting response of glaciers contained on steep, south-facing slopes (Cluster 2) indicates they are likely to survive this demise longer.
Lastly, the positioning of the GoEA across the Main Divide of the Southern Alps is particularly important from a management perspective. Because the icefields straddle the Main Divide, they contribute meltwater to the catchments of the Rangitata River 710 on the east coast, and the Wanganui and Whataroa rivers to the west. Given the important role of snow and ice in New Zealand's water resources, changes in the volume of these icefields have the potential to impact the hydrology of these rivers in the future.

Limitations and developments of the albedo method
The results presented in this study rely on the inherent accuracy of MODImLab. The specific error associated with the 715 MODImLab retrieval over the GoEA is uncertain due to the lack of in situ data. However, Dumont et al. (2012) quantified the relative error between field measurements and the MODImLab 250 m broadband albedo to be approximately ± 10% (RMSE = 0.052), confirmed by Sirguey et al. (2016). This value is an average estimate, with error recognised to be up to twice as high around the mixed-pixel margins of glaciers compared to the clear snow/ice pixels near the centre (Dumont et al., 2012). As a result, it is expected the uncertainty will be variable between glaciers, where estimates of ̅( ) on large glaciers (i.e. Lambert 720 Glacier) may be more reliable than those of small glaciers with high edge effects (i.e. East Lambert Icefall) (Brun et al., 2015).
Supporting Rabatel et al. (2017) and Davaze et al. (2018), we find that the artefacterror in the albedo retrieval caused by surface shading does not affect the identification of the summer ̅ yr min . On glaciers where shading is most prominent (steep sloping, south-facing, topographically incised), we show the ̅ yr min is typically reached much earlier than early-April (when the shadinginduced decrease in ̅( ) begins). Importantly, thisThis means that while New Zealand has a large number of glaciers that 725 exhibit similar hypsometry to those glaciers in Cluster 2, in most cases, the albedo retrieval will still yield reliable values of ̅ yr min .
The masking technique proposed in this study should simplify the process of objectively converting glacier outlines to glacier masks. As we show in Sect. 3.2 and 4.1, the new approach to masking glacier boundaries over the GoEA yields values of ̅( ) over the important ablation period within 3% of the more subjective method employed by previous studies. The success of this 730 technique may support the rapid application of the albedo method to new glaciers using existing outlines, regardless of whether they are made up of mixed-pixels and/or debris-cover. In addition, the snow-covered-area filter may help to reduce some of the problems of using a static glacier mask (e.g. the glacier extent changing significantlysubstantially over the observation period).
Finally, the use of MODIS imagery to apply the albedo method remains applicable only for glaciers that are large enough to 735 allow albedo to be retrieved and averaged over a number of pixels. Davaze et al. (2018) show that ̅ yr min can be successfully retrieved on glaciers covering as little as 0.5 km 2 with good correlation to annual mass balance. This is consistent with the agreement we obtained between ̅ yr min and EOSS on Vertebrae Col 25 (0.7 km 2 or 11 pixels), although it is believed to stretch the reliance on MODIS for the albedo method. The temporal resolution achievable by combining multi-spectral imaging from higher resolution sensors such as Sentinel-2 and Landsat 8 OLI justifies the development of snow and ice albedo products to 740 resolve ̅ yr min on smaller glaciers. Such development is desirable to advance the albedo method and support the widespread monitoring of New Zealand glaciers.

Conclusion
The results from this study represent the next meaningful step towards the use of the albedo method for widespread monitoring of glaciers in the Southern Alps. Following the successful application to Brewster Glacier and Park Pass Glacier, we have 745 produced ana 19-year long coherent seasonal signal of glacier-wide albedo on the previously unstudied Gardens of Eden and Allah (GoEA). These results have supported and advanced key aspects of the methodology that will be important forinform future applications in the Southern Alps and beyond. The key findings can be summarised as: 1. A new objective glacier masking approach has been developed that compares favourably to a more traditional manual method of identifying suitable pixels to calculate glacier-wide albedo. 750 2. The annual minimum glacier-wide albedo ( ̅ yr min ) for individual glaciers ranges between 0.42 and 0.70, and can occur as early as mid-January and as late as the end of April. The timing of ̅ yr min appears to have shifted to later in the year over the 19-year period on glaciers that are not steep or south-facing.all glaciers other than those located on steep, south-facing slopes. However, there is only a weak relationship between the delay in timing and the magnitude of ̅ yr min , which suggests that glacier-wide albedo is not necessarily lower if it is delayed. 755 3. The glacier-wide surface albedo anomaly for the 19-year period explains 55% of the variability in the average annual departure of the 50 index glaciers from the EOSS programme (EOSSAlps). The largest negative departures in ̅ yr min (lower than average albedo) are consistent with high snowlines, with the 2018 departure the most negative on record. This is consistent with Brewster Glacier and was caused by a marine and terrestrial heatwave over New Zealand (Salinger et al., 2019a). 760 4. The MODIS record of ̅ yr min for Vertebrae Col 25 explains less than half of the variability observed in the EOSS SLA record (R 2 = 0.43, p = 0.003). However, the relationship is stronger when compared to other GoEA glaciers, with Angel Glacier having the strongest relationship with EOSS observations at Vertebrae Col 25, accounting for 69% of its variance (p = 1.8×10 -5 < 0.001). Lambert, East Lambert, and Eve glaciers each capture half or more of the variance (52%, 50%, and 55%, respectively). The relationship between EOSS observations at Vertebrae Col 25 and ̅ yr min of 765 each glacier is related in order of importance to topographic shading, slope and aspect.
5. The EOSS programme has reported on how strongly each of the 50 index glaciers behaviour is related to the mean of all remaining glaciers (known as EOSSAlps), and pair-wise regression shows there is high intra-correlation between glaciers. However, the albedo method enables the variability in response of individual glaciers to be explored in more detail, revealing that topographic setting plays an importanta second order rolecontrol in addition to the regional 770 climate signal (first order control). The albedo method captures enough individual glacier variability on the GoEA to firmly question the validity of the hypothesis that glaciers in the Southern Alps behave as a single climatic unit. 6. For the first time, MODIS imagery acquired by the Aqua platform (MODIS A ) has been used successfully to increase the temporal resolution of albedo monitoring using the MODImLab algorithm. There is some evidence to suggest it is capable of capturing diurnal variability in albedo as controlled by changes to snow and ice properties during 775 daytime. Despite cloud being more frequent during the afternoon, especially in summer, there are advantages in using MODIS A due to higher incident radiation (less shading) on some slopes at certain times of the year. 7. Cloud cover results from MODIS imagery acquired by the Terra (MODIS T ) and Aqua (MODIS A ) platforms show the spatial and temporal variability of clouds. The frequency of cloud in pixels west of the Main Divide is as high as 90% during summer months, and reaches a minimum of 35% in some locations in winter. There is a strong gradient in 780 cloud cover frequency between regions west and east of the Main Divide, and specific areas appear to be consistently more cloudycloudier than others. These complex cloud interactions deserve further attention as they are likely to play an importanta considerable role in glacier surface energy and mass balance.
The key findings presented in this research have provided a platform to further develop and extend the application of the albedo method to monitor glacier behaviour in the Southern Alps. The next logical step will be to attempt to resolve spatial and 785 temporal patterns in glacier behaviour inan assessment of all of the main glaciated areas of the Southern Alps together at a high temporal resolution, which will complement observations made by the EOSS programme, and expand our understanding of the linkages between glaciers and the climate system. There is some urgency to do this as our observations of the GoEA, and those obtained from traditional glaciological observations elsewhere, suggest that glaciers in the Southern Alps are undergoing an unprecedented decline at present. 790 Data and code availability. MODIS data used in this research are freely available from the Level 1 and Atmosphere Archive and Distribution System (LAADS) Web. NZSoSDEM is freely available from the koordinates.com geographical data repository. The MODImLab software is available upon request from PS.

795
Author contributions. PS and NC initiated and coordinated the study. NC and PS obtained funding for the research. AD and PS processed and analysed the MODIS data. AD ledwrote the writingfirst draft of the manuscript with inputs, while PS and support from co-authors.NC were responsible for the submission and revision of the final research.
Competing interests. The authors declare that they have no conflict of interest. 800 Tables Table 3: Correlation between the median timing and magnitude of the ̅ yr min and glacier hypsometry using the Spearman's rank coefficient (r) for the 12 outlet glaciers of the GoEA. 2-tailed correlation is significant at the 0.05 level (*) and the 0.01 level (**).