We appreciate the reviewer's valuable comments, which have helped us improve the clarity and quality of the manuscript.

Reviewer, Comment 1:

Figure 1: Abbreviation AM is not explained/introduced.

Response:  

Thank you for pointing this out. We have now introduced and explained the abbreviation “AM” (which stands for Apuseni Mountains) in the figure caption and in the main text we added the abreviation. 

Reviewer, Comment 2:

4. Analysis of the GNSS position (coordinates) time series: It seems that the estimation of velocities, jumps, annual and semi-annual parameters is based on a white-noise assumption, i.e. assuming no time correlation. However such time series usually contain significant stochastic time-correlation. Neglecting this might lead to over-optimistic confidence regions (as maybe the case in figure 3) and implausible statements about the significance of parameters. Have the authors considered this aspect?

Response: 

Thank you for raising this point. In the current analysis, the method is quite straightforward. Our algorithm estimates the velocity, seasonal signals and jumps in a single batch solution for each time series. The accuracies (sigmas) are derived from the WRMS of the residuals. They are computed as 2/3*WMRS divided by the length of the time series. This avoids overly optimistic estimates of the accuracy. It has proven to be a reliable and robust method, yielding realistic velocity and accuracy estimates (sigmas), comparable to literature solutions using different algorithms.

Reviewer, Comment 3:

4.1 For the modeled jumps in the time series: Have the authors tried to identify the reasons behind (antenna change or some local effect, ...)?

Response:

We thank the reviewer for this observation. The possible causes of the modeled jumps, including antenna changes and local site effects, are discussed in the manuscript (see Section 4.1, lines 168–173). We have now revised the text slightly to make this discussion more prominent and easier to locate. 

Reviewer, Comment 4:

4.2: A reference to table S3 in the supplementary material is missing.

Response:

Thank you for your observation. We have revised the manuscript to include an appropriate in-text reference to Table S3 in the Supplementary Material. This should now clarify the connection between the main text and the supplementary content.

Reviewer, Comment 5:

4.4: The reference mentioned for the STRAINTOOL is not given in the References of the manuscript.

Response:

Thank you for pointing out the omission. We have now included the full references for STRAINTOOL in the References section as shown below.

Anastasiou D.G., Papanikolaou X., Ganas A., and Paradissis D.: StrainTool: A software package to estimate strain tensor parameters. Zenodo; Version 1.1, https://doi.org/10.5281/zenodo.5501234, 2021

Shen, Z. K., Wang, M., Zeng, Y., and Wang, F.: Strain determination using spatially discrete geodetic data, Bull. Seismol. Soc. Am., 105, 2117-2127, https://doi.org/10.1785/0120140247, 2015.

Reviewer, Comment 6:

Figure 3: Some large velocity vectors seem to deviate significantly from their neighbors (their motion pattern), e.g. ORAD, CLU2, TGMS. Shouldn't these deviations be addressed?

Response:

We appreciate the reviewer’s observation. As discussed in Section 4.2, we applied an outlier detection method to identify stations exhibiting significant deviations from the regional trend. This includes stations such as ORAD, CLU2, and TGMS. These stations are now represented as black vectors in an updated version of Figure 3.

As mentioned, partially based on your comments, we made many changes in the text of the final manuscript to give better explanations. This version is not online (yet).

We appreciate the reviewer's valuable comments, which have helped us improve the clarity and quality of the manuscript.

Reviewer, Comment 1:

Figure 1: Abbreviation AM is not explained/introduced.

Response:  

Thank you for pointing this out. We have now introduced and explained the abbreviation “AM” (which stands for Apuseni Mountains) in the figure caption and in the main text we added the abreviation. 

Reviewer, Comment 2:

4. Analysis of the GNSS position (coordinates) time series: It seems that the estimation of velocities, jumps, annual and semi-annual parameters is based on a white-noise assumption, i.e. assuming no time correlation. However such time series usually contain significant stochastic time-correlation. Neglecting this might lead to over-optimistic confidence regions (as maybe the case in figure 3) and implausible statements about the significance of parameters. Have the authors considered this aspect?

Response: 

Thank you for raising this point. In the current analysis, the method is quite straightforward. Our algorithm estimates the velocity, seasonal signals and jumps in a single batch solution for each time series. The accuracies (sigmas) are derived from the WRMS of the residuals. They are computed as 2/3*WMRS divided by the length of the time series. This avoids overly optimistic estimates of the accuracy. It has proven to be a reliable and robust method, yielding realistic velocity and accuracy estimates (sigmas), comparable to literature solutions using different algorithms.

Reviewer, Comment 3:

4.1 For the modeled jumps in the time series: Have the authors tried to identify the reasons behind (antenna change or some local effect, ...)?

Response:

We thank the reviewer for this observation. The possible causes of the modeled jumps, including antenna changes and local site effects, are discussed in the manuscript (see Section 4.1, lines 168–173). We have now revised the text slightly to make this discussion more prominent and easier to locate. 

Reviewer, Comment 4:

4.2: A reference to table S3 in the supplementary material is missing.

Response:

Thank you for your observation. We have revised the manuscript to include an appropriate in-text reference to Table S3 in the Supplementary Material. This should now clarify the connection between the main text and the supplementary content.

Reviewer, Comment 5:

4.4: The reference mentioned for the STRAINTOOL is not given in the References of the manuscript.

Response:

Thank you for pointing out the omission. We have now included the full references for STRAINTOOL in the References section as shown below.

Anastasiou D.G., Papanikolaou X., Ganas A., and Paradissis D.: StrainTool: A software package to estimate strain tensor parameters. Zenodo; Version 1.1, https://doi.org/10.5281/zenodo.5501234, 2021

Shen, Z. K., Wang, M., Zeng, Y., and Wang, F.: Strain determination using spatially discrete geodetic data, Bull. Seismol. Soc. Am., 105, 2117-2127, https://doi.org/10.1785/0120140247, 2015.

Reviewer, Comment 6:

Figure 3: Some large velocity vectors seem to deviate significantly from their neighbors (their motion pattern), e.g. ORAD, CLU2, TGMS. Shouldn't these deviations be addressed?

Response:

We appreciate the reviewer’s observation. As discussed in Section 4.2, we applied an outlier detection method to identify stations exhibiting significant deviations from the regional trend. This includes stations such as ORAD, CLU2, and TGMS. These stations are now represented as black vectors in an updated version of Figure 3.

As mentioned, partially based on your comments, we made many changes in the text of the final manuscript to give better explanations. This version is not online (yet).

Dear reviewer,

Thank you for your constructive comments and suggestions.

We have made the necessary modifications throughout the manuscript based on your further suggestions.

We have updated all figures in the manuscript except for Figures 6 and 8, which remain unchanged as they already accurately represent the intended data and findings. You can view all these changes in the updated version of the manuscript, which will be uploaded to the journal platform. 

Point-to-point answer to your questions/suggestions:

Q1. This well-written paper describes results from 20 years of continuous GPS measurements in Romania. This is a useful addition to the existing literature. Some improvements can be made to the reference frame definition and the uncertainty estimation. The comparison with previous results is there, but could be expanded. I am not fully convinced by the vertical velocities and their interpolation. The interpretation is fine, with a range of hypothesis provided in terms of the underlying dynamics. 

It may be useful to show seismicity on Figure 1, together with the location of GPS stations but without colors and site names. This would be particularly useful to accompany the introduction and tectonic setting sections of the paper. The details of the cGPS site names and network affiliations can be shown in the supplements, as they are not essential to the understanding of the paper.

A1. Thank you very much for your positive feedback! We have replaced Figure 1 with a map showing station locations and seismicity, as suggested. 

Q2. My first significant comment concerns the definition of the Eurasia-fixed reference frame, briefly described in lines 148-150. First, the publication from which ``the most recently published rotation pole solution for that plate'' is used should be cited. Second, one needs to know how well the velocities presented here fit a Eurasia-fixed frame. In other words, do stations allegedly on stable Eurasia have velocities that do not statistically differ from zero? This is particularly important since the ``residual'' velocities discussed in the rest of the paper are very small, on the order of 1~mm/yr.

For instance, I am surprised that sites located on the East European Platform -- supposed to be stable Eurasia? or perhaps not? -- show a consistent NW-direct residual on the order of 1~mm/yr.

The authors may be using the Eurasia-ITRF14 angular velocities from Altamimi et al. 2017 (ITRF2014 plate motion model), but the consistency between Altamimi's frame definition and the authors is unlikely to be accurate at the 1~mm/yr level since they are (most probably) not using the same models.

The most rigorous way to determine such a frame would be for the authors to include in the analysis a number of sites located on stable Eurasia, for instance those used by Altamimi et al. 2017. Then the authors could compute their own Eurasia-ITRF14 angular velocity and use it to defined a Eurasia-fixed frame consistent with their analysis. But if the authors choose to use an independently-published angular velocity, then they must shown that residuals are close to zero at sites located on stable Eurasia.

A2. We converted the reference plate of our solutions to the Eurasian tectonic reference plate using the ITRF14 rotation parameters of that plate (Altamimi et al., 2017). The reference will be included in the final document. This plate model has proven to be a reliable model for our analyses since it also works perfectly for other areas in Europe. So, there is no need to compute our own plate model.

We were also surprised by the general NE motion of the EEP. But our results are consistent with Piña-Valdés et al. (2022), showing the same trend for a larger area around Romania. However, that is not a subject of our study.

In response to your question, we made a careful comparison of our results versus Piña‐Valdés (87 stations with common solutions). For a few stations, we find differences of up to 2.0 mm/yr. These are probably due to the differences in the length of the time series. But on average, the differences are less than 0.08 mm/yr.

We also note that the residual velocities on nominally stable Eurasian sites are on the order of ~1 mm/yr with formal uncertainties of ~0.2 mm/yr, making them statistically significant at approximately the 5σ level. These velocities are derived from multi-year (>4 years) GPS time series, which substantially reduce random noise and provide robust linear trends. Even in a globally defined Eurasia-fixed frame, small non-zero residuals can naturally arise from minor intraplate deformation, local geophysical effects, or slight differences between our network realization and the published rotation parameters. The observed residuals are therefore consistent with expected variations within the Eurasia-fixed frame and support the broader patterns of intraplate motion discussed in the manuscript.

Q3. My second significant comment concerns the estimation of velocity uncertainties, not explicitly described in the paper. One can assume from lines 156-158 that the uncertainty is that of a least-squares fit to the time series, assuming linear and seasonal terms plus offsets. If it is the case, then how do the authors deal with the well-known issue of time-correlated noise in cGPS time series and its impact on velocity uncertainty estimates? I suspect that the uncertainties estimated here actually underestimate more realistic ones that would account for time-correlated noise. That is an important issue in a region where residual velocities are very small.

Several published options exist to fit GPS time series using a noise model that accounts for time-correlated noise, for instance with CATS (Williams GPS Solutions 2008) or Bos et al. J. of Geodesy 2008.

As a side note, the authors should quote the confidence interval they choose instead of using "sigma". Note that a 1-sigma uncertainty in 2 dimensions (east-north velocities for instance) correspond to a 39% confidence interval. A much safer 95% confidence interval would correspond approximately to 2.4-sigmasi in 2 dimensions.

A3. We responded to the first reviewer who had a similar question about our method of calculating uncertainties. Our method is quite straightforward. Our algorithm estimates the velocity, seasonal signals, and jumps in a single batch solution for each time series. The accuracies (sigmas) are derived from the WRMS of the residuals. They are computed as 2*WRMS/2.4/Tspan, where Tspan is the length of the time series in years, resulting in 95% confidence level accuracies. This avoids overly optimistic estimates of the accuracy. It has proven to be a reliable and robust method, yielding realistic velocity and accuracy estimates (sigmas), comparable to literature solutions using different algorithms. We clarified this in section 4.1 of the manuscript.

Q4. Section "Estimating a gridded, smoother horizontal velocity field": how is the uncertainty on actual GPS velocities propagated to the gridded velocities? By the way, Figure 3 should show uncertainty ellipses for the gridded velocities.

A4. The formal uncertainties of our GPS horizontal velocities are on the order of 0.02 mm yr⁻¹, whereas the velocities themselves are typically 1-2 mm yr⁻¹. The propagation of such small uncertainties would have a negligible effect on the resulting velocities. For this reason, we did not explicitly compute or display uncertainty ellipses for the gridded field. We will, however, clarify in the text that the station-level uncertainties are small relative to the observed signal and therefore do not significantly influence the derived deformation pattern.

Q5. I understand in section 4.4 that the strain rates are computed from the raw GPS velocities, not the gridded ones. Is that right? If so, what is the point of calculating gridded velocities? By the way, I assume that STRAINTOOL estimates velocities on a regular grid as well. Perhaps those should be shown on Figure 3? But I may be misunderstanding section 4.4 -- which perhaps calls for some clarification in the text

A5. In section 4.4, we mentioned: “we further estimate strain rate from the interpolated horizontal vector field of GPS velocities [...] This algorithm interpolates our gridded solution to derive horizontal strain rates across the region, using a weighted least squares approach on a more dense regular grid.”

The gridded velocity field was used as input for the STRAINTOOL. We modified the text to clarify this.

Q6. Section 5.1 and throughout the paper: I suggest replacing "appear to show" by simply "show". The agreement between solutions seems ok in the horizontal, but there seems to be some scatter in the vertical?

A6..We have revised the text in Section 5.1 and throughout the paper to replace "appear to show" with "show", where appropriate. This change improves clarity and strengthens the wording.

Q7. The comparison of the velocities in this work with those from Pina-Valdes et al. 2022 and Serpelloni et al. 2022 on figure 7 must show the 3 different solutions in different colors.

A7. As for Figure 7: It does not make sense to use 3 different colors. All vectors around and south of Romania are Piña-Valdés, D’Agostino, and Serpelloni solutions. We did not analyze those stations ourselves because we don’t have access to the original GPS data. So, we see no need to plot Piña-Valdés, D’Agostino, and Serpelloni in different colors. In the final figure, all vectors will be plotted in black. The purpose of this figure is only to show our Romanian solutions in a larger geodetic setting.

Q8. I suggest the authors also review and cite D'Agostino et al. Earth and Planetary Science Letters 539 (2020), which intersects the study area presented here.

A8. For D’Agostino et al. (2020, Earth and Planetary Science Letters), which indeed covers a region overlapping with our study area, we found the study to be relevant and informative. We did include these solutions in our analyses of the regions around Romania. Inside the Romanian territory, there are only 11 common solutions, with most of them with large uncertainties due to the short time series.

Q9.Line 285: I see no complexity in the strain rate field, which is rather smooth, with a consistent NS extension throughout and significant EW compression in the SW corner. I would actually like to read a complete interpretation of Figure 5. Is there consistency between the earthquake focal mechanisms and the GPS-derived strain rate principal axes? What do we learn about the sense of slip on the main active faults? Are compressional events in Vrancea consistent with the principal axes of the geodetic strain rates?.

A9. We described this relationship in the manuscript in chapter 6.2 “Correlation with fault systems and seismicity”. While the GPS-derived strain field is indeed relatively smooth at the regional scale, earthquake focal mechanisms and previously published stress inversion results (Petrescu et al., 2021) reveal local pockets of compression and complex faulting. These focal mechanisms often occur on secondary or inherited faults and accommodate transitions between compressional, strike-slip, and extensional regimes away from Vrancea. The main active faults primarily accommodate the regional N-S extension through normal and oblique-normal slip, while smaller structures adjust locally to accommodate deviations from the regional strain.  It is a complex challenge to jointly interpret them, mostly because short-wavelength deformation is superimposed on the smooth, long-wavelength strain captured by GPS, yet there is no other attempt to do this in a region where destructive earthquakes occur in a seemingly stable continental setting, where previous GPS studies relied on short-term campaign data.

Q10. Figure 5 could be improved by making the "beachballs" smaller, as they are quite cluttered. Also, it would be useful to show on the figure the depth of these events, as I assume that some are crustal but others much deeper. Perhaps color-coding the black quadrants of the focal mechanisms as a function of hypocenter depth would help?

A.10.  We have now coloured the focal mechanisms by depth, as suggested, and scaled them according to magnitude. 

Q11. I am a bit puzzled by the vertical velocities shown on Figure 4 and its interpolation. For instance, the uplift observed on the East Carpathians essentially depend on 2 sites with vertical velocity significantly different from the neighboring sites. Can we be certain that one is not looking at local anomalies due to monument stability or other local process? It is a bit troubling that this does not appear in Pina-Valdez' solution, for instance. The same holds for the SE Carpathians foredeep.

A11.  We are not certain which specific stations the reviewer is referring to. The uplift observed in the Eastern Carpathians is supported by a cluster of stations (not just isolated ones), including ARYB, VAMO, VATR, TPLT, BICA, PTNT, and BIST, all located within a 2°×2° area. South of this cluster, another group of stations (GHE, CIUC, TGTS, BACA, SFGH) shows a consistent downward trend. 

To ensure data reliability, we excluded all stations displaying anomalous time series or evidence of monument disturbance. Additionally, we only retained stations with at least four years of stable observations, ensuring that the derived vertical velocities reflect long-term signals. These quality control criteria are detailed in Sections 4 and 4.1. 

From a geological perspective, the observed uplift in the Vrancea region, along the Putna River, is consistent with independent thermochronological evidence reported by Necea et al. (2005, 2021). The river profile shows that there is a difference between the knee point and river water divide; hence, the river is now trying to reach equilibrium, hence the observed uplift can be explained. This geomorphic signal supports an interpretation of ongoing regional uplift, in agreement with the geodetic results presented here.

Necea, D., Fielitz, W., and Matenc,o L.: Late Pliocene–Quaternary tectonics in the frontal part of the SE Carpathians: Insights from tectonic geomorphology, Tectonophysics, 410, 137-156, https://doi.org/101016/j.tecto.2005.05.047, 2005.

Necea, D.,  Juez-Larré, J.,  Matenco, L.,  Andriessen, P. A.M., and  Dinu C.: Foreland migration of orogenic exhumation during nappe stacking: Inferences from a high-resolution thermochronological profile over the Southeast Carpathians, Global Planet. Change, 200, 103457, https://doi.org/10.1016/j.gloplacha.2021.103457, 2021.

Q12. It seems to me that the authors are placing too much confidence in their interpretation on local vertical velocities, especially at local scale. This comment applies to the South Carpathians (lines 415-416), but also to the SE Carpathian foredeep.

A12. Our interpretation of localized uplift in the South Carpathians is based on multiple independent datasets, including Bouguer gravity anomalies, thermochronology, and recent studies on neotectonics and geomorphology (Necea et al., 2005, 2013, and 2021; Merten et al., 2005, 2010; Poncos et al., 2022; Gailleton et al., 2021). We have carefully considered uncertainties in these observations, and the combined evidence consistently points toward differential isostatic compensation and/or lithospheric-scale processes as plausible mechanisms. To ensure clarity, we have slightly reworded the text to emphasize that these interpretations are well-supported but not uniquely constrained

Gailleton, B., Sinclair, H. D., Mudd, S. M., Graf, E. L. S., & Mațenco, L.C.: Isolating lithologic versus tectonic signals of river profiles to test orogenic models for the Eastern and Southeastern Carpathians. J. Geophys. Res.Earth Surface,126, 8, https://doi.org/10.1029/2020JF005970, 2021.

Necea D., Fielitz W.,  Kadereit A., P.A.M. Andriessen P.A.M., and Dinu C.: Middle Pleistocene to Holocene fluvial terrace development and uplift-driven valley incision in the SE Carpathians, Romania, Tectonophysics, 602, 332-354, https://doi.org/10.1016/j.tecto.2013.02.039, 2013.

Merten, S., Andriessen, P. A. M., Juez-Larré, J., Bertotti, G. V., & Dunai, T. J.: Dating the exhumation of the Romanian Carpathians: first results from apatite (U-Th)/He thermochronology, Abstract from Geophysical Research Abstracts, 7, 08138, https://meetings.copernicus.org/www.cosis.net/abstracts/EGU05/08138/EGU05-J-08138.pdf, 2005.

Merten S., Matenco L., Foeken J. P. T., Stuart F. M., Andriessen P. A. M.: From nappe stacking to out-of-sequence postcollisional deformations: Cretaceous to Quaternary exhumation history of the SE Carpathians assessed by low-temperature thermochronology, Tectonics, 29, https://doi.org/10.1029/2009TC002550, 2010

Poncos, V., Stanciu, I., Teleaga, D., Maţenco, L.,  Bozsó, I.,  Szakács, A., Birtas, D., Toma, S.A.,  Stanica, A., and  Radulescu, V.: An Integrated Platform for Ground-Motion Mapping, Local to Regional Scale; Examples from SE Europe. Remote Sensing, 14, 1046, 10.3390/rs14041046, 2022.

Q13. Conclusion:

- The authors write "Our results mark a significant improvement in spatial coverage and resolution", "Our extended and more reliable data" -- I am wondering compared to which previous work?

A13. The earlier geodetic work in the region, based on campaign-style GNSS measurements, is described in the Introduction (lines 62 - 68). However, we agree that the connection between that background and our claims of improvement was not clearly stated in the Conclusions section. A new sentence has been added in the Conclusions section to make this comparison explicit: “Earlier studies in the region relied on campaign-style GNSS observations, Van der Hoeven et al., 2005. In contrast, our dataset includes continuous GNSS data from 133 stations spanning 20 years, providing improved spatial density and temporal resolution.”

Q14. The authors write "The relative motions between these regions generate a complex strain rate pattern with zones of extension, compression, and shear, which closely align with observed regional seismic activity." What is meant by "align", especially since there is no seismicity map in the paper? For instance, do the authors mean "coincide" or also "share the same mechanisms"?

A14. The seismicity map used to support this statement was included in the supplementary material (Figure S1). Still, we acknowledge that this may not have been sufficiently clear (included) in the main text. As stated earlier, we have incorporated your suggestion and included the seismicity in the first figure.

Q15. Throughout the paper: I suggest toning down (or removing) terms such as "fundamental", "vital", "crucial", "complex", "fresh", "significant". After all, this is only geodesy!

A15. We appreciate the reviewer’s suggestion to tone down terms such as “complex” or “significant.” However, we respectfully maintain these descriptors because the region genuinely presents a unique and scientifically challenging setting. Despite being far from plate boundaries and considered part of the stable Eurasian continent, this area exhibits persistent deformation detectable with modern geodetic methods and hosts destructive earthquakes, including in the Vrancea slab region. The last geodetic studies relied on short-term campaign data with large oscillations, published more than 20 years ago, leaving substantial gaps in our understanding of the ongoing strain field. Moreover, the interactions between the subcrustal Vrancea Slab, the overlying crust, and regional plate motions generate stress and deformation patterns that are far from trivial to interpret. Given these factors, we believe it is appropriate to describe the region and our study as complex and significant, and we prefer to retain this terminology to accurately reflect the scientific context and novelty of our work.