Can shrubs help to reconstruct historical glacier retreats?

In the 21st century, most of the world’s glaciers are expected to retreat due to further global warming. The range of this predicted retreat varies widely as a result of uncertainties in climate and glacier models. To calibrate and validate glacier models, past records of glacier mass balance are necessary, which often only span several decades. Long-term reconstructions of glacier mass balance could increase the precision of glacier models by providing the required calibration data. Here we show the possibility of applying shrub growth increments as an on-site proxy for glacier summer mass balance, exemplified by Salix shrubs in Finse, Norway. We further discuss the challenges which this method needs to meet and address the high potential of shrub growth increments for reconstructing glacier summer mass balance in remote areas.


Introduction
In terms of climate change research, glacier dynamics have experienced a strong focus throughout the past decade (e.g. Nordli et al 2005, Oerlemans 2005, Rosenzweig et al 2008. This is because most of the glaciers world wide have been retreating since the end of the climatic depression in the 18th century, known as the Little Ice Age (e.g. Overpeck et al 1997, Maisch 2000, Mann 2002. For the 21st century, further glacier retreat, resulting from global warming, is predicted (e.g. Maisch 2000, Hall andFagre 2003). Yet, glacier dynamic models are getting increasingly complex and long-term glacier mass balance records are needed for the validation of such models (Marshall 2005).
However, most of available glacier mass balance records only span several decades and show data gaps, due to the location of most glaciers in remote Arctic or Alpine areas (e.g. Oerlemans 2005, Kjøllmoen 2011). Proxies are Content from this work may be used under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. frequently used to close such gaps and prolong these records (e.g. Holzhauser et al 2005, Loso andDoak 2006). The two main components of a glacier's mass balance are winter mass balance-when a glacier accumulates mass due to snowfall-and summer mass balance-when a glacier loses mass due to melting (e.g. Nordli, 2005). The age of terminal moraines can be dated through lichenological studies (Beschel 1973, Bradwell 2001, Reyes et al 2006, Bradwell 2007 and winter mass balances can be reconstructed by ice-core analyses or by ground penetrating radar sounding measurements above the equilibrium line (Kohler et al 1997). Past glacial summer mass balance (GSM) has earlier been derived from tree-rings growing in the vicinity of the glaciers (e.g. Watson and Luckman 2004, Laroque and Smith 2005, Wood et al 2011. One advantage of using tree-rings for the approximation of GSM is their annual resolution. Trees are, however, often absent near glaciers, in particular in the Arctic (Weijers 2012). Shrubs are, in contrast, known to grow far beyond the tree-line and could serve as a substitute for trees. Several studies have shown that shrubs act as good proxies for the regional climate, in particular summer temperatures  summer temperature is the main driver of GSM Olesen 1989, Giesen andOerlemans 2009), it therefore seems likely that shrubs may act as on-site proxy for the reconstruction of GSM. The aims of our study were to investigate whether shrubs can serve as a proxy for past GSM and to discuss the challenges and requirements of this approach.

Material and methods
The study site is located in the Scandic Mountains of southern Norway at the south facing slope of the mountain 'Lille Finsenuten', roughly 3 km north of the Hardangerjøkulen icecap (lat: 60 • 35 59 N; long: 7 • 28 32 E; figure 1). The site was chosen for the length of the mass balance record of this local icecap (1963-now) and the availability of several shrub species close to the glacier.
In total, 24 shrub samples were taken near the local shrub line, between 1310 and 1360 m above sea level (asl), 12 of each of two investigated shrub species (Salix lapponum L. and S. glauca L.). We chose Salix shrubs as they were the dominant shrub species growing at the site and are abundant throughout the Arctic (which would lead to a widespread applicability of the method, if proven successful). The sampled shrubs were harvested in August 2011 and treated according to the serial sectioning technique (Kolishchuk 1990, Hallinger et al 2010. The basal shrub sections (i.e. the above-ground stem section nearest to the roots) were prepared as colored microtome thin sections and ring widths were measured along two radii at a resolution of 0.01 mm using a ring-width measuring stage (LINTABLE 5, RINNTECH, Heidelberg, Germany). For two individuals (one for each species) all sampled stem discs (i.e. five and seven, respectively) were measured to locate possible missing outer rings which enables correct dating of ring-width series. Radial measurements were cross-dated visually and averaged using TsapWin professional (version 0.59, RINNTECH, Heidelberg, Germany). Average Gleichlaeufigkeit (glk; i.e. a measure of the similarity between two chronologies), average inter-serial correlation (rbar, i.e. average of all possible correlation coefficients), expressed population signal (EPS, i.e. a probability, that all shrubs contributing to a master chronology may represent the total population), and mean sensitivity (i.e. a measure of the year-to-year variability of a chronology) were calculated (see also Schweingruber (1988) for further details). In addition, running EPS (rEPS) values over a moving window spanning 15 yr were computed for each single year. We did that to highlight the periods of the respective chronology which express an EPS above 0.85, a value that is considered as a threshold of acceptable statistical quality for dendroclimatological analyses (Wigley et al 1984). Three master chronologies (i.e. the average ring-width series, derived from all series included) were built, one for each of the two species and one including both species with prior detrending of individual series using a cubic smoothing spline. This detrending method was chosen in order to remove the low-frequency variability of radial growth that may be the result of biological or stand effects (dplR package in 'R'; Bunn 2008Bunn , 2012. Age related trends which could have justified the application of a negative exponential detrending were not observed in the raw data. For further explanation of the dendrochronological methodologies, we refer to Cook and Kairiukstis (1990).
Monthly precipitation and temperature data (1961-2012) of the climate station 'Finsevatnet' (1210 m asl; figure 1) was retrieved from the Norwegian Meteorological Office (http:// eklima.met.no). This data series is, however, not continuous: it has gaps between 1994 and 2002. We therefore estimated monthly climate data for this period through linear regression between the Finsevatnet and Geilo (841 m asl) climate stations for temperature (r 2 = 0.98, p 0.001) and the Finsevatnet and Eikemo (178 m asl) climate stations for precipitation (r 2 = 0.75, p 0.001). Finally, on-site monthly mean temperatures were approximated by subtracting 0.65 K from the 'Finsevatnet' temperature for each 100 m of altitudinal difference.
Glacial net-, winter-, and summer mass balance of the Hardangerjøkulen icecap between 1963 and 2010 were extracted from Elvehøj (2011). There, annual net balances had been obtained by application of the 'stratigraphic method' (Østrem and Brugman 1991), which determines the altitudinal difference of glacier surface between two consecutive years. Winter and summer mass balances had been determined separately by stake measurements (Kjøllmoen 2011).
The detrended ring-width series (RWI), glacier mass balance series and monthly climate data series were tested for correlations (Pearson correlation) between each other. In addition, Gleichlaeufigkeit (glk) was calculated among these variables. To verify the reconstruction of GSM from RWI, two models calculating GSM from RWI were calibrated for the periods from 1975-1992 and 1993-2010. These periods were chosen, as rEPS revealed to be higher than 0.85 for the period from 1975-2010. The resulting model parameters from the respective calibration period were subsequently used to reconstruct GSM of the respective other (verification) period from RWI. Finally, the complete reconstructed and recorded GSM series combined from the two verification periods were compared to each other visually and statistically using the squared Pearsons' correlation coefficient (r 2 ), reduction of error (RE), sign test, and product means (PM) test (for detailed descriptions on the calibration/verification procedure we refer to Fritts (1976) and Fritts et al (1990)). Further, correlation coefficients between RWI and summer temperatures and among RWI and GSM were computed over different periods with different sample depths (i.e. the number of samples included in the RWI varies throughout the time series as shrubs were not even aged) and then these period and sample dependent correlation coefficients were tested for Spearman's rank correlation (accounting for non-normal distribution of the data, determined with Shapiro-test) with the number of samples representing this particular period. We did so to analyze the influence of sample depth on the model error. This is of higher interest, if the described method is applied to reconstruct GSM beyond instrumental records, as the sample depth of the chronology can often be regarded as a measure of reliability. All statistical analyses were carried out in 'R' (version 2.12.0, R Foundation for Statistical Computing, Vienna).

Results
GSM showed strong and significant negative correlations with average June, July, August and September temperatures (r = −0.49, −0.40, −0.66 and −0.40 respectively, all p-values < 0.01) and the strongest and most significant correlation (r = −0.81, p 0.001) if these were averaged to 'melting season temperatures' (i.e. average temperature from June to September; figure 2). Glacial winter mass balance showed the strongest correlation with winter precipitation sums (i.e. precipitation sum from January to March; r = 0.85, p 0.001).
The maximum measured shrub age was 77 yr (Salix glauca), and average shrub age 28 yr. Wedging rings were observed frequently. Only a few missing outer rings were observed in the basal stem discs, which lead to a decrease in sample depth of the chronologies from 2006 onward. Basic statistics of the chronologies are given in table 1. rEPS of the master chronology was higher than 0.85 from 1975 to 2010, i.e. for the period with a sample depth larger than 5. The master chronologies of both species were significantly   correlated with each other (r = 0.67, p 0.001, glk = 0.70) and showed significant positive correlations with average peak growing season temperatures (i.e. July-August temperatures-the time of the growing season with the strongest impact on shrub growth according to correlation analyses, r = 0.36, p < 0.01, glk = 0.67 for S. glauca and r = 0.61, p 0.001, glk = 0.65 for S. lapponum). The master chronology consisting of both species was positively correlated (r = 0.45, p < 0.001) with peak season temperatures as well (figure 3). Sample depth and correlation coefficients among RWI and peak season temperature showed a strong positive correlation (r = 0.91, p 0.001; figure 4), i.e. the higher the sample depth (and thus in general, the shorter the time series), the higher the correlation among RWI and temperature data. In the same way, the correlation coefficients among RWI and GSM were positively correlated with the sample depth (r = 0.49, p = 0.01), but the correlation coefficients only improved slightly above a sample size of 4 ( figure 4).
GSM was negatively correlated with each chronology (r = −0.45, p < 0.01, glk = 0.35 for S. glauca and r = −0.52, p < 0.001, glk = 0.33 for S. lapponum chronologies,   respectively, and r = −0.48, p < 0.001, glk = 0.28 for the combined master chronology, see also figure 5). If only the period 1975-2010 (i.e. with EPS > 0.85 of the combined master chronology) was considered the correlation coefficient among GSM and RWI of the combined chronology was −0.67 (p 0.001, glk = 0.28) and even −0.80 (p 0.001, glk = 0.18) if the maximum available sample depth (i.e. 24 shrubs, from 1995 to 2006) was considered (see also table 2 for correlation coefficients of the other chronologies). Please note that the mentioned glk values among GSM and RWI are lower than 0.5 due to the negative correlation among these two variables. Comparison of the complete (1975-2010) reconstructed and measured GSM revealed a significant and high correlation (r = 0.67, p < 0.001) and a glk of 0.71 (see also figure 6). The verification statistics further reveal that the reconstruction of GSM from RWI is robust for the considered

The potential of shrub RWI as proxy for GSM
The results indicate that RWI of S. glauca and S. lapponum are suitable proxies for GSM. Both, the RWI chronologies and the GSM of the Hardangerjøkulen icecap were significantly correlated with summer temperatures. Summer temperature is known to influence GSM (e.g. Braithwaite and Olesen 1989, Hagen and Liestøl 1990, Giesen and Oerlemans 2009 and several investigations have shown that shrub growth is strongly affected by summer temperatures in cold environments (Baer et  I.e. in summers warmer than average, shrubs will generally express above average growth increments while GSM will be lower than average (i.e. values will be more negative and more melting will occur) and vice versa for summers cooler than average. Due to summer temperature as the common driver, it is reasonable to use RWI as proxy for GSM, despite the lack of a mechanistic link between these components. As RWI of the investigated shrub species only correlated with July and August temperatures, while GSM correlated with June-September temperatures, the explanatory power of RWI on GSM was lower compared to that of the climate data. Further, other environmental parameters than temperatures also affect shrub growth, contributing additional noise to the RWI-GSM relationship.
The reconstruction of GSM from RWI (with model parameters derived from two different calibration periods and applied to the respective other period) revealed fairly high correlation coefficients and glk with recorded GSM over the whole calibration/verification period (i.e. 1975(i.e. -2010. In addition, the calibration/verification statistics (table 3) give confidence in the reconstructed GSM. Therefore, we consider the observed relationships among RWI and GSM as robust and thus reliable in terms of reconstructions.
Our findings are in line with Xiao et al (2007), who found a similar indirect relationship in the Qilian Mountains, north-western China. However, these authors reported a weaker correlation between RWI and the glacier equilibrium line altitude. This is likely caused by too short a calibration period (only 9 yr) but possibly also due to a stronger noise in their data related to micro-topographical conditions of the shrub habitats (Xiao et al 2007).
We explain the discrepancy between the effective period of correlations among temperatures and RWI (July-August) and temperatures and GSM (June-September), by the onset of ice melting as soon as air temperatures exceed 0 • C. In contrast, cambial activity in woody plant species commonly initiates when air temperatures have exceeded 5 • C for several days and ceases as soon as daylight duration surpasses a specific threshold or first remarkable night-time freezing temperatures occur (e.g. Körner 2003, Rossi et al 2007, Xiao et al 2007. Therefore, the GSM variability of a certain part of the melting season (i.e. the spring period when temperatures are between 0 and 5 • C and the autumn period after the first frost and/or daytime threshold but with temperatures still above 0 • C) can most likely not be approximated by RWI. The lower correlation coefficients of the S. glauca chronology with GSM when compared to the S. lapponum chronology can be explained by the fact that at its beginning the S. glauca chronology consists of only one individual over a period of 14 yr (instead of only 1 yr with a sample depth of 1 in the S. lapponum chronology). When looking at correlation coefficients for periods with higher sample depths in the S. glauca chronology, this discrepancy was no longer observed. Therefore, parts of the chronologies with such a low sample depth should not be considered as already suggested by rEPS-values lower than 0.85 in the respective period.
Thus, as long as monthly, 'on-site' climate data are available, they have to be considered as more precise proxies for GSM. However, in many remote areas-such as Arctic or Alpine areas-the density of climate stations is low, while the spatial variability of climate can be rather high (especially in mountainous regions, e.g. Körner 2003). Approximating GSM by 'off-site' climate station data in these areas would likely produce a remarkable error. In addition, only few climate records are longer than hundred years and those that are, are often not situated in the vicinity of glaciers. Watson and Luckman (2004) and Laroque and Smith (2005) successfully reconstructed GSM of selected North American glaciers over the past 300 yr from tree-ring records, explaining GSM variances between 26% and 51%. If only considering the period with EPS > 0. 85 (i.e. 1975-2010), models derived from the data in our study would be able to explain 45% and even 64% of GSM variation for the period with the maximum sample depth (i.e. 24 from 1995 to 2006). Recalling that trees in many cases are absent in the vicinity of glaciers (in particular in the Arctic), RWI of selected shrub species (e.g. Juniperus with known ages of >200 yr; Hallinger et al 2010) could thus serve as the most suitable proxy known for the reconstruction of GSM, in case climate data are absent, and as long as growth of these shrubs is mainly determined by summer temperatures. Giesen and Oerlemans (2009) computed a twodimensional ice-flow model for the Hardangerjøkulen icecap and compared their model output with the same glacier mass balance data as used in our study (Kjøllmoen 2011). The winter mass balance was modeled with a high precision (r = 0.92), while modeled GSM showed a lower correlation with measured data (r = 0.75). Thus, regarding GSM, the explanatory power of the RWI investigated in our study for the same data set is comparably high (r = −0.67 for the part of the chronology with EPS > 0. 85, i.e. 1975-2010).
Weijers et al (2010) found up to 183 yr old Cassiope tetragona shrubs growing in Endalen on Svalbard, of which the growth variability was well explained (r 2 = 0.41) by July temperatures in the period between 1912 and 2008. Several shrub species are known to reach ages of more than 70 yr (Schweingruber and Poschlod 2005). Currently, several shrub chronologies from all over the Arctic are available; many of them positively correlated to summer temperatures (Myers-Smith et al 2011). Thus, it seems likely that more long-lived shrub species are suitable for the reconstruction of GSM. Investigations of the available data sets within a network could allow for large-scale GSM reconstructions in areas where no GSM data is available.
On the one hand, such reconstructions would provide detailed information on historical land-surface albedo changes in Arctic and Alpine areas, which would increase our understanding of the observed increase of global air temperatures throughout the 20th century. On the other hand, such reconstructions could provide an opportunity to better reconstruct the contribution of glacier ablation in Arctic and Alpine environments to the observed global sea level rise in the past century (e.g. Shepherd andWingman 2007, Radic andHock 2009).

Challenges and potential causes of errors
Concerning the increasing variance that was explained by our models with increasing sample depth, we want to stress that the rEPS value should ideally be higher than 0.85 for the reliable reconstruction of GSM as also stated by Wigley et al (1984). This means, that a sufficient amount of shrub individuals needs to be sampled-in the case of our study at least five. In addition, the reliability of a master chronology also increases with higher sample depths, as it enables the detection of (partly) missing rings.
A possible explanation for the increase in correlation strength (among climate and growth) with higher sample depths is the elimination of micro-site effects from individual shrubs through the incorporation of more shrubs in the master chronology (e.g. Cook and Kairiukstis 1990). Thus, individuals reflecting specific micro-site conditions will have less impact on the master chronology if other samples even out this noise. A possibility to lower model errors related to such kind of unrepresentative individuals would be to eliminate shrubs from the master chronology, which express low correlation and Gleichlaeufigkeit coefficients (e.g. Tape et al 2012). Another explanation for the observed increase in explained variance with increasing sample depth (and thus a shorter period of model calibration) could be, that the climate has changed throughout the investigated time period in a way that it becomes more important for shrub growth, leading to a successively better correlation among RWI and GSM. However, a time-series analysis of the main drivers of glacier dynamics (i.e. summer temperature and winter precipitation) did not indicate such a change.
In our study we detrended the ring-width series by application of a smoothing spline, which removes lowfrequency variability in time series that is assumed to be caused by biological or stand effects . We also investigated the correlations among raw ring-width series and GSM and found similar but slightly lower correlation coefficients (r = −0.6 instead of −0.7). This is likely related to some unknown biological effects in the sampled shrubs that were not related to summer temperatures and therefore added an additional noise to shrub RWI, thus decreasing the explained variance of the computed models. However, as soon as long-term reconstructions of GSM by RWI are aimed at, it seems possible that the application of a smoothing spline may remove a part of the desired signal in shrub RWI, in particular the observed increase in air temperatures since the Little Ice Age. To test the validity of the applied detrending, the investigation of a longer time series would be helpful. Possibilities to solve this problem (if it occurs) could be (I) the application of a regional curve standardization (e.g. Esper et al 2003), (II) a power transformation of ring widths with subsequent standardization and variance stabilization (Tape et al 2012), or (III) standardization by dividing ring-width series with their horizontal mean to conserve the temperature signal (Weijers et al 2012).
One major challenge of the described reconstruction is to find a sufficiently large number of high-aged and suitable (i.e. sensitive) shrubs. For instance, in our investigation the oldest shrub had an age of 77 yr, while the second oldest was only 44 yr old. In our study, a sample depth of five was necessary to achieve rEPS-values above 0.85. Thus our data would only allow for reliable GSM approximations throughout the past 36 yr, i.e. from 1975 to 2010, which would be no gain to the local GSM record of the Hardangerjøkulen icecap . However, there are several monitored glaciers in Norway (and most likely also in other parts of the world) where GSM measurements (if any) last back for shorter time periods (Kjøllmoen 2011), and there, the development and application of the described method with a larger sample size of comparably old shrubs could achieve valuable reconstructions. In addition, other long-lived shrub species across Arctic-Alpine habitats may provide longer time series for the reconstruction of GSM.
Possible causes of error in the described method are missing and wedging rings, commonly observed in shrubs growing at high altitudes/latitudes (Myers-Smith et al 2011, Wilmking et al 2012. To minimize this, the application of the serial sectioning technique (Kolishchuk 1990) and measurements along at least two radii of colored thin sections is strongly recommended. Another, more time consuming but likely more precise measurement technique would be the direct measurement of basal area increments as suggested by Baer et al (2006).
An additional requirement of the described method is calibration data of GSM. Many glaciers in remote areas have not been investigated for their mass balance. In that case, it seems possible to derive GSM approximations either from remote sensing images or from physical models, linking summer temperatures (which also can be reconstructed from RWI) with GSM (e.g. Olesen 1989, Braithwaite andZhang 2000). In the latter case, summer temperatures would be reconstructed from shrub RWI and these summer temperatures would then allow for the calculation of GSM with physical/glaciological models.

Conclusion
Our study shows that RWI from shrubs growing near glaciers can serve as proxies for GSM. If successfully applied to respectively longer chronologies, this method allows for extending GSM time series beyond instrumental records. Extended GSM records may bear the potential to (I) assess the impact of glacier melt within a particular area on the observed sea level rise of the 20th century if the method is applied to a network of shrub RWI (e.g. reconstructing the summer ablation of the Greenlandic Ice Shield), (II) increase the reliability of glacier models through extended calibration/verification periods and (III) reconstruct historical glacier retreats if successfully combined with lichenological and glaciological investigations. However, to examine these potentials in detail further investigations are needed.