We thank Dr. Ferré very much for his appreciation of our work! We agree that a more detailed presentation and discussion of the results of the Bayesian analysis would be helpful.

We suggest adding few figures in order to (1) show marginal posterior distributions of SHPs (for one exemplary site) to further discuss how parameter uncertainty was constrained by the observations; and (2) visualize the results of the uncertainty propagation for cumulative recharge and magnitudes of recharge peaks at all sites. In the revised manuscript, we will discuss the propagated parameter uncertainty in context with parameter sensitivities which we have established in Schübl et al., 2022 https://doi.org/10.1016/j.jhydrol.2022.128429 with similar HYDRUS models and site conditions. Results in the study here showed, that the prediction of cumulative recharge was affected by less uncertainty than predicting magnitude in recharge peaks (related to the higher uncertainties in the soil hydraulic conductivity parameter Ks, and lower uncertainties in VGM shape parameter n).

To find out how site conditions affect uncertainties, we performed a systematic analysis in Schübl et al. (2022) (see above). In the study here, we obtained a similar result in that uncertainties in SHPs were higher in coarser soils. Our results here also suggest, that parameter uncertainties were higher when the estimated measurement error was high relative to the temporal variability of the measurements. We would like to expand the discussion on this.

The second part of the paper was aimed at characterizing the Austrian sites based on measured and modeled hydrological variables and site properties. Since uncertainties in long-term estimates used for characterization were generally low, and since the methods we used for that purpose required unique values, we did not include the uncertainty ranges in this part of the analysis. We agree, that in the revisions, we should state more explicitly where and why we proceeded using the ML estimates.

Finally, we can draw more conclusions from discussing uncertainties in SHPs and long-term recharge estimates resulting from the Bayesian analysis: It has often been found that in-situ field measurements of soil water content are not sufficient for the accurate and precise estimation of SHPs (Jacques et al., 2002 https://doi.org/10.1016/S0022-1694(01)00591-1; Ritter et al., 2003 https://doi.org/10.1016/S0378-3774(02)00160-9; Scharnagl et al., 2011 https://doi.org/10.5194/hess-15-3043-2011). However, in this study we have found that while SHPs were partially affected by considerable uncertainties (especially soil hydraulic conductivity and residual water content parameters), the propagated uncertainty in groundwater recharge estimation was low. This was related to the fact that the VGM shape parameters, especially n, were associated with small uncertainties. At locations with very wet/humid climate conditions, the uncertainty was small in proportion to estimated recharge sums. At locations in a dry climate on the other hand, uncertainties were considerable in relation to absolute rates. There especially, uncertainties in SHP estimation should be reduced, for example by improving the data quality and resolution or by including other observation types, such as soil pressure heads.

We think that by adding these points to the manuscript we will give more insights and get more value out of the Bayesian analysis.

We are very happy about the positive and motivating feedback on our work from Reviewer#2! We thank him/her for reviewing our manuscript.

We thank Reviewer#3 for his/her appreciation of our paper and the helpful feedback. We address the main comments and questions here and provide answers to the detailed comments in the attached PDF file.

The Reviewer raises an important point of discussion to better clarify our approach.

The Bayesian inference can be applied directly to obtain Soil Hydraulic Parameters (SHPs), if data include pressure head, water content, and conductivity (e.g., laboratory measurements derived from the simplified evaporative method). Instead, our study uses field scale observations of volumetric water content to inversely estimate the most probable distribution of SHPs that generated observations. Therefore, the assumption of homoscedastic and normal errors (reflected in the likelihood function) refers to TDR measurements, not to SHPs. The posterior distribution of SHPs is inferred by combining two components: 1) the prior distribution, which reflects the modeler’s believe about Soil Hydraulic Parameters (SHPs) before considering measurements (in our case, soil water content), 2) the likelihood, which describes the probability that a parameter set drawn from the prior has generated the observations. By combining the likelihood and the prior with Nested sampling  (or Markov-chain Monte Carlo) and HYDRUS, we obtain a posterior distribution of the most probable SHP values, which reflects the parameters’ uncertainty:

1. Prior: As it is frequently the case in vadose zone hydrology (e.g., Brunetti et al. 2020 https://doi.org/10.1016/j.jhydrol.2020.124681, Wöhling et al. 2015 https://doi.org/10.1002/2014WR016292), we assumed bounded uniform priors to avoid making important assumptions on the shape of the posterior, and let the data tell us what is its shape. However, at the same time, we imposed hard boundaries on the parameters to avoid the possibility to obtain physically unrealistic values.
2. Likelihood: We assumed that sensor errors are normally distributed. This is a widely established approach in inverse vadose zone modeling (e.g., Schelle et al. 2012 https://doi:10.2136/vzj2011.0169, Gao et al. 2019 https://doi.org/10.2136/vzj2019.03.0029).

We also thank the Reviewer for pointing out the problem of equifinality and non-uniquess, which gives us the opportunity to better clarify the utility of the Bayesian approach to address this aspect, and to discuss the limitations of the sensitivity analysis compared to the Bayesian probabilistic approach.

The Bayesian inference is precisely conceived to have a statistical rigorous appraisal of the “equifinality and non-uniqueness”. The Bayesian approach infers a marginal posterior distribution that exposes the parameters’ uncertainty, and their interaction (e.g., correlation). If the resulting uncertainty is high (i.e., wide posterior), then data are not informative for that parameter. In this case, the modeler has two choices: 1) ask for other measurements (e.g., disk infiltrometer for Ks) to have more informative priors and run again the Bayesian analysis; 2) honestly communicate what is the parameters’ uncertainty with the data available, and more importantly, how the estimated uncertainty propagates in the model simulations. This is what we precisely did in our study.

The sensitivity analysis is certainly a valuable tool, which we use frequently in our research. However, it will not provide any more meaningful information compared to the Bayesian analysis for this work. A global sensitivity analysis will sample the parameters’ space (frequently ineffectively as MCMC techniques are much better in finding high-probability regions), and then decompose the variance to identify influential and uninfluential factors. But this is already better targeted in a Bayesian analysis: influential parameters are those that exhibit leptokurtic posteriors, while uninfluential factors are those that have similar prior-posterior distributions (in our case flat). The sensitivity analysis might have some utility in high dimensions for numerical sampling reasons, but this is not the case and is beyond the purpose of the present study.

In general, we agree that it is helpful to include pressure head data, as it can help to identify SHPs with even less uncertainty than with soil water content data. However, at the sites in this study we had some issues with soil pressure head measurements: (1) they were not available for all sites which would have impaired the comparability of results between locations,  (2) they were composite from different measurement techniques (tensiometers and gypsum blocks) and included sudden shifts and large gaps in time series. Altogether, we found the measurements to not be reliable enough to be used in this study. We agree that results for SHP estimates might change when including the available soil pressure head data, therefore we will include this in the discussion. We also agree that this would be a very interesting topic for further studies with improved field measurements of soil pressure heads.

We used the PCA Biplot in this study to visualize the two clusters of hydrologically similar sites in context with the variables according to which they have been characterized. This kind of visualization, with individuals (sites or samples) and variables in the same plot, has been used in several studies, for example by  Rodríguez et al. (2020) https://doi.org/10.1007/s10750-020-04201-5 (Figure 3), Luna et al. (2018) https://doi.org/https://doi.org/10.1002/eco.1896 (Figure 7), Gibson et al. (2019) https://doi.org/https://doi.org/10.1016/j.ejrh.2019.100643 (Figure 7).

R tutorials showing this kind of plot are also available, e.g.,

[1]https://f0nzie.github.io/machine_learning_compilation/detailed-study-of-principal-component-analysis.html (See 4.20 Biplot)

[2] https://www.datacamp.com/tutorial/pca-analysis-r

[3] https://finnstats.com/index.php/2021/05/07/pca/

Citing from [1] for the use of Biplots with variables and individuals, the focus is “…on the direction of variables but not on their absolute positions on the plot. Roughly speaking, a biplot can be interpreted as follows: an individual that is on the same side of a given variable has a high value for this variable; an individual that is on the opposite side of a given variable has a low value for this variable.” We wrote the code for our plot in Python using the sklearn module (Pedregosa et al., 2011 https://arxiv.org/abs/1201.0490) which we should cite properly (the module is cited in the manuscript for the site clustering but not the PCA plot). We agree with the comment in the manuscript that the data is not well represented in the Biplot of Figure 4(b) and it does not add further insight, we therefore suggest deleting it.