Unravelling the long-term, locally-heterogenous response of Greenland glaciers observed in archival photography

. We present an approach for extracting quantiﬁable information from archival aerial photographs to extend the temporal record of change over a region of the central eastern Greenland Ice Sheet. The photographs we use were gathered in the 1930s as part of a surveying expedition, and so they were not acquired with photogrammetric analysis in mind. Nevertheless, we are able to make opportunistic use of this imagery, as well as additional, novel data-sets, to explore changes at ice margins well before the advent of conventional satellite technology. The insights that a longer record of ice margin change bring is 5 crucial for improving our understanding of how glaciers are responding to the changing climate. In addition, our work focuses on a series of relatively small and little studied outlet glaciers from the eastern margin of the Ice Sheet. We show that whilst air and sea surface temperatures are important controls on the rates at which these ice masses change, there is also signiﬁcant heterogeneity in their responses, with non-climatic controls (such as the role of bathymetry in front of calving margins) being extremely important. In general, there is often a tendency to focus either on changes of the Greenland Ice Sheet as a whole, 10 or on regional variations. Here, we suggest that even this approach masks important variability, and full understanding of the behaviour and response of the Ice Sheet requires us to consider changes that are taking place at the scale of individual glaciers.

in the cryosphere. Such abilities are important in our efforts to investigate the links between climate and glacier change, and variations in how glaciers respond to climate change across the world. Despite the undeniable value of satellite observations, 60 and the insights they have afforded into cryospheric change, the period prior to the launch of ERTS1 is characterised by relative data sparseness ::::::::::::::::: (Goliber et al., 2021). In light of this, significant advances would be gained from further expansion of the record of glacier change into the past.
Here, we exploit a series of images gathered obliquely along the east coast of Greenland between Kangerdlugssuak and Umivik, along a c. 260 km-long section of coastline between 66.3°and 68.4°N, taken for surveying purposes. These images 65 were gathered between 1930 and 1931 as part of the British Arctic Air Route Expedition (BAARE), which was carried out in an effort to discover the possibility of a new and shorter transit route between the UK and Canada. This route, in part, passes over Greenland, and one of the mission's aims was to survey the eastern and central parts of Greenland -the section of the proposed route that was least well known. The BAARE survey team did this using a ship and sea-plane, in an effort to photograph and map the coastline (The SPRI Picture Library, 1999). 70 The imagery that we utilise ::::: utilize : here provides oblique views of a number of outlet glaciers in two nearby regions. The focus of our study is the opportunistic 'snapshot' that these images provide of the state of these glaciers during the BAARE survey in 1930 (Shepherd and IMBIE Team, 2020). In this paper, we utilise ::::: utilize Structure from Motion (SfM) approaches to build georeferenced orthophotos of the terrain within the BAARE imagery, so as to extract information about glacier extent from over 90 years ago.

75
This provides an exceptional and important additional key to understanding ice mass change in this region way before the advent of satellite technology. To further supplement our investigations, we also add in additional steps between 1930 and the start of the Landsat record with orthophotos that we generate from imagery of the region from the now-declassified 1960s Corona ::::::::: CORONA satellite mission (Shin, 2003), and similarly from aerial photography from the 1980s (Bjørk et al., 2012).
Finally, we also explore the Landsat record from 1985 up to the present day. Overall, this suite of data provides unprecedented 80 insights into the changes that have taken place over >90 years in this relatively poorly understudied ::::: studied : region of East Greenland, where important changes are nevertheless known to have taken place. To investigate this further, we also explore changes in air temperature, sea surface temperature and mass balance in an effort to identify the drivers of glacier change here.

Study area
Some of the imagery acquired during the BAARE expedition covered a part of Eastern Greenland that approximately corre-85 sponds to the 'Central Eastern Region' (CE) as defined by Mouginot et al. (2019) in their delineation of Greenland into 7 discrete regions. Other, more recent work, by King et al. (2020) divided the GrIS up into just 4 regions, with the BAARE sector under investigation here corresponding to the 'Southeast Region' (SE). This SE region is of great interest because whilst in all other parts of Greenland, glacier thinning is due (at least in part) to glacier discharge being greater than the balance flux (indicating dynamic disequilibrium), in this SE region, the primary cause of thinning prior to 2000 was increased surface melt, 90 indicating that this region responds more rapidly to climatic forcing, reinforcing observations made by Hugonnet et al. (2021)  of glacier fluctuations in agreement with precipitation and temperature fluctuations. Dividing the ice sheet into discrete regions like this has proved to be a powerful approach for exploring broad-scale patterns. The biggest glaciers in the area are Hutchinson Gletscher, Polaric Gletscher and KVJ Steenstruup Nordre Brae, but the region is a mountainous and dense fjord and valley system and so glacier outflow is primarily dominated by relatively small outlet glaciers. The vast majority of the glaciers in our study are outlets of the GrIS but as shown in Figure 1, a small number are identified as being smaller local glaciers and ice caps (GICs) peripheral to the margins of the main GrIS (Rastner et al., 2012).

100
As a consequence of the mountainous/alpine terrain, as well as its climatology, this area has much higher accumulation rates compared to the rest of the GrIS (and surface mass balance remains positive, which is discussed subsequently, cf. Figure 7c)).

Archival photography
In order to retrieve historical information on the geometry of the 24 East Greenland glaciers we used three sets of photographic 105 data. Each of these sets was obtained using a significantly different imaging set-up, thus requiring customisation of the required processing methods.
British Arctic Air Route Expedition, aerial oblique images -1930-31 The BAARE took place between July 1930 and August 1931. It was a privately funded expedition to investigate the feasibility of a new and shorter air passage between England and Canada. Part of the survey involved photogrammetric reconnaissance.

110
This was done with the use of two De Havilland DH.60 Moth planes with Gipsy 1 engines (Stephenson, 1932;Aviation Safety Network, 1999). One of the planes was equipped for taking vertical and oblique photographs (Watkins et al., 1932). A Williamson P14 camera with a lens of known focal length of 2209.8 mm (7.25'), and 127 mm × 101.6 mm (5' × 4') glass plates with envelope adaptors for changing slides in daylight were used. Each flight took approximately 90 minutes and the plates were changed every 30 seconds. The time interval allowed for about 65% overlap on the photographs (Watkins, 1930).

115
During the summer of 1930 a total of 9 photographic flights (18h 20m) producing 450 plates were carried out. This covered the area from Bjorne Bugt up to and including Kangerdlugsuak Fjord, and also some parts of Sermilik Fjord and Angmagssalik Island. In the summer of 1931, due to poor weather conditions and the subsequent required dismantling and then rebuilding of each of the aircraft over winter and early spring, only two flights of 7 hours were carried out covering the area of Sermilik Fjord from Sermilik up to Umivik (Watkins et al., 1932). From all of this work, only 248 photographic plates remain, with scanned 120 versions held in the Picture Library of the Scott Polar Research Institute (The SPRI Picture Library, 1999). Unfortunately, the remaining images were lost after the expedition due to poor operational logistics by the returning party, as well as technical problems with the processing of the plates. Also some batches of plates were deemed to be unusable by cursory inspection and were subsequently destroyed. For our study we used 73 images obtained during the summers of 1930 and 1931.

125
The CORONA satellite mission was a clandestine surveillance program lead ::: led by the CIA (Central Intelligence Agency) of the United States of America, and the U.S. Air Force, aimed at gathering spatial data for the creation of maps of vast remote areas for intelligence purposes (Goossens et al., 2006). Its existence was not acknowledged until the data entered the public domain in 1995. The Corona :::::::: CORONA : data can currently be obtained (as digital high-resolution scans; 7 µm) from the EarthExplorer website (U.S. Geological Survey USGS, 1995;Shin, 2003).

130
The CORONA mission was a vanguard in the early days of satellite surveillance. As such it piloted the use of sophisticated methods of shutter and camera construction, as well as an innovative and sometimes unreliable means of data retrieval. As a result, each mission had different operating specifications and a large amount of the data collected was not successfully retrieved. Generally speaking, for each mission the plan was to launch the satellite to a predetermined height, capture images on photographic film and then allow the satellite to return through the atmosphere and disintegrate on entry. Prior to this, a capsule 135 5 was jettisoned from the satellite and parachuted towards Earth containing the exposed film. This capsule was intercepted on descent by a plane (Galiatsatos, 2005). Since this was an intelligence collecting mission the capsule was designed to selfdestruct if it was not intercepted before it reached a critical height. There was only one known occurrence when the selfdestruction mechanism did not work as intended, and the capsule landed on the surface of the earth (Pieczonka et al., 2011).
In our study we used images taken on 24th September 1966. For this mission the KH-4M camera was used. This was 140 an upgrade version of the KH-4 camera -the first stereoscopic camera used in space, providing 75% overlap. The KH-4A (Keyhole-4A) carried two J-1 (in earlier missions KH-3 cameras of 3.66 m resolution) panoramic cameras, with a focal length of 61 cm, and a ground resolution of 2.7 m to 7.6 m. It also carried a 4 cm index camera, with a focal length of 38 mm, a ground resolution of 162 m, and frame coverage of 308 km × 308 km. The J-1 cameras were placed on an M (Mural) mount, one pointing 15°aft from the vertical and the other 15°forward (Galiatsatos, 2005). The minimal flight height was 180 km and 145 the duration of each mission was 14-15 days. Additional meta-data, such as ephemeris, ground velocity of the platform and the scan rate, the photographic coordinates of the principal points and the fiducial marks are not available (Shin, 2003).
The images obtained with the CORONA cameras have a complex image geometry (Casana and Cothren, 2015). The panoramic cameras used (also for aerial photography) work on the general principle that during the scanning process the lens and the scan arm moves while the film remains stationary. In this case the lens rotates around the second nodal point allowing 150 the cylindrical focal plane to keep the image of distant objects sharp. As a result a 'bow-tie' shaped region is photographed and becomes compressed into a rectangular image. This effect creates significant panoramic image distortions. Additional significant imaging issues associated with those pictures are scan position distortion resulting from motion during the scanning process, image motion compensation distortion, tipped panoramic distortion, and geometric distortions resulting from roll, pitch, yaw and altitude instability. Many of these effects could be rectified with rigorous geometric distortions corrections, as 155 is done with current satellite imaging systems. However, the lack of available meta-data makes such an approach intractable, thus necessitating a more customised approach (Shin, 2003;Galiatsatos, 2005).
Greenland 1:15000 scale, vertical aerial images -1978-87 Aerial photographic missions were carried out between 1978-87 by the Geodetic Institute, the National Cadastre and Survey of Denmark, and the Danish Geodata Agency. More recently, these organisations were merged and renamed as the Agency for 160 Data Supply and Efficiency (SDFE), which holds the records of the survey including the original photographs, scans, flight plans, and calibration data for both cameras. There is also GCP (Ground Control Point) data (obtained via triangulation, aerotriangulation and Doppler measurements) and this was used for the creation of a DEM model in the early 2010s (Bjørk et al., 2012). The photographic data covers all of Greenland together with the surrounding smaller islands, but excludes the interior of the ice sheet (Korsgaard et al., 2016). A WILD RC10 camera with a nominal focal length of 87.72 mm was used to collect 165 super-wide-angle photographs at planned flying heights of 13 000 m. The images were captured on photographic film, in black and white and with 8 fiducial marks on each image. For our study we used 58 images obtained on 30th of July and 14th of August 1981. All were captured in favourable weather conditions and with at least 66% overlap between frames. used. This is a relatively new product provided by the National Science Foundation (NSF) and the US National Geospatial-Intelligence Agency (NGA), which has been produced since 2015. The dataset is constructed by combining in-track and crosstrack high-resolution (about 0.5 m) imagery acquired using the DigitalGlobe constellation of stereoscopic optical imaging satellites, and including WorldView-1, WorldView-2, WorldView-3, and GeoEye (Meddens et al., 2018). It is created using Surface Extraction with the TIN-based Search-space Minimization (SETSM) algorithm (Noh and Howat, 2015). ArcticDEM 175 raw products are additionally georeferenced by alignment to ICESat point cloud that has high 0.01 m ± 0.07 m accuracy but coarse measurement footprint of 70 m (Morin et al., 2017). The current version of the ArcticDEM offers a highest resolution of 2 m raw data strips (day stamped) and a DSM mosaic averaged over time and area with a resolution of 2 m.
There has been little research aimed at establishing the defined accuracy or consistency of these models. However, our own experiments and the results reported in Błaszczyk et al. (2019) both suggest that the 2 m day stamped strips are often 180 wrongly aligned and prone to artefacts. These artefacts are usually 3D representations of cloud cover or random 'tower'shaped elements (Crosby, 2016;Meddens et al., 2018). Moreover, the artefacts are hard to recognise in raster format due to a lack of corresponding texture, but become obvious after export to a point cloud format. Also, many of the strips are incomplete, having empty pixels. Thus it would be only possible to use them if combined with additional strips of the same area obtained on a different date (Barr et al., 2018). Lastly, the strips of adjoined areas have not all been captured at the same time. ArcticDEM : quality, the accuracy for our GCP was assumed to be the same as the pixel size of the ArcticDEM 2 m mosaic. It is also important to mention that we encountered areas on the mosaic that were clearly artefacts and we removed them from further analysis.

Structure-from-Motion based orthomosaics
Structure-from-Motion (SfM) has rapidly become one of the most popular means of obtaining 3D data from image sequences.

195
Most SfM algorithms seek to simultaneously estimate a 3D scene model (sparse point clouds), camera intrinsic parameters (focal length, centre of projection etc.) and camera extrinsic parameters (3D pose, translation, rotation) from a set of overlapping images. This is aided, if needed, by geographical localisation information provided by GCP or camera path data. In  point cloud with 3D points consisting of matched 2D features (Ryan et al., 2015). Using fixed camera parameters so-obtained, a dense point cloud can be estimated using a process called Multi-View Stereo (MVS). This process is often based on performing dense binocular stereo between pairs of images with a large overlap. As a result multiple depth maps are effectively combined.
A dense point cloud can be transformed into a mesh and then rendered with textures extracted from the original images (Park 205 and Lee, 2019; Yurtseven et al., 2019).
The outputs of SfM can be used in a number of ways. One of the most relevant to the work reported here is the creation of an orthomosaic or orthophotomap. An orthomosaic is a combined image created by the seamless or near seamless merging of the original images projected onto the plane (or DEM/mesh model) and then transformed to the required projection. During this process the images are ortho-rectified (geometrically corrected) such that the scale is uniform and so that a photo or image 210 adheres to a given map projection (Lamsters et al., 2020;Agisoft Metashape, 2020). In the work reported here, we created orthophotomaps based on the images described in Section 2.1. This allowed us to create archival orthophotomaps and compare them to current Landsat-based satellite images (described in Section 2.2) thus providing detailed information on glacier front movements or overall glacier movement in the region studied. All of the datasets used the same projection, namely WGS84/NSIDC Sea Ice Polar Stereographic North (EPSG:3413). In order to georeference the archival orthophotomaps, a 215 number of GCPs obtained from the ArcitcDEM :::::::::: ArcticDEM were used. This number varied with the size of the area constituting each orthophotomap, but generally the 1960s dataset required the largest number of points due to its unique distortion properties. For each orthophotomap between 35 and 100 GCPs were used (Table 1) .

2.2 Archival orthophotomaps
The archival orthophotomaps were produced with the use of Agisoft Metashape (Agisoft Metashape, 2020). Since we did 220 not have sufficient data to select the images with the best overlap or with the best light conditions, it was decided to use all the available images for this procedure. Initially we divided the 1930s data set into 5 regions for the production of 5 orthophotomaps. However, the 1960s and 1980s datasets were significantly different in terms of their extent and overlap. This forced us to divide our area into different sub-regions (Table 1) in order to produce mosaics of the same glaciers as covered by the 1930s images. In Figure 2 examples of typical source imagery are provided so that these can be compared. as well as 225 derived SfM-based orthophotomaps. Table 1 describes the data used and the accuracy of the results obtained from them. For all of the images we found the corresponding areas on the ArcticDEM model and created GCPs on stable, non-ice covered bare ground, which we assume to be fixed over the time period covered. The GCP placement accuracy calculated during the processing of the mosaics was around 1 pixel and in 90% of cases was smaller than 0.8 pixel. The spatial accuracy in metres varied, but in most cases was 230 less than 20 m, and did not exceed 15 m in the X and Y directions separately. This result can be considered satisfactory when taking into account the quality of the ground truth model, problems with the definition of the stable areas for GCPs, and the age and type of the archival images. We also considered isostatic uplift of GCPs over the 90 year period of our investigation.
Shepherd and IMBIE Team (2020) explored rates of isostatic uplift rates over Greenland via a number of GIA models. On average, across our region of investigation, these models suggest uplift rates are approximately 0 mm a −1 to 2 mm a −1 , which 235 equates to a maximum potential mismatch of 18 cm over 90 years. In light of other much larger uncertainties, we do not consider this potential error source to be of significance. More detailed information on the creation of the orthophotomaps can be found in the supplementary materials.

Satellite imagery
In order to complete our record and extend the period covered by our study to the present day, we use satellite records to 240 explore glacier change from 1985 to the present day. Imagery were downloaded from the USGS website EarthExplorer (U.S. Geological Survey USGS, 1995). We sought images with minimal cloud-cover that were gathered in July/August, when we expect surface snow cover to be at a minimum. We sourced imagery from Landsat-5 for the years 1985 and 1995, from Landsat-

Quantifying changing glacier margins
After the image processing was complete, seven discrete time steps of glacier extent were generated, covering a period of 89 years. For each year for which we had data, the margins were delineated manually. Prior to doing this, we investigated the use of semi-automatic margin detection approaches but found that a manual approach was more suitable and more accurate (Paul et al., 2017;Rippin et al., 2020). These delineated margins were then analysed to investigate glacier change over time. Specifically, Our derived orthophotomaps are shown in the third row with the approximate correspondence to source images shown in red.

Additional datasets
In addition to our focus on imagery as outlined above, subsequent investigations and analyses also make use of a range of 260 additional environmental datasets. Here we briefly outline these and their sources.

Surface mass balance (SMB)
Direct measurements of surface mass balance data from this part of Greenland are very rare, and in fact are only available for 275 any significant duration from a single glacier -Mittivakkat Gletscher (Mernild et al., 2011;Bjørk et al., 2012). It is this lack of mass balance data that makes exploring frontal change as a proxy for mass balance so important (Bjørk et al., 2012). Here we utilise ::::: utilize modelled surface mass balance data following the approach of Wake et al. (2009) and Box and Colgan (2013) for South East Greenland and have explored how this varied over our study period. As with SST, mass balance fluctuations are displayed with reference to a baseline defined as the mean SMB over the period 1951 to 1980 ).  Figure 3 shows the overall net glacier terminus change that has taken place over the 89 year period of investigation. This figure gives unique insights into glaciological changes. This is unique to our study because this area is relatively little studied, and also unique because here we extend the record of change back beyond the era of imagery available from the satellite record alone.

290
This gives important insights into glacier extent and change in the pre-satellite era. Figure 3 shows that over this time period all glaciers in our study area experienced overall retreat. Significantly, those glaciers in the southern study region retreated by up to 3 km. In the northern study region, all glaciers again showed retreat. However, there was more variation with some glaciers in this relatively small area experiencing total retreats of several kilometres. Others showed much smaller amounts of retreat. Figures   295 4 and 5 break these frontal changes down into the individual time-steps available to us from our suite of imagery. These data are also summarised in Table 2. Figure 4 shows changes over each time step in the northern region (cf. Figure 1). In the earliest period  all glaciers appear to be retreating. By contrast, in subsequent time periods some glaciers appear to show In the southern region ( Figure 5), there are fewer glaciers to consider. Indeed we do not have any data for the period of 1930-1985for glaciers 20 . :: It :: is ::: also ::::::::: important :: to ::: note :::: that ::: we :: do ::: not :::: have :::: data ::: for ::: the ::::: 1960s ::: for ::: this :::::: region and 21. We are also 305 missing the 1966 time step for glaciers 22-23 wheres we do have data from 1930s and 1985 for those glaciers, and so parts A and B of Figure 5 are blank. As with the northern region, the 1985-1995 timestep was one in which advance of some glacier terminii :::::: termini also took place, while retreat was experienced by others. Moving beyond 1995, as with the northern region the retreat of most glaciers resumed at an increased rate in the most recent period (2015-2019).
In Figure 6, we divide the glaciers in our study regions according to type, and use box plots to visualise the range of variation 310 in response of different glaciers. These plots reveal that there is substantial variability between individual glaciers, which is not unexpected due to the complexity of individual glacier response. We can identify a general trend from retreat during the early parts of our study period through to slight advance in the middle part of our study and then increasing rates of retreat in more recent decades. The marine terminating glaciers appear to be more dynamic, showing greater diversity of change between individual glaciers than land-terminating glaciers. However, both types of glacier show the same overall pattern, as described 315 above. Taking ice-sheet outlets as a whole, we see that period of 1960s-1990s is one characterised by advance. Both before this period and since, retreat has been more dominant. In particular, more retreat is apparent since the turn of the 21st century.
Local glaciers meanwhile, have displayed retreat throughout (with possible equilibrium in the 1990s) and perhaps a slight trend towards increasing rate of retreat in more recent years. Notably, land-terminating glaciers and local GICs that are peripheral to the GrIS (which also tend to terminate on land) show much less variability (i.e. shorter whiskers in the box plots). This may 320 indicate a more consistent response to climatic drivers of change.
It is also worth noting that amongst the wide variability in behaviour, some glaciers demonstrated marked stability throughout the period of our study. The precise details pertaining to this will be dealt with in depth in the discussions.

Sea surface temperatures
Mean annual sea surface temperatures (SST) used in this study were determined from the Hadley Centre Sea Ice and Sea Surface Temperature data set (which is itself taken from the Met Office Marine Data Bank (MDB)). Detailed information on how this data set was created can be found in //www.metoffice.gov.uk/hadobs/hadisst/index.html and Rayner et al. (2003).
In Figure 7b), white circles represent measurements of SST while the red line is a 10 year rolling mean. We compare these 345 values with a baseline (black horizontal line) calculated as the mean SST of the 1951-1980 period . Triangles indicate timing of our glacier front observations from either aerial imagery (purple) or satellite imagery (green).
We are able to explore SST over a period of more than one hundred years. During this period it fluctuates around the baseline (calculated as the mean of the SST over the period 1951-1980; Figure 7b)). In the first time-step  lower than average SST was initially apparent. This is followed by a rise above the mean later on. Immediately prior to this time-step,

350
SSTs were below the mean. In the second time-step , SSTs dropped below the baseline. After 1985, SSTs start to increase and continued to do so, reaching values higher than at any preceding time. SSTs very closely track changes in annual air temperatures (Figure 7a), showing the same broadscale variation.

Surface mass balance
Surface mass balance (SMB) in South East Greenland (Wake et al., 2009;Box and Colgan, 2013) can be seen in 7c). Black 355 circles represent measurements of SMB while the red line is a 10 year rolling mean. We compare these values with a baseline (black horizontal line) calculated as the mean SMB of the 1951-1980 period . Triangles indicate timing of our glacier front observations from either aerial imagery (purple) or satellite imagery (green). There is a significant amount of variability in SMB over the 20th century, reflecting the importance of the 10-year running mean in order to discern longer term patterns in behaviour (Figure 7c)). It is important to note that over the entire 20th century, SMB was positive in this 3.5 Summarising frontal advance/retreat, SMB, SST and annual air temperature changes Our investigation of SST, annual air temperatures and indeed SMB changes over both the 20th century and the early part of the 21st century is in an effort to explore likely drivers for the changes in ice front positions that we observe in the archival imagery.
A glacier's terminus position is controlled by the balance between a) the amount of ice being added to the parent ice mass, b) 370 the flow of this ice towards the terminus, and c) losses at the terminus induced by melt (either atmospheric or marine) as well as potential iceberg calving. We are interested in investigating patterns between these controlling environmental variables. We also make links between these variables and the behaviour of the ice fronts explored in both our northern and southern regions.
Both air temperature and SST show the same broad trends, dividing the 20th and early part of the 21st century that is covered by our investigation into:

375
(i) an early period (up until~1930) of cooler than baseline temperatures; (ii) a sustained period of warmer than baseline temperatures (up to~1965); (iii) a shorter period of cooler temperatures (up to~1990/1995); (iv) a period of warming that continues up to the present day but which shows some flattening in recent years.
Of note, however, is that the SST fluctuation is much more subdued than that in the air temperature. This is particularly the 380 case in period (ii), when SST is only very marginally above the baseline. There is also a lag, such that SST variations are not only subdued but also lag several years behind air temperatures. This is, of course, not surprising, since ocean temperatures rise as a consequence of atmospheric warming or cooling.
Variations in SMB do not track changes in air temperature or SST in a simple or direct way, but clear trends are apparent.
Period (i) is one of fluctuating SMB around the baseline. Period (ii) is one in which SMB becomes increasingly positive before 385 declining again. This is followed by continuing declining SMB and then further fluctuations around the baseline in period (iii) and (iv). Superimposed on this is a higher frequency variability in SMB. A simple assumption that SMB responds directly to time-integrated air and/or SST changes is therefore not apparent, nor would it be expected. The associated complexity is a consequence of how the controls on energy inputs into large ice masses vary on a range of timescales, and also the temporal lag between these inputs and an ice mass responding. It is also a consequence of other controls on ice mass response, such 390 as changing oceanic circulation patterns and geomorphological controls. Explanations for the complexity of the response are considered in detail in the discussions. It is interesting to observe that there are clear and broadscale patterns in the SMB response that could be attributed to variations in air temperature and SST.
We also observe these external forcings playing out in the changing extent of the outlets, but with a degree of complexity possibly reflecting a lag in the response of the ice masses. For instance, the cooler period (i) does not immediately manifest 395 itself as ice front advance during this same period. Rather it is some years later (most notably Figures 4b and c) where we see glacier advance. This is also reflected in the positive SMB in period (ii) which clearly manifests itself as advance of many outlets. Sustained and widespread retreat at a growing rate commences just as period (iv) begins, continuing up to the present day. There is thus general and broad scale manifestation of these changing parameters in both SMB and glacier frontal response.
However, not all outlets respond in the same way, and so these too are considered in the discussions.

Outlet response to regional climatic trends
Our work with archival imagery has enabled us to extend the record of glacier frontal change beyond the limits of the satellite record. We have also been able to do this in a region of Greenland that is relatively poorly studied. We have shown that glacier frontal positions varied over this time period, alongside limited measurements of surface mass balance. These varying 405 glacier extents occur in response to changes in both air temperatures and SSTs, which fluctuate between cool/warm/cool/warm conditions (around our baseline). This suggests that the underlying drivers of these changes are air and ocean temperatures.
In general, in our data we see that the overall (regional) trends in glacier change (as observed in the box plots of Figure 6) do track the prevailing climatic forcing. Greater rates of retreat take place during the warmer period (approximately 1925-1964), with a more subtle slowing of this response during the cooler periods (approximately 1905-1925 and 1964-1996), and isolation, although there is some clear variability, our data also shows that the overall trend is flat or at least subdued. However, considering the contemporary period as a whole, we believe that temperature trends do show an overall increasing trend. This is particularly so in the record of minimum air temperature, which may be significant when considering the role of elevated 415 minimum temperatures on the net amount of melt that takes place. In the contemporary period, we also see warmer seas, as well as a larger increase in positive degree days. There is considerable variability from year to year in the positive degree days during this period, which perhaps reflects the compensating short term warming and cooling events referred to by Hanna et al. (2021).
However, glacier frontal response to these climatic drivers is more complex and time-lagged. Of course not all glaciers respond in the same way, with the same magnitude or at the same rate. This indicates that there are additional controls too. Figure 6 demonstrates this significant heterogeneity. Here we have subdivided the glaciers in our study area according to type.
We see that glaciers of a different type respond differently to external drivers. The marine terminating glaciers in our study 430 region are a) more dynamic, b) show more retreat and c) show more varied behaviour than land terminating glaciers. Such behaviour is well-documented (Moon and Joughin, 2008;Murray et al., 2015), and highlights that the oceans (currents, tides and bathymetry) and SST changes :: (as :::: well :: as :::::::::: subsurface :::::::::: temperature ::::::: changes) : have a vital role in the stability of these ice masses.
Such complexity of response, and variability amongst marine terminating glaciers is also an observation reported recently by 435 Wood et al. (2021), and which is discussed in more detail below. For local glaciers and land terminating glaciers, we observe that whilst in earlier periods, changes in these glaciers were relatively small, larger changes have become more apparent recently. We propose that this 'switch' could be representative of SMB becoming an increasingly important driver of change in recent years, as has been documented elsewhere (cf. Wood et al. (2021)).
The study attempts to categorise glaciers according to their geometry and water-depth. In the area covered by our study, the vast majority of glaciers (in fact all but one) are described by Wood et al. (2021) as 'noncategorized,' which means that the bathymetry and water properties are unknown. Our long-term investigation of these glaciers and the observation of their ap-460 parent stability, suggests these glaciers may sit in relatively shallow water on shallow ridges. This prevents the intrusion of warm deep water which would further enhance mass loss . Although their bathymetry is currently entirely unknown , it is possible that in the future these glaciers may pass a tipping point when they retreat off their pinning ridge into deeper waters. This would see a switch from their current status of having little frontal change, to a phase  1966,1985,1995,2005,2015  with much more rapid retreat. At present this is very much speculative, but ongoing monitoring of these outlets is therefore of 465 great importance.
In addition to these previously discussed outlets in which frontal change is minimal, there are three outlets that in particular show 2-3 times more retreat than these. Two of these are part of Glacier 8 (see Figure 8) which has several outlets. The two southernmost ones showed rapid and large-scale retreat between 1930 and 1966, but then displayed very little change over the years since then. In contrast, Glacier 10 showed relatively modest retreat from 1930 to 2015, but then large-scale and rapid 470 retreat in the four years to 2019 (Figure 8). Figure 9 shows similar behaviour. Here, Glacier 14 appears to be very stable, with minimal fluctuation around its terminus over the duration of the study period, albeit with an overall trend towards modest retreat. However, there is some complexity within this glacier alone. This is due to the fact that the eastern side of this very wide marine-terminating outlet shows more consistent and substantial retreat. Glacier 15 shows much more retreat with several large retreat 'steps', but with the most significant retreat step being between 2015 and 2019.  1966,1985,1995,2005,2015 and 2019. for glaciers 14 and 15 as shown in Figure 1A. (B) Frontal positions are again shown, but with colouration indicating ice velocity as well. Ice velocities derived from Joughin et al. (2010) and the background is the 80s mosaic. Wood et al. (2021) similarly reported that many of Greenland's marine-terminating glaciers have sped-up and lost mass as a consequence of warming ocean waters, but that there are some glaciers that have exhibited small or no retreat. The explanation presented by them for this minimal retreat is that this is a result of water being shallow or outlets resting on shallow ridges. It may well be that this also helps to explain the diversity of behaviour we identify. Many of our study glaciers show little retreat over the 90 year study period and although we do not have bathymetry data, we propose that these understudied glaciers also 480 sit on ridges and/or in shallow water. Where our glaciers have shown periods of more significant retreat for some part of the 90 year investigation, we propose that these periods of change indicate when glaciers move off pinning ridges into deeper water, even though they then may subsequently become grounded again and thus their retreat slows.
With regards to the differing behaviour of the two parts of the front of Glacier 14 (see Figure 9), we are fortunate to have bathymetry data (A. Bjork :: Bj : ø :: rk, personal communication, January 2021) for the region directly abutting the ice front ( Figure   485 9B). This reveals starkly different topography, with that in front of the more stable region being significantly deeper than that in front of the more changeable region. However, on closer inspection of the bathymetry data it is apparent that directly in front of the western part of this glacier, there is a subtle shallowing of the bed. This could suggest the presence of a ridge which pins the glacier and thus explains why it appears to have a stable front. We do see higher ice velocities here (see Figure 9B) and so it is also possible that the apparent relative stability of of the western outlet arises because comparatively high calving rates are 490 offset by higher ice velocities delivering ice more rapidly to the ice front. Conversely, the eastern outlet lies in shallower water but retreat is nevertheless more substantial. Ice velocities are lower here and so calving and/or melting is not countered by ice flow from inland (i.e. lower velocities than in the west ( Figure 9B)).
Finally, Glacier 15 shows significant frontal retreat and high surface velocities. Following the thinking described above, we propose that this suggests that despite the delivery of large amounts of ice to the calving front from inland, significant retreat 495 is still occurring, and so this glacier may be losing the greatest amount of mass overall.

Conclusions
Our investigation has shown the potential of archival imagery that was not originally (and thus not optimally) collected for the purpose of photogrammetric investigations of glacier change. It is thus an important demonstration of the powerful quantitative data that resides in such imagery. This archival imagery has enabled us to extend the record of change of a number of little 520 studied glaciers that reside in the central-eastern part of Greenland (Mouginot et al., 2019), back by several decades beyond the beginnings of the satellite record. Being able to do this is of great benefit since a longer time series of glacier change enables a better understanding of how ice masses have responded to climate to be developed (Dyurgerov and Meier, 2000).
Our focus here has been on a series of outlet glaciers from the Greenland Ice Sheet, and an investigation of how these have varied alongside a number of other controlling environmental parameters for which we also have long-term records. Our study 525 covers~90 years and is the first such dedicated study in this region and over this duration. It deals with changes of a number of previously poorly studied glaciers that have perhaps been largely overlooked. One of our key findings is that climate forcing exhibits strong controls on glaciers in the region generally, and that there is a very close link between air temperatures and SSTs. Arguably, SSTs are more important as we see larger scale significant retreat of outlets terminating in water as the oceans have warmed.

530
However, our study region contains a number of different types of glacier. We observe that it is the marine terminating glaciers that show the greatest mass loss, particularly in the more recent period. Aside from our observation of the importance of climatic forcing, we also highlight significant local variations and the potential importance of non-climate-related factors.
Above all, one of our primary conclusions is that there is enormous variability in how glaciers respond to the climatic and non-climatic drivers. In particular, we propose that the great variability in the retreat of marine terminating glaciers (both in 535 terms of the magnitude and timing of retreat) is strongly ::: may :: be : controlled by the presence or lack of shallow ridges which act to pin glaciers as they retreat. We :: In ::: our :::::::::::: interpretation, :: we : envisage an undulating submarine/subglacial topography which has meant that some glaciers have showed periods of much greater or lesser retreat, and some are apparently stable in their position. Such a situationlends , :: if :::::::: accurate, ::::: would :::: lend : itself to the possibility of future periods of comparatively rapid retreat of glaciers that appear to be stable, and likewise future stabilisation of other glaciers that may currently (or in the past) have 540 shown more significant retreat. Catania et al. (2018) showed similar results ::: also ::::::: provided ::::: such :::::: insights : for western Greenland which gives us :: and ::: so ::: we :::: have : greater confidence in our interpretation here. The novelty of our investigation is not only that we have shown such behaviour in a previously unstudied reagion ::::: region : of Eastern Greenland, but also that our use of archival imagery allows us to identify that such behaviour has been occurring over a longer time period than it has been previously able to show. This helps to demonstrate the rich insights that can be gained from the processing pipeline we demonstrate 545 here. In the past, regional investigations across the Greenland Ice Sheet have been key (e.g. Mouginot et al. (2019); King et al. (2020)). This has been important for exploring broad scale regional behaviour and responses. However, our work here, in which we have focused on glacier-to-glacier heterogeneity, shows that within regions there is great complexity, with even adjacent glaciers behaving very differently. In our efforts to better understand the complexity of the response of the Greenland Ice Sheet to a warming climate, we propose that it is increasingly important to consider the variability between outlet glaciers because 550 of the variation in responses that we have identified here. :: We :::: also ::::::: support, :::: and :::::: further :::::: stress, ::: the :::: need ::: for ::::: much ::::::::: improved ::::::::: knowledge :: of ::::: fjord :::::::: geometry, :: as ::::::: initially :::::: called ::: for :: by ::::::::::::::::: Porter et al. (2018) :::::: because :: of ::: its ::::::: probable :::::::::: importance :: in :::::::::: controlling :: the :::::::::::: heterogeneity :: in :::::: glacier ::::::::: behaviour. : In addition, our work has also highlighted how difficult it is to analyse overall glacier response from investigations of frontal variations alone. An important future direction would be to focus on surface elevation change and also to explore the subglacial topography of these outlets to predict likely future 'jumping' periods of retreat, or 555 indeed stabilisation.
Code availability. There is no specific code to be made available, but we will happily discuss our approach on request.

570
Competing interests. There are no competing interests are present.