Pintori et al. 2021

The manuscript “Common mode signals and vertical velocities in the great Alpine area from GNSS data” by Pintori et al. 2021 utilizes variational Bayesian Independent Component Analysis to separate loading signals due to atmospheric and hydrologic sources from vertical GNSS time series to further isolate uplift rates due to mountain building in the greater alpine region. This also serves to provide additional information on the source of the common mode signals in GNSS time series. In this study, they also compare the vbICA method to other methods of constructing the common mode. Finally, they remove the non-mountain building signals from the data and estimate uplift rates. Overall, the article is well written and is a useful contribution. However, there are a few areas that would benefit from additional clarification and added detail.

General comments:

ICA method in general

To my understanding, an ICA is meant to separate the different sources out of a given signal. This leads me to interpret that each component would thus represent a different source. However, throughout the paper it seems to me that each component is not necessarily a singular source (eg Figure 2). This component is obviously due to some combination of tectonic trend as well as some seasonal signal likely due to hydrology or atmospheric loading. If the point of using an ICA is to separate out different sources so that you can further isolate a specific source, how can you be sure that you are fully capturing the signal you think you are (in this case the tectonic signal or later on NTAL/HYDL)?

To this end, why are you decomposing the “source” signals (NTAL and HYDL)? Wouldn’t a component from the GNSS decomposition represent the NTAL signal? Or the HYDL signal? And then why do you combine NTAL IC1 and HYDL IC1 and compare them to the GNSS IC1? This implies this component is a portion (and only a portion) of two very different sources. How do you know that’s all that’s in there? I suppose, what I’m asking is some further clarification in the text about (1) what the different components actually mean in terms of “sources” insofar as are they “sources”? or just spatially independent signals/temporally independent and thus could be heavily influenced by certain things but not necessarily the entire signal (2) further explanation for the motivation behind decomposing the source signals (NTAL/HYDL) and why it’s necessary. I realize some of this is not specific to this paper but ICA in these applications in general but I think the text would greatly benefit from further explanation.

Decomposition of the NTAL/HYDL

The spatial pattern of the different components from the NTAL and HYDL are incredibly similar. Is this due to how the algorithm works or are these signals just by coincidence showing very similar spatial patterns. How much of the variance due each of these components represent? I think including that, maybe even just in the figures would be helpful for interpretation of the different components.

Temperature

I agree that the fourth component is well correlated with temperature. However, temperature is just a strong seasonal signal so couldn’t this signal be something else? In lines 369-370, you mention that when temperature increases the stations in the mountains subside. I’m just confused by what physical mechanism would cause this. The two mechanisms that you list for temperature in lines 505, don’t explain why the mountains would experience downward deflection during warm periods. Can you provide further explanation for the physical cause of this? I think in the paper you indicate too heavily that this component is due to temperature fluctuations (especially in Figure 8 and the associated text, the conclusion and abstract) and don’t necessarily support this. Correlation does not always indicate causation. I think further data and text is needed to support this finding. Especially since this is mentioned in the abstract (line 16) as well as the conclusion (line 586/593).

Application of vbICA for removing NTAL and HYDL

Martens et al 2020 (J. of Geodesy) highlighted the importance of removing NTAL and NTOL signals from GNSS timeseries to reduce scatter/dispersion. In lines 351-353, you mention, vbICA may not be able to separate the NTAL vs HYDL signals. Why not just remove the signals using the GFZ products instead of using the ICA method? Does removing the ICA reduce the scatter more than just removing the signals to begin with? -

Minor comments:

GNSS processing - Do remove signals due to earthquakes? In the supplement you mention removing offsets due to equipment changes but don’t mention offsets or post seismic signal removal. Does the ICA capture earthquake signals? Wouldn’t this be a good signal to remove to better isolate the uplift?

Line 123: grammatical issue - “Since they allow to account”

Line 254: grammatical issue – “its temporal evolution has not a domination frequency”

Lines 230: What reference frames are you using for the NTAL and HYDL models?

Lines 250-251: There are no units for y-axes on the temporal portion of the components. What are the units? Are there any? To construct the signal at a given spot, do you multiply the temporal by the spatial displacement for that point? It would be helpful for understanding the figures.

Lines 340: For the second component, you list Pearson’s correlation coefficient in addition to the Lin’s. Can you list the Pearson’s for component 1? And in the third component, is this the Pearson’s or the Lin’s coefficient?

Line 338: How many stations have displacements above 3mm?

Line 389: Is the k value -2 for both?

Line 404-408: I think there’s a typo here.

Consequently, the unfiltered time series are modeled only with the linear trend plus the temporal correlated noise, while the unfiltered time series modeling annual and semi-annual terms are also included.

Are both unfiltered? I think the first one should be filtered, yes?

Section 5.3: If you are removing the linear trend, then are your uplift rates non-tectonic uplift? Or are you adding that back in? Just confusing since in the introduction it seemed like you were settling up to better estimate uplift rates due to tectonics? I think it’s fine to remove the linear trend for comparison of stacking methods ect but for Figure 13 and discussion in 5.3 is this with the linear trend removed or included? I the nontectonic uplift? Or are you adding the linear trend back in? Can you clarify?

Many of the figures appear blurry. Additionally, the font on the axes of many of the figures is incredibly difficult to read (eg Figure 3) and would benefit from larger font size.

Pintori et al. use a version of the ICA method (called variational Bayesian ICA) to decompose vertical GPS position time series and hydrology/atmospheric predicted loading time series around the European Alps. They study the agreement between the ICs extracted from the GPS series and from the loading models for the period from 2010 to 2020. Their main conclusions are that 1) the vertical GPS series can be separated in a tectonic linear motion and variations caused by temperature and atmospheric/hydrology loading; and that 2) improved tectonic velocities are obtained by correcting the GPS series using ICs obtained from the GPS series themselves.

While the volume of work is of note, especially concerning the GPS data processing, I do not think the conclusions are supported by the data and methods used by the authors. It is reasonable to say that temperature variations, atmospheric pressure variations and hydrology load variations contribute to the variations observed in vertical GPS time series, especially at the annual period, as GPS positions react to these and many other phenomena together. A completely different thing is to say that the observed GPS variations of vertical position *are* originated or explained by these processes, as the authors repeatedly state in the manuscript. This is a clear misinterpretation of their analysis and I develop my reasoning in the paragraphs below.

Before that, and assuming conclusion 1 is right, it’s very surprising that the authors do not try to remove the modeled loadings from the GPS series to test the impact on the estimated velocities. Instead, conclusion 2 is based on removing the GPS ICs from the GPS series, i.e., conclusions 1 and 2 are totally unrelated. The GPS ICs were obtained from GPS series that were previously detrended, explaining the small change of the estimated velocities from the filtered series. The ICA filtering also explains the reduction of the noise in the series and, therefore, of the estimated velocity uncertainty from the filtered series. Where I think this approach fails is that the raw series (used to estimate the velocity, the filtered velocity being very similar) and the filtered series (used to re-estimate the velocity uncertainty) are not consistent and therefore the velocity and its “improved” uncertainty are not consistent either. The authors could have tried a more aggressive filtering, like a band-pass filter leaving the trend and high-frequency noise only, or could have not consider colored noise in the velocity estimation (both ways are equivalent) and they will get even smaller velocity uncertainties. Unfortunately, this will not give any valuable information on the quality of the velocity and your ability to extrapolate it to understand tectonic physical processes. The only way to improve velocity estimates is to understand and reduce variability in the GPS series with proven corrections and models. If the white noise is more visible in the filtered series is probably because the GPS ICs absorb together a significant portion of the power-law noise that typically dominates the variance of the detrended GPS series, though this is not very clear from the IC PSDs in Fig. 3. Precisely, the power-law noise in the GPS series is only mentioned briefly and its influence on the GPS ICs and on the correlation with the loading ICs is not discussed at all.

With respect to the GPS ICs and their attribution of a geophysical origin, I enumerate below several points raising concerns on the authors’ approach. Generally, many past publications have shown than GPS series and loading models do not see the same thing, except partly for the annual variation. Most of the variance in the loading model series is concentrated at the annual period. Compared to the PSD of the loading models, the GPS series contain a relatively higher variance at long periods with a distinct PSD slope and a PSD much richer in periodic artifacts at short periods. The authors briefly comment on the systematic errors that are present in the GPS series, but they do not try to make the GPS series more consistent with the model series. For instance, it is known the annual draconitic variation could significantly affect the comparison to the solar annual variation of the loading models. The results obtained by the authors are confusing (see points below) and do not refute findings from past publications, contrary to their claims to successfully separate geophysical signals from the GPS series. For instance, authors show no evidence that the HYDL series significantly explain variations in their GPS series. The GPS and NTAL annual seem to partly agree (see points below), so the authors introduce a thermal annual component in the discussion without providing strong evidence nor explanation of its spatial pattern. It is also probably worth mentioning that, if the GPS series were effectively explained by the combination of atmospheric/hydrology loading and temperature variations, as the authors claim, we should get the same GPS series out of the same GPS data when using different software, different strategies and different corrections. However, this is often not the case, especially when comparing global and regional GPS solutions.

Other general points:

1. While I understand the objective of the ICA applied to the GPS series is to separate the variability into independent processes, I cannot understand the rationale for applying ICA to the NTAL and HYDL series. What are the independent processes to be separated in the atmospheric pressure loading or water loading? Even more confusing are the results from the comparison of a single GPS IC to a single NTAL/HYDL IC and the claim that the GPS series are explained by both. The ICA analysis is forcing the NTAL/HYDL series into non-gaussian independent components, even if they do not exist physically. This probably explains why the total NTAL annual is split across ICs with spatial patterns as orthogonal as possible. The same spatial patterns are found for the GPS series, probably because once the trend, offsets and annual are removed from the GPS series, what is left is a Gaussian or near Gaussian series with temporal & spatially correlated noise and also the above-mentioned systematic periodic errors. It may be that the easiest way for the ICA to force the separation of these residual series into ICs is by making their spatial patterns orthogonal (see another possible explanation in point 5 below). The authors’ conclusion that GPS and loading see the same spatial patterns is therefore not very solid.
2. The GPS and NTAL/HYDL series have different spatial samplings, which must complicate the interpretation of their comparison. Also related to the spatial sampling, it must be difficult to extract accurate NTAL values in the Alps due to the pressure model resolution and the short-scale changes in topographic gradient, making its comparison to the GPS series even less trustworthy. I suspect similar limitations exist when comparing GPS and HYDL model series in a mountain range.
3. Each dataset used by the authors is decomposed in different numbers of ICs: 7 for GPS, although only 4 are discussed, and 3 for the model loadings. Then they compare the first 3 individual ICs and find weak correlations between them. The authors conclude on the origin of the individual GPS ICs based on their correlation to the individual loading ICs. However, this criterion is very weak, especially with correlation values around 0.6. As an example, similar (Pearson’s) correlation values would be obtained between a pure sinusoidal and the same sinusoidal delayed almost pi/3, which is roughly two months if the sinusoidal has a period of one year. When subtracting one sinusoidal from the other, it is clear that we are not correcting much. The ratio of explained variance between the different ICs would have been more appealing, but, it is not clear that the individual ICs from different datasets correspond to the same fraction of the total signal (see point 1). So maybe the ICA method is not well adapted to this problem or should not be applied to the NTAL/HYDL series (see point 1). A band-pass filtered comparison of GPS and loading series would probably be more informative here. Also rather than filtering the GPS series, I think it would have been better if the authors had shown how the loading models change the variance of the GPS series, as it is done in many other publications. The loading would need to be computed at the station locations. It would have been even better to show how the GPS variance changes (not necessarily reducing) all along its power spectrum when correcting the loads.
4. The authors are processing a regional network and aligning it to a global linear frame (IGb14) that does not include seasonal variations. The frame alignment of the daily solutions from regional networks acts as another CME-like filtering of the series, not discussed by the authors, but probably similar to the SFM method. The filtering is more efficient as the network size is smaller, but the authors do not provide enough information on this point. It is then difficult to interpret the common network-wide annual signal shown by the GPS IC1. I would expect the regional frame alignment would absorb part of this common GPS annual signal, making it difficult to compare to the loading model and also leaving an amplitude much smaller than the residual station-dependent annual signal that is probably captured by the IC4. However, the numbers in table 1 indicate the opposite, assuming the average “of the amplitude of the maximum displacement” is somehow related to the annual amplitude, which is not clear either. The annual variation is the most prominent signal in NTAL with amplitudes typically of a few mm, less than 1 cm at the center of large continental masses. So it’s not clear what the authors mean with atmospheric loading amplitudes larger than 2 cm. It is also not mentioned which frame was used to create the loading series and whether they were detrended like the GPS series, especially the HYDL series.
5. The 2nd and 3rd GPS ICs are particularly interesting. These represent daily E/W and N/S network tilts with a rather flat spectrum. The NTAL and HYDL show similar spatial tilts, but their physical meaning is dubious (see point 1) and their spectral content is completely different: mostly seasonal for NTAL and mostly interannual for HYDL. The origin of these network tilts is very likely not the same among the datasets, as stated by the authors. In addition, if the whole GPS network is truly moving like these two ICs and it is not an artifact of the ICA separation, I would first think of a problem with the reference frame alignment. As said in point 4, network-wide common mode signals, including daily tilts and annual up & downs, should be at least partly (if not totally) absorbed by the frame alignment as these signals are not included in the linear reference frame and the network size is probably not large enough. Figure 7b must be wrong as there is no annual variation in the GPS IC2.

The manuscript "Common mode signals and vertical velocities in the great Alpine area from GNSS data" by Francesco Pintori et al. presents

how ICA decomposition of GNSS time series in the alpine area allows to separate sources of deformation and then retrieve with a better

uncertainty the velocity field in Europe. The authors process the daily  GPS observations with GAMIT/GLOBK software, using subnetworks later tied to IGb14 reference frame. The obtained 2010-2020 time series have then been analysed in order to explore the origin of the common modes, and the potential of Independant Component analysis to extract these modes with a more "physical" basis and filter the time series. The ICA method used in the paper is the vbICA, a bayesian multivariate source separation method. The ICA analysis conducted here is performed in two steps, one, with 8 components, allows to extract and correct the trend (the velocity), the other, with detrended GNSS data as input, contains seven components. In parallel, hydrological and atmospheric loading predictions from two institutes are also analysed with vbICA with three components. These three components corresponds mostly to a uniform spatial pattern, an E-W trend and a N-S trend. The GNSS components appear well correlated to the hydrological plus atmospheric loads components, proving the loading origin of these components. A last component is clearly seasonal and presents spatial variation at small wavelength, in phase with temperature variations. The four vbICA components are used to correct the GNSS time series, which allow a new estimation of the velocity, in very good agreement with the first estimation but with a much smaller error estimation. The authors also compare different methods for common mode estimation, the stacking Filtering method, or weighted stacking filtering method to the filtering obtained by an Independant component analysis.

Overall I found the manuscript interesting and worth of publication, as it shows a convincing correspondance between what is referred as "common modes" and the atmospheric and hydrological loading. However, I think that the paper, although well written, is quite hard to follow, with numerous abbreviations, and comparisons which could be better presented and illustrated. I have also a few scientific comments that can be adressed. I suggest a major revision.

Here are my suggestions:

* I find intriguing that the main three components that are discussed here correspond to a uniform pattern, an E-W tilt and a N-S tilt.

These three components correspond to the largest perpendicular spatially correlated signals possible.

(1) Can you change the color scale of all panels of IC1, to show how uniform it really is ? For example GNSS IC1 should be plotted with a

20-32 scale.

(2) For IC2 and IC3, how significantly different from a tilt the components are ?

(3) the loading models appear to predict mainly very long wavelength features, corresponding to the first three components.

Is this true ? Can you show an example of the predicted load-induced displacement map ? The percentage of the variance do the three components is indicated to be > 97%. For atmosphere, I guess pressure variations are large-scale such that the earth response is also at large-scale. But I would have thought that hydrological loading should be more local. Can you comment on that ?

* The seasonal contribution should not be named temperature contribution. This would suggest a thermal contraction effect which is far from being proven. A lot of signals could be seasonal. Unless you prove that there is a strong correlation between the IC4 and temperature beyond the seasonal term (ie at higher frequency) the correlation appears fortuitous. Fig 8 shows that temperature seems to have higher frequency fluctuations not observed in IC4, but it s hard to tell from the figure only.

I suggest to rewrite the paragraphs and sentences related to this seasonal contribution of unknown origin everywhere in text.

* The statistics shown (mean, median, standard deviation) in tables and discussed in text are not well presented. I suggest to move S4

in the main text, it is quite graphical and shows better the agreement in terms of distribution than Tables 1 and 2, that could be moved to supplementary material. lines 289 to 292 could be replaced by a more readable text.

* The part on correlation coefficients is confusing where it should not.

If you consider that your signal is a sum of IC like Xi(x,y)*Ti(t), then we expect to provide the correlation coefficient

between Ti-GNSS and Ti-HYDR for example, or Ti-GNSS and Ti-ATM, and of Xi-GNSS with Xi-HYDR or XI-ATM. Only two values describing the temporal and spatial correlations would be sufficient. Here, it took me time to understand that, because you add Xi_ATM(x,y)*Ti_ATM(t) and Xi_HYDR(x,y)*Ti_HYDR(t), your spatial and temporal correlations stop being independant from each other. This is why I guess you provide ion Fig6 a spatial map of the temporal correlation of the GNSS and HYDR+ATM.  Could you please clarify for the reader why you end up with such a plot ?

In fact, if you had made and ICA on (ATM+HYDR) directly, may be you would have obtained a similar result but easier to compare (ie an independent comparison in space and time).  The "blue points" on fig. 6 in the middle of the tilt, in opposite phase, have no real significance, as the spatial patterns of ICs do not exactly correspond to each other. I find more significant the peak in the ditribution, of 0.65 for IC2 and of 0.55 for IC3 which are significant numbers although the PSDs of the Ti do not really match.

* Once ATM and HYDR loads are proven to be good estimators of the common modes, why not use them to correct the time series ?

The advantage is that you can then anticipate that possible decadal trends of ATM and HYDR would then be removed from the time series and thus provide a better displacement rate due to tectonics. Here, the trend is first estimated from a first ICA, removed from GNSS time series, and then a new ICA is performed to extract ICs, that will correct the raw GNSS data, before a new trend estimation. How can you be sure that the last estimation will not be "by construction" biased towards the first ? On the other hand line 219-220 of 3.1 suggests that the separation of tectonics trend from other potential non tectonic trends is already done by the first ICA. Can you clarify this point ?

Figures :

ICA figures:

- change color scales of IC1 for all plots to show lateral variations

- temporal vector: normalisation should be made by variance and not by min/max (if I understood correctly) for the reader to visualie the relative amplitude of each term. Min/max can be outliers.

Figure 6: change colorscale to see changes in correlation coefficient for IC1 (the colorscale is completely saturated in the red).

Don't use "Lin" abbreviation but linear

Figure 7: panel b is identical to panel a

Abstract:

First sentence : too complicated. Simplify and clarify

line 10: associated with : modeled from

line 11: processes: drop

line 16-17 : Atmospheric .... gradients: rewrite

Introduction

First sentence: "active geophysical processes on land, ice and atomosphere": ground displacement on atmosphere. Rewrite.

In general : a lot of references are missing on mountain uplift, both observations and mechanisms. Please provide some refs outside Italy.

Id. for lines 68-80

line 117: give principle of CMC Imaging

line 190: pdfs --> PDFs (and elsewhere)

line 192: drop "that"

line 216: a priori any temporal : rewrite

line 389: k=-2 for both noise and flicker : correct text

line 391: avoid + in text

line 506: elastic hydrological load ---> elastic response to hydrological load

* Don't use "lin" abbrevation but replace by linear correlation coefficient.