Fluctuations of a Greenlandic tidewater glacier driven by changes in atmospheric forcing: observations and modelling of Kangiata

. Many tidewater glaciers in Greenland are known to have undergone signiﬁcant retreat during the last century following their Little Ice Age maxima. Where it is possible to reconstruct glacier change over this period, they provide excellent records for comparison to climate records, as well as calibration/validation for numerical models. These glacier change records therefore allow for tests of numerical models that seek to simulate tidewater glacier behaviour over multi-decadal to centennial timescales. Here we present a detailed record of behaviour from Kangiata Nunaata Sermia (KNS), SW Greenland, between 1859 and 2012, and compare it against available oceanographic and atmospheric temperature data between 1871 and 2012. We also use these records to evaluate the ability of a well-established one-dimensional ﬂow-band model to replicate behaviour for the observation period. The record of terminus change demonstrates that KNS has advanced/retreated in phase with at-mosphere and ocean climate anomalies averaged over multi-annual to decadal timescales. Results from an ensemble of model runs demonstrate that observed dynamics can be replicated. Model runs that provide a reasonable match to observations always require a signiﬁcant atmospheric forcing component, but do not necessarily require an oceanic forcing component. Although the importance of oceanic forcing cannot be discounted, these results demonstrate that changes in atmospheric forcing are likely to be a primary driver of the terminus ﬂuctuations of KNS from 1859 to 2012. We pro-pose that the detail and length of the record presented makes KNS an ideal site for model validation exercises investigating links between climate, calving rates, and tidewater glacier dynamics.


Introduction
Calving from tidewater glaciers (TWGs) presently accounts for up to 50 % of the mass loss from the Greenland Ice Sheet ( Van den Broeke et al., 2009). Determining controls on tidewater glacier dynamics over decadal to centennial timescales is crucial to understanding their contribution to sea level in a warming climate (Alley et al., 2010;Vieli and Nick, 2011). The ability to achieve this in Greenland has been restricted in part by the relative lack of TWG terminus observations prior to the satellite age, as well as evidence of terminus locations being spread across a disparate array of sources. However, the synthesis of these sources has previously allowed for multi-decadal to centennial records of TWG glacier behaviour to be reconstructed (e.g. Csatho et al., 2008;Bjørk et al., 2012;Weidick et al., 2012).
Such records provide potentially excellent calibration and validation records for numerical modelling efforts (Vieli and Nick, 2011). That is to say, numerical models that are capable of replicating observed terminus behaviour over decadal to centennial timescales will be better placed to predict the future behaviour of a TWG over similar timescales. Despite this, there remain few examples of modelling efforts that have attempted to calibrate their results against multi-decadal observational records (e.g. Colgan et al., 2012). The ability of most numerical models to replicate dynamics over such timescales using realistic inputs therefore remains largely untested.
By undertaking calibration/validation exercises, the sensitivity of the terminus to different climatic forcing can also be evaluated (e.g. Nick et al., 2013;Cook et al., 2014;Lea et al., 2014a). This is achieved by comparing the sensitivity of a modelled glacier to climate forcing against observations . With a knowledge of realistic ranges of forcing, this allows for evaluation of the relative importance of each in contributing to the observed TWG behaviour.
Changes in oceanic forcing are significant drivers of TWG retreat in Greenland (Murray et al., 2010;Straneo et al., 2010;, but their relative importance between glaciers appears to be dependent on geographical location, glacier geometry , and potentially fjord connectivity with the open ocean (Straneo et al., 2012). Model-based studies have also helped to demonstrate the sensitivity of some major outlet glaciers to air temperature changes (via enhanced runoff increasing crevasse water depth; Nick et al., 2013;Cook et al., 2014).
Where multi-decadal to centennial timescale climate data exist alongside records of terminus position, these provide the potential for robust evaluation of both numerical models and the importance of different drivers of TWG terminus change. In this study we aim to (1) reconstruct the fluctuations of Kangiata Nunaata Sermia (KNS), SW Greenland, from 1859 to present, coinciding with the availability of climate records; (2) use these data to evaluate the ability of a well-established climate-driven numerical ice-flow model in order to replicate its dynamics; and (3), in conjunction with fjord topography data, assess controls on the terminus stability of KNS over multi-decadal to centennial timescales.

Field site and climate data
KNS is the largest TWG on the west coast of Greenland, south of Jakobshavn Isbrae ( Fig. 1; . It is known to have undergone significant retreat since its Little Ice Age maximum (Weidick et al., 2012), retreating a total of 22.6 km, with at least 12 km of this retreat occurring prior to 1859, when climate forcing data are unavailable (Lea et al., 2014a). It is situated ∼ 100 km inland from Nuuk at the head of Godthåbsfjord, and currently has a calving flux of  A continuous record of mean monthly air temperature is available at Nuuk from 1866 to present (Vinther et al., 2006;Cappelen, 2012). Temperatures at Nuuk are known to be strongly correlated with those near to the terminus region of KNS throughout the year (Taurisano et al., 2004). For this reason, we take the Nuuk record as an indicator of the atmospheric forcing at KNS.
As with all TWGs around Greenland, there are no long observational records of fjord water temperatures adjacent to KNS, though detailed hydrographic studies of the fjord have been undertaken recently (Mortensen et al., 2011(Mortensen et al., , 2013. A shallow ∼ 80 m sill at the entrance to Godthåbsfjord at Nuuk has been suggested to limit the connectivity of the fjord to warm ocean waters at depth. In fjords where shallow sills do not exist, the incursion of these warm ocean waters is thought to have significantly affected the stability of TWGs Straneo et al., 2012). The presence of the shallow sill in Godthåbsfjord also results in significant tidal mixing at the fjord entrance, allowing for sea surface waters to be incorporated at depth, which are then advected into the fjord (Mortensen et al., 2011). These intermediate-level mixed waters have been proposed to significantly influence the energy available for submarine melting at the termini of the TWGs in Godthåbsfjord (Mortensen et al., 2013).
Due to the impact of surface waters near the fjord entrance on the energy balance of the fjord (Mortensen et al., 2011(Mortensen et al., , 2013, and the potentially restricted influence of warm coastal currents at depth (Straneo et al., 2012), we suggest that sea surface temperatures (SSTs) provide a good indicator of the relative oceanographic forcing affecting KNS. Such data have also been used previously to interrogate the role of oceanographic forcing on TWG stability where observations at depth are unavailable (e.g. McFadden et al., 2011;Bevan et al., 2012). The HadISST1 1 • × 1 • data set provides SST estimates for the period 1871-present (Rayner et al., 2003), with annual averages for the area immediately offshore from Nuuk (62 to 64 • N, 51 to 53 • W) used as an indicator of oceanographic conditions affecting Godthåbsfjord. Although the data used will in part be based on interpolation of observations (especially in the earlier part of the record), the data have been validated for west Greenland against independent records back to 1875 (Hanna et al., 2009). This therefore provides confidence in the results obtained from the HadISST1 data set.

Glacier reconstruction data
A combination of geomorphology, maps, photography (ground-based, oblique-aerial, and vertical aerial images), and satellite imagery are used to reconstruct the terminus dynamics of KNS. By 1859 KNS is known to have retreated between 12 and 15 km from its Little Ice Age (LIA) maximum extent (Lea et al., 2014a). The post-LIA maximum glacial geomorphology of KNS has been mapped, while previous analysis of a photograph taken in the 1850s, and a map published in 1859 places the terminus position somewhere inside the limit of a significant glacier readvance/stillstand (Lea et al., 2014a). We refer to this as the Akullersuaq Stade (after the headland that its maximum extent adjoins), previously referred to as the "1920 Stade" (Weidick et al., 2012). This event is renamed due to the uncertainty of the exact timing of the glacier maximum.
Where the full terminus cannot be observed in photographs, terminus position is determined indirectly using the GIS-based analyses described below, in conjunction with evidence from maps (e.g. Lea et al., 2014a). Subsequent to 1921, intermittent direct observations of the terminus are available, enabling mapping of terminus positions from imagery (list of sources in Table 1). Terrestrial photo J. Møller in Bruun (1917Bruun ( ) 1921 Terrestrial photo A. Nissen in Weidick et al. (2012Weidick et al. ( ) 1932 Terrestrial photo A. Roussell in Roussell (1941)  Landsat panchromatic band imagery was used to map terminus positions for 1987-2012. Cloud-free Landsat scenes were selected for analysis, acquired as late in the melt season as possible, or just after its end. The start of November was used as the latest date from which images could be selected, since mélange in the fjord has been observed to freeze beyond this, causing the terminus to advance (Mortensen et al., 2011;Sole et al., 2011). The majority of images were acquired during September or October, though cloud-free images for 1993 and 2003 were only available for dates in August (30 August 1993 and 9 August 2003 respectively). No suitable images were available for the years 1988-1991 and 1998, meaning that annual resolution rates of terminus change were acquired for 1992-1997 and 1999-2012 (Table 1).
For the entire length of the record, where more than 1 year separated terminus observations, annually averaged rates of change were calculated. This provides a continuous record of the trends in behaviour, as well as inter-annual variability of KNS for the period spanning 1859-2012. This behaviour could then be directly compared to atmospheric and oceanic climate data.
Each terminus position was quantified using an adaptation of the box method (Moon and Joughin, 2008;Howat and Eddy, 2011), called the curvilinear box method (CBM; see Lea et al., 2014b, for details). This has a marked advantage over the centreline tracking or standard box methods, as it is capable of accounting for changes in terminus geometry while also accurately tracking changes in fjord orientation (Lea et al, 2014b). Furthermore, the box used to calculate terminus change is always centred on the glacier/fjord centreline, which is also the flow line used for the numerical model. Consequently, terminus positions and observed changes in position derived using the CBM can be compared directly to model output.

The TWG model
The numerical model used is specifically designed to simulate the dynamics of TWGs along a flow band (Nick et al., 2010). It has been successful in replicating the dynamics of marine terminating outlets in both Greenland (e.g. Vieli and Nick, 2011;Nick et al., 2012;Lea et al., 2014a) and Antarctica (Jamieson et al., 2012(Jamieson et al., , 2014, and has also been used to make centennial timescale projections of the future contribution of Greenland's major TWG outlets to global sea level . The model uses a stretched grid, allowing for a robust treatment of grounding line dynamics (Pattyn et al., 2012), while basal, lateral, and longitudinal shear stresses are accounted for. Bed topography data for the majority of the catchment are provided by Bamber et al. (2001), though the lower 40 km is generated using a mass-continuity-based bed reconstruction (Morlighem et al., 2011), validated against available OIB/CReSIS flight lines . Fjord width (Fig. 3c) is defined as the sum of the minimum linear distances from a point on the flow line to either side of the fjord (Fig. 3a). Where available, fjord bathymetry data are also used where KNS has retreated following its LIA maximum ( Fig. 3c; Weidick et al., 2012). Sensitivity analyses conducted by Lea et al. (2014a;their Fig. 10) for this bed configuration demonstrated that the model exhibits broadly comparable patterns of retreat behaviour where bed elevation is varied within an uncertainty of ±50 m.
A constant height versus surface mass balance (SMB) relation is used to calculate SMB for the ablation zone of KNS (Eq. 1). This is derived from the average RACMO SMB model output for 1958(Van Angelen et al., 2013. Table 2. List of parameters and constants used for running the model.

Parameter/constant Value
Ice densityρ i 900 kg m −3 Meltwater densityρ w 1000 kg m −3 Proglacial water body densityρ p 1028 kg m −3 Gravitational accelerationg 9.8 m s −2 Friction exponentm 3 Friction parametersµ and λ 1 Glen's flow law exponentn 3 4.5 × 10 −17 where b(x) is the SMB for position x (the along-flow coordinate) on the model flow line and h(x) is the glacier elevation for position x on the flow line. Due to the tendency for overestimation of accumulation in RACMO in this region (Van As et al., 2014), positive SMB values in the upstream section of the modelled glacier are prescribed to be lower than the RACMO output, allowing for the glacier to maintain its contemporary elevation profile. Irrespective of this, SMB forcing has previously been demonstrated to be of minimal importance to results of modelled TWG dynamics over the timescales that are being investigated (Lea et al., 2014a). The model is initialised using a glacier geometry approximating that of the Akullersuaq Stade maximum (ASM), derived from geomorphological mapping of associated trimlines ( Fig. 1). Constants and parameter values used are summarised in Table 2, while the initial tuning procedure followed for this configuration is the same as that used by Lea et al. (2014a). Surface runoff (Van As et al., 2014), air temperature (JJA average), and SST (annual average) data are used to drive changes in crevasse water depth (d w ) and submarine melting (M) respectively. Although seasonal cycles in velocity are observed at KNS within 20 km of its terminus , at locations >35 km from the terminus these have been demonstrated to have negligible effect (∼ 1 %) on net annual motion (Sole et al., 2011). Given that the timescales of interest are annual to decadal, seasonal variability in basal and lateral sliding is therefore not included within the model experiments. The model uses an effective pressure sliding law, allowing for it to replicate a typical tidewater glacier velocity profile, accelerating towards its terminus. Two zones of constant basal roughness (upstream and downstream) are prescribed to allow for the model to replicate observed elevation and velocity profiles (Lea et al., 2014a). This also ensures that particular areas of the fjord are not biased towards advance/retreat behaviour. All parameters which control the

Relating crevasse water depth to air temperature
Changes in the value of d w have previously been related to runoff variability (e.g. Nick et al., 2010;Cook et al., 2012Cook et al., , 2014, and have been successfully used as a climate-linked forcing directly affecting terminus change . This is achieved via a physically based crevasse water depth calving criterion (Benn et al., 2007), where crevasse penetration depth, as well as potential for calving, is enhanced by d w . However, the only previously used scaling of surface runoff to d w requires a baseline d w value to be prescribed, which it cannot fall below (Nick et al., 2013, their Eq. S3).
To remove the need to define a minimum d w value at the beginning of each model run, we present a new, unrestricted parameterisation that relates seasonal changes in monthly surface runoff to d w , and allows for d w to freely evolve due to changes in annual runoff (Eq. 2).
where d wNew is the new crevasse water depth for a particular month; d wPrev is the crevasse water depth from the previous month; α 1 is the coefficient relating crevasse water depth sensitivity to changes in runoff; R year represents total runoff for a given year (Gt yr −1 ); β month is the fraction of annual runoff occurring in a particular month; and R base is a baseline/long-term average annual runoff total (Gt yr −1 ), equivalent to the annual volume of water that is either refrozen within the glacier or drains from the crevasse to the bed. This assumes that the rate of refreezing/drainage of water from crevasses is constant from year to year. Where annual runoff exceeds R base , the average annual d w will therefore increase, and where runoff falls below R base , the average annual d w will decrease. Dividing annual runoff into each month's contribution also allows for the direct incorporation of d w 's seasonal variability. The value of d w will therefore reach its annual minimum prior to the onset of the melt season, and peak in August. The coefficient α 1 allows for the sensitivity of d w to changes in runoff to be adjusted, and is used as a tuning parameter.

Definition of β month
The fraction of annual runoff occurring in each month, β month , is derived from analysis of each month's average runoff from the catchments of both KNS and Akullersuaq Sermia (AS) over the period 1960-2012, as given by highresolution SMB modelling of the region . The runoff values for KNS and AS are summed since the glaciers were confluent for much of the time since their LIA maximum, including a significant portion of the period of interest of this study (see below, and Wedick et al., Runoff values prior to 1960 are therefore estimated using the relation that exists between average June-July-August (JJA) air temperatures (A JJA ) from Nuuk for 1960-2012 (Cappelen, 2012) and the modelled runoff values (r = 0.75). A regression equation is generated from this (Eq. 3), allowing for runoff estimates (Gt a −1 ) for the period 1866-1959 to be made from the Nuuk air temperature ( • C) record (Vinther et al., 2006;Cappelen et al., 2012).
Combined with the 1960-2012 modelled values, this produces a continuous record of estimated annual runoff for 1871-2012. Average monthly variability in runoff is superimposed on this record using the β month term.

Confluence with AS: adjustments to d w and ice flux
While KNS and AS are confluent in model simulations, variability in d w at the terminus is driven by total runoff values from both catchments. The confluence area of the two glaciers is defined on the model flow line as being 5 km across, lying between 4 and 9 km from the 2012 terminus position. However, as KNS retreats through the confluence with AS, the runoff contribution from AS to the terminus is removed, meaning that d w needs to be scaled to reflect this. Modelled annual runoff totals for each catchment show that KNS and AS respond directly in phase with one another (r = 0.99), with KNS accounting for 70.3 % (MARv3.2) or 74.6 % (RACMO2) of total runoff (Van As et al., 2014). To allow for this reduction in runoff as KNS retreats through the confluence, the value of d w is multiplied by a scale factor, γ , that will have a fixed value for each model run of between α 2 (a confluence scaling factor) and 1, such that Because AS and KNS will at times be partially confluent, the value of γ is also scaled linearly with respect to the relative position of the terminus through the confluence, such that γ = 1 when they are fully confluent, and γ = α 2 when fully diffluent. Values are varied linearly between α 2 and 1 for terminus positions within the confluence according to where x conf is the distance of the terminus through the confluence and X conf is the total flow-line distance over which the confluence occurs. Due to uncertainty regarding the precise scaling of runoff to d w as KNS retreats through its confluence with AS, as well as other confluence effects, α 2 is used as a tuning parameter within the model. The extra ice flux contribution from AS when confluent with KNS is estimated to be approximately one-sixth of that of KNS, based on the contemporary across-glacier velocity profiles  and terminus widths of AS and KNS. This extra flux is added to the modelled glacier as positive SMB at the confluence of KNS and AS, distributed along the flow line proportionate to the contemporary AS acrossglacier velocity profile (Lea et al., 2014a).

Relating submarine melt to sea surface temperature
Submarine melt rate (M) has previously been linearly related to deep-ocean temperature (DOT) variability using a scaling coefficient (Nick et al., 2013, their Eq. S2). Using this parameterisation, the highest values of M (expressed in this study in km 3 a −1 ) are associated with the highest scaling coefficients. Therefore high scaling factor values would also be linked to the highest inter-annual variability of M. This study takes a slightly different approach in that (1) M is scaled to SST rather than DOT, for reasons relating to fjord circulation explained above, and (2) we introduce a constant (minimum) baseline M rate, M base , which is added to the linear relation with SST. We therefore calculate M (km 3 a −1 ) according to where α 3 is a submarine melt rate scaling coefficient, and T year is the annual average SST. This allows for multiple minimum background rates M base to be tested for different model runs, with various sensitivities of M to changes in SST superimposed upon this using α 3 .

Model experiments and evaluation
Tuning parameters α 1 , α 2 , α 3 , and M base were varied randomly within prescribed limits for a total of 1500 Monte Carlo-style model runs. These were defined at the start of each run's spin-up period and held constant throughout. The limits for each of the tuning parameters were (1) α 1 , between 0 and 1.5; (2) α 2 , between 0.3 and 0.8; (3) M base , between 0 and 0.7 km 3 a −1 ; and (4) α 3 , between 0 and 0.3. These ranges of α 1 and α 2 were chosen to reflect a wide range of potential forcing scenarios, while the values of M base and α 3 were chosen so total submarine melt rates could potentially range from 0 km 3 a −1 to values that exceed those estimated for other TWGs in western Greenland Enderlin and Howat, 2013). This allowed for the different potential drivers of the observed terminus change to be comprehensively assessed. Runs were conducted for the period 1871-2012, given that this is the period when both atmospheric and oceanic climate records are available. The model was initialised at approximately the ASM profile and terminus position, as defined by the geomorphology, and given the duration of the spin-up period to stabilise for the given forcing scenario. During spin-up, d w was allowed to freely evolve by up to ±3 m a −1 to allow for the terminus to stabilise at the ASM, with R base and T year held constant. These were defined as the 1871-1920 runoff average (3.107 Gt yr −1 ) and SST average (2.605 • C) respectively. These values were used for spin-up as it is known the ASM was attained at some point within this time window. Model results were evaluated against their ability to replicate observed terminus dynamics, where absolute terminus positions are known (i.e. 1921 to 2012). The period from 1871 to 1920 therefore effectively becomes a transient spinup period, where the model is driven using real climate data, though terminus position is only known within a range. The ability of each model run to replicate observed dynamics was determined using a weighted regression (R 2 ) calculation, with the weighting of each terminus observation calculated according to where w is the observation weighting in the regression calculation, n is the terminus observation, k is the total number of terminus observations, and D is the date of the terminus observation. Each terminus observation is therefore temporally weighted according to the median length of time elapsed between the terminus observations that occur before and after observation n. This ensures that the evaluation of model performance is not biased towards the last ∼ 20 years, where there is a comparatively high density of observations. Model runs were counted as successful where (1) the difference between the modelled and observed 1921 position was < 500 m, (2) the weighted R 2 was > 0.85, and (3) the gradient of the resulting line of regression was > 0.85.

Glacier reconstruction results
The geomorphology shows distinct upper and lower sets of lateral moraines on both sides of the fjord, with fluted moraines occupying the intervening space (Fig. 1a). The upper set are associated with the LIA maximum (Lea et al, 2014a), while the lower set were formed during the Akullersuaq Stade. Fridtjof Nansen's (1890) account of the first traverse of Greenland in 1888 includes a drawing from a photograph showing AS and KNS to be confluent, though the terminus position itself is not visible. Although the original image could not be traced or an exact date of acquisition determined, it is likely to have been taken some time near to the publication date of 1890. Maps from 1859, 1860, 1866, and 1885 all show the terminus of KNS to be adjoining Akullersuaq and fully confluent with AS (Kleinschmidt, 1859;Poulsen, 1860;Brede, 1866;Rink, 1866;Jensen, 1885). While it is possible that some details on the maps were copied following Kleinschmidt (1859), the addition of detail such as lakes on plateaus near to KNS by Jensen (1885) provides confidence that this map faithfully records the contemporary terminus position. There is nothing to suggest that KNS became diffluent from AS at any time from 1859 to 1885. However, due to a lack of map detail and the Nansen (1890) drawing not including the terminus, these sources cannot be used to provide absolute terminus positions.
The earliest images of KNS are from the 1850s and 1903. Both are taken from approximately the same position, with the terminus partially obscured by foreground topography (Weidick et al., 2012). The presence of medial moraines in each image demonstrates that KNS was confluent with AS. Lea et al. (2014a) quantified the terminus position uncertainty for the 1850s photograph using viewshed analysis. Similar analysis has been undertaken for the 1903 image, showing that the uncertainty in terminus position is the same as for the 1850s image (Fig. 3). The maximum terminus extents for both images are therefore located behind a headland corresponding to the ASM on the eastern side of the fjord (Figs. 1a, 3).
It is not currently possible to say from any observational evidence when the ASM was attained, only that it occurred sometime between 1859 and 1920. The climate anomalies for the period (compared to    (Cappelen et al., 2012;Vinther et al., 2006). (d) Annual SST anomaly for the area 61 to 6 • N, 51 to 56 • W at annual resolution (white bars) with red line showing the averaged SST anomaly between terminus observations (Rayner et al., 2003).
air temperature (AT) and SST anomalies were, on average, anti-phased for the period 1871-1903 (Fig. 4c, d), though AT and SST anomalies are in phase (negative/near-baseline) for 1903-1920. Conditions are therefore more likely to have been conducive for glacier advance during the latter period. Terminus position was mapped directly for the remaining images, providing a record of 29 terminus positions spanning the period 1921-2012 ( Figs. 1 and 4). The first direct terminus observation (1921) shows a slight retreat from the ASM. Subsequent to this, KNS retreated a total of 9.7 km at a nonuniform rate up to 2012, interrupted by short periods of readvance (Fig. 4a, b). Averaged retreat rates of −116 m a −1 are observed between 1921 and 1946, before a rapid retreat of 3.9 km within the 2-year period from 1946 to 1948 (Figs. 1a,  4). Between 1948 and 1968 KNS retreated on average by −97 m a −1 , before readvancing by +60 m a −1 up to 1979 (Fig. 4b). A terrestrial photograph taken in 1965 with the majority of the terminus obscured shows the termini of KNS and AS to be fully diffluent.
The 1921-1968 period of sustained retreat was accompanied by positive average AT and SST anomalies (Fig. 4c, d). The highest AT anomalies occurred during the period 1928-1941, though the largest retreat (between 1946 and 1948) occurred during a comparatively less extreme period of positive AT and SST (Fig. 4).
From 1979 to 1987 KNS retreated by −658 m in total (−82 m a −1 ), before readvancing by +758 m from 1987 to 1992 (+152 m a −1 ). Using the near-complete 20-year annual record of terminus fluctuations from 1992 to 2012, KNS advanced for 4 out of 5 years between 1992 and 1997, followed by retreat in 11 out of 13 years from 1999 to 2012 at an average rate of −103 m a −1 . The latter included eight annual retreats of > 100 m, with the largest retreats occurring in 2004 (−438 m) and 2005 (−316 m). These periods of advance and retreat behaviour occurred during periods of inphase negative and positive climate anomalies respectively.
Where temporal density of observations was high, terminus behaviour that was anti-phased with the prevailing climate anomalies was also observed (i.e. advancing during positive temperature anomalies, or retreating during negative temperature anomalies). Examples of this include a retreat of −626 m in 1995, when both climate anomalies were negative, while terminus advances occur in 2008 and 2009 despite markedly positive AT and SST anomalies (Fig. 4). At annual resolution, the magnitude of terminus retreat/advance was also found to be unrelated to the magnitude of either climate anomaly for each particular year. Based on interpolated terminus positions between observations, terminus widths were consistent at ∼ 3.5 km from 1932 to 1946, and ∼ 4.2 km from 1968 to 2012, when terminus change was comparatively slow (Fig. 5b). Although fjord depths at the terminus for these periods were more variable, they did not exceed a range of ±22 m. Fjord width and depth at the terminus displayed two step changes during the retreats between 1921 and 1932 and between 1946 and 1948 (black lines, Fig. 5a-c). During the first of these, both width and depth increased (by ∼ 550 and 44 m respectively), whereas width increased but depth decreased during the second (by ∼ 700m and 146 m respectively).

Model results
From a total of 1500 model runs conducted, 29 runs (1.9 %) successfully replicated the observed dynamics of KNS according to the criteria outlined above (Fig. 5). Following the initiation of climate forcing in 1871 (Fig. 5d, e), the results of each run are highly comparable up to 1884, with little modelled terminus change observed. Following this, for the period 1884 to ∼ 1910, 6 of the 29 runs (21 %) show evidence of multi-annual terminus retreats and equivalent readvances of > 750 m with periodicities of 2-4 years. A further seven runs (24 %) show evidence of at least one short-lived (< 5year) oscillation in terminus position of > 750 m between 1884 and 1920. None of these model runs significantly exceed the ASM position, and thus they are in agreement with the geomorphological evidence presented, as well as the position of the 1921 terminus observation.
All model runs retreat to the observed 1932 position between modelled years 1929 and 1936, via a single retreat event of ∼ 1 km. Subsequent to this, modelled retreat to the observed 1946 position is gradual, before the model successfully replicates a large topographically controlled retreat from the 1946 position. There was varying success in modelling the exact timing of this retreat (observed between 1946 and 1948), with the model ensemble predicting it to occur anywhere between 1943 and 1962. The position where the modelled terminus restabilises following the retreat through the AS confluence is generally too far advanced by ∼ 1 km compared to the position following the 1946-1948 retreat. All model runs then go on to over-predict terminus extent for the 1968 observation by between 0.35 and 1.59 km.
Though no model runs exactly match the precise interannual terminus fluctuations from 1968 to 2012, they do capture the general multi-annual to decadal pattern of retreat observed. This is characterised by general terminus stability within a range of ±500 m for the period 1968 to ∼ 1999, before the terminus begins to retreat ∼ 2 km towards the 2012 position. All of the successful model runs identified predict KNS to be in a more retreated position in 2012 than observed by a range of 0.32 to 5.04 km. Where a significant difference between observed and modelled terminus positions has occurred by the end of the model run in 2012, the divergence begins in 2010 at the earliest. This coincides with a widening of the modelled fjord associated with the uncertainty in fjord topography upstream of the contemporary terminus (Fig. 5b).
The distributions of tuning parameters for successful runs are shown in Fig. 6, with the distribution of all histograms shown to be non-normal. Submarine-melting-related tuning parameters α 3 and M base , tended towards the mid-to lower ends of the ranges tested (Fig. 6c, d). Values of α 3 peak between 0.075 and 0.1, though there is no clearly defined peak in the distribution of M base values.
In contrast, none of the d w -related tuning parameters (α 1 and α 2 ) approach 0 ( Fig. 6a, b), with the lowest values being 0.412 and 0.389 respectively. Construction of a correlation matrix comparing all tuning parameter values for all successful runs also demonstrates a significant inverse relationship between the value of α 1 and the AS confluence parameter, α 2 (r = −0.92). While other significant correlations are observed (Table 3), these are not of sufficient strength to allow for confident conclusions to be drawn.

Observed terminus behaviour
From 1903 to 2012 AT and SST anomalies covaried, with the terminus generally undergoing retreat during periods of positive anomalies and advancing/stabilising when near/below baseline climate (Fig. 4). Exceptions to this  1987-1997) and positive (1998-2012) climate anomalies, the terminus responds in phase with the climate anomalies. This demonstrates the risks of using short data sets (2-5 years) to determine how a TWG is responding to climate forcing, highlighting the inherent noisiness, potential importance of antecedence, and the non-linearity of TWG response to climate.
A notable caveat to this occurs where significant topographically controlled glacier retreats occur (i.e. those driven by changes in fjord width and/or depth). These events could potentially skew annually averaged terminus change rates when attempting to characterise terminus response to climate forcing. The relative importance of this will be entirely dependent on the magnitude of individual events, and most significant where there is potential for multi-kilometre topographically controlled retreat. For example, if the 1946-1948 retreat event was not temporally well constrained, it could have significantly biased the terminus change rate values between 1936 and 1968 (Fig. 4b).
The 1946-1948 retreat occurs where the fjord widens and shallows at the terminus, while the 1921-1932 retreat is associated with a fjord widening and deepening (Fig. 5a-c). The 1946-1948 retreat is therefore likely to have been controlled by changes in lateral topography rather than basal topography, whereas the 1921-1932 retreat (if it occurred rapidly, e.g. in 1-2 years) likely resulted from a combination of both. In the periods between these > 1 km retreats, both fjord width and depth at the terminus remained largely consistent (Fig. 5b, c). While kilometre-scale, rapid retreat of KNS is likely due to a combination of retreat into fjord widenings or deepenings (e.g. Mercer, 1961;Carr et al., 2013Carr et al., , 2014Porter et al., 2014), the 1946-1948 retreat helps to demonstrate that destabilising changes in one aspect of fjord topography can dominate stabilising changes in the other, until a new equilibrium is reached.
Since TWGs exhibit varying degrees of non-linearity in response to climate forcing, the identification of where and when these rapid multi-kilometre retreat events occur is crucial for interpreting the causes of terminus fluctuations. Where comparatively smaller (i.e. < 500 m) climatically The Cryosphere, 8, 2031-2045, 2014 www.the-cryosphere.net/8/2031/2014/ anti-phased advance/retreat events occur, their effect on average terminus change rates can be mitigated by averaging change over timescales up to or greater than a decade. For example, extending the 1992-1997 average (51 m a −1 retreat) to cover the period 1987-1997 (91 m a −1 advance) provides a more representative impression of multi-annual terminus behaviour, since five out of the six observations available show terminus advance. Where observations are separated by > 1 year, interpreting the absolute values of terminus change rates should therefore be done with caution. In most cases these values will be more representative of the average direction (i.e. advance/retreat), rather than the average distance of terminus change.
With uncertainties due to topographic controls on terminus stability taken into account, observations of terminus change over a period of several years provide a better indication of a TWG's response to climate forcing. However, for this study, deconvolving the relative importance of AT versus SST in driving terminus change is difficult using observations alone, given that both climate drivers vary in phase for 1903-present. It could potentially be argued that AT is the primary driver of change, since the 33-year period of positive anomaly SST from 1871 to 1903 had relatively little impact on the terminus stability of KNS. However, a narrow and relatively shallow fjord geometry in this region could also have been a significant factor in stabilising the terminus during this time (Fig. 3c). Arguably this becomes less likely when it is considered that, while SST was similar for the period 1921-1948, positive AT allowed for KNS to retreat through the same section of fjord and through its confluence with AS within 26 ± 1 years (Fig. 4). However, given the lack of certainty in terminus position between 1871 and 1920, it is not possible to robustly verify these arguments.

Implications of modelling
The observed terminus behaviour of KNS from 1921 to 2012 was successfully replicated by 29 of 1500 model runs using surface runoff and SST records as drivers of terminus change. This demonstrates that the parameterisations used to scale these climate records to d w and M respectively can successfully be used to simulate the observed pattern of tidewater glacier behaviour over centennial timescales. Where the observational record is of sufficient detail to resolve interannual terminus fluctuations , the model does not replicate these. This is to be expected given (1) the flowband nature of the model and associated depth and width integrations over each grid cell, meaning that fluctuations of terminus configurations such as the creation of calving bays cannot be replicated (e.g. Fig. 1b); (2) the uncertainty in fjord bathymetry and geometry potentially affecting relative terminus stability; and (3) the use of single terminus observations as notionally definitive indicators of annual terminus change, where the stochastic nature of calving and associated sub-annual terminus fluctuations make any direct one-to-one comparisons to modelled results inappropriate. Valid comparison of model results to observations should therefore only be attempted over multi-annual timescales where terminus dynamics within calving bays, sub-annual calving events and fine-scale uncertainties in fjords, and basal topography become comparatively less significant.
For successful model runs, the interrelationships between the parameter values that determine d w and M sensitivity to the climate records also inform the relative importance of changes in atmospheric and oceanic forcing in driving terminus change. The lack of any significant relationship between α 1 and α 3 demonstrates that a change in model sensitivity to surface runoff is not offset by any change in model sensitivity to SST (e.g. a higher α 1 would not need to be offset by a lower α 3 for the model run to match observations). Taken alone, this evidence indicates that either atmospheric forcing (via surface runoff) dominates oceanic forcing (via SST) or vice versa. However, the occurrence of runs where α 3 does not significantly exceed 0 (i.e. where runs experience negligible M variability) demonstrates that the model can successfully reproduce observed behaviour with nearly no changes in oceanic forcing from year to year. Although some successful model runs did have significant inter-annual M variability (e.g. the maximum range of M values for an entire 141-year model run was 0.76 km 3 a −1 ), each model run always requires significant atmospheric forcing variability to allow for it to replicate observations. The importance of oceanic forcing variability can therefore not be entirely discounted.
The model demonstrates that knowledge of atmospheric forcing (via runoff), without needing to vary oceanic forcing, can be sufficient to reproduce realistic patterns of observed glacier behaviour at KNS over the last century. However, the precise physical mechanism by which air temperature could drive observed change requires further investigation. For example, though a combination of modelled and empirically estimated runoff values has been used to drive changes in d w to force the model, subglacial runoff variability is also known to drive rates of submarine melting at the terminus (Jenkins, 2011;Xu et al., 2012;Sciascia et al., 2013). Therefore we do not rule out that the behaviour observed could also be explained by calving driven by seasonal changes in submarine melt rates, which are in turn a function of subglacial runoff (e.g. Sciascia et al., 2013).
The relative insensitivity to changes in oceanic forcing is not necessarily surprising given the hydrographic setting of KNS -located at the end of a > 100 km long fjord system that is thought to be largely insulated from changes in ocean conditions due to the presence of a shallow sill at its entrance (Mortensen et al., 2011(Mortensen et al., , 2013. This has previously been used to suggest that recent changes in ocean conditions (e.g. Straneo and Heimbach, 2013) have not affected the dynamics of KNS significantly (Straneo et al., 2012). The results presented here are therefore compatible with this argument.
The overestimation of terminus retreat by 2012 of every successful run is thought to result from the poor knowledge of fjord width geometry beyond the contemporary glacier terminus. Upstream of the 2012 terminus, the lateral ice margins are used to define model glacier width, leading to a likely overestimation of the prescribed fjord width. The divergence between the actual and prescribed fjord width is likely to increase upglacier, increasing the likelihood of model error in this area. This explains why significant divergence from the observational record only occurs once the modelled terminus has retreated ∼ 1.5 km beyond the 2012 terminus (Fig. 5a-c). Given the shallowing of the fjord bathymetry upstream of the 2012 terminus (Fig. 5c), fjord width uncertainty is likely to be the major cause of the model overestimating retreat (Fig. 5b). This also substantiates observations that destabilising changes in fjord width can dominate stabilising changes in fjord depth. Any attempt at modelling the future fluctuations of KNS will therefore require both improvements to subglacial topography estimates and comprehensive assessments of fjord width uncertainties as part of any predictions.

Conclusions
Utilising multiple lines of evidence, it has been possible to reconstruct terminus fluctuations of KNS from 1859 to 2012. This study therefore completes the record of terminus fluctuations of KNS from its LIA maximum, in 1761 (Lea et al., 2014a), up to the present, providing one of the longest and most detailed records of observed TWG change in Greenland. The length and detail of this record, in conjunction with existing data sets providing boundary conditions, therefore make KNS an ideal validation site for models aiming to simulate outlet glacier retreat and/or the impact of calving on tidewater glacier dynamics. At present the major boundary condition uncertainty is fjord topography, though what is known is sufficient for the model used in this study to replicate observed dynamics over multi-decadal to centennial timescales.
Results from numerical modelling show that the fluctuations of KNS can be simulated through parameterisations that link surface runoff to a crevasse-water-depth-based calving criterion. Changes in crevasse water depth and/or runoffdriven rates of submarine melt are therefore suggested as potential drivers of observed change. Although ocean-driven changes in submarine melt rates are not always required for the model to replicate the observed length variations of KNS, results do not allow for their importance to be discounted entirely.
Observations of KNS show it to respond in phase with AT and SST anomalies over multi-annual to decadal timescales from at least 1921 to 2012 (i.e. retreating during positive temperature anomalies, and advancing during negative temperature anomalies). However, where inter-annual comparisons to AT and SST are possible , climatically anti-phased terminus fluctuations are observed. This highlights the inherent noisiness of terminus response over short timescales, the potential importance of antecedence, and the dangers of using similarly short calibration periods for predictive modelling efforts.
Results from numerical modelling successfully capture the terminus dynamics of KNS over multi-annual to decadal timescales, though not precise inter-annual fluctuations. This is due to a combination of uncertainties in fjord topography as well as the approximations inherent to the depth and width integrations associated with using a one-dimensional flowband model.
Nevertheless, this study demonstrates that simple flowband numerical models of tidewater glaciers can be used to capture TWG dynamics over multi-annual to centennial timescales. This provides validation that these models can be useful tools for palaeo-, contemporary, and prognostic modelling efforts. However, the primary challenge to their use as predictive tools remains the accurate definition of subglacial topography and fjord width, which exert dominant controls on glacier stability. Any future efforts at prognostic modelling of TWGs should therefore seek to account for these uncertainties in addition to those associated with sensitivity to climate forcing.