The manuscript by Piana Agostinetti and Sgattoni presents a novel scheme for identifying times in an earthquake arrival time data set in which a change in the pattern of differential times occurred (presumably due to a change in the velocity structure). This can be useful for the purpose of 4-D tomography (3-D plus temporal change) when data from different, relatively stable time periods should be separated. My first suggestion is that the rather vague title be modified to better reflect this key theme of the manuscript.

The authors need to be more careful about the use of the term's error (errors) versus uncertainty (uncertainties). I believe that for the majority of their uses of error (errors), uncertainty (uncertainties) is the correct word to use.

I have difficulty understanding their constructions of the covariance matrices C_e, C_e*, and C_e,w. They refer to Piana Agostinetti and Malinvervo (2018) for the explanation of C_e,w, but in my opinion, this needs to be explained here. Where C_e comes from is not at all clear, and then C_e* is the final version of C_e, which comes from Figure 4b, but what is Figure 4b based on? Further confusion is added because Figure 4b is not the least bit diagonal in appearance, but C_e* is assumed to be diagonal. The entire development section needs to be made crystal clear and thereby be reproducible.

I also found the presentation of the weights on page 11 to be overly confusing. One would expect that for a DD datum for which the events cross a change point, the weight would be small. However, the w_ij factor that enters the exponent in Equation 9 is a sum of "the weights associated to (sic) each changepoint", giving the impression that the corresponding weight W_ij(m) can be large. It confuses things even more by having the variable m be used for two different things in Equation 10.

The discussion of the MCMC sampling could be improved. The four "moves" have specified probabilities (0.4, 0.4, 0.1, and 0.1, respectively). How is the move to be made actually chosen? Is it the case that some random number between 0 and 1 is drawn, and if it is between 0 and 0.4, move 1 is selected, etc.? Is model fitness used as a factor, as it usually is in MCMC? If so, how? Clarity would be appreciated. I also recommend that the synthetic tests referred to on page 12 be presented to provide confidence in the claimed parsimonious character of the method.

A recent paper by Roecker and coworkers (Double Differencing by Demeaning: Applications to Hypocenter Location and Wavespeed Tomography, BSSA, 2021) shows that a demeaning approach performs as well as, and is mathematically equivalent to, the DD approach. This paper should be cited and the potential application of the new method presented here to their demeaning approach should be addressed, if possible.

In general, a methods paper is more valuable if it is actually applied to something tangible. Here, the authors stop at determining changepoints, which by themselves don't really teach us anything about Katla volcano. Isn't there something worthwhile that can be done with these data and the new method?

My final point is that the manuscript needs major improvements in the writing. I have tried my best to provide suggestions for such improvements in the annotated scan of the paper.

The authors provided an automatical and data-driven approach to pre-process seismic data and applied their method to the Katla volcano data to define the changepoints. The aim is clear and the method is suitable. My major suggestion is that the way how they present their method and its application is hard to follow and can be re-organized so that the idea can be clearly conveyed and understood. I have the following suggestions and comments:

1. The authors give a clear basic background introduction to 'double-difference data' and? 'Bayes theorem'. However, besides introducing the background, this part should also serve the following sections or the whole paper. For example, one of the main focuses of the paper is to construct the covariance matrix of the error, which is only simply mentioned in the introduction part. If the introduction part can be more specific, say focus on the methods used in the paper, the readers are more clear how these methods fit in the current Bayes theorem, and the following parts can also be much easier to be understood.

2. Following #1, I also found that the Data and Method parts seem to be somehow mixed. Maybe first clarify the methods and then talk about their application would be better, especially for a method-focused paper.

3. As for the inconsistency when compared with the cross-correlation (CC), my understanding is that this paper's method can reflect the similarities/differences between any traces (because the covariance matrix is used), however, the CC is calculated according to one major event. I wonder, what if you calculate the CC within each window, instead of with the largest event? Whether the comparison will be more consistent or not?

4. As for the inconsistency when compared with the seismicity rate, I think it might be dangerous to draw the conclusion that 'the time-history of the seismicity rate should be carefully evaluated', based on the current comparison and evidence. From my understanding, the data the authors used are the individual traces, and the time-related information, also in the data, is not considered (e.g., what is the time interval between two traces). If this is the case, in my opinion, the authors are looking at 'different parts' of the data from the people who study seismicity rates. For example, assuming a synthetic case (two time periods) with a constant subsurface elasticity but with a different stress/strain field. The waveforms in these two time periods should be exactly the same (therefore no changepoint is defined from this paper's method), but the seismicity rate can be different. Therefore, the discussion on seismicity rate needs more considerations.

5. A minor question about the terminology. From my understanding, when there is one event and one station, and we compare the differences between the observed data and the synthetic data, this can be valued as 'single-difference'. And when there is a station- or event-pair, this is called 'double-difference'. But in this paper, only the observed data (observed travel time) is mentioned. I wonder if it is still suitable to call it double-difference data when the synthetic data part is not involved in?

6. Some of the typos:

Line 100: Such solution have (has)

Line 108: where p(m|d) represent (represents)

Line 173: we compute compute (delete one of the compute)

To conclude, I think the aim of this paper is clear: how to find changepoints of the data, and their method should be suitable. However, the current structure prevents readers from fully understanding what is and how to exactly use their method. I think a major revision, especially on the structure (how to better present the method), is needed to make it a better paper. The method should be easier to be understood, rather than harder, through the authors' re-construction. Then, the following-up applications and the scientific findings behand have the chance to be addressed, as it is cleared mentioned at the beginning.

Please see attached for my review