Proper orthogonal decomposition of ice velocity identifies drivers of flow variability at Sermeq Kujalleq (Jakobshavn Isbræ)

The increasing volume and spatio-temporal resolution of satellite-derived ice velocity data has created new exploratory opportunities for the quantitative analysis of glacier dynamics. One potential technique, Proper Orthogonal 10 Decomposition (POD), also known as Empirical Orthogonal Functions, has proven to be a powerful and flexible technique for revealing coherent structures in a wide variety of environmental flows. In this study we investigate the applicability of POD to an openly available TanDEM-X/TerraSAR-X derived ice velocity dataset from Sermeq Kujalleq (Jakobshavn Isbræ), Greenland. We find three dominant modes with annual periodicity that we argue are explained by glaciological processes. Mode 1 is interpreted as relating to the stress-reconfiguration at the glacier terminus, known to be an important control on the 15 glacier’s dynamics. Modes 2 and 3 together relate to the development of the spatially heterogenous glacier hydrological system and are primarily driven by the pressurisation and efficiency of the subglacial hydrological system. During the melt season, variations in the velocity shown in Modes 2 and 3 are explained by the drainage of nearby supraglacial melt ponds, as identified with a Google Earth Engine MODIS dynamic thresholding technique. By isolating statistical structures within velocity datasets, and through their comparison to glaciological theory and complementary datasets POD indicates which glaciological 20 processes are responsible for the changing bulk velocity signal, as observed from space. With the proliferation of optical and radar derived velocity products (e.g. MEaSUREs/ESA CCI/PROMICE) we suggest POD, and potentially other modal decomposition techniques, will become increasingly useful in future studies of ice dynamics.


Introduction
The surface flow of glaciers is one of the most important measurements to assess the health of the cryosphere: critical to 25 understand its mass balance and changing dynamics. It is an emergent property resulting from the complex interaction of numerous processes occurring on a variety of spatial and temporal scales. Prising apart the influence of the physical processes responsible for ice dynamic patterns has been a key aim in glaciology for decades and remains challenging for researchers.
Technological advances and a growing appreciation of the need to resolve a wider range of glacier behaviour has led to satellite-derived ice velocity images being acquired with increasingly regularity (Joughin et al., 2020a;Vijay et al., 2019;Lemos et al., 30 2018;Fahnestock et al., 2016), which, in turn, has allowed for the development of new analytical techniques and treatments of datasets (Riel et al., 2021;Greene et al., 2020).
One potential technique, Proper Orthogonal Decomposition (POD), also known as Empirical Orthogonal Functions (EOFs) within meteorology and oceanography, has proven to be a powerful and flexible technique for revealing coherent structures in 35 a wide variety of environmental flows. POD may have value in glaciology because it exactly describes a series of snapshots from a flow field with the product of ranked spatially orthogonal "modes" of spatial weighting, and one-dimensional "temporal" coefficients (eigenvectors). Importantly, in many cases the variance of the flow field is well described by just a few dominant orthogonal modes which, potentially, relate to distinct glaciological flow mechanisms. In this paper we explore the applicability of this eigen-decomposition technique to a series of colocated satellite-derived ice velocity images for the 40 first time. We use openly available ice velocity datasets and focus on the exceptionally well-sampled trunk of a major marineterminating glacier in West Greenland.
Greenlandic marine-terminating ("tidewater") glaciers are an intriguing test area for our analysis because they account for the majority of the ice sheet's contribution to global sea level rise , and anticipating their future behaviour 45 is critical to improving our understanding of the health of the ice sheet (Choi et al., 2021;Catania et al., 2019). However, understanding the spatiotemporal dynamics of marine-terminating glaciers remains challenging; in part due to the competing and potentially interrelated atmospheric, geological, oceanographic, and glaciological factors affecting glacier geometry, flow and frontal ablation (e.g. Fahrner et al., 2021;King et al., 2020;Slater et al., 2019;Brough et al., 2019, Bevan et al., 2019Felikson et al., 2017;Straneo et al., 2010;Holland et al., 2008). In addition to the widespread, averaged multi-year trend of 50 terminus retreat, the seasonal velocity patterns of these glaciers are often classified as being dominated by one of three types (Vijay et al., 2019;Moon et al., 2014), each relating to a specific dominant mechanism (e.g. Howat et al., 2010): (1) systems that are dominated by frontal stress-reconfigurations, and summer velocity is strongly linked to terminus position; (2) systems where the supply of meltwater to the ice-bed interface increases subglacial water pressure and enhances basal motion, and peak velocity occurs close to peak melt; and (3) where the subglacial drainage system adapts to the influx of subglacial water 55 and develops from an inefficient, high-pressure, low-transit-time drainage system to an efficient, channelised system, and velocity tends to decrease throughout the melt season. This exemplifies the challenge of separating the relative importance of different flow mechanisms and forcings from the observed bulk signal: usually examined as points on a centreline or averaged over a region of interest.

60
Through the analysis presented in this paper we find that distinct spatial and seasonal patterns of ice velocity exist within the ice velocity dataset. We interpret the most dominant mode as being primarily due to glacier stress reconfiguration due to variations of frontal ablation. We then hypothesise that the smaller, subtler second modal signal is due to glacier hydrological forcing. By analysing timeseries of nearby supraglacial melt ponds, several published glaciological and environmental datasets, and errors associated with the input data we build a case for this interpretation and, therefore, the broader applicability of POD 65 to ice velocity datasets.

Proper Orthogonal Decomposition
POD is based on the decomposition of a geophysical field into a series of ranked eigenvectors, spatially orthogonal modes, and singular values (related to eigenvalues) that each partially explain the dataset variance, analogous to principal component 70 analysis. The technique was originally derived to identify coherent structures, i.e. zones in which velocity fluctuations are correlated, in turbulent flow (Lumley, 1967). More recently, however, the ability of a POD to extract spatial pictures of regions of coherence has resulted in successes in relatable fields, e.g., fluvial hydrodynamics (Higham, 2017), granular rheology  and solar flare magnetohydrodynamics (Albidah, 2020).

75
The principle of a POD is that a series of snapshots of a centred (i.e. temporal mean removed) flow field can be exactly reconstructed by the product of spatially orthogonal "modes" of spatial weighting by their singular values, and one-dimensional "temporal" coefficients (eigenvectors). A POD places no restriction on the temporal coefficients, such that these can be periodic, aperiodic, and with amplitudes which vary over time. Typically, in a fully converged dataset, a large proportion of the variance of the flow field can be explained by a few dominant modes. The orthogonal nature of each mode means that 80 spatially coherent forcing mechanisms on glacier motion may be identifiable by examining each mode individually. Strongly recurrent spatial structures relating to a consistent forcing typically manifest clearly in POD modes, however more spatially variable dynamics are commonly captured with POD but "split" over two or more modes. If the time period over which the POD is performed is long, then these signals are likely to be uninterpretable due to being spread over many modes that individually appear incoherent. If the process is split over just two modes these are described as "paired" quasi-conjugate 85 modes and are identifiable by displaying common but phase-shifted features within their eigenvectors, and similarly-valued eigenvalues that disrupt the otherwise typical, smooth decay of eigenvalues with increasing mode number. Incoherent lowmagnitude signals due to, for example, measurement or environmental noise are accommodated for in higher order modes. For a more detailed discussion readers are directed to several detailed reviews of POD and related techniques: Weiss (2019), Jolliffe & Cadima (2016), Navarra & Simoncini (2010). 90 Previously, POD has been applied to a variety of polar science problems (Datta et al., 2019;Bian et al., 2019;Mernild et al., 2017;Nield et al, 2012). In ice dynamic studies, however, the technique has been underutilised. To our knowledge its only published application to a time series of ice-surface velocity is by Mair et al. (2002) who decomposed measurements of ice velocity from a network of survey stakes on Haut Glacier d'Arolla during the 1995 melt season. Their analysis clearly indicated 4 that although variations in the local driving stress dominated spatial variation in surface velocity, hydrologically-induced flow patterns were identified that defined the declining spatial extent of enhanced basal sliding in response to meltwater inputs as the subglacial drainage system's hydraulic efficiency evolved. In particular, they found that Mode 2 corresponded to the known location of an incipient subglacial drainage channel. The timing of a change in sign of the coefficient of this mode corresponded to the timing of a switch from inefficient, distributed drainage to efficient, channelised drainage, as deduced from results of 100 multiple dye tracing experiments (following Nienow et al, 1998). Mode 2 was therefore interpreted to represent the changing basal sliding response to meltwater inputs to an evolving subglacial drainage system, switching from a high pressure, high basal sliding response, to a low pressure, low sliding response. This early work provides a strong scientific motivation to explore the applicability of POD to modern, high-resolution satellite derived ice velocity datasets.

105
To determine the POD modes we calculate the eigenvectors and singular values (related to the eigenvalues) using an economy singular value decomposition (SVD) (i.e. only solving the number of timesteps, T, eigenvectors). We create a collection of column vectors (n, t) by column vectorising each of the individual i × j image snapshots where n = 1,2, … N = ij and t = 1,2 … T are the number of equally spaced time-ordered snapshots. The SVD of the matrix is given by: where, * = ( * ) * , (2) * = ( * ) * , and * = * = , (4) 115 where * is a conjugate transpose and is the identity matrix. From this decomposition ∈ ℝ N×T relates to the orthogonal spatial modes, each of ∈ ℝ T×T (the eigenvectors) relate to the temporal evolution of the spatial modes, and both the spatial modes and temporal coefficients are ordered by the ordered singular values contained by the diagonal of ∈ ℝ T×T . The eigenvalues associated with each mode are ( ) 2 . As such we describe the relative importance of each mode as the percentage variance explained by each mode: 120 If we now transform each column of back to image coordinates; these can be interpreted as the spatial patterns of uncorrelated weighting which partially explain the geophysical field. We provide a visual representation of POD for a series of flow snap shots in Figure 1 and outline in Section 2.2.2 how POD is applied in this study to ice velocity data.

Study area and ice velocity data 130
For this exploratory study, the selection of a study site and the availability of ice velocity data are interlinked. To increase the likelihood of revealing statistical structure through POD the study area must show clear seasonal variability, and its velocity be well sampled. With these criteria we selected Sermeq Kujalleq in West Greenland, also known as Jakobshavn Isbrae ( Figure   2). Herein, we refer to this glacier with the initialism "SK-JI" to distinguish it from other Greenlandic glaciers with the same official name. SK-JI is covered by the "W69.10N" sub-region of the "MEaSUREs Greenland Ice Velocity: Selected Glacier 135 Site Velocity Map Version 3" archive (Joughin et al., 2020a, Joughin et al., 2020bJoughin et al., 2010). The primary data for this study are colocated ice-surface velocity magnitude rasters from this sub-region. The dataset has a nominal 100m xyresolution, is derived from SAR imagery from 2008-2019, and hosted by the National Snow and Ice Data Center.
The ice velocity of this MEaSUREs sub-region is frequently measured with small image-pair separations and is therefore likely 140 to capture the most complete range of glacier behaviour possible. Error associated with each pixel is approximately ± 10-20 m a -1 but may be subject to geolocation errors up to 100 m, as outlined by Joughin et al. (2020a). We assign each velocity image the centre date from the image pairs used to derive it and record the image pair separation. In total 562 velocity magnitude images, and corresponding error maps, with centre-dates from 17-Jun-2008 to 20-Dec-2019 were accessed.

145
From this dataset we select spatial and temporal limits in which the applicability of POD can be best assessed. Spatially, this selection is driven by data coverage as expressed on a pixel-by-pixel basis as a percentage of those 562 velocity images which contain a measurement (Figure 2d). Riel et al. (2021) identified travelling waves in the trunk of SK-JI with a propagation speed of up to 400 m day -1 . To minimise the potential effect of this manifesting in POD modes and potentially limiting the identification of signals due to subglacial hydrology we further limit our study to a 4 km long section, such that any travelling 150 wave within 11-day image-pair velocity image will not be resolved. Furthermore, the traveling wave amplitude decreases sharply at 10 km from the terminus position (Riel et al., 2020, their Fig. 2), and so we seek to identify a region upstream of this. Temporally, this selection is driven by our requirement to identify a multi-year period where: fluctuations in velocity are primarily annual; where no, or very little, multi-year trend is observed, and thus where stationarity of the dataset is approached.
We therefore select an exceptionally well sampled region of the SK-JI trunk ( Figure 2) and an analysis period of summer 2012 155 (01-June) to end of winter 2016 (01-March). The inter-annual velocity variability of SK-JI is strongly related to its terminus position (Lemos et al., 2018;Bondzio et al., 2017) and this period corresponds to a time of relatively stable frontal position with only a minor multi-year slowdown, contrasting to after 2016 when the rate of slowdown increases markedly (Lemos et al. 2018). Our study area contrasts to previous studies which have investigated the velocity of SK-JI with high temporal resolution at the near-terminus (Cassotto et al., 2015;Sundal et al., 2013). 232 velocity images are available covering the study 160 period and area, 229 of which have an image-pair separation of 11 days; two (centre dates: 7-Aug-2013; 9-June-2015) have a separation of 22 days; and one (centre date: 1-Aug-2013) 33 days. The median difference between velocity image centre-dates is 5 days, the 25 th percentile 1 day and 75 th percentile 10 days.

Data processing and POD
We prepare the dataset by reading in each ice velocity raster and saving a "stack" of colocated images, and a table of relevant metadata including raw filename, image pair dates, separation, and centre-dates. In MATLAB R2020a the i-j-t stack is then rearranged into the form of (Eq. 1). Some small gaps exist in the satellite-derived velocity data, in the very southwest corner 175 of our study area 4 measurements are missing in the pixel time series', and over the vast majority of the study area only 1 measurement is missing. To fill-in these minor gaps we interpolate temporally the satellite-derived ice velocity measurements using D'Errico, (2020). We then resample the velocity measurements using a piecewise, non-overshooting cubic interpolation (Akima spline) to be regularly spaced according to the median measurement gap (5 days). A SVD of this matrix is then calculated and a three-value moving average taken of V to dampen very high frequency fluctuations, following the best practice 180 of Higham et al. (2016). These fluctuations are potentially due to some measured velocities being derived from image pairs closely spaced in time, but collected with ascending and descending orbits, noted to result in velocity discrepancies in some areas by Riel et al. (2021). The velocity matrix is then reconstructed (i.e. the inverse of Eq. 1), and the mean of each pixel timeseries subtracted to find the velocity anomaly of the whole dataset. The POD is then performed as in Section 2.1 and Figure 1. A MATLAB script for this procedure is provided with this manuscript. 185

Characterisation of surface hydrology
To test whether melt-season velocity fluctuations can be attributed to surface-to-bed glacial hydrology, as in Mair et al. (2002), we use satellite imagery within Google Earth Engine to characterise daily supraglacial meltwater pond behaviour. We focus on two immediately upstream sub-regions ( Figure 2) where ponds form quasi-annually. These correspond to the "saturated 190 crevasses" CV2 and CV5 identified by Lampkin et al. (2013) as forming in 2007 and delivering enough meltwater during drainage to enhance basal sliding in a distributed subglacial hydrological system. We assume that rapid changes in ponded water extent corresponds to spatially and temporally focussed delivery of meltwater to the ice bed.
To begin, we derive outlines of the maximum annual supraglacial melt pond extent for the melt season (1 May to 30 September,195 2009-2019) from daily Moderate Resolution Imaging Spectroradiometer (MODIS) Terra images within Google Earth Engine (Gorelick et al., 2017). Level-2 MOD09GQ band 1 surface-reflectance data (250 m resolution) were cloud masked where coincident band 7 (500 m resolution resampled to 250 m resolution) reflectance values exceeded 0.09. Supraglacial melt ponds were defined using a dynamic threshold approach (Selmes et al., 2011;Williamson et al., 2017), where a central pixel's value had a reflectance of less than 0.65 against the mean reflectance within a moving 21-by-21 pixel window. 200 These daily supraglacial melt pond area masks were subsequently summed to create annual composites of maximum melt pond extent and false positives were removed where the ponded water either: did not appear twice within six successive days from first appearance; was not observed during a period of five cloud-free days before the third occurrence; did not appear on three or more occasions over the melt season; or was observed on more than 50 % of cloud free days, acting to remove land 205 or protruding bedrock (Cooley and Christoffersen, 2017;Liang et al., 2012;).
We then create the two study areas around the maximum melt pond extents 2009-2019 for high-resolution temporal analysis ( Figure 2b). Here we analyse daily MODIS imagery within these two subregions using the thresholding method above, recording the cloud cover and the "melt pond" area within each box. Importantly, we are unable to perform the same stringent due to the rough crevassed terrain. We reject all measurements when cloud cover over each subregion is >10 %. We identify potential drainages qualitatively when no melt pond area is consistently detected after a period of consistent melt pond detection. Comparison of these dates with available corresponding 30 m spatial resolution Landsat-7 and -8 imagery demonstrates we accurately bracket major melt pond drainage events (Appendix A). 215 In addition to the remote sensing approach outlined above, we use runoff output from the Modèle Atmosphérique Régional (MAR) v3.11 (ftp://ftp.climato.be/fettweis/MAR/) extracted from the SK-JI catchment  and portioned into 500m elevation bands using the dataset's surface height layer to characterise the seasonal runoff in, and upstream, of our study region. 220

Ancillary data
We draw upon several open datasets to aid interpretation of POD results: dates rigid melange is present in front of SK-JI and relative terminus position from Joughin et al. (2020b), and topographic datasets from BedMachine v3 (Morlighem et al., 2017).
From these topographic datasets we estimate the subglacial water flowpaths by calculating subglacial flow accumulations 225 (using Schwanghart & Scherler, 2014) from the Shreve (1972) hydraulic potential surface where ice-bed effective pressure is zero. In our study area subglacial flowpaths are strongly constrained by subglacial topography, nevertheless we also calculate Shreve flowpaths with varying subglacial water pressures (ρw), expressed as a fraction of ice overburden (ρi) following Wright et al. (2016). We find subglacial water confluences and flowpaths do not show major changes unless ρw/ρi <0.5: such extreme low values are not observed in the boreholes of Wright et al. (2016), and as such consider our simple modelling of flowpaths 230 robust and appropriate for this area.

POD
The decomposed dataset is shown in Figure 3 and Mode 1 dominates the POD results. Its spatial weighting broadly reflects the mean ice velocity pattern ( Figure 3b) and its temporal coefficient peaks annually in late summer and reaches its minima in 235 spring. This is confirmed by the power spectral density (Figure 3c) showing an almost exclusive annual periodicity. The autocorrelation of Mode 1 results in an integral length scale of 20, confirming that for our purposes this dataset is fully converged following the criterion of Tropea et al. (2007). There is no negative spatial weighting, meaning the entire study area increases (spring to summer) and decreases (autumn and winter) in-phase, with a higher seasonal amplitude in the downstream of the study area. 240 Modes 2 and 3 show similar spatial structures with positive weighting to the west (down glacier) of the study area and negative weighting in the east of the study area (up glacier). This indicates a switching velocity differential across our study area, where one area shows above average velocities due to this mode while another shows below average velocities. Their temporal coefficients show annual periodicity and, generally, peak temporally in summer and have troughs during winter. Importantly, 245 the temporal coefficients resemble quasi-conjugate pairs, in that they have similar Var(%) values and higher-frequency components appear mirrored in one another. For example, during the series of peaks in winter-2012 to spring-2013, and the major signal excursion autumn-2014. Higham et al. (2018) showed that similar modal pairs in decomposed particle image velocimetry datasets relate to the formation of the von Karmen street in the lee of a flow obstacle. This indicates a common, but spatially variable or moving, process is responsible for Modes 2 and 3 in our dataset. When summed these quasi-conjugate 250 pairs may provide a picture of the overall distribution of the process responsible, we investigate this further during POD interpretation (Section 3.2).
Modes 4 and 5 show some common spatial structures, but no clear periodicity. Broadly spatial weightings are negative in the central east-west axis of the study area and negative towards the study area margins. The most prominent feature is the N-S 255 narrow band of very high absolute values to the west of the study area. Mode 6 and lower modes show no discernible spatial structure either temporally or spatially, and account for a very small amounts of the dataset variance.
In Figure 4 we plot the percentage of variance explained by each mode, calculated from the decomposition singular values ( ). As expected, the variance is strongly dominated by Mode 1 (Var(%) = 98.1). When plotted on a logarithmic axis 260 we see a distinct "kink" at Modes 2 (Var(%) = 0.9) and 3 (Var(%) = 0.5), which disrupts the otherwise smooth decay of eigenvalues characteristic of incoherent signals and noise in the higher modes. This is evidence that Modes 2 and 3 contain some physically-meaningful information, that they are paired in some manner and warrant further investigation.

270
(d) The temporal coefficients, V, for Modes 1 -6. To aid the interpretation of seasonality the sign (+/-) of Modes 1-3 are reversed such that V peaks in summer, because this is applied to both U and V this has no effect on the results.

Interpretation
The modes extracted from the decomposed dataset are purely statistical features and do not necessarily bear any relation to glaciological processes. In this subsection we compare the spatial modes and temporal coefficients with complementary glaciological and environmental datasets to assess whether the forcing from particular glaciological and environmental 280 processes may contribute to the patterns revealed by the POD analyses. With this approach we explicitly do not assign glaciological meaning to the short wavelength features in the spatial modal weighting, U. This is important following our initial observation that Modes 2 and 3 appear paired, and as such must be combined before any inference of broad forcing mechanism can be made. We also discuss the published errors of the original ice velocity dataset to ensure that we do not attach unwarranted glaciological significance to features in the decomposed dataset. 285

Glaciological and environmental factors
In Figure 5 we compare the temporal coefficients (V; Figure 5a) with a variety of potential glaciological forcing mechanisms.
Mode 1 has a strong seasonal cycle that tracks the advance and retreat of the terminus position (Figure 5e dynamic enhancement and adjustment due to the changes in terminus stress associated with variations in calving or grounding line position. This physical interpretation of Mode 1 explains its dominance (Var(%) = 98.1) and the increase in the positive spatial weighting in the direction of the calving terminus. This is expected as several other authors have also noted that stress reconfiguration due to frontal ablation dominates the bulk velocity signal of SK-JI (e.g. Joughin et al., 2020a;Khazendar et al., 2019;Lemos et al., 2018;Bondzio et al., 2017). Figure 6 shows two distinct, secondary, early-summer peaks in the velocity 300 at points I-VI (most obvious in 2012 and 2014) also reflected in Mode 1. The cause of these peaks are not simply attributable to terminus position change or any forcing mechanism. Both occur in the highest two melt years and so may be due to meltwater input. That these are in Mode 1 indicates a widespread decoupling of the bed due to early season melting and subsequent pressurisation of the subglacial hydrological system, similar to "spring events" first seen on valley glaciers (Mair et al., 2003;Iken et al., 1983). 305

grey bars denote rejected data impacted by cloud cover (see Section 2.3). (e) SK-JI terminus position from Joughin et al. (2020b) with higher (lower) values indicating a more retreated (advanced) position, for illustrative purposes Mode 1 V is scaled and overlain.
Following from our interpretation that Mode 2 and 3 are quasi-conjugate pairs, i.e., they are strongly related by a common process but with features shifted in phase, in Figure 5b we sum of Modes 2 and 3 (herein referred to as Modes 2+3) to gain an overall picture of the timing of the process responsible. Higham (2018)   One powerful use of POD is the ability to reconstruct the dataset using only selected modes: allowing us to see the combined 350 effect of Modes 2 and 3 on ice velocity. In Figure 7a-c we compare Modes 2 and 3 U, and the locations of persistent surface meltwater ponds and simple subglacial flowpaths after Shreve (1972). In Figure 7e the velocity due to Modes 2 and 3 is reconstructed at representative Points I-VI and compared to the temporal variation in supraglacial hydrology from Figure 5. This approach offers a manageable way of interpreting the large POD results dataset in detail. Points I-VI were chosen to sample to the domain approximately evenly in space and prior to subsequent analysis. The significance of this lies not in 355 explaining the detail of temporal and spatial glacier dynamics within a relatively small area, but rather lies in assessing the extent to which the POD modes are capable of discerning hydrological forcings, even on a glacier where glacier geometry, mass balance regime and calving processes are more dominant.
We firstly discuss how different modes of delivery of water to the bed may influence the dynamic patterns revealed by Modes 360 2 and 3. Figure 7 reveals contrasting dynamic patterns between points I-III, which display high frequency variations superimposed on a broad annual cycle, and points IV-VI, which show fewer high frequency variations but a clearer annual cycle (Figure 7e). Points I, II and III are located down glacier of all the supra-glacial melt ponds and are most proximal to the northern area, N, where the data indicate that ponded areas are larger and drainage events more frequent than the southern area, S (Figure 5d, A1). Consequently, these locations will be most influenced by hydrological forcing as described in type 1a 365 above which would explain pronounced fluctuations in velocity anomalies some of which are clearly coincident with identifiable meltwater drainage events. Points IV -VI display less short-term velocity variability as they sit downstream of the less active area S. The different locations also vary in the sign of their dynamic response to potential melt pond drainage events, as identified in Figure 5d. For example, points I and III show an increase in velocity coincident with potential drainage events whilst point II shows a decrease in velocity coincident with these events. To explain this and other patterns we secondly discuss 370 how subglacial drainage system structure affects dynamics, invoking conceptual models of their spatial and temporal evolution.
Shreve's subglacial flow accumulations (Figure 7c), derived from hydraulic equipotential mapping, indicate the likely preferential drainage axes (PDAs) within the study area. These are the broad axes along which temporally dampened, catchment-wide melt from up-glacier will be directed during summer (i.e. type 1b above), and where base-level water flow, 375 generated from frictional melting, will flow throughout the year. Bearing in mind that, in this region, subglacial flowpaths are strongly constrained by topography, (Section 2.4) we are cautious to not over-interpret their precise location, and refer to reconstructed PDAs in conjunction with other independent data sources . Furthermore, uncertainties in basal topography create uncertainties in the exact locations of PDAs (Mackie et al., 2021): hence in areas where multiple PDAs are confluent, we have higher confidence that these are areas of high hydrological activity. Shorter term dynamism and transients will undoubtably 380 result in subglacial flowpaths outside of these axes (Stevens et al., 2018;Wright et al., 2016), but given that the velocity dataset analysed here are 11-day image-pairs we would be unlikely to be able to resolve any such higher temporal and spatial resolution hydro-dynamic coupling. The dynamic response to hydrological forcing will vary through time as the drainage system evolves.
Following Rothlisberger (1972), when discharge through these regions is consistently high (as would be likely during the summer melt-season), turbulent flow will develop which may result in the development of hydraulically efficient, possibly 385 channelised, low-pressure subglacial flow (i.e. type 2a above) with minimal impact on basal sliding and ice velocity. This would explain why Points IV, V and VI, which are located close to major confluences in flow accumulation, show negative velocity anomalies during the summer melt-seasons. The developed of efficient, drainage at Point II (located directly over a PDA) also explains the negative velocity anomalies during drainage lake drainage events, as noted above. Conversely, during the rest of the year, all drainage axes will have closed by ice deformation and the glacier will likely be underlain by an 390 inefficient, distributed drainage system (i.e. type 2b above). During these periods, the areas with confluent PDAs will become the most likely areas of the bed to experience low magnitude, base-level water flows which will promote higher water pressures and enhanced basal sliding. Consequently, points IV, V and VI show positive velocity anomalies during the winter. The annual cycle of points I and III shows negative velocity anomalies during the winter which likely reflects their isolation from flow confluences and hence that they are likely to be located over unconnected regions of the glacier bed where there is minimal 395 hydrological forcing in the winter (following Hoffman et al., 2016). In summer, the forcings from high magnitude melt pond drainage events discussed above create positive velocity anomalies at points I and III since short-term, high magnitude discharge is forced into inefficient drainage systems, hence provoking high basal water pressures. Notable fluctuations in velocity seen at points I-III in December 2012 may be due to release of stored water from englacial or subglacial sources (Pitcher et al., 2020). 400 Despite the complex temporal and spatial modal patterns described by Mode 2 and Mode 3 we propose that together these modes do define the component of dynamic variability that is strongly influenced by changing patterns of subglacial hydrological forcings. The seemingly complex behaviour can be explained with recourse to long-standing hydrological theory and simple conceptual models for meltwater delivery mechanisms and drainage system structures. If the POD analyses could 405 be performed over a larger spatial scale, or across a glacier that was less directly dominated by calving dynamics, then POD modes might be defined that more obviously reflect subglacial hydrological forcing, similar to Mair et al. (2002). This will become possible as the temporal resolution and spatial extent of satellite derived glacier and ice sheet dynamics continues to improve.

SAR-derived velocity error analysis 415
In this subsection we consider the published errors of the SAR velocity product of Joughin et al. (2020b) itself, and their likely effect, if any, on our POD mode results. In our domain and study period the median error for individual pixels is 11-19 m a -1 , with maximum (minimum) values 18-21 (5-10) m a -1 . Importantly, these values relate to the bulk, non-decomposed velocity record and not individual modal components of the velocity as identified with POD. The velocity fluctuations in the reconstructed velocity series using only Modes 2 and 3 (Figure 7) vary by an order of magnitude more: typically, from -100 420 to +100 m a -1 with extreme values almost twice this. As such, even if we were to make the simplistic assumption that Modes 2 and 3 experience all the error, rather than some portion of it, our interpretation would remain unchanged. The time series of velocity error magnitudes at Points I-VI are highly variable measurement-to-measurement and do not show a seasonal pattern, unlike Modes 2 and 3. A relevant study of POD and measurement errors comes from Wang et al. (2016) who showed that errors in a particle image velocimetry dataset are accommodated in the higher order modes, not the lower modes due to physical 425 processes. This study, along with the above discussion, lead us to conclude similar behaviour in this study. Nevertheless, to check the robustness of our decomposition, we create an "error envelope" for the Mode 2 and 3 reconstructed velocities shown in Figure 7. We do this by repeating the POD with the published error added to the velocity field, and then with the error subtracted. This results in variations to the velocity anomalies shown in Figure 7 of around ±1 m a -1 which are negligible in the context of the velocity fluctuations due to Modes 2 and 3, and our interpretation that they are dominated by glacio-430 hydrological induced ice motion.

Interpretation Summary
To summarise, Mode 1 relates strongly to the terminus driven variability that is known to strongly dictate the bulk observed velocity signal of SK-JI. Modes 2 and 3 together describe the seasonal development of the spatially variable subglacial 435 hydrological system influenced by episodic input of meltwater and the changing efficiency of subglacial hydrological drainage axes. Mode 4 and lower modes show no interpretable spatial or temporal patterns and therefore likely do not contain useful statistical information.

Discussion
The work presented here is an important contribution to efforts to explain how patterns of ice motion, observed from space, 440 can reveal different glaciological forcing mechanisms. As the outlet glaciers of Greenland retreat, and their geometry and flow evolves, it is likely that the relative importance of processes governing their ice dynamics will also evolve, though this is less easily assessed. This study acts as a proof of concept, building on the work described by Mair et al. (2002), that the eigendecomposition of ice velocity can detect the influence of changes in glacial hydrology on ice dynamics. Figure 7 shows that hydrologically-forced flow contributes to velocity variabilities of approximately ±100 ma -1 to the ice flow in our study area of 445 SK-JI. This contribution varies seasonally but also in response to sudden inputs of meltwater.
Our study area within SK-JI was chosen owing to the excellent and openly available ice velocity dataset, but it is arguably not the ideal place to test the applicability of POD to ice velocity due to the dominance of stress patterns created by changing terminus geometry on the overall seasonal velocity, as outlined by previous work (Joughin et al., 2020a;Lemos et al., 2018;450 Bondzio et al., 2017). The application of POD analyses to a velocity record from a marine-terminating glacier with seasonal behaviour more influenced by subglacial hydrology such as Kangiata Nunaata Sermia (e.g. Davison et al., 2020, Sole et al., 2011 or a land terminating section of the ice sheet, especially one with complementary discharge measurements such as Russell and Leverett Glacier (e.g. Davison et al., 2019;Chandler et al., 2013) might yield more insights into their hydrodynamic coupling behaviour. This underlines the importance of continued production of high-resolution, all-year SAR datasets 455 of ice velocity in key, dynamic areas. Practically this involves the continued exploitation of Sentinel 1 and TerraSAR-X/TanDEM-X for ice velocity and, looking forward, to future missions such as NASA-ISRO SAR (NISAR).
In this analysis, although we only characterise supraglacial melt pond drainages in the immediate locality of our study area, these appear to explain the melt season major Modes 2 and 3 reconstructed excursions in ice velocity. At Sermeq Kujalleq 460 (Store Gletsjer), the drainage of supraglacial lakes has been shown to rapidly increase the efficiency of the subglacial drainage system and decrease ice velocity at a point near the terminus , strengthening our interpretation in this study that the drainage of supraglacially ponded water can both slow and enhance ice velocity. The interpretation of decomposed surface velocity fields presents a potential method for comparing field data to model output in a more developed manner than straightforward correlations. We propose comparing low-dimensional representations of measured ice velocity (i.e. separated 465 POD Modes) with ice sheet model outputs calculated with differing parameterisations of subglacial hydrology and basal motion. This would provide a means of assessing the efficacy of models to account for spatially and temporally variable hydrodynamic coupling, ultimately leading to more physically robust predictions of future ice sheet behaviour.
The POD technique applied here is complementary to that of Riel et al. (2021), who also decomposed the ice velocity record 470 of SK-JI as imaged by TanDEM-X and TerraSAR-X from Joughin et al. (2020b;2010). In that study the velocity time series of a wider area is decomposed into the seasonal cycle and transient multi-annual parts to show the propagation of a kinematic wave emanating from the glacier terminus. In this study we have isolated years where data stationarity can be reasonably argued, essentially removing the multi-annual transient. Consequently, our Mode 1, primarily due to seasonal terminus fluctuations, is broadly analogous to their seasonal velocity signal. Our Modes 2 and 3, that apparently relate to 475 glaciohydrological processes, occur on shorter timescales than the seasonal signal extracted by Riel et al. (2021). Riel et al. (2021) show that their B-spline approach can be used to detrend multi-annual velocity variations. This raises the prospect of applying POD to the reconstructed seasonal velocity only. This could address the issue of stationarity that limits the applicability of POD on datasets with long-term trends: a common feature of ice velocity time series. One consideration that might complicate this, however, is that in the four years analysed in this study, the influence of subglacial hydrology appeared the influence of hydrological forcing on spatial patterns of ice motion may be divided across many more separate modes which would make physical interpretations more difficult. One factor to be considered in further use of POD is the length scales over which basal variability is transferred to the ice surface velocity, and modal components of that surface velocity. This study, and Mair et al. (2002), offers empirical evidence that changes in basal hydrology are reflected in modal components of the 485 velocity field on length scales close to the limit where "bridging stresses" might obscure the velocity response of basal hydrology had the bulk velocity signal been interpreted alone. Future work could investigate these effects further, potentially using POD to focus on key glaciological features (e.g. supraglacial lakes, subglacial lakes, patterns of basal drag and/or ice stream sticky spots) over timescales and length scales that capture their dynamism.

Conclusions and outlook 490
In this study we have presented the first application of proper orthogonal decomposition (POD) to a series of satellite-derived maps of ice velocity. By focussing on an exceptionally well-sampled area of the SK-JI trunk from 2012-2016 (Joughin et al., 2020b) we significantly develop the POD/EOF analysis using survey stake observations of ice velocity on Haut Glacier d'Arolla by Mair et al. (2002) and apply it to a Greenlandic marine-terminating glacier.

495
We explored statistical structures or "modes" in the data that partially explain dataset variance. The first Mode 1 clearly relates to the bulk velocity signal that previous researchers have shown links strongly to terminus-driven stress reconfiguration, transferred upstream rapidly due to low bed and lateral resistance (Joughin et al., 2020a;Bondzio et al., 2017). However, a second seasonal pattern exists which we hypothesise is due to glacio-hydrological processes. We present independent evidence from remote sensing and several published datasets to support this interpretation and explain the observations in terms of 500 proximity to the confluences of subglacial preferential drainage axes and the input of meltwater to the subglacial environment from two nearby supraglacial melt ponds. Short-term dynamism of glaciers is often explained in the context of evolving subglacial drainage systems and with respect to the manner of water delivery to the bed. Our POD analysis offers a powerful means to disentangle such complex influences from other higher-order drivers of ice dynamics. The decomposed modes may contain some small artefacts that relate to the original dataset and its errors. Modes do not necessarily relate to a single 505 glaciological process and signals due to spatially variable processes may be split between modes. However, we have shown that with careful interpretation these effects can be accommodated for. As ice velocity measurements continue to increase in spatiotemporal resolution and continuity, these will likely reduce in the future. POD, and potentially related decomposition techniques such as dynamic mode decomposition (Schmid., 2010), are flexible, 510 powerful and quick techniques that can be applied to other colocated measurements of ice surface velocity, derived from satellites, time-lapse imagery or GPS stake networks. The next step is to apply POD to other high-resolution records of ice surface velocity from glacier systems with a different balance of processes with contrasting propagation speeds. Where hydrology is a more important control on the overall velocity pattern, its associated modal expression would be more dominant.
We have applied POD only to the ice velocity magnitude, however, extending the analysis to polarstereographic grid x and y 515 components is relatively simple and may provide insights to explain ice velocity patterns in the vicinity of key glaciohydrological features. Furthermore, POD may offer a low-dimensional method of comparing surface velocity measurements with coupled glacier hydrology and ice dynamics models that simulate hydrological and ice dynamic processes explicitly and coincidently. We assert that computationally lightweight modal decomposition techniques such as POD can become a key method in illuminating the effects of individual glacial processes on combined emergent patterns of ice dynamics. 520

Appendix A
To evaluate the results of the daily MODIS dynamic thresholding method as outlined in Section 2.3, we compare U.S.
Geological Survey Landsat-7 ETM+ band 4 and Landsat-8 OLI band 5 scenes with < 30 % cloud cover either side, and close to, the drainage dates detected using the daily MODIS thresholding method as shown in Figure 5. This acts as a validation of 525 the automated Google Earth Engine approach, and that the daily, but relatively low resolution, MODIS can detect melt pond drainages in this setting. It is important to note that the time between identified drainage and a suitable Landsat-7/8 image acquisition varies considerably, and this must be considered when comparing the two datasets. Furthermore, due to the large contrast in spatial resolution between MODIS (250 m), and Landsat-7 and -8 (30 m) we expect some disparity between the two approaches. The overall aim here is to increase confidence that the timing of the delivery of meltwater from the surface, 530 to englacial, and ultimately the subglacial environments, is correctly bracketed.
Results are shown in Figure A1. In the less active, southern melt pond area one drainage event is identified per year, occurring in July. For 2012 and 2013, drainages S1 and S2 in Figure A1, the available images show a slightly smaller ponded area indicating that only a partial drainage took place or refilling partially occurred in the days after drainage and before the Landsat 535 acquisition. The 2013 drainage, S3, however is well identified by our technique. In 2012 and 2013 a secondary drainage ( Figure   A1c-d and A1j-k) takes place. These are visible in Figure 5d as a short return to zero melt-pond area conditions after short periods of consistent meltwater detection, but these were not judged to be clear enough for identification as a potential drainage event during analysis (Section 2.3). In both these cases however, drainage corresponds closely to the drainage of the larger, more active, northern melt pond area which is, in all cases, well-identified with the MODIS-derived record. In totality, this 540 comparison increases our confidence in our approach correctly identifies melt water leaving the ice surface and, in turn, our interpretation that Modes 2 and 3 are strongly influenced glacier hydrology.

Code availability 550
An annotated MATLAB script to perform the POD from the stack of ice velocity magnitude rasters, and recreate the results of this study is available at https://doi.org/10.5281/zenodo.4699392.

Data availability
Data generated and presented in this manuscript is available at https://doi.org/10.5281/zenodo.4699392. All secondary data used in this manuscript is provided from the relevant citation within the reference list. 555

Author contribution:
DWA and DWFM conceptualised the study and led the interpretation. DWA led the POD analysis, data synthesis and manuscript writing. JEH provided guidance on the POD analysis. SB and JML performed the MODIS Google Earth Engine processing and advised on the interpretation of supraglacial hydrology. SB performed MAR data analysis. IJN provided guidance on the implications of results for ice dynamics and modelling. All authors contributed to and commented on drafts 560 of the manuscript.