Review for:

Proper orthogonal decomposition of ice velocity identifies drivers of flow variability at Sermeq Kujalleq (Jakobshavn Isbrea) Ashmore et al. 2021

Overview

This paper presents a novel analysis of a velocity timeseries maps measured using TanDEM-X/TerraSAR-X  and analyzed using proper orthogonal decomposition. POD is a formal method of timeseries into different modes which have correlated velocity structure. By examining the time content and spatial distribution of the modes (i.e. spatially correlated time series components that contribute to the time variability of the velocity signal)  interpretations are made to gain additional insight seasonal dynamics in the study area that would not be apparent without the use of POD.

Overall, the paper is very well written, and the methodology and application of POD to the study area is clearly explained. The authors present this study as a proof of concept, citing that the application of the technique would be much more useful where the seasonal signal due to hydrology was more pronounced. I agree that the location (Jakobshavn) is less than ideal given the dominance of velocity changes due to terminus variability and the small hydrologic component of the signal. Nonetheless, the technique could be quite interesting if applied to a different region and ideally over a much more expansive area.  

My main concerns which have to do with error and improving the discussion are outlined below.

Major comments

1. The study area of examined is only 2.5 x 4 km and the identified modes 2 – 6 show complicated spatial structures on the scale of a few hundred meters. Given that the ice thickness in the Jakobshavn trunk in this region is ~ 1.5 km thick, the authors should consider whether or not observing well resolved spatial structures at this scale is physically realistic. The length scale of basal variability that can be resolved at the surface is explored thoroughly by several papers by Gudmundsson.
2. The error analysis is important but unconvincing. The authors analyze whether or not spatial pattern of errors correlates the spatial pattern of modes. However, the modes identified are low amplitude and contribute less than < 5% of the total velocity variance. Are these spatial and temporal variations above the error threshold? To assess this it might be worthwhile putting the error as an envelope on the modeled mode velocities in Figure 6 and maybe add an additional figure using the modeled velocity for specific snap shots in time where the spatial patterns are strong (i.e. during lake drainages) to determine whether spatial variability in velocity are above the error threshold.
3. According to Figure 3 panel A the six modes identified only account for a combined ~ 75% of the time series data variance. Is this correct? Is the other 25% of the data variance noise? If so, to me it seems like modes-2-6 are below the error threshold. This ties into the previous question, but if this is the case, an explanation of why this doesn’t matter would be useful.

Line by line comments:

80: Would you consider this data set fully converged? Is there a variance threshold to be achieved?

200: Did you visually check to determine weather the 5% threshold is rejects questionable images?

215: Does using a different assumption on the effective pressure (i.e. effective pressure = some fraction of overburden) change the flow paths significantly from what is presented here?

General comment: Maybe change Points I-VI to alphabetic identifiers. Even though it the Modes use numeric identifiers and the Points use roman numerals it can still be a bit confusing in the text.  

340: How robust are the PDA’s if the pressure assumptions in Shreve’s calculation are changed?

340-342: Not necessarily true. Water routing pathways can end up being really complex and counterintuitive because they ultimately reflect spatio-temporal changes in the basal pressure field and thus are sensitive to transients. For an example of this see Stevens et al. 2018.

343: Uncertainties in water pressures can also make a large difference…see Wright et al. 2015.

346: and elsewhere: The use of provoke seems a bit strange. Consider changing.

Paragraph starting at 339: I find the interpretation in this paragraph quite specific and speculative for the evidence given. While the author’s use the PDA as a reference, the location of the PDA is uncertain. All the points besides III are look similarly close to the PDA (and generally closer than an ice thickness to each other), making the spatial explanations unconvincing. As an example, you claim that points 4-6 are under the influence of channelized drainage because of their inferred proximity to drainage pathways and negative velocity anomalies, but point 2 which is inferred to be directly over a drainage pathway has does not show a negative velocity anomaly? How is this consistent?

I would advise focusing the discussion more on the processes which could result in the fluctuations observed at the points and less on the proximity to inferred drainage pathways. This will keep you from forcing the velocity timeseries into a “box” based on our conceptual understanding of subglacial hydrology and the location of the inferred PDA (which is poorly constrained). By doing this, you could use the measurements to either support or refute the hypothesis that a PDA calculated using Shreve’s assumption likely controls hydrology in the region. Focusing on process will also allow for more discussion on some of the more interesting aspects (i.e. the large winter fluctuations).

General comment: Why did you chose these point locations?  

General comment: The fact that you infer such a complex velocity structure (if over error threshold) over short length scales imprinted on the main signal is a really interesting result.  It might be worth discussing this, as it seems more fundamentally interesting than the relatively unconstrained PDA analysis.  

418-419: Without mapping upgradient lake drainages and determining if there is a velocity response the analysis present does not imply this.

Overview:

The manuscript "Proper orthogonal decomposition of ice velocity identifies drivers of flow variability at Sermeq Kujalleq (Jakobshavn Isbræ)" explores the application of POD to ice velocity time series in order to efficiently decompose time series into orthogonal spatial and temporal modes. This decomposition could potentially isolate velocity signals originating from distinct forcing mechanisms, such as changes in terminus position or enhanced basal sliding due to subglacial hydrology. These signals may or may not be readily apparent from the raw velocity time series, so POD could potentially reveal key signals to enhance our understanding of ice dynamics.

Overall, the manuscript was well-written, the methods were clearly described, and the results and supporting data are quite intriguing. POD and its variants have always been powerful analysis tools for spatiotemporal data, and this paper demonstrates the utility of POD to the increasing availability of velocity data over the ice sheets. However, POD also has well known limitations, primarily due to the assumption of orthogonality of the modes. While orthogonality doesn't necessarily affect lower-order reconstruction of the data, it does leave users vulnerable to over-interpretation of the separate modes, and I believe the authors need to explore these vulnerabilities a bit further. I describe in detail the implications of orthogonality below. My other main comment is concerned with the error analysis, which I also describe below.

I do believe that POD will be an important tool for glaciologists. I also believe that potential users should be well aware of the expected behavior of POD-type algorithms on different classes of dynamical signals. If the authors can sufficiently address and discuss these issues, then this work will indeed be a valuable contribution.

Major Comments:

1) Currently, my biggest concern is the application of POD to signals that propagate in space, which encompasses most of the signals related to ice dynamics. Taking the Sermeq Kujalleq (SK) data presented here as an example, the largest signal (Mode 1) is clearly related to the dynamic response of the glacier to stress perturbations at the terminus.  However, previous work has shown that this response has a finite phase velocity, i.e. it takes some time for the signal to propagate upstream (~400 m/day; Riel et al., 2021, as a recent example). The fact that POD enforces every point in the modeling domain to have identical temporal eigenvectors (scaled by the spatial coefficients) means that every point has the same phase for Mode 1, which was stated by the authors. The consequence of this is that POD will end up splitting the spatially-propagating signal into multiple modes (with the same periodicity), just as the authors hypothesized for the conjugate Mode 2+3 pair. It's entirely possible that Modes 2 and 3 are partially "conjugates" for Mode 1, i.e. they are compensating for the error of modeling a finite phase velocity signal with an infinite phase velocity signal. Additionally, the spectral content of an upstream-propagating signal can change, e.g. pulse broadening. All in all, any interpretation of the weaker modes for subglacial hydrology is likely to be corrupted by the propagation of the compensation signal for Mode 1.

In general, the splitting of a spatially-propagating signal into multiple modes will be dependent on the phase velocity. When the phase velocity is very high, as would be the case for a dynamic wave resulting from redistribution of longitudinal stresses, then the splitting will be minimized. However, for slower moving waves, as would be the case for kinematic waves or subglacial hydrological effects, then one would expect the signal to be split into more modes, which would hinder any subsequent interpretation of these modes.

With that said, the supporting data for the hydrological interpretation for Modes 2+3 are compelling (and very cool). If the authors can find a way to account for the leakage of Mode 1 into the other modes to isolate the hydrology signal, then the interpretation would be more convincing. Therefore, I would suggest the following initial diagnostic steps and improvements:

a) Expand the spatial extent of the POD analysis (e.g., Figure 3) so the readers can see a wider picture of how the spatial coefficients change in space. I think the spatial extent shown in Figures 2b-d would be sufficient. This way, we can get a better sense of any mode splitting due to signal propagation.

b) Include a plot of the singular values for the decomposition. I know these values were included as labels in Figure 3a, but I think a plot can give a better sense of any potential mode splitting. Slowly propagating signals will be split into a higher number of modes, and one would see a more gradual decay of the singular values. Rapidly propagating signals will have a much sharper spectral decay.

c) For the purposes of a methods paper, it would be hugely beneficial to compare the decomposition between SK and a glacier with known hydrological effects that dominate its dynamics. I would expect quite different spectral decay. Since that would entail a lot of additional work, alternatively, the authors should add some discussion on how different classes of dynamic signals will have different phase velocities and how those phase velocities would affect the POD results.

2) The error analysis performed here is a bit too simplistic. Mainly, the errors are analyzed as if they are independent from pixel-to-pixel. However, it is possible for the errors to have some non-zero spatial correlation (off-diagonal terms in the scene-wide data covariance matrix), which could be quantified empirically via a variogram over areas that are known to have minimal ice velocity. A more robust analysis would be to generate random realizations of the velocity data (perhaps with covariance matrices with non-zero off-diagonal terms) and re-do the POD in a Monte Carlo sense. How would the spatial coefficients change with the different data realizations? Since POD is computationally efficient, this likely wouldn't take too much compute time.

Line-by-line Comments:

- Line 39: Might be too soon to state the modes correspond to "distinct flow features".

- Line 77: The distinction between eigenfunctions and eigenvectors can be a little misleading since an eigenfunction is instantiated as an (eigen)vector. In other words, they are both vectors, and my intuition is that the temporal vectors are closer to being instantiations of actual eigenfunctions. I would suggest to either replace the use of eigenfunctions with eigenvectors or to add a bit more detail into why the use of eigenfunctions is more appropriate.

- Line 79: What is meant by a "fully converged dataset"?

- Line 81: Orthogonal does not by definition mean independent. Achieving statistical independence (e.g., minimizing mutual information between modes) tends to be the objective of Independent Component Analysis and can result in modes with very different spatial and temporal eigenvectors. More importantly, orthogonality does not guarantee that independent forcing mechanisms will be expressed in the different modes.

- Line 99: Is there any indication of the numerical algorithm being performed by the "economy" SVD? Is it performing a full SVD and just removing a certain number of singular values? Or is it performing some form of randomized SVD? This information may be useful for reproducibility purposes.

- Line 124: I think it's a little confusing to state "the study area must be dynamic" since that may be interpreted as only dynamic time series can be decomposed with POD. I think what's being stated here is that dynamic glaciers tend to be well-observed by satellite platforms, so I would suggest some slight re-wording.

- Line 141: Is there really a requirement for the time series to be stationary in order to decompose them? What happens if non-stationary time series are used? I ask mainly because multi-year signals can be quite important for understanding the full spectrum of ice dynamics, and it would be interesting to see how POD handles those.

- Line 145: A bit awkward to describe a slowdown with an "increased rate". It would be clearer to just state the decrease in velocity values after 2016.

- Line 165-166: Are the formal data errors accounted for during the inpainting and interpolation procedure? If not, then every point is assigned equal weight, which could affect final interpolated signal.

- Line 167: Is there any difference between post-smoothing the columns of V and pre-smoothing the time series of X?

- Line 195: What is the resolution for the "high-resolution" MODIS analysis?

- Line 220: Here's the main issue: the results are stating that the entire area moves in-phase, but previous results show the phase changes with upstream distance, with a finite phase velocity of ~400 m/day.

- Figure 3: It would be useful to include a length scale bar in the images.

- Line 286: "we sum of Modes 2 and 3" -> "we sum Modes 2 and 3"

- Line 289: What are possible artefacts? If one were to interpret a splitting of a spatially-moving, seasonal signal into multiple modes as artefacts, then those would also exhibit strong seasonality.

- Line 425: Conversely, it would also be useful to apply POD to ice flow model outputs (forced with distinct mechanisms) to determine how well the different signals can be separated.

- Line 440: Again, this is a key point: spatial patterns of ice flow will split into multiple modes based on the propagation speed of the signal. Slower propagation speeds will be split into more modes.

- Line 460: Consider re-use of the word artefacts here as well. Mixing of signals into multiple modes could be considered an artefact but is definitely not a minor effect.

- Line 467: My intuition is that hydrological effects will have slower spatial propagation speeds than redistribution of membrane stresses (such as the response to terminus forcing), which would cause more modal mixing.