Surges of Harald Moltke Bræ, north-west Greenland: Seasonal
modulation and initiation at the terminus

Abstract. Harald Moltke Bræ, a marine-terminating glacier in north-west Greenland, shows episodic surges. A recent surge from 2013 to 2019 lasted significantly longer (6 years) than previously observed surges (2–4 years) and exhibits a pronounced seasonality with flow velocities varying by one order of magnitude (between about 0.5 and 10 m/day) in the course of a year. During this six-year period, the velocity always peaked in the early melt season and decreased abruptly when meltwater runoff was maximum. Our data suggest that the seasonality has been similar during previous surges, and, to a much lesser extent, during the intermediate quiescent phases. It is peculiar to Harald Moltke Bræ that the seasonal amplitude is amplified episodically to constitute glacier surges. The surge from 2013 to 2019 was preceded by a rapid frontal retreat and a pronounced thinning at the glacier front (30 m within 3 years). We discuss possible causal mechanisms of the seasonally modulated surge behaviour by involving various system inherent factors (e.g. glacier geometry) and external factors (e.g. surface mass balance). The seasonality may be caused by a transition of an inefficient subglacial system to an efficient one, as known for many glaciers in Greenland. The patterns of flow velocity and ice thickness variations indicate that the surges are initiated at the terminus and develop through an up-glacier propagation of ice flow acceleration. Possibly, this is facilitated by a simultaneous up-glacier spreading of surface crevasses and weakening of subglacial till. Once a large part of the ablation zone has accelerated, conditions may favour substantial seasonal flow acceleration through seasonally changing meltwater availability. Thus the seasonal amplitude remains high for two or more years until the fast ice flow has flattened the ice surface and the glacier stabilizes again.


. Flow velocity derived from Landsat, Sentinel-1, and TerraSAR-X at a point close to the terminus of Harald Moltke Brae. For Sentinel-1 we use two different velocity data sets: one processed at Denmark's National Space Institute (DTU Space), referred to as Sentinel-1 (D), and another one provided by the Programme for Monitoring of the Greenland Ice Sheet (PROMICE), referred to as Sentinel-1 (P).
Harald Moltke Brae is located in north-west Greenland at about 76.5°N (Fig. 2). It is one of three glaciers that terminate in the Wolstenholme Fjord (Mock, 1966). At its present (2020) front the glacier is about 5 km wide. Starting at the present terminus its main stream can be tracked about 65 km upstream (Fig. 2, red line). At a distance of about 10 km from the 30 present glacier front the Blue Ice Valley Glacier (Mock, 1966) flows into Harald Moltke Brae. The boundary between these two streams is clearly visible by a medial moraine along the entire distance from their confluence to the terminus. Based on the Arctic Digital Elevation Model (ArcticDEM) and flow lines identified in Landsat images, we estimate the size of the overall drainage basin to be about 1500 km 2 consisting of Harald Moltke Brae (1200 km 2 ) and Blue Ice Valley Glacier (300 km 2 ). As another remarkable feature, a 3 km long and 1 km wide lake abuts to the northern side of Harald Moltke Brae at about 20 km. 35 It might be an additional source of freshwater influx into the glacier system. Already prior to the era of satellite remote sensing, significant changes in the dynamic behaviour of Harald Moltke Brae were reported. Wright (1939) observed an exceptional advance of the glacier front by about 2 km between 1926 and 1928 and inferred that the average surface flow velocity in this interval was at least 1000 m/year (2.7 m/day). Mock (1966) used the displacement of ice-surface features visible in aerial and terrestrial photographs to show that between 1954 and 1956 the 40 average velocity was about 1 m/day, ten times higher than the average velocity between 1937 and 1938.
Based on its periods of clearly accelerated ice flow, Harald Moltke Brae has been assumed to be a surge-type glacier (Rignot et al., 2001;Rosenau et al., 2015;Joughin et al., 2010;Hill et al., 2017Hill et al., , 2018. In general, the flow behaviour of a surge-type glacier is characterized by an alternation of long periods of low velocity (3-100 years, quiescent phases) and comparably short periods with velocities increased by at least one order of magnitude (1-10 years, active phases or surge events) (Jiskoot, 1999; The black triangle marks the point (67.628°W, 76.588°N) to which all velocity time series in this paper refer. Benn and Evans, 1998). During a quiescent phase, the lower part of the glacier flows slower than the upper part. Therefore, the transition part, called reservoir area, thickens dynamically while the lower part, called receiving area, thins (Jiskoot, 1999).
This pattern reverses in the active phase. When ice flow in the receiving area is enhanced during the active phase the glacier front advances rapidly. Most known glacier surges started with an increase of flow velocities in the upper part of the glacier and propagated down-glacier (Raymond, 1987;Solgaard et al., 2020;Wendt et al., 2017). However, some glaciers, such as the 50 Aventsmarksbrae in Svalbard (Sevestre et al., 2018), exhibit an opposite sequence with an initiation of the surge at the glacier front and its subsequent propagation up-glacier.
Only a small proportion of about 1 % of the glaciers worldwide is assumed to exhibit a surge-type behaviour (Jiskoot, 1999). Such glaciers normally cluster in certain regions (such as Alaska, Svalbard and the Karakoram) where the conditions are favourable for a surge behaviour (e.g. a soft glacier bed and an at least partially temperate regime) (Benn et al., 2019). 55 Some surge glaciers are also known in Greenland, such as the Hagen Brae in north-east Greenland (Solgaard et al., 2020). In north-west Greenland, however, no surge glaciers are known except for Harald Moltke Brae (Hill et al., 2017).
Two different mechanisms, both involving an internal instability, have been considered as possible causes of glacier surges: (A) thermal and (B) hydrological (Murray et al., 2003). (A) The base of a polythermal glacier can be partly frozen to the glacier bed such that the ice flow is in part restricted during the quiescent phase. The increasing shear stress due to the steepening of 60 the reservoir area during the quiescent phase can subsequently trigger a feedback with an initially slow movement producing friction heat. This leads to a further acceleration and enables basal sliding over large parts of the glacier base (Murray et al., 2003;Jiskoot, 1999). (B) A hydrologically driven surge event occurs when the subglacial drainage system closes after having been efficient over several years (quiescent phase), and becomes inefficient. When an increasing amount of meltwater meets this inefficient drainage system, subglacial water pressure increases and induces basal sliding (Murray et al., 2003;Jiskoot, 65 1999). Previous studies tended to explain the surge behaviour of Harald Moltke Brae by a feedback mechanism associated with weakening of soft subglacial sediments, that is, by a thermally driven mechanism (Hill et al., 2018;Joughin et al., 2010;Mock, 1966).
While system inherent drivers are responsible for cyclic surge behaviour, the seasonality of the flow velocities is mostly ascribed to external forcing changing in the course of a year. According to Moon et al. (2014) seasonal velocity variations of 70 a large number of Greenlandic glaciers result from either the seasonally changing ice front position or the varying meltwater availability.
In this paper, we focus on the three active phases of Harald Moltke Brae detected by satellite remote sensing, 1999-2000, 2005-2006 and 2013-2019. The observation of the surge 2013-2019 is unprecedented in terms of the simultaneous strong seasonality with velocities decreasing to the level of the quiescent phase always in summer. In addition, the active phase 75 2013-2019 is significantly longer (6 years) than the previous two active phases , 1999-2000 and 2005-2006 (always 2 years).
In contrast to most known surge glaciers, the surges at Harald Moltke Brae are shown to be initiated at the glacier front and to develop up-glacier. We analyse and interpret this extraordinary flow behaviour in terms of possible causal processes. For this, we correlate spatio-temporal flow patterns to several system-inherent parameters (e.g. glacier geometry) and external parameters (e.g. surface mass balance). We support the interpretation by analysing dynamic ice-height changes predicted from 80 observed velocity variations.

Ice flow velocity data sets
To determine flow velocity fields for outlet glaciers in Greenland, Rosenau et al. (2015) developed a processing scheme based on the feature tracking method using suited Landsat images available since 1972. Resulting velocity fields are used in the 85 present study. In addition, we included three velocity datasets derived from SAR offset tracking: Sentinel-1 (P) and Sentinel-1 (D) processed by PROMICE  and DTU Space, respectively, and TerraSAR-X provided by MEa-SUREs (Joughin et al., 2020). Tab. 1 provides technical details about these four data sets. All are combined to derive a time series of monthly averaged velocity fields with a spatial and temporal coverage as high as possible (Appendix A). Thereby, an almost seamless time series can be inferred for the period after 2013. Gaps in the Landsat data set caused by polar night are 90 filled by data from SAR based techniques. The accuracy of the joint velocity time series is shown to be better than 0.5 m/day over most of the time, but can exceed 1 m/day in single months with rapidly changing ice flow (Appendix A). Table 1. Overview over the technical details of the velocity data sets from Landsat, Sentinel-1 (D), Sentinel-1 (P) and TerraSAR-X. The velocity fields are given in form of regular spatial grids. The time basis refers to the acquisition time difference of the image pairs used for the velocity determination. The temporal resolution is the time difference between consecutive velocity fields.

Ice front position
For the period between 1916 and 1960 ice front positions are taken from the previous studies (Koch, 1928;Wright, 1939;Davies and Krinsley, 1962). For the period after 1972, ice front positions are digitized on the basis of Landsat images. The 95 variation of the front position is measured by the average distance from its position in 1916 applying a method proposed by Moon et al. (2008) (Appendix B).
In addition to the measured surface elevation changes, we computed the monthly dynamic ice-height change rates to be 105 expected due to the flow velocity distribution. This provides information on geometric glacier changes with a high sampling rate, and enables a better understanding of how these changes are related to the flow velocity changes. To do so, we consider only the horizontal components of ice flow and assume parallel ice flow as well as a constant density of the glacier ice. Then, the relationship between ice flow and the dynamic height change at a given point is (1)

110
H denotes the glacier thickness, t denotes time, x denotes the position along the flow line and v is the velocity averaged over the vertical column of the glacier. We evaluate v by multiplying the observed surface velocities with a factor of 0.9 This factor adopted in the absence of more specific information is a rough approximation and a potential source of error particularly 5 https://doi.org/10.5194/tc-2020-256 Preprint. Discussion started: 6 November 2020 c Author(s) 2020. CC BY 4.0 License.
for surge-type glaciers (Wu and Jezek, 2004;Andersen et al., 2015). Note that the total surface height change sums up from this dynamic height change and the height change due to surface mass balance (SMB).

115
We implement Eq. 1 by using velocity fields and ice thickness data (see next section) interpolated to the main flow line. To suppress noise, we subsequently fit a 3rd degree polynomial to the results.

Bedrock topography
BedMachine (Morlighem et al., 2017) provides both ocean floor and bedrock topography in a gridded format. Additional bedrock data are available along profiles of airborne ground penetrating radar measurements by the Center for Remote Sensing 120 of Ice Sheets (CReSIS). Based on a comparative analysis of BedMachine and CreSIS data (Section 3.7 and Appendix C) we use BedMachine data only for areas above 16 km.

Parameters of external forcing
The Regional Climate Model (RACMO2.3p2) provides surface mass balance (SMB), air temperature, precipitation and snow melt on a monthly basis. We infer the monthly mean SMB, precipitation and snow melt averaged over the overall drainage 125 basin.
Additionally, we compute the ice-mass flux through a cross section orthogonal to the flow lines located close to the glacier front at about 18 km. To determine the cross-sectional area, we use the ArcticDEM and BedMachine.

Visible features
We also assess changes traceable by visual inspection of the Landsat scenes. We focus on four different features: lakes on 130 the glacier surface, outflow of subglacial meltwater at the glacier front (meltwater plumes), sea ice coverage in the fjord, and calving events (Appendix E). We assess these features by distinguishing between three states: not visible, moderate and strong occurrence.
3 Characterization of spatio-temporal patterns of glacier geometry, flow velocity and external influences 3.1 Front position 135 Fig. 3 shows that most of the active phases between 1916 and 2020 were accompanied by an advance of the terminus interrupting its long-term retreat (Mock, 1966). From 1926 to 1928 the frontal advance was more pronounced than in later active phases. During the surge event in 1954-1956, the glacier front retreated slightly.  The average SMB was exceptionally high (10 kg m −2 day −1 ) in 2013, the year of surge initiation. In contrast to the SMB maximum in 2004, the SMB maximum in 2013 is due to a low meltwater runoff rather than a high accumulation rate.
Directly before (September 2013) and after a phase of accelerated ice flow (August 2014), the velocities were about 1 m/day 175 or less over almost the entire glacier (Fig. 7). A slight acceleration already occurred in July and August 2013 (Fig. 6). From September to December 2013 the velocities increased rapidly in the part below 20 km, whereas the parts further upstream were 8 https://doi.org/10.5194/tc-2020-256 Preprint. Discussion started: 6 November 2020 c Author(s) 2020. CC BY 4.0 License. less affected or not affected at all. Subsequently, also the upper parts of the glacier accelerated significantly leading to velocities of more than 5 m/day within most of the glacier area below 25 km. Fig. 7e show that a small part at the glacier front with an extension of 2 −3 km had slightly

Ice-mass balance
The time series of the mass flow through a cross-sectional area close to the glacier front (Fig. 12)

Bedrock and ice surface geometry
We examined the bedrock topography from BedMachine (Morlighem et al., 2017) and from CReSIS, both along the same CReSIS flight path (Fig. 13)

Seasonal mechanism
Types (a) and (b) correspond to the seasonality of the glacier types (2) and (3) According to Moon et al. (2014), type (3) differs from (2) in autumn as there is still enough meltwater available to raise the water pressure significantly and, thus, to accelerate the ice flow. At Harald Moltke Brae, we find a rather converse relationship: Compared to years with low meltwater runoff (e.g. 2017), in years with larger meltwater runoff the ice flow slows down a bit earlier (already in July instead of August), the velocities decrease to a lower level and the deceleration is not directly followed 290 by a rapid velocity increase in autumn (leading to type (b) in the following year) (see Fig. 8). Potentially, more meltwater leads to a more efficient subglacial drainage system that prevails for a longer time in the year so that less meltwater is trapped at the be more similar than suggested at a first glance in Fig. 1 and Fig. 6. The low velocities at the beginning of 2020 indicate a transition to a new quiescent phase and, thus, a continuation of the surge behaviour.
As the observed accelerations 1926-1928 and 1954-1956 refer to the average over a period of 2 years, the true maximum velocities are probably significantly higher than the documented velocities of 3.6 and 1 m/day, respectively. In addition, the exact timing of the surges is unknown. They may have lasted longer or shorter than 2 years. In summary, there is no argument 310 that the surge behaviour before 2000 was significantly different from that thereafter.

Classification of flow patterns
We distinguish between four different flow patterns (A-D) based on the spatial distribution of flow velocities. Each pattern is associated with a certain profile of dynamic ice-height changes and, thus, has certain consequences for the stresses and the mass redistribution in the glacier system.

315
(A) A first pattern is characterized by low velocities (< 0.3 m/day) in the lower part of the glacier and moderately higher velocities further upstream. This causes only minor dynamically induced changes in the receiving area and a dynamic thickening in the reservoir area. This pattern can be found for much of the quiescent phases, e.g. in March 2012 (see Fig. 10). (B) The second pattern equals (A), except for an area near the glacier front where the velocities exceed 1 m/day. This pattern arises shortly before the surge initiation in September 2013 (see Fig. 7). Thus, additionally to the thickening in the reservoir 320 area a pronounced thinning occurs directly at the glacier front. (C) A third pattern has a steeply sloped velocity profile (with a maximum at the glacier front) in the lower 10 km of the glacier, whereas larger parts further upstream remain at low velocities similar to patterns (A) and (B). This was the case in December 2013 (see Fig. 7), for instance. This pattern may be associated with only moderate height changes at the glacier front and in the upper parts of the glacier, whereas the middle part of the glacier is rapidly thinning (see Fig. 14). This results in a decrease of surface slope at the terminus and a significant increase of 325 surface slope in the middle part. (D) The fourth pattern mostly reverses pattern (A) so that the glacier thickens in the receiving area and thins in the reservoir area (see Fig. 14). Thus, the effect of the ice dynamics on the glacier geometry during the quiescent phase is compensated. This pattern is shown, e.g., by the velocities in spring and early summer 2014 (see Fig. 7).
During initiating all surges. This sequence reflects an evolution of dynamic changes within the glacier where first, through pattern (B), mass from the middle part of the glacier is mobilized to compensate the mass deficit at the glacier front leading to (C).
Subsequently, the surplus of mass in the upper part of the glacier is set into motion to compensate the mass deficit in both the middle and the lower part resulting in pattern (D).

335
In the years after the surge initiation, the sequence appears to be different. There is a rather direct switch between the (A) and (D), whereas the patterns (B) and (C) are largely absent (see Fig. 15).

Glacier geometry during the quiescent phases
Similar to other surge-type glaciers, the quiescent phase of Harald Moltke Brae is characterized by slower flow in the lower part und faster flow in the upper part. However, the transition between slower and faster flow is rather smooth which implies just a slight dynamic thickening (few metres over three years, see Fig. 13). We therefore hypothesize that for the initiation of the surge, the increase of driving stresses in this transition zone is less important than the decrease of resisting stress 350 at the glacier front implied by the frontal retreat. The findings are opposite down-glacier propagating surges e.g. observed at Variegated Glacier (Raymond, 1987) and at Bivachny Glacier (Wendt et al., 2017). The pronounced thinning at the glacier front is compatible with an up-glacier propagating surge observed at Aventsmarksbrae (Sevestre et al., 2018). Thus, a significantly larger thinning at the glacier front compared to a lesser thickening in the upper part of the glacier could be a key factor for the surge initiation at the terminus and an up-glacier propagation.  (Wendt et al., 2017;Raymond, 1987). This could not be observed at Harald Moltke Brae. However, the presence of a seasonality parallel to surges has been taken as a clear indicator for a hydrologically driven mechanism at other glaciers (Wendt et al., 2017;Raymond, 1987). Mayer et al. (2011) explained the seasonal modulation of a 375 surge by the Kamb drainage-switching theory. We conclude that the hydrological mechanism provides a overall more plausible explanation for the surges at Harald Moltke Brae than a thermally driven mechanism. Additionally, at a thick glacier, the steepening of glacier due to thinning at the terminus could increase the driving stress (Sevestre et al., 2018). The resulting large net stress in flow direction may cause a small area at the glacier front to accelerate (flow pattern B). This induces different effects that facilitate a further acceleration and the upward propagation of the surge: crevassing of the ice surface enabling more meltwater to reach the glacier base (Sevestre et al., 2018), weakening of subglacial till (deformation of the glacier bed), longitudinal stresses (tension) and a further surface thinning/steepening. As a result, the Regarding temporal velocity variability, we identified types of seasonality which indicate a hydrological control of the seasonal velocity changes. The time of a rapid deceleration in July or August suggests a switch between an inefficient and an efficient subglacial drainage system due to the changing amount of meltwater availability. This, however, does not provide an explanation for the significant increase of seasonal amplitude during the surges.

420
Over the past two decades, the overall discharge of the glacier was larger than the overall ice mass accumulation. Thus, the overall thickness of the glacier was decreasing which could be a reaction to the observed long-term retreat of Harald Moltke Brae. However, the glacier does not react instantaneously to the terminus retreat by dynamic thinning. The reason might be a restricted ice flow during the quiescent phases due to internal factors, probably related to the glacial hydrology. This possibly causes an increasing instability and results in an alternation of two reverse feedback mechanisms which involve pronounced 425 acceleration and deceleration, respectively.
By distinguishing between different patterns of spatial distribution of flow velocities we could demonstrate that the surge develops first at the glacier front, and that it is propagating rapidly up-glacier within a few months thereafter. Therefore, we reject those velocity fields with a time basis longer than 60 days in order to avoid strong smoothing effects.
Further, we do not use Landsat ice flow velocity fields with a time basis shorter than 25 days as these are assumed to be of comparatively lower accuracy.

455
Velocities are analysed in form of spatial fields, time series with respect to a point close to the glacier front or velocity profiles along a flow line.

Estimation of monthly velocity fields 460
To ensure consistency when determining a joint time series, all given data grids are transformed to the coordinate system of the Landsat data (Polar Stereographic Grid North). Subsequently, monthly medians are calculated separately for each data set.
The resulting monthly time series are merged by computing the monthly mean for each cell. A similar approach was applied to determine semimonthly velocity fields. and data availability (Morlighem et al., 2017). The mentioned distinct topographic features coincide with transitions between 485 different methods and different associated uncertainty levels (Fig. C1a,c). In particular, different interpolation methods prevail