University of Birmingham Recent changes in area and thickness of Torngat Mountain glaciers (northern Labrador, Canada)

The Torngat Mountains National Park, northern Labrador, Canada, contains more than 120 small glaciers: the only remaining glaciers in continental northeast North America. These small cirque glaciers exist in a unique topoclimatic setting, experiencing temperate maritime summer conditions yet very cold and dry winters, and may provide insights into the deglaciation dynamics of similar small glaciers in temperate mountain settings. Due to their size and remote location, very little information exists regarding the health of these glaciers. Just a single study has been published on the contemporary glaciology of the Torngat Mountains, focusing on net mass balances from 1981 to 1984. This paper addresses the extent to which glaciologically relevant climate variables have changed in northern Labrador in concert with 20th-century Arctic warming, and how these changes have affected Torngat Mountain glaciers. Field surveys and remote-sensing analyses were used to measure regional glacier area loss of 27 % from 1950 to 2005, substantial rates of ice surface thinning (up to 6 m yr−1) and volume losses at Abraham, Hidden, and Minaret glaciers, between 2005 and 2011. Glacier mass balances appear to be controlled by variations in winter precipitation and, increasingly, by strong summer and autumn atmospheric warming since the early 1990s, though further observations are required to fully understand mass balance sensitivities. This study provides the first comprehensive contemporary assessment of Labrador glaciers and will inform both regional impact assessments and syntheses of global glacier mass balance.


Introduction
The glaciers of the Torngat Mountains, northern Labrador, Canada, occupy a unique physiographic and climatic setting at the southern limit of the eastern Canadian Arctic. Their proximity to the Labrador Current provides temperate, maritime summer conditions yet very cold and dry winters. Examination of Torngat glacier change may provide insights into the deglaciation dynamics of other similarly situated small-glacier populations (e.g. around the south Greenland ice sheet peripheries and Canadian Arctic Archipelago; Nilsson et al., 2015;Gardner et al., 2012) and other small glaciers retreating into cirque basins in temperate high-mountain settings (e.g. European Alps and Pyrenees, North American Rockies and Cascade Mountains). In addition, the glaciers of the Torngats contribute to the total land ice mass of very small glaciers (e.g. Bahr and Radic, 2012), are an important component of the local Arctic tundra and fjord ecosystem (Brown et al., 2012), and form part of the cultural landscape of the Labrador Inuit. Whilst early exploration and reconnaissance survey flights led to classification and mapping of some Labrador glaciers (e.g. Forbes, 1938;Henoch and Stanley, 1968), a complete inventory and analysis of the glaciers of the Torngat Mountains, covering a total active glacier area in 2005 of 22.5 ± 1.8 km 2 (not including relict ice), has only recently been completed and incorporated into global glacier datasets (Way et al., 2014;Pfeffer et al., 2014;Arendt et al., 2015). Combining the new regional glacier inventory (Way et al., 2014) with field-and remote-sensing-based analyses, Way et al. (2015) showed a consistent glacial response to centennial-scale regional climate warming by documenting a 52.5 % reduction in ice extent since glaciers reached their Little Ice Age (LIA) maxima (dated to between 1581 and 1673).
The only other work published on the contemporary glaciology of the Torngat Mountains was based on field research expeditions conducted by Memorial University, Canada, during the early 1980s (Rogerson, 1986;Rogerson et al., 1986b). Between 1981 and1984, these authors measured net mass balances using the glaciological method, from centreline and distributed transects of ablation stakes and snow pits, on clean and debris-covered ice at four cirque glaciers in the Selamiut Range, south of Nachvak Fjord, in the central part of the national park (Fig. 1). In all years between 1981 and 1984, the Abraham and Hidden glaciers in the adjacent McCornick River valley experienced negative net balances, ranging from −0.21 to −1.28 m yr −1 water equivalent (w.e.). Nearby higher-elevation Minaret Glacier and Superguksoak Glacier (at 1.37 km 2 , the largest glacier in the region) had near-zero balances in 1981, positive balances in 1982 (both 0.28 m yr −1 w.e.), and slightly negative balances (both −0.16 m yr −1 w.e.) in 1983 (uncertainties are unreported). These results were linked to annual variations in winter precipitation (Rogerson, 1986). Aside from the 2005 inventory (Way et al., 2014), no more recent information is available regarding the health of these four glaciers or of the remaining glacial ice in Labrador. A key question is the extent to which glaciologically relevant climate variables in northern Labrador have changed in concert with Arctic warming during the 20th century (Kaufman et al., 2009), and how these changes have affected Labrador glaciers.
This study uses field, remote-sensing, and homogenized climate and reanalysis datasets to examine the response of glaciers in northern Labrador to prevailing climatic conditions since the mid-20th century. The 2005-dated inventory of Way et al. (2014) (derived from very high resolution aerial orthophotos) is used as a basis for historical change assessment since the mid-20th century by comparing to glacier outlines derived from archival aerial survey photography (1950). Recent glacier changes are calculated using glacier areas derived from high-resolution spaceborne electro-optical imagery (2008). Regional volume and volume change estimates are calculated using volume-area scaling relations. Field and remote-sensing surveys of ice surface topography at three of the four glaciers visited in the early 1980s are used to calculate rates of surface elevation change between 2005 and 2011, and to calculate geodetic mass balances. Finally, climate and reanalysis datasets are selected to view both 20thcentury area changes and 1980s and late-2000s glacier mass balances in the context of prevailing climatic trends. This study builds on earlier work (Brown et al., 2012;Way et al., 2014;Way et al., 2015) and provides the first comprehensive, regional-scale remote-sensing and field-based assessment of the contemporary state of Labrador glaciers and their likely climate sensitivities, the results of which may be incorpo-rated into both regional impact assessments and syntheses of global glacier mass balance.

Study area and methods
The Torngat Mountains are the southernmost mountain range in the eastern Canadian Arctic, occupying a position towards the north of the Labrador Peninsula, south of Baffin Island (Fig. 1). The synoptic climatology of the region is strongly influenced by the Labrador Current, which transports cold Arctic water into the Labrador Sea and southward along the Labrador coast (Myers and Donnelly, 2008). The influence of the Labrador Current produces maritime summer conditions and colder, drier winters due to extensive seasonal sea-ice formation (e.g. Kvamstø et al., 2004). Precipitation during winter is generated by storms that are directed by the influence of the Canadian Polar Trough on the position of Arctic and polar fronts (Moore et al., 2011). Glaciers in the Torngat Mountains are preserved due to their proximity to this coastal moisture source, as well as their topographic setting (reducing incoming solar radiation receipt) and, in some cases, the insulating effect of supraglacial debris (Rogerson, 1986;Rogerson et al., 1986b;Way et al., 2014). The four glaciers measured in the early 1980s are influenced by a combination of these factors. Abraham and Hidden glaciers ( Fig. 2; labelled A and H) are less than 20 km from the coast, have partly debris-covered termini, and are located below ∼ 500 m high backwalls comprising Cirque Mountain (1568 m above sea level (a.s.l.)). Superguksoak Glacier (Fig. 2, labelled S) is similarly situated and is debris-covered over ∼ 40 % of its surface area. Minaret Glacier (Fig. 2, labelled M) does not have extensive supraglacial debris cover but is located ∼ 350 m higher than the other glaciers and is shaded by a large backwall comprising Mount Caubvick (1652 m a.s.l.). The interplay of altitude, supraglacial debris, slope, aspect, and height of (and proximity to) cirque backwalls has been shown to influence the spatial variability of surface mass balance (Rogerson, 1986).

Derivation of area changes
A glacier inventory derived from late-summer (August) 2005 Parks Canada 1 : 40 000-scale colour aerial photographs was used as a reference for documenting historical and late-2000s glacier area changes. Full details of the 2005 inventory data collection, image processing, ice area delineation, and error assessment are provided by Way et al. (2014). Historical glacier outlines were manually digitized from 107 scanned diapositives of Royal Canadian Air Force survey photographs, acquired during a 2-week period in August 1950 (Labrador (LAB) series 47,48,50,79,[91][92][93][94][95]. Oneto-fifty-thousand-scale diapositives scanned at 1200 dpi provided digital files with a corresponding ground pixel resolution of ∼ 1 m. Image orthorectification in the absence of accurate camera calibration information was undertaken using a rational functions model in PCI Geomatica software (v.10.3). The model calculates correlations between image pixels and ground locations using a ratio of polynomial functions with coefficients calculated based on the number and quality of ground control points (GCPs). Between 15 and 20 GCPs per image frame were collected from 2005 orthophotos and an undated 18 m resolution Parks Canada digital elevation model (DEM), to correct for topographic distortion due to terrain relief (e.g. Kääb, 2005). The resulting root mean square (rms) error between the 1950 orthophotos and GCPs was < 5 m across all images.
To examine mid-to late-2000s area changes, scenes from the SPOT5-HRG (Satellite Pour l'Observation de la Terre High-Resolution Geometric) image archive were analysed ( Table 1). The very small size of Labrador glaciers and the need to obtain both cloud and snow-free scenes during summer precluded the use of medium-resolution spaceborne visible and near-infrared sensor archives such as ASTER (Advanced Spaceborne Thermal Emission and Reflection Ra- diometer) and Landsat. SPOT5-HRG images at 10 m ground resolution in multispectral mode (three bands between 0.50 and 0.89 µm) and 5 m in monospectral (panchromatic) mode (0.51-0.73 µm), providing sufficient spatial resolution to accurately delineate Labrador glacier outlines following pan sharpening (e.g. Paul and Kääb, 2005;. Cloud-and snow-free scenes used in this work were acquired during late July and early August 2008 (Table 1). Glacier extents were delineated by overlaying 2005 outlines onto 1950 and 2008 imagery sets and identifying where detectable changes had occurred (e.g. Fig. 2). Due to missing or poor-quality 1950 photo frames, outlines of 101 glaciers were digitized (out of a total of 124 active glaciers in the 2005 inventory). The areas of these 101 glaciers represent 81 % of the 2005 total ice area. Due to the large swath coverage of SPOT5-HRG scenes (60 × 60 km, continuous tiles), all 124 active glaciers in the 2005 inventory were successfully identified in the 2008 imagery and redigitized. Glacier outline adjustment (mapping) took place following the protocols of the Global Land Ice Measurements from Space (GLIMS) initiative (Racoviteanu et al., 2009). Indistinct and debris-covered ice margins were mapped using morphological features such as lateral meltwater streams, debris flow bands, and changes in surface slope (facilitated by three-dimensional visualization of orthophotos draped over the Parks Canada regional DEM) (e.g. Paul and Kääb, 2005;Racoviteanu et al., 2009;Way et al., 2014).
Following the approach of Paul et al. (2013), the error (precision) of glacier area delineations was determined by independent remapping of ice margins. Three separate operators mapped the extent of 10 of the largest glaciers in the region (∼ 15 % of the total ice area), selected to include a representative range of margin characteristics (bare ice, snow-covered, debris-covered, shadowed, and unshadowed). The rms difference between operator extents (0.009 and 0.019 km 2 for 1950 and 2008 imagery) corresponds to mapping errors of ∼ 2 and ∼ 5 m (equivalent buffer widths) along respective glacier margins. Buffers of ±5 and ±10 m were subsequently assigned at 1950 and 2008 ice margins to account for cumulative errors in mapping, image planimetry, and orthorectification (e.g. Bolch et al., 2010;Way et al., 2014).

Regional ice volume estimation
Volumes of individual glaciers in 1950, 2005, and 2008 were calculated and summed regionally, using the volume-area (V -A) power law relation: V = c a A γ (e.g. Chen and Ohmura, 1990;Bahr et al., 1997). Aggregate volumes were calculated for ensembles of glaciers, as errors can be large when reported at the individual glacier scale (e.g. Chen and Ohmura, 1990). Bahr et al. (2015) note that a 1σ error around a mean of the multiplicative scaling parameter c a would produce a 34 % error in calculated ice volume. However, volume errors calculated by V -A scaling for regional-scale aggregate ensembles of glaciers are typically < 25 % (Meier et al., 2007). In the absence of independent measurements of glacier volume from which to derive local scaling coefficients, c a was assigned equal to 0.191 m (3−2γ ) , and the scaling exponent was fixed to the theoretical constant γ = 1.375. These values are based on worldwide means of scaling parameters and a theoretical derivation of the scaling exponent (Bahr et al., 1997(Bahr et al., , 2015. To account for variations in the scaling parameters due to region, glacier, or other factors (spatial or temporal), uncertainties in volume (errors) were estimated by calculating volume sensitivities to the choice of scaling coefficients (e.g. DeBeer and Sharp, 2007; . These coefficients were derived empirically from direct volume measurements of glaciers in the Alps, Cascades, and similar regions (Chen and Ohmura, 1990); of small glaciers in the Canadian Cordillera (DeBeer and Sharp, 2007); and from purely physical considerations (Bahr et al., 1997). Total volumes of all glaciers were calculated with the scaling parameters (units m (3−2γ ) ) / exponents (dimensionless): 0.206/1.360, 0.115/1.405, and 0.210/1.360. The maximum difference between calculated volumes was then taken as the measure of uncertainty (error) in total volume calculations (e.g. DeBeer and Sharp, 2007; . This approach is considered preferable to a component error propagation calculation, as uncertainty in volume calculations due to the choice of scaling coefficients is likely to be an order of magnitude larger than that from the principal individual error sources (e.g. area measurement errors).

Ice surface elevation changes
Field surveys of glacier centreline ice surface topography were undertaken at Abraham and Hidden glaciers during August 2008, and at Minaret Glacier during August 2008, using dual-frequency (L1 and L2) Global Positioning System (GPS) instruments. Base-station data were collected at 1 Hz using tripod-mounted Leica GPS 500 receivers fixed over two survey markers drilled into exposed bedrock in the forefields of Hidden/Abraham and Minaret glaciers (for locations see Fig. 2). Ice surface topography was measured at accessible points along glacier centrelines at ∼ 5 m horizontal spacing using an additional GPS 500 receiver in rover mode connected to a survey-polemounted antenna. Simultaneous base-station positions were recorded during surveying to provide differential processing capability. Rover survey points were differentially corrected (including for pole height) using post-processing tools in Leica GeoOffice software. Following differential correction, positions with three-dimensional accuracy < 10 cm were excluded from further analysis. This procedure resulted in elimination of ∼ 15 % of survey points at higher elevations proximal to cirque backwalls due to the poorer configuration of available satellites during surveying in those areas.
Elevation changes (δh) were calculated by differencing individual crossover spot heights (h) within a 1 m radius of the preceding year's survey position. Corrections for along-and across-track slope were deemed unnecessary due to the small crossover window size and shallow surface slopes of the accessible (and thus surveyed) parts of the ice surface (typically < 5 • ). Whilst some surface elevations were excluded from the difference analysis due to positional deviations from the centreline track, sufficient crossovers remained to adequately characterize centreline elevation change without the need for interpolation or plane fitting (e.g. . Errors in individual crossover elevation changes were determined empirically by calculating the root sum of squares (RSS) of errors in receiver position. In addition, three-dimensional GCPs collected in 2008 over bare rock surfaces around each glacier were used to re-process 10 (Abraham and Hidden) and 5 m resolution (Minaret) DEMs, from 2005 Parks Canada air photos. Ice surface elevation changes were then calculated between 2005 DEM cell pixels and 2008 GPS survey positions.

Geodetic mass balance
Centreline ice surface elevation changes were extrapolated using glacier hyposometies to calculate geodetic mass balances from 2005 to 2008, from 2008 to 2009 (all), and from 2009 to 2011 (Abraham and Hidden glaciers only). Hypsometries were extracted from 10 m (Hidden and Abraham) and 5 m (Minaret) -posting digital elevation models of each glacier surface. DEMs were derived from stereophotogrammetric processing of 2005 aerial survey images (Way et al., 2014). Individual centreline elevation changes (δh) were averaged into 10 m elevation bins and applied to all pixels in each bin across the glacier surface (e.g. Nuth et al., 2010). In the event that no δh measurements were available (due to safe-access issues during surveying or culling of low-accuracy points; see Sect. 2.3), that bin was assigned the averaged value of the next adjacent bin. Total volume changes were then calculated by multiplying δh by the pixel area and summing across the glacier surface (e.g. Barrand et al., 2009;Nuth et al., 2010). Volume changes were converted to geodetic mass balances in water equivalent units (m yr −1 w.e.) by dividing by the glacier area, assuming an ice density of 850 ± 60 kg m −3 (Huss, 2013), and in the case of 2005-2008 and 2009-2011 measurements, dividing by the time between epochs.
Uncertainties in volume change and mass balance were calculated following an approach similar to Nuth et al. (2010). Individual point elevation change errors (E PT(δh) ) were calculated empirically from errors in DEM vertical accuracy or GPS receiver position (Sect. 2.3). Total elevation change errors, E Z , were then calculated by combining E PT(δh) in quadrature with extrapolation errors (E EXT ), reducing by the square root of the total number of independent measurements (e.g. Arendt et al., 2006;Nuth et al., 2010), considered to be one measurement per 10 m elevation bin, given the presence of weakly positive spatial autocorrelation. E EXT quantifies the uncertainty in applying centreline δh measurements to entire elevation bins. Following Nuth et al. (2010), E EXT was approximated by the standard deviation of the averaged δh in each elevation bin. Where bins had fewer than five δh measurements, E EXT was set to twice the glacier average, while bins with extrapolated measurements were assigned 3 times the average. Total volume change errors, E VOL , were calculated as the RSS of E Z multiplied by the 2005 glacier area. As all elevation measurement surveys were conducted in the same month of each year, no additional uncertainties were assigned as a result of seasonal changes.

Results and discussion
3.1 Regional area and volume changes Between these two most recent years, uncertainties from the maximum difference between scaling coefficients (Sect. 2.2) were larger than scaled volumes as a result of the relatively small changes in glacier area. Using the average of the old and new glacier areas, the total volume change between 1950 and 2005 corresponds to an area-averaged long-term thinning rate of −0.16 ± 0.05 m yr −1 (from V -A scaling). The scaled ice volumes from 124 glaciers of the complete inventories of 2005 and 2008 were respectively 0.48 ± 0.06 km 3 and 0.47 ± 0.05 km 3 .

Selamiut Range/Cirque Mountain glaciers
Labrador's largest and only previously studied glaciers are situated in the Selamiut Range and proximal to Cirque Mountain, south of Nachvak Fjord (Figs. 1, 2). Focusing on the four sites visited in the early 1980s, Superguksoak Glacier (the largest glacier in Labrador) decreased in size from 1.47 ± 0.08 km 2 in 1950 to 1.37 ± 0.03 km 2 in 2005 (−7 %) and then further to 1.35 ± 0.11 km 2 in  (Fig. 2). Although no area increases were observed, it is important to note that individual year measurements may mask considerable interannual variability in glacier area change (including the possibility of glacier growth resulting from positive-mass-balance years) (e.g. Rogerson, 1986).
Hypsometries derived from 2005 stereo-photogrammetric DEMs reveal the area-altitude distributions of Abraham, Hidden, and Minaret glaciers (Fig. 4). Abraham and Hidden have bimodal hypsometric distributions, with peaks at ∼ 800 and 900 m a.s.l. at Abraham, and ∼ 875 and 975 m a.s.l. at Hidden Glacier, and no sustained accumulation zone. In contrast, all of the ice of Minaret Glacier exists at higher elevations than Hidden and Abraham, predominantly between 1250 and 1350 m a.s.l. Rates of ice surface elevation change varied between glaciers, and particularly between measurement intervals (Fig. 4). From 2005 to 2008, Abraham, Hidden, and Minaret glaciers thinned at average rates of −0.98, −2.32, and −0.93 m yr −1 , respectively. Between 2008 and 2009, the largest thinning rates were measured at Abraham Glacier (up to 6 m yr −1 at 850 m a.s.l.). These very large thinning rates (compared to resultant area changes) decreased further up-glacier yet remained around 1 m yr −1 at higher elevations. Over the same time period, Hidden Glacier, which has a similar area-altitude distribution to Abraham yet is more extensively shaded and debris-covered, experienced slightly lower thinning rates, up to 4 m yr −1 at lower elevations and 1-2 m yr −1 above ∼ 860 m a.s.l. This elevation-dependant control on valley glacier elevation change was also observed between 2005 and 2008 at all glaciers, and between 2008 and 2009 at Minaret Glacier, with substantial thinning (∼ 4 m yr −1 ) again at lower elevations (1200 m a.s.l.), decreasing to 1-2 m yr −1 above ∼ 1250 m a.s.l. Thinning rates at Abraham and Hidden glaciers were greatly reduced during the period 2009-2011. While relatively few measurements of elevation change were available at Abraham Glacier, they provided an average thinning rate of 0.47 ± 0.22 m yr −1 . An elevation-dependent thinning rate was observed at Hidden Glacier between 2009 and 2011, decreasing from 1 m yr −1 at lower elevations to 0.2 m yr −1 at 1000 m a.s.l. (Fig. 4). The large thinning rates measured at Abraham, Hidden, and Minaret glaciers between 2008 and 2011 far exceed the area-averaged longterm    (2005-2008 and 2008-2009 only). Considering uncertainties, small mass balance values from 2010 to 2011 may be considered as close to zero rather than negative. The larger uncertainties associated with Abraham and Minaret Glacier balances reflect the relative paucity of available elevation change measurements at these sites, largely a result of safe-access issues during surveying at lower elevations and poor GPS satellite coverage within the upper cirque basins.
While geodetic balances are available for years between 2005 and 2011 only, they may be compared with net balances from field measurements obtained in the early 1980s. Net balances at Abraham Glacier in the balance years 1981, 1982, and 1983 were −0.27, −0.46, and −1.28 m yr −1 w.e. (Rogerson, 1986). The geodetic balance measured between 2008 and 2009 (−1.72 ± 0.41 m yr −1 w.e.) is thus the most negative mass balance yet recorded at this glacier. The 2005-2008 geodetic balance at neighbouring Hidden Glacier (−1.83 ± 0.19 m yr −1 w.e.) is more negative than any balance year during the early 1980s (1981, 1982, and 1983: −0.24, −0.21, and −0.81 m yr −1 w.e., respectively) and is the most negative balance yet recorded at this or any other glacier in Labrador. This strongly negative post-2005 mass balance is similar in magnitude to that measured at nearby Terra Nivea ice cap on southern Baffin Island (Papasodoro et al., 2015). While these 6 years of discontinuous measurements represent the only geodetic mass balance observations yet recorded at Labrador glaciers, it is important to note that they may mask considerable interannual variability in accumulation and ablation at the sampled glaciers and within the wider regional population.

Climatological drivers
Two gridded observation-based climate datasets and two reanalysis products were selected to provide a northern Labrador climatological context for the late-20th-and 21stcentury glaciological observations. Seasonal 2 m air temperatures and precipitation totals (where available) were extracted from the homogenized observational datasets Can-GRID and Berkeley Earth (Zhang et al., 2010;Rohde et al., 2013), and from the NCEP/NCAR R1 and JRA-55 reanalysis products (Kalnay et al., 1996;Kobayashi et al., 2015). An intercomparison of these observational data and reanalyses is worthwhile as both approaches and products have data coverage over northern Labrador yet utilize different input data and interpolation schemes. They may also vary in time period, domain, and spatial resolution. A more detailed overview of widely used surface-based and reanalysis datasets for temperature and precipitation trend analysis in the wider Canadian Arctic is provided by Rapaić et al. (2015). Following Rapaić et al. (2015), Berkeley Earth, NCEP R1, and JRA-55 temperature anomalies were evaluated for the Torngat Mountains region (58-60 • N, 295-297 • E) using the CanGRID dataset as a reference. The strongest regression model fit was found between CanGRID and observation-based Berkeley Earth temperature anomalies (r 2 = 0.96, Fig. 5a). The best model fit between Can-GRID and reanalysis temperature anomalies was JRA-55 (r 2 = 0.89), with noticeably less spread and a better model fit than NCEP R1 (r 2 = 0.73, Fig. 5b, c). A time-series plot of summer (JJA) 2 m air temperature anomalies in the Torngat Mountains from all four products confirmed the 20thcentury Arctic summer warming trend (mean of 0.85 • C between 1990 and 2011) (e.g. Rapaić et al., 2015) yet showed a clear cold-temperature bias in NCEP R1 since the mid-1990s (Fig. 5d). Due to this bias and their stronger regression model fits to CanGRID reference data, Berkeley Earth (global availability) and JRA-55 (including precipitation) were selected for further examination of seasonal climate indices.
Time-series plots of seasonal 2 m air temperatures (JJA, SON, DJF, MAM) and total accumulation season precipitation (SON+DJF+MAM) were created to examine seasonal climatic conditions in the Torngat Mountains for the period associated with glacier change measurements . Seasonal temperatures were most similar between Berkeley Earth and JRA-55 datasets in autumn (Fig. 6b) and showed similar interannual variability yet 1-2 • C magnitude temperature difference for summer, winter, and spring seasons (Fig. 6a, c-d). Each of the time series showed long-term warming trends, superimposed with interannual and pentadscale variability, consistent with observations throughout the rest of the province of Labrador (Finnis and Bell, 2015;Way and Viau, 2015). Warming trends were most pronounced during summer and autumn after 1990 ( Fig. 6a-b). Winter temperatures decreased between 1950 and 1990, before rising again post-1990 (Fig. 6c). The post-1990 warming may in part explain individual years of negative mass balance in the late 2000s. Total ablation season precipitation from JRA-55 reanalysis shows 8 years of above-average precipitation between 1974 and 1986, and then a slightly increasing precipitation trend since the early 1990s (Fig. 6e). Year-to-year variability in climate-mass balance linkages are presented in Fig. 7.
While these climate data are both temporally and spatially averaged and thus may not be representative of climate conditions at individual glaciers, it is possible to examine seasonal temperature and precipitation trends in order to suggest potential drivers of glacier change. Glaciers in Labrador lost 27 % of their surface area between 1950 and 2005. As both Berkeley Earth and JRA-55 show largely stable seasonal temperatures between 1950 and the early 1990s, it may be suggested that a majority of this area loss occurred after 1990. Pentadal means of seasonal air temperatures from both records show summer and autumn warming of around 2 • C by 2010-2015, compared to pre-1995 levels (Fig. 6a,  b). This summer warming suggests greater energy available for ablation season melt, while warmer autumn temperatures (which include several seasonal average temperatures above 0 • C) may also indicate an increasing proportion of winter precipitation falling as rain.
The field mass balance measurements of Rogerson (1986), which included years of positive mass balance, were explained by the author in terms of variability in winter snowfall -a finding supported by the precipitation data presented here. JRA-55 precipitation totals during 1976JRA-55 precipitation totals during -1985JRA-55 precipitation totals during pentads were ∼ 8 and 12 % greater (in 1976JRA-55 precipitation totals during -1980JRA-55 precipitation totals during and 1981JRA-55 precipitation totals during -1985, respectively) than the preceding pentadal period (1971)(1972)(1973)(1974)(1975) (Fig. 6e). The only year of negative mass balance at all four glaciers (1983, mean net balance of −0.60 m yr −1 w.e.) coincided with a very low precipitation total during the accumulation season (462 mm).  glaciers between 2005 and 2008, following a high precipitation total of 726 mm in 2006. Larger negative balances at Hidden glacier during this time period may suggest greater insensitivity to precipitation. Geodetic balances more negative than −1 m yr −1 w.e. at all glaciers between 2008 and 2009 followed a low precipitation total of 452 mm during the 2007-2008 winter (which may explain the 2008 snowfree summer ice conditions), and summer temperatures more than 1 • C warmer than the 2005-2010 pentadal mean (itself 1 • C higher than the 1950-2015 mean) (Fig. 6). The lower thinning rates and less negative geodetic balances at Abraham and Hidden glaciers between 2009 and 2011 follow precipitation totals some 17 and 22 % greater than the preceding year (2008) accumulation season total (Fig. 6e), and an extended ablation season with melting conditions extending into the September-November period during 2010 (Fig. 6b). While only tentative climate-mass balance linkages may be drawn from these limited data, interannual variability in winter precipitation and post-1995 climate warming may play important roles in contemporary Labrador glacier mass changes (similar to many other Northern Hemisphere mountain glacier systems). Longer mass balance records with improved temporal resolution and more local-scale climatological observations are required to confirm these findings.

Summary and conclusions
This study presents the first comprehensive analysis of the contemporary state of small glaciers in the remote Torngat Mountains of northern Labrador, Canada. These glaciers occupy a unique position as the most southerly ice masses in the eastern Canadian Arctic. Very little was previously known about the recent status of Labrador glaciers, with a single study published on the mass balance of four glaciers during the early 1980s. Due to their very small sizes and remote locations, analysis of Labrador glacier change required an experimental design comprising both detailed field observations and high-resolution remote-sensing image analysis. Glaciers in Labrador lost 27 % of their surface area between 1950 and 2005, a larger proportional loss than that from similarly sized glaciers in nearby southern Baffin Island over the same period. They then shrunk a further 3 % to 21.80 ± 2.2 km 2 by 2008, with an associated (scaled) total ice volume of 0.47 ± 0.05 km 3 . Three cirque glaciers measured in the early 1980s and revisited from 2008 to 2011 thinned at rates between 0.5 and 6 m yr −1 at lower elevations at measurement periods between 2005 and 2011, mostly in excess of the long-term  thinning rate derived from volume-area scaling (−0.16 ± 0.05 m yr −1 ). These very high thinning rates resulted in consistently negative geodetic balances at all three glaciers between 2005 and 2011, with a balance of −1.83 ± 0.19 m yr −1 w.e. at Hidden Glacier between 2005 and 2008 being the most negative mass balance yet recorded at this or any other glacier in Labrador, and similar in magnitude to that measured at nearby Terra Nivea ice cap. Both geodetic balances between 2005 and 2011 and net balances from early-1980s field measurements were linked to prevailing climatic conditions from gridded observational and reanalysis-based climate datasets. These findings suggested that Labrador glacier mass trends may be controlled by variability in winter precipitation, and increasingly by strong summer and autumn atmospheric warming since the early 1990s, though further observations are required to confirm these linkages.

Data availability
Datasets are available from the corresponding author on request.