General Comments

The manuscript presents a thoughtful and accessible methodology for performing clustering analysis for glacial seismology. Given increased interest in continuous seismic monitoring of the cryosphere, the study is timely and instructive and would be of value to members of both the seismological and cryosphere communities. Furthermore, the manuscript is well written and strikes an appropriate balance between introducing basic theory and demonstrating careful application and thorough analysis.

One of the chief concerns I have is the utility of k-means on a parameter space with 30 dimensions. Though the authors took care to appropriately select features for clustering, the procedure was not conceived with dimensionality reduction in mind, nor is dimensionality reduction mentioned once in the literature review, although several papers were cited that make use of it (e.g., autoencoders, PCA). Given the dubious utility of distance-based clustering metrics in high dimensions (Aggarwal et al. 2001; Aggarwal & Reddy 2014), I would at least expect acknowledgement of this limitation. I hesitate to require re-performance of what is already a substantial analysis presented in this manuscript, but a good candidate for future work would be repeating the analysis but with reduced dimensionality, e.g., using the top 5 or 10 components from PCA.

My other main critique is the use of the k-means(++) algorithm. K-means has several important limitations, namely that it performs poorly with overlapping clusters and that it assumes equal variance in all dimensions. It is sensitive to outliers, which is exacerbated by the high dimensionality of the feature space. I would instead suggest using an expectation-maximization (EM) approach to refine the cluster definitions, as provided in Gaussian mixture model (GMM) clustering. GMM handles multivariate distributions and overlapping clusters better than k-means, and is trivial to implement with scikit-learn. However, GMM is still subject to the curse of dimensionality.

In summary, the authors should at a minimum revise the manuscript to acknowledge the curse of dimensionality and the limitations of k-means clustering. Reworking the entire analysis to reduce dimensionality and implement GMM clustering, while desirable, is likely beyond the scope of this work.

Specific Comments

[Line 62] I would be careful with your use of "high-dimensional" here. As you correctly identify, *k*-means clustering relies on Euclidean distances; however, Euclidean distance becomes a less useful metric for clustering since, with increasing dimensionality, data become sparse and distance less meaningful (the curse of dimensionality). It thus becomes increasingly difficult to distinguish clusters, since all the distances between points appear the same. An informative exploration of this phenomenon is given by Aggarwal et al. (2001, http://link.springer.com/10.1007/3-540-44503-X_27).

[Line 96] I would suggest citing the paper (Jenkins et al. 2021), not the AGU abstract.

[Line 102] You cite another AGU abstract by Sawi et al. - if they have a related paper, you should cite it instead.

[Lines 170-176] I appreciate the thorough explanation, but feel it can be said more concisely.

[Line 175] "...the features demonstrate comparable distributions." Do they? Perhaps you mean to say they are distributed over comparable scales? Because the distributions themselves are quite unique: some are normal, some bimodal, etc.

[Line 178-179] Let me first say that I am pleased to see this careful treatment of your input features, acknowledging that not all features are useful for the clustering analysis. However, I think you should explain your reasoning further or at least provide a citation to a useful reference, for the benefit of your readers. What type of bias are you referring to? Is it distinct from the bias you strive to eliminate by standardizing your features? Why is the inclusion of too many features bad?

[Sec. 3.2] Doesn't k-means as implemented by scikit-learn use random restarts, as well? If so, you should include this in your description and discuss how the centroids are determined accordingly.

[Lines 207-208] The silhouette score is tricky to describe, and I'm afraid your explanation leaves me confused. Particularly confusing is the phrase, "each cluster set of means." Please clarify the explanation.

[Sec. 3.1, Sec. 3.2, Fig. 3] You index clusters by r, but also use "r" as the correlation coefficient. I would suggest de-conflicting your notation, e.g., indexing clusters as k=1,...,K and reserving r for the correlation coefficient. Furthermore, I would suggest italicizing the correlation coefficient so it is abundantly clear you are referring to a parameter. In the caption of Fig. 3, I initially interpreted "(r)" as "right-hand."

[Line 262 & Suppl. Line 46] Why is k=14 a "realistic" maximum?

[Sec 3.3.2, Sec. S3.3, Fig. 5, Fig. S6] I applaud your effort to thoroughly examine the "evolution" of clusters. However, as currently written, the first paragraph of Sec 3.3.2 leaves me with a sense of confusion at best, or at worst that the analysis may be flawed. How are you able to track the progeny of k clusters from from the previous k-1 clusters? You mention that k-means has built-in pseudo-randomness (and random restarts), but between Section S3.1 and this section, I fail to understand how (and why) you overcame this. Are you setting the seeds manually, and preventing restarts? Though you refer to Section S3.3 in a previous section (Line 209), I would refer to it again here, as there are crucial details in it that you should relate to this part of the text.

The last thing I'll say about this section & related supplement/figures is that I had to re-read them several times to become convinced that you are not proposing that the composition of clusters at k+1 is dependent on the composition at k. The discussion in paragraph 2 (Line 261) ultimately settled the question, but I think if you clean up the first paragraph, you can alleviate the confusion. With the question settled, Fig. 5 is a nice analysis of how cluster composition changes.

[Suppl. Line 64] This is not a priori; it is a posteriori since you must run experiments to determine the optimal number of clusters.

[Suppl. Line 66]  "that is, ..." This is self-evident and redundant.

[Suppl. Lines 65-71, Fig. S3b] What justifies ignoring the outliers?

[Fig. 7] The top row is an interesting analysis, but I'm not convinced of the value of the bottom row. To me, it seems you are comparing apples to oranges. Because the datasets are essentially different due to the inclusion/exclusion of additional features, there are too many factors that affect the cluster assignment of data, including pseudo-random seeding, restarts, and, even if those are controlled, the vastly different dimensionality of the parameter space. All this is to say, clusters in one dimensionality do not look like clusters in another dimensionality, as indicated by the top row of subplots.

Technical Corrections

[Line 101] A space is missing between "fracture" and the citation.

[Line 182] Just use "variance."

[Fig. 2] This is a rather trivial comment, but your 4th column is missing x-axis tick labels. My preference (certainly not a prescription) for this type of figure is to have the same y-scale for all plots, and print the axis labels on just one subplot.

[Suppl. Line 58] Change to "where n=100 is the number of bins...", etc.

Dear authors, dear editors,

I have read with great interest the scientific article titled "Towards the systematic reconnaissance of seismic signals from glaciers and ice sheets - Part B: Unsupervised learning for source process characterisation" submitted by Latto et al., for publication in The Cryosphere. In this article, the authors propose to systematically explore seismic data acquired by a network of stations deployed on the Ross Ice Shelf during the austral summer of 2010-2011. The data processing pipeline, based on a priori detection (detailed in a companion paper), relies on the clustering of seismic events through the deployment of the K-means clustering method on a feature space computed from curated seismic signals features. The authors discuss the influence of feature selection and the number of clusters (one of the hyperparameters of the K-means method) on their results. They demonstrate that this approach is at least capable of revealing relatively pure clusters containing microseisms generated by stick-slip phenomena associated with the dynamics of the ice shelf. This clustering also allows for the identification of new microseismic events associated with tidal forcings.

The paper proposed by Latto et al. is remarkable for several reasons. Firstly, it is very well written and easy to follow. The literature review is particularly relevant while remaining concise. All critical information is contained within the paper, but the authors also provide a significant amount of supplementary results that address questions the reader may have. The division between the main content and supplementary content is particularly relevant. The results are convincing, and the discussion on clusters, methodology, and especially the choices of hyperparameters and features is comprehensive and very interesting. The figures are of good quality, although a few minor improvements could be made (see my comments below). Overall, this is an excellent contribution that will undoubtedly have a significant impact on communities interested in cryo-seismicity and other applications in environmental seismology. Therefore, I strongly recommend the publication of this article. However, I do have some minor comments that I will detail below.

General Comments :

The choice of the K-means clustering method appears appropriate for the present study, and the following suggestions are by no means an invitation to completely revise this paper. I believe it is already comprehensive and insightful enough for publication as it is. However, I would like to draw the authors' attention to the fact that this clustering method may not be the most relevant for working with the features proposed by Provost et al. (2017). These features often have values distribution overlapping  between each class of events, but K-means is not able to consider the "fuzzy" boundaries between different cluster. Methods such as Gaussian Mixture Models seem to be more suitable, and I suggest that the authors at least explore this family of methods in future work.

The choice of the clustering method used, however, is secondary compared to the more important question of feature selection. In this article, the authors propose to reduce the number of features by calculating pairwise correlation coefficients between features. This step is absolutely necessary for K-means to work well (as demonstrated in this study). However, this requirement limits the robustness and versatility of the approach proposed. The chosen set of features may work well for this dataset, but the correlations may be different for another dataset. Furthermore, even correlated features can carry complementary and, more importantly, relevant information for discriminating certain events (in supervised approaches like Random Forest or Gradient Boosting, removing correlated features usually reduces precision scores). An alternative approach to reducing the number of features would be to reduce the dimensionality of the feature space using dimensionality reduction methods, which would eliminate the cross-correlation selection step and retain some of the information carried by each feature. I encourage the authors to consider the possibility of using methods like PCA (although it rarely works on seismic data), or even better, t-SNE (Van der Maaten and Hinton, 2008) or UMAP (McInnes, Healy, & Melville, 2018) for future work.

The question of feature normalization is also critical. The authors chose to normalize by the standard deviation. This normalization helps to homogenize the overall distribution in the feature values space, but it sacrifices some level of information about each feature. The absolute value of the feature is probably as important as the position of the feature value in the overall distribution. Figure 4 presented in this paper is a clear example of this. What would your clustering result look like if you applied K-means (or another method) with non-normalized values in a two-dimensional space (Characteristic frequency - Duration) or three-dimensional space (Characteristic frequency - Duration - Peak Amplitude)? Have the authors tested their approach without normalization?

Finally, why was standard deviation chosen for normalization? If the goal is to preserve the properties of the distributions of each feature, there are other normalization approaches that may be more relevant (Yeo-Johnson transform, Quantile Transform, Unit Vector Scaling, Sigmoid scaling?). I suggest testing these in future work.

Minor Comments : 

L.102 : A missing space “fracture(Hammer) [...]”

L.163-164 : This needs some clarification : you computed the median, the mean or something else of the values of the features of the seismic signals recorded at each station for a given event ? Have you also considered adding in your feature arrays the standard deviation of the values of the features of the seismic signals recorded at each station ? This can provide some valuable information on the location of the source and on other geometrical properties of the event (e.g. near of far field origin, dip, etc.). 

L.175 : If for each feature you have the same distribution then those features become useless for any identification methods. Please rephrase and clarify.

L.207 : There are other estimators used to try to determine the relevance of a clustering result besides the Silhouette test (e.g. Davies-Bouldin score, Elbow score, Calinski-Harabasz index, Dunn index). It would have been interesting to see how these factors evolve in relation to the number of clusters chosen and in comparison to the ideal number decided by you.

L.231-232 : You find 136 stick-slip events identified as such by Pratt et al. (2014), 4 new ones, but how many did you miss from the Pratt et al. (2014) catalog ?

L.241-243 : Indeed, but don't you lose this information by your normalization of the features (see general comment)? The duration/peak amplitude correlation also seems quite strong (as is often observed for microseismic sources, which are sometimes assigned a "duration magnitude"). Is one of these two features therefore excluded from your analysis?

Figure 4 and Figure 7 miss subplot labels “a”, ”b”, ”c”, etc.

Figure 6 (b) : The X-axis label “Days after December 14, 2010” is not convenient for interpretation I think. I suggest giving the real dates. 

L.317 : Define “less-defined”. Maybe quantify this (percentage of events from a given class?)