Subglacial Drainage Evolution Modulates Seasonal Ice Flow Variability of Three Tidewater Glaciers in Southwest Greenland

Surface‐derived meltwater can access the bed of the Greenland ice sheet, causing seasonal velocity variations. The magnitude, timing, and net impact on annual average ice flow of these seasonal perturbations depend on the hydraulic efficiency of the subglacial drainage system. We examine the relationships between drainage system efficiency and ice velocity, at three contrasting tidewater glaciers in southwest Greenland during 2014–2019, using high‐resolution remotely sensed ice velocities, modeled surface melting, subglacial discharge at the terminus, and results from buoyant plume modeling. All glaciers underwent a seasonal speed‐up, which usually coincided with surface melt onset, and subsequent slow‐down, which usually followed inferred subglacial channelization. The amplitude and timing of these speed variations differed between glaciers, with the speed‐up being larger and more prolonged at our fastest study glacier. At all glaciers, however, the seasonal variations in ice flow are consistent with inferred changes in hydraulic efficiency of the subglacial drainage system and qualitatively indicative of a flow regime in which annually averaged ice velocity is relatively insensitive to interannual variations in meltwater supply—so‐called “ice flow self‐regulation.” These findings suggest that subglacial channel formation may exert a strong control on seasonal ice flow variations, even at fast‐flowing tidewater glaciers.


Introduction
The dynamic response of Greenland's tidewater glaciers to changes in environmental conditions remains a key uncertainty in predictions of future sea level rise . Each summer, meltwater produced at the ice sheet surface reaches the ice sheet base, increasing basal water pressure and causing seasonal speed-ups of both land-terminating glacier margins (Davison et al., 2019) and tidewater glaciers (Moon et al., 2014;Vijay et al., 2019). At land-terminating glacier margins, continual subglacial water flow during the summer months causes the formation of hydraulically efficient subglacial channels. These enable the rapid evacuation of meltwater, decreasing basal water pressure, and ultimately cause the overlying ice to decelerate in late summer to speeds slower than those prior to the melt season (an "extra slow-down") (e.g., Sole et al., 2013). In addition, this late-summer extra slow-down scales with meltwater supply such that annually averaged ice velocity is insensitive to interannual variations in meltwater supply-the so-called "ice flow self-regulation" van de Wal et al., 2015).
At tidewater glaciers, however, the relationship between meltwater supply and ice velocity appears to be more complicated. While the majority of tidewater glaciers undergo a seasonal meltwater-induced speedup (an "early-summer speed-up"), only~40% of them (Vijay et al., 2019) experience a seasonal extra slowdown, similar to that of land-terminating margins (so-called "type 3" glaciers; Moon et al., 2014). In contrast, other tidewater glaciers do not undergo an extra slow-down and instead decelerate back to (but not below) pre-melt season speeds ("type 2"; Moon et al., 2014). It has been widely hypothesized (e.g., Bevan et al., 2015;Kehrl et al., 2017;Moon et al., 2014;Sundal et al., 2011;Vijay et al., 2019) that the extra slow-down occurs at type 3 glaciers because they develop efficient subglacial channels during summer (as has been inferred at land-terminating margins). Due to the absence of this compensatory extra slow-down, it has been suggested that type 2 glaciers may accelerate on annual timescales as meltwater supply and the early-summer speed-up increase in a warming climate. However, the difficulty in measuring ice velocity close to tidewater glacier termini at sufficiently high temporal and spatial resolution, and in observing tidewater glacier subglacial drainage systems, has meant that these hypotheses have not been thoroughly tested, leaving a gap in our understanding of how tidewater glaciers may respond to seasonal and longer-term variations in meltwater supply.
Given the difficulty of observing the necessary components of tidewater glacier systems, much has been inferred through comparison with land-terminating sectors of the GrIS, which are more accessible. The similarity in seasonal dynamic behavior of type 3 glaciers to that of land-terminating glaciers (e.g., Bartholomew et al., 2010Bartholomew et al., , 2011Sundal et al., 2011) led Vijay et al. (2019) to suggest that the underlying processes controlling dynamics may be the same. Specifically, the development of hydraulically efficient channels during the melt season is thought to reduce water pressure across large areas of the bed (Hoffman et al., 2016;Sole et al., 2013), leading to the extra slow-down. It is not clear, however, to what extent subglacial channels can form beneath tidewater glaciers. Theoretically, fast-flowing ice and small hydraulic potential gradients expedite subglacial channel closure and promote only slow subglacial channel growth (Röthlisberger, 1972). Therefore, channel formation may be subdued at fast-flowing tidewater glaciers (Doyle et al., 2018;Kamb, 1987), especially where the bed deepens inland. Indeed, distributed near-terminus subglacial drainage systems have been inferred at fast-flowing tidewater glaciers (Bartholomaus et al., 2016;Chauché et al., 2014;Fried et al., 2015;Jackson et al., 2017;Rignot et al., 2015;Slater et al., 2017). On the other hand, runoff-driven buoyant plumes are often visible adjacent to tidewater glacier termini (e.g., Bartholomaus et al., 2016;Jackson et al., 2017;Slater et al., 2017), indicating the efflux of large volumes of meltwater from one or more discrete sources, indicative of efficient subglacial channels. Therefore, it seems that large, efficient channels can form at times beneath some tidewater glacier margins. However, no studies have yet linked the occurrence (or not) of these channels to seasonal changes in ice velocity, which would provide valuable insight into interactions between hydrology and ice dynamics at tidewater glaciers.
Here, we investigate the extent to which seasonal changes in ice velocity are controlled by evolution in the hydraulic efficiency of subglacial drainage at three tidewater glaciers in southwest Greenland (Figure 1), Kangiata Nunaata Sermia (KNS), Narsap Sermia (NS), and Akullersuup Sermia (AS), during 2014-2019. To do this, we derive high-resolution velocity estimates close to the termini of these tidewater glaciers by feature tracking of optical and radar satellite imagery. We compare these time series to observations and modeling of environmental forcings and simple inference of subglacial hydraulic efficiency. The three study glaciers, which are exposed to similar climatic variability, are representative of a spectrum of medium-sized outlet glaciers, with grounding line depths ranging from 60 m at AS to 250 m at KNS and speeds from~2 m d −1 at AS to over 15 m d −1 at KNS. We might therefore expect a range of hydrology-dynamic responses to similar meltwater supply variability, which should be applicable to many other Greenlandic tidewater glaciers.

Ice Velocity 2.1.1. Offset Tracking Procedure
We estimated ice velocity primarily from feature and speckle tracking of Sentinel-1a and Sentinel-1b Interferometric Wide swath mode Single-Look Complex Synthetic Aperture Radar (SAR) amplitude images, acquired using Terrain Observation with Progressive Scans (TOPS). We utilized over 350 Sentinel-1a/b repeat-pass image pairs between January 2015 and June 2019. The majority of pairs acquired since the launch of Sentinel-1b in April 2016 had a 6-day separation period (baseline), while older pairs had a 12-day baseline. In the following, the first and second images acquired in a given pair are referred to as the "master" and "slave" images, respectively. These data were supplemented, and the time series extended back to January 2014 by our own feature tracking of Landsat-8 imagery and the GoLive data set (Fahnestock et al., 2015). Below, we describe the Sentinel-1 and Landsat-8 processing chains.
Prior to tracking, Sentinel image pairs were focused and co-located to within 0.1 pixels in the Generic Mapping Toolbox for SAR imagery (GMTSAR; Sandwell et al., 2011aSandwell et al., , 2011b. Traditional SAR image alignment using image cross-correlation or enhanced spectral diversity fails with TOPS-mode data and is not appropriate over fast-flowing ice, due to coherence loss between images (Nagler et al., 2015;Sandwell et al., 2011b). Image co-location instead utilized precise orbit ephemerides (3-5 cm accuracy) prior to August 2018 and restituted orbit data (10 cm accuracy) afterwards (Fernández et al., 2015; https://qc.senti-nel1.eo.esa.int/) and the Greenland Ice Mapping Project (GIMP) Digital Elevation Model (DEM; Howat et al., 2014), to interpolate the slave image to the grid of the master image. For both Landsat and Sentinel image pairs, each image was split into overlapping image patches. Cross-correlation of corresponding image patches in co-located master and slave images was used to determine ice displacement during the baseline. Due to the different nominal range and azimuth resolution of the Sentinel imagery (~2.3 and~14.1 m, respectively), the Sentinel images were oversampled in the azimuth direction by a factor of two prior to cross-correlation (Khvorostovsky et al., 2018). To minimize information loss over image patch edges, image patches should be approximately four times larger than the maximum expected displacement (Thielicke & Stamhuis, 2014). Therefore, Sentinel-1 image patches were approximately 1 × 1 km (400 × 100 Single-Look oversampled range and azimuth pixels), with approximately 800-m overlap in both directions (300 × 126 pixels). The smaller computational demands imposed by the 15-m Landsat imagery allowed us to perform multiple cross-correlation passes on each image pair, resizing and deforming image patches according the previous pass (Adrian & Westerweel, 2011;Thielicke & Stamhuis, 2014). Image patch size and spacing therefore varied from 480-1,920 and 180-960 m, respectively, equal in both directions. This enabled us to track features close to the terminus without sacrificing either resolution or accuracy further upglacier, where flow speeds are lower. To improve feature identification, we pre-processed images using contrast-limited adaptive histogram equalization and a Butterworth high-pass spatial frequency filter. The latter removed image brightness variations with a wavelength greater than approximately 1 km (de Lange et al., 2007), ensuring that only smaller, moveable surface features were tracked.
Tracking of the co-located and filtered images was undertaken in MATLAB, within a version of PIVsuite (Thielicke & Stamhuis, 2014; https://uk.mathworks.com/matlabcentral/fileexchange/45028pivsuite) adapted for ice flow. Computationally efficient subpixel displacement estimates for each image patch were made by obtaining an initial estimate of the cross-correlation peak using a fast Fourier transform and then upsampling by a factor of 100 the discrete Fourier transform using matrix multiplication of a small neighborhood (1.5 × 1.5 pixels) around the original peak location (Guizar-Sicairos et al., 2008). The resulting velocity estimates were filtered in several stages ( Figure S1). Correlations with a signal-to-noise ratio (defined as the ratio of the primary cross-correlation peak to the average of the remaining cross-correlation field) less than 5.8 were removed (de Lange et al., 2007). Remaining spurious estimates were removed primarily using an image segmentation filter, a threshold strain filter (Rosenau et al., 2015), and a kernel density filter based on the paired displacements in the range and azimuth directions for each image patch (Adrian & Westerweel, 2011). Additional filtering based on velocity magnitude and flow direction removed remaining spurious estimates. The filtered velocity fields derived from Sentinel-1 imagery were transformed from radar to map coordinates using the GIMP DEM (Howat et al., 2014) and were posted on a 150 × 150 m grid.
For the analysis presented here, we averaged ice velocity within 1 × 1 km fixed-position regions of interest, located close to the calving front of each glacier (locations in Figure 1). Only dates on which more than 50% of the region of interest contained data were sampled. A median velocity error of 0.06 m d −1 was estimated by sampling velocity over bedrock areas.

Validation of Velocity Estimates
To validate our method of velocity estimation, we compared our ice velocity estimates derived from Sentinel-1 imagery to those from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) Sentinel-1 data set (https://www.promice.dk/) over identical time periods ( Figure S2). We find that our estimates differed from those of PROMICE by 0.25% on average, with a standard deviation of 4.4% and that this difference was independent of ice velocity. The small, non-systematic difference between our estimates and those of PROMICE arise due to (i) the different spacing at which the data are binned (150 m for our data and 500 m for PROMICE); (ii) differences in the size and spacing of image patches used to estimate ice velocity, which will affect the locations over which velocity is estimated; and (iii) differences in the degree of post-processing (Merryman Boncori et al., 2018).

Terminus Position and Ice Mélange
For each glacier, we digitized terminus positions from Landsat-8 and Sentinel-2 imagery during 2014-2019 using the Google Earth Engine Digitisation Tool (Lea, 2018). Terminus position change was calculated using the multicenterline method in the Margin change Quantification Tool (Lea, 2018). Ice mélange, and (at KNS) a seasonal ice tongue (Motyka et al., 2017;Moyer et al., 2017), hindered accurate terminus identification typically during November-May each year.
The timings of ice mélange (and, at KNS, ice tongue) weakening, disintegration, and reformation were identified based on our ice velocity estimates over the mélange and using visual assessment of satellite imagery. We identified mélange weakening from crack formation, changes in color of the mélange (which we interpreted as thinning), and mobility (based on velocity estimates or loss of coherence). Mélange disintegration was defined as complete clearing of the mélange from the fjord adjacent to the calving front. These events likely bracket the period during which buttressing forces provided by ice mélange decreased each year. During times of cloud cover (when Sentinel-1 images were not available) or image sparsity, we recorded the date of the first clear image.

Surface Melt and Subglacial Discharge
We extracted the time series of average RACMO2.3p2 modeled surface runoff rates (Noël et al., 2016(Noël et al., , 2018 within the regions of interest on each of our study glaciers (black boxes in Figure 1). Using these, we defined the start of the melt season for each glacier as the first day of a period of at least three consecutive days when the runoff rate was greater than 1-mm water equivalent per day (i.e., following the Danish Meteorological Institute definition). To gain further insight into surface meltwater generation, we also analyzed air temperature data acquired at PROMICE weather station NUK_L (https://www.promice.dk/WeatherStations. html) and defined the onset of positive temperatures as the first day of a period of at least three consecutive days when temperatures exceeded 0°C. We also estimated subglacial discharge at the terminus of each glacier using RACMO2.3p2 modeled surface runoff, which was spatially and temporally integrated over each glacier's subglacial catchment, delineated using hydropotential analyses (Shreve, 1972) bounded by BedMachine v3 topographic data (Morlighem et al., 2017). For simplicity, surface runoff was assumed to access the bed immediately, was routed to the terminus at 1 m s −1 (Chandler et al., 2013;Cowton et al., 2013), and used to estimate mean daily subglacial discharge rates.
Our time series of subglacial discharge at the glacier terminus was derived by making several simplifying assumptions. To summarize, we assume (i) that meltwater accesses the bed immediately (no supraglacial or englacial storage), (ii) that meltwater accesses the bed where it is generated (no supraglacial routing), and (iii) a subglacial transit velocity of 1 m s −1 . We quantify and discuss the impact of these assumptions on our conclusions in Text S1.

Subglacial Hydraulic Efficiency
We used the visible presence or absence of plumes at the fjord surface adjacent to each glacier, combined with simple plume modeling, as an indicator of near-terminus subglacial hydraulic efficiency. Plumes are often, but not always, visible at the fjord surface during summer. Previous modeling work  demonstrates that typically, relatively little subglacial discharge (<50 m 3 s −1 ) from a single channel is required to cause plume surfacing at these glaciers. Therefore, when modeled subglacial discharge is high, yet no plume is observed, one possible explanation is that subglacial discharge emerged from multiple points along the terminus, such that the discharge at each outlet was less than~50 m 3 s −1 . Although this does not provide direct information about the efficiency of the near-terminus subglacial drainage system, the spatial distribution of water efflux across the grounding line is suggestive of an "inefficient" near-terminus subglacial drainage system (e.g., Slater et al., 2017).
We recorded the presence of subglacial discharge plumes at the fjord surface using all available Landsat-8, Sentinel-2, and Sentinel-1 satellite imagery during 2014-2019 ( Figure S3). For simplicity, we adopted a binary system to classify near-terminus subglacial hydraulic efficiency. When plumes were visible at the fjord surface, we assumed that the subglacial drainage system was efficient. When plumes were not visible, we used buoyant plume theory (Jenkins, 2011;Morton et al., 1956;Slater et al., 2016) to estimate the minimum number of outlets required to prevent plume surfacing (assuming discharge was split evenly between outlets as in Slater et al., 2017). If two or more outlets were required to prevent plume surfacing, we assumed an "inefficient" drainage system. We emphasize that we cannot provide more specific information on the likely morphology of these "inefficient" systems and that, under the terms of our classification, "inefficient" does not preclude the existence of subglacial drainage channels, so long as not more than 50% of the total discharge (and normally much less) is carried by a single channel. Very little is known about drainage system morphology near the termini of tidewater glaciers, but we speculate that this definition of inefficient could include anything from a linked cavity network or porous till through to a network of multiple transient channels.
We forced the plume model (Slater et al., 2016) with our time series of grounding line subglacial discharge, while 28 conductivity, temperature, and depth (CTD) casts, acquired 32-90 km from the KNS terminus during 2014-2016 (http://ocean.ices.dk), were used as ocean boundary conditions (Figure 1; Mortensen et al., 2013Mortensen et al., , 2014Mortensen et al., , 2018. The results presented below assume a half-conical plume geometry (consistent with the geometry of plume surface expression). Changes to model boundary conditions (including subglacial discharge) and parameters, within parameter uncertainty, resulted in only minor changes to the timing of periods of inferred efficient and inefficient subglacial drainage (Text S1; Figures S4 and S5).
On days when it was not possible to make an inference about the efficiency of the subglacial drainage system, we recorded data gaps. These gaps are due to (i) an absence of satellite images, (ii) mélange or cloud cover, or (iii) insufficient subglacial discharge for a single modeled plume to reach the fjord surface. In addition, there may be periods erroneously classified as "inefficient" during times when (i) plumes did reach the fjord surface, but were not visible in satellite imagery (e.g., due to the plumes being below the resolution of the satellite image), or (ii) there was some other reason for lack of plume surfacing, such as freshening of the fjord surface due to iceberg melt and surface runoff (De Andrés et al., 2020). We further discuss the assumptions and sensitivities of this method in Text S1.

Results
We observed seasonal ice velocity variations at each glacier that were qualitatively similar to the type 2 and type 3 patterns identified by Moon et al. (2014) and Vijay et al. (2019). We compared these ice flow variations to changes in terminus position, ice mélange characteristics, runoff, and inferred subglacial hydraulic efficiency. Overlaid is a 24-day Butterworth low-pass filtered velocity time series, colored red (blue) when it is higher (lower) than the pre-melt season speed. The vertical violet (gray) lines indicate times of inferred inefficient (efficient) subglacial drainage. Areas with no vertical lines indicate that either no images were available or that modeled subglacial discharge was zero. The vertical dash-dot gray lines delimit each calendar year and are for visual guidance only. (b) Width-averaged terminus position in red crosses (lower values indicate a more retreated position). The horizontal bars indicate mélange presence: cyan sections indicate a "strong" mélange while blue sections indicate transition periods, when the mélange was present but appeared to be weakening or reforming. (c) Modeled subglacial discharge at the terminus (red) and modeled surface runoff (gray). (d) Air temperature from PROMICE station NUK_L and the timing of key events (black cross ¼ acceleration; red bar ¼ positive temperature onset; red downward pointing triangle ¼ melt onset; green triangle ¼ terminus retreat onset; blue cross ¼ mélange weakening; blue circle ¼ mélange breakup).

Akullersuup Sermia
AS began to accelerate between mid-March and early May each year (henceforth the "early-summer acceleration"), reaching peak speeds within about a month (Figure 2a). Peak speeds were on average 18% greater than those prior to the early-summer acceleration (henceforth "pre-acceleration" speeds). Ice flow then typically decelerated rapidly, falling to approximately 40% below pre-acceleration speeds by midsummer and then remained low until early September each year, before gradually accelerating over winter (henceforth "recovery"). This seasonal pattern was observed throughout the topographically constrained (outlet) part of the glacier (Figures 3 and S6-S8). Short-lived speed-ups, coincident with spikes in modeled surface melt, were superimposed on this seasonal pattern (e.g., November 2017). Average ice velocity between April and March 2016-2017 was 4.7% greater than during 2017-2018 and 0.2% lower than during 2018-2019 (Table S1).
The early-summer acceleration was usually difficult to distinguish from the gradual acceleration over the preceding winter, and so it was difficult to associate the onset of acceleration with a particular forcing ( Figure 2). In most years, the early-summer acceleration seemed to begin before any of our defined forcings; however, there were short-lived periods of above-zero air temperatures during or just prior to the onset of acceleration in every year ( Figure 2). Peak speeds occurred within 1-2 weeks of the first observation of plume surfacing during every year except 2015, when plume surfacing occurred~1 month later. Ice velocity subsequently decreased and usually remained low until plume surfacing ceased around September each year. The uniformity of the seasonal velocity variations across the outlet part of the glacier (Figures 3 and  S6-S8) suggests a spatially consistent control that was not affected by proximity to the glacier terminus. We found no correlation (R 2 ¼ 0.006, p ¼ 0.19) between ice velocity and terminus position during our study period (Figure 4a).

Narsap Sermia
NS displayed qualitatively similar behavior to AS, but the relative magnitude and duration of the early-summer acceleration and the late-summer slow-down differed markedly (Figure 5a). In every year except 2017, the early-summer acceleration began around the time of both runoff onset and terminus retreat, but before visible ice mélange weakening. In 2017, the early-summer acceleration is indistinguishable from the 2016/2017 winter recovery, with ice velocity steadily increasing from around early March 2017, a time with frequent excursions to positive temperatures (Figure 5d). Peak speeds were up to 25% greater than pre-acceleration values, and velocity remained elevated relative to pre-acceleration values for 2-3 months.
Beginning in midsummer, ice flow decelerated toward a minimum in early September that was approximately 10% below pre-acceleration velocities. The seasonal transition from accelerating to decelerating ice flow coincided closely with a switch to inferred efficient drainage (signaled by plume surfacing) in all years except 2017. In every year, this summertime deceleration occurred despite continued terminus retreat. After reaching a velocity minimum at the end of the melt season, ice velocity gradually recovered each winter ( Figure 5) but remained below pre-acceleration speeds for the majority of each winter. As with AS, this seasonal pattern was similar throughout the outlet part of NS (Figures 6 and S6-S8). An exception to this gradual winter recovery occurred during January 2019, when we observed a short-lived (1-2 weeks) and high-magnitude (~25%) speed-up throughout the outlet part of the glacier.  (Table S1). We observe moderate positive correlations between glacier terminus position and ice velocity during individual years (R 2 ¼ 0.13-0.68, p < 0.007), but a weak positive correlation when considering multiple years (R 2 ¼ 0.1, p < 0.001) (Figure 4b).

Kangiata Nunaata Sermia
KNS was characterized by seasonal velocity variations that resembled the type 2 behavior described by Moon et al. (2014) and Vijay et al. (2019). The early-summer acceleration coincided closely with surface runoff and/or positive temperature onset, which was usually several weeks prior to observed terminus retreat and visible mélange weakening. The summer speed-up was both more pronounced (peaking up to 40% above pre-acceleration values) and more sustained (ice velocity remained elevated relative to pre-acceleration speeds for at least the entire melt season each year) than at either AS or NS. During 2016-2018, when our velocity and plume data are most complete, the transition from accelerating to decelerating flow occurred with, or shortly after, an increase in inferred drainage system efficiency ( Figure 7a) and regardless of whether the terminus was still retreating. In contrast to AS and NS, there was little or no extra slow-down in any year. An exception to this occurred in October 2016, when ice velocity in the lower 6 km of KNS briefly dipped below pre-acceleration speeds (i.e., there was a minor extra slow-down) following the drainage of the large ice-dammed lake Isvand, which produced a plume adjacent to KNS (Figures 1, 7a, 8, and S9). Like AS and similar to NS, ice velocity was unrelated (R 2 ¼ 0.002, p ¼ 0.7) to terminus position during the study period ( Figure 4c

Discussion
We observed pronounced and differing seasonal velocity variations at three contrasting tidewater glaciers exposed to similar climatic variability. At AS and NS each year, we observed an early-summer acceleration,  subsequent deceleration to below pre-acceleration speeds, and gradual acceleration over winter (Figures 2, 3, 5, and 6). Although there were differences between these glaciers (discussed below), the seasonal ice velocity pattern of both AS and NS fits the "type 3" pattern identified in Moon et al. (2014). In contrast, KNS did not usually undergo an extra slow-down and so broadly fits the type 2 classification, with perhaps some type 3 behavior (Moon et al., 2014). Our more detailed observations therefore support the classification of AS and KNS by Moon et al. (2014). At most tidewater glaciers, changes in terminus position or subglacial hydrology are thought to be the dominant drivers of seasonal dynamics (e.g., Moon et al., 2015), but disentangling these contrasting processes is difficult (e.g., Fried et al., 2018). In the discussion below, we argue that evolution in subglacial hydraulic efficiency can explain the key features of the seasonal ice flow variations at these tidewater glaciers. We therefore build on previous studies (Moon et al., 2014;Vijay et al., 2019) by providing additional evidence to support the hypothesis that subglacial drainage evolution both occurs and exerts an important control on ice dynamics at tidewater glaciers. Furthermore, our time series enable a more detailed description of the subtle variations in these temporal velocity patterns both between glaciers and between years, allowing a more robust interpretation of their controls.
Each of the three glaciers studied underwent a temporary speed-up each year, usually commencing between mid-March and mid-May (Figures 2, 5, and 7). This acceleration always began before any visible mélange weakening, indicating that changes in buttressing by ice mélange do not serve as a key control on the seasonal dynamics of these glaciers. In many cases (especially at KNS but also at NS), the onset of acceleration coincided approximately (at the resolution of the data) with the onset of surface runoff and/or positive temperatures. Sometimes, this also coincided with terminus retreat, which may have contributed to the observed acceleration (Fried et al., 2018). At AS, acceleration usually occurred prior to any obvious forcing (Figure 2). We suggest that this is partly due to the difficulty in distinguishing the early-summer acceleration from the ice flow recovery over the preceding winter. In addition, it is possible that surface melting on the lower part of the glacier, caused by brief excursions to positive temperatures, did occur in sufficient volume to affect ice dynamics, but was not captured by RACMO2.3p2.
At all three glaciers, ice velocity was generally greatest near the beginning of the melt season, when meltwater runoff was rising rapidly. This behavior resembles that of land-terminating glaciers and occurs when the drainage system is continually challenged by rapidly increasing meltwater inputs, causing frequent spikes in water pressure (Bartholomew et al., 2012;Harper et al., 2007;Schoof, 2010), cavity expansion Iken, 1981;Kamb, 1987), and/or sediment deformation (Iverson et al., 1999). Our observations of seasonal meltwater-induced speed-ups were relatively small (16-40%) compared to land-terminating glaciers (180-400%; Sole et al., 2013;van de Wal et al., 2008van de Wal et al., , 2015, though the maximum speed-ups we observe are likely reduced by smoothing over the 6-to 12-day image baseline. Similarly, modest seasonal meltwater-induced speed-ups (typically less than 15%) have been observed at several other Greenlandic tidewater glaciers (Ahlstrøm et al., 2013;Andersen et al., 2010;Bevan et  and may be subdued relative to land-terminating glaciers because of the already low basal resistance at tidewater glaciers (Shapero et al., 2016).
As subglacial channels become larger and more efficient, further increases in meltwater input have a more limited impact on basal water pressure, so velocity stabilizes then falls (Röthlisberger, 1972). Such behavior has been observed across land-terminating sectors of the ice sheet (Bartholomew et al., 2012;Cowton et al., 2016). Although direct evidence of subglacial drainage evolution is even more challenging to obtain at tidewater glaciers than at their land-terminating counterparts, the velocity and plume observations and plume modeling presented here provide indirect support for a similar process at our study glaciers. At all glaciers, velocity decreased part way through the melt season, but this occurred later at the faster-flowing KNS and NS than at AS. In addition, there was in most years a seasonal progression toward more efficient drainage, as inferred from more frequent plume surfacing. In some cases (e.g., NS in 2016 and KNS in 2017), the appearance of plumes at the fjord surface coincided closely with the transition from acceleration to deceleration. This provides further support for the role of changes in subglacial hydraulic efficiency in modulating the seasonal velocity patterns at these glaciers, though this evidence is treated with caution given the approximate nature of the method used to infer hydraulic efficiency and the interannual variation in the timing of observed plume surfacing with respect to peak velocities.
After peak velocities were reached at AS and NS each year, ice velocity fell to below pre-acceleration values, despite continued surface melting, and was followed by recovery over winter. The deceleration occurred earlier and seasonal velocity minima were lower relative to pre-acceleration speeds at AS than at NS. This pattern of flow variability provides further support for the hypothesis that subglacial drainage evolution modulated the observed seasonal flow variations at these glaciers. At land-terminating margins, the extra slow-down is thought to be caused by drainage of water from weakly connected areas of the bed (Andrews et al., 2014;Hoffman et al., 2016) toward persistent and efficient subglacial channels . The extra slow-down is therefore more pronounced where the late-summer subglacial channels are more hydraulically efficient . Based on this understanding, we propose that the seasonal subglacial drainage system at AS became more efficient than at NS and KNS because the extra slow-down was most pronounced at the former. Furthermore, flow recovery at both AS and NS did not begin until the end of each melt season, presumably only after subglacial channels had closed, thereby allowing re-pressurization of the subglacial drainage system by basal melting and rain events (e.g., November 2017).
KNS displayed little or no seasonal extra slow-down near the terminus, suggesting that the subglacial drainage system was not usually efficient enough to induce a widespread reduction in basal water pressure to below pre-acceleration values. However, two sets of observations indicate that seasonal increases in subglacial hydraulic efficiency still acted to dampen the seasonal speed-up each year, thereby limiting interannual velocity. First, the transition from acceleration to deceleration each summer usually closely followed a switch to predominately channelized drainage (as inferred from the plume observations and modeling). Moreover, this transition began earlier and the subsequent deceleration was faster in 2016, when the melt season began earlier and modeled meltwater discharge was greater. Thus, although there was little seasonal meltwater-induced extra slow-down (as occurred at AS and NS), annual average ice velocity at KNS was still lower during 2016, the year with the greatest runoff (Table S1). Second, following the drainage of Isvand in October 2016 ( Figure S9), we observed a plume at the fjord surface and a short-lived velocity perturbation characterized by a significant speed-up followed (crucially) by an extra slow-down and subsequent recovery (Figures 7 and 8). These observations suggest that even short-lived pulses in runoff supply can potentially form channels efficient enough to induce a compensatory slow-down. This raises the possibility that as runoff supply increases in the future, current type 2 glaciers may transition toward type 3 behavior.
There are differences in the seasonal velocity patterns of each glacier that we argue provide further insight into the evolution of the hydrological system beneath tidewater glaciers. For example, summertime peak speeds at AS were much smaller relative to pre-acceleration speeds and occurred much earlier, than at NS. Those at NS were in turn smaller and occurred earlier than those at KNS. One possible explanation for this is that subglacial channels developed earliest (and grew fastest) at AS and latest (and grew slowest) at KNS. If our assumption that plume surfacing indicates the presence of an efficient subglacial drainage system is correct, then our time series of plume surfacing supports this explanation because plumes appeared first and were most persistent at AS and appeared last at KNS. Given the difficulty of directly observing 10.1029/2019JF005492

Journal of Geophysical Research: Earth Surface
subglacial channel development at tidewater glaciers, we cannot prove that efficient subglacial channels do develop more rapidly at AS. However, it is worth briefly considering the likely theoretical conditions for channel development at each glacier ( Figure 9). AS is grounded throughout, with a height above flotation of~100 m near the terminus. In contrast, the lower 1.5 km of NS is potentially at flotation, which would reduce hydraulic potential gradients and therefore reduce subglacial channel growth rates (Röthlisberger, 1972). This suggests that readily available metrics such as height above flotation may provide some insight into the influence of hydrology on the dynamics of other tidewater glaciers. Under this reasoning, however, it is not clear why the seasonal velocity patterns of KNS and NS differ so markedly, since their flow speed and height above flotation are similar.
As AS is well grounded and flows only 2-3 times faster than is typical of land-terminating sectors, it is perhaps unsurprising that seasonal velocity patterns, and inferred drainage system evolution, resemble those at land-terminating glaciers. However, it is notable that qualitatively similar seasonal velocity variations occurred at NS, which at~12 m d −1 is flowing 3-4 times faster than AS and an order of magnitude faster than is typical of land-terminating glaciers. Therefore, although channel formation is theoretically hindered by fast ice flow and weak hydraulic potential gradients (Kamb, 1987), our observations suggest that efficient subglacial drainage configurations are able to form and persist long enough to modulate the relationship between meltwater runoff and ice velocity even at fast-flowing tidewater glaciers. Furthermore, our data indicate that increasing hydraulic efficiency during the melt season can dampen or offset meltwater-induced speed-ups even at fast-flowing tidewater glaciers like NS and KNS, regardless of whether the seasonal velocity pattern is "type 2" or "type 3." If we accept that the dynamic behavior of these glaciers is qualitatively consistent with ice flow selfregulation, controlled by changes in subglacial hydraulic efficiency van de Wal et al., 2015), we would expect that there would be limited sensitivity in annually averaged ice motion to interannual variations in runoff. Unfortunately, however, our time series of unbroken velocity estimates is not long enough to allow us to confirm this behavior (Table S1). Nevertheless, meltwater runoff is believed to be a significant influence on tidewater glaciers via other mechanisms, particularly as a driver of submarine melting at their termini Straneo et al., 2011;Xu et al., 2012). Over longer timescales, increased runoff may therefore lead to increased submarine melting, tidewater glacier retreat (Cowton et al., 2018; Slater  (Nick et al., 2009;Pfeffer, 2007). We therefore emphasize that our findings do not preclude a potentially important role for increased meltwater runoff in the dynamic evolution of tidewater glaciers and it remains a subject requiring further investigation to determine more precisely the extent and rate of subglacial channel development, and the efficacy of ice flow selfregulation, at contrasting tidewater glaciers.
Our results contrast with studies that have identified terminus retreat and mélange disintegration as key drivers of seasonal flow variations (Howat et al., 2010;Moon et al., 2015). While we argue that the evidence presented here strongly supports the hypothesis that the seasonal velocity variations at our study glaciers are modulated by subglacial hydrology, it is possible that changes in terminus position and mélange buttressing also play a role. For example, a proportion of the early-summer acceleration at each glacier (especially KNS) occurred while the terminus was retreating each year (similar to Fried et al., 2018). However, despite the partial temporal overlap between acceleration and retreat, there are several lines of evidence indicating that variations in terminus position and mélange buttressing did not exert a dominant control on either the magnitude or timing of the velocity variations observed here: (1) the onset of acceleration always occurred several weeks prior to observed terminus retreat or visible mélange weakening; (2) peak velocity always occurred before the most retreated terminus position; (3) gradual winter acceleration began before mélange reformation and continued after mélange reformation; and (4) over interannual timescales, we find little relationship between ice velocity and terminus position at any of our study glaciers (Figure 4). At NS, we observed deceleration during terminus retreat, resulting in a positive correlation between ice velocity and terminus position during individual years (Figure 4b)-the opposite to that expected if terminus position was driving the changes in ice velocity.
It is likely that the drivers of seasonal ice flow variations differ between glaciers, and we note that while KNS is the largest tidewater glacier in southwest Greenland, our study glaciers do not represent the full range of Greenland's tidewater glaciers. At glaciers that are faster, thicker, and more lightly grounded grounded than KNS, basal friction may be lower (Shapero et al., 2016) and calving events tend to be larger, thereby potentially increasing the importance of terminus position changes over subglacial hydrology as a control on ice dynamics (e.g., Bevan et al., 2015;Kehrl et al., 2017). Nevertheless, an examination of a broader sample of glaciers (Vijay et al., 2019) showed that over 50% were characterized by seasonal ice flow variations similar to those observed here and which we argue are controlled primarily by subglacial hydrology. Therefore, while our study has focused in detail on a few glaciers within a single fjord system, we expect our findings concerning the role of subglacial hydrology in driving seasonal ice flow variability to be relevant to a large number of small to medium-sized tidewater glaciers around Greenland.

Conclusion
We use high-resolution ice velocity estimates, observations of terminus position and ice mélange, modeled subglacial meltwater discharge, and inferred subglacial hydraulic efficiency to investigate drivers of seasonal ice flow variability of three contrasting tidewater glaciers, with similar climatic forcing, in southwest Greenland. At all three glaciers, we find little relationship between ice velocity and variations in terminus position or ice mélange occurrence. Instead, we infer that surface-derived meltwater inputs drive pronounced seasonal changes in ice velocity characterized by early-summer flow acceleration followed by deceleration either to, or below, pre-acceleration speeds. We argue that the amplitude and longevity of the seasonal acceleration and deceleration are controlled by the development of hydraulically efficient subglacial channels. We suggest that this behavior is qualitatively consistent with ice flow self-regulation (where in warmer years, faster summer ice flow is balanced by slower winter motion, resulting in limited net annual differences in ice motion), which has been observed over extensive land-terminating sectors of the GrIS but not near the termini of tidewater glaciers. Therefore, changes in subglacial hydraulic efficiency likely exert a strong control on the seasonal dynamics of many of Greenland's small to medium-sized tidewater glaciers. The net impact of this hydrodynamic coupling on annual and interannual timescales nevertheless remains uncertain and requires further investigation.