We thank Reviewer 1 for their thorough and positive feedback (e.g., "Overall, this is one of the best scientific papers I have ever read, especially considering its highly technical nature”). We appreciate their encouraging comments on the manuscript’s presentation, including the clarity and conciseness of the writing, the precision of the mathematical formulae and the efficacy of the figures (e.g., “All aspects of the paper appear to be very carefully prepared, including citation of relevant literature, readability of the text by a broad audience, presentation of methods and results, and technical precision in presentation of the mathematics and statistics”). We are also grateful they recognised our efforts towards contextualising this research within the existing literature and giving a balanced discussion about the advantages and disadvantages of the proposed method.

We thank Reviewer 2 for their insightful comments and suggestions. We address the reviewer's concerns point by point below.

Comment 1

“First, ECMWF also provides high-resolution reanalysis precipitation data (ERA5 Land, hourly, 0.1 degree, ~ 9 km), which is not considered in the manuscript. How does the generated MFGP precipitation estimates compare with ERA5 Land precipitation data?”

ERA5-Land is a reanalysis dataset that provides a consistent view of the evolution of land variables at an enhanced spatial resolution of 0.1° x 0.1° (9 km) compared to ERA5’s resolution of 0.25° x 0.25° (31 km). This is produced by running a land surface model to regenerate some of the land components of the ERA5 climate reanalysis. For atmospheric forcing, it uses ERA5 atmospheric variables such as air temperature and precipitation at a 0.1° resolution by linearly interpolating the driving variables to the ERA5-Land grid. Although other forcing variables are corrected, this is not the case for precipitation. For further details please see Muñoz-Sabater et al. (2021) and ECMWF (2024). Precipitation characteristics from ERA5-Land are therefore very similar to ERA5 (Gomis-Cebolla et al., 2023; Xu et al., 2022; Xin et al., 2022). They also should theoretically perform worse than the linear regression models presented in our paper (Table 1-3) which also include elevation as a predictor.

Moreover, we would like to stress that we have carefully chosen four gridded precipitation datasets for the High Mountain Asia region to evaluate our model against (based on existing literature). These are:

• APHRODITE, a gridded rain-gauge interpolated dataset for Asia considered the gold standard for precipitation in High Mountain Asia,
• CRU-TS4, a global gridded rain-gauge interpolated dataset. 
• A bias-corrected high-resolution regional climate model simulation, which used the Weather Research and Forecast (WRF) model at a spatial resolution of 5 km, with precipitation output corrected using local rain-gauge data for the region investigated in this manuscript.
• TRMM, a satellite-based precipitation dataset designed to improve our understanding of precipitation in the current climate.

To address the reviewer’s concern, we will make the connection between ERA5-Land and our linear regression models clearer in Section 4.

Gomis-Cebolla, J., et al. (2023), Evaluation of ERA5 and ERA5-Land reanalysis precipitation datasets over Spain (1951–2020), Atmos. Res., 284 https://doi.org/10.1016/j.atmosres.2023.106606

Muñoz-Sabater, J., et al. (2021), ERA5-Land: A state-of-the-art global reanalysis dataset for land applications, Earth Syst.Sci. Data,13, 4349–4383, https://doi.org/10.5194/essd-13-4349-2021

ECMWF (2024). ERA5-Land: data documentation. Available at: https://confluence.ecmwf.int/display/CKB/ERA5-Land%3A+data+documentation 

Xu, J., et al. (2022), Do ERA5 and ERA5-land precipitation estimates outperform satellite-based precipitation products? A comprehensive comparison between state-of-the-art model-based and satellite-based precipitation products over mainland China, J. Hydrol., 605, 127353, https://doi.org/10.1016/j.jhydrol.2021.127353.

Xin, Y., et al. (2022), Evaluation of IMERG and ERA5 precipitation products over the Mongolian Plateau, Sci. Rep., 12, 21776, https://doi.org/10.1038/s41598-022-26047-8.

Comment 2

“Second, the authors only consider a very simple machine learning model, i.e., linear regression, and complex deep learning models that require a lot of training data, including Convolutional Conditional Neural Processes (ConvCNP) and Convolutional Gaussian Neural Processes (ConvGNP). They neglect simple machine learning methods that do not need many data, such as random forest and support vector machines.”

The main goal of the paper was to provide calibrated uncertainty distributions for precipitation in this area. This information allows hydrologists to quantify the probabilities of extreme events and policymakers to make better decisions with limited resources as highlighted in Section 1. Although there are many possible machine learning benchmarks thatwork on smaller datasets, Random Forests and Support Vector Regression are not inherently probabilistic models.

However, we appreciate that the case for the performance of Gaussian processes could be better contextualised by including these models. To address the reviewer’s concern we implement Random Forests and Support Vector Regression for the validation experiments which we present in a new appendix. These new results are shown below in Tables A-C. The Random Forests perform similarly to the GPs doing better in the Beas and Sutlej Basins but worse over Europe. The Support Vector Regression models perform similarly to the linear regression models. Nevertheless, only the Gaussian Processes output calibrated uncertainty estimates.

Table A: Comparison of model performance metrics trained on ERA5 data for the 'data rich' setup over Europe. We include a linear regression model, a Random Forest, a Support Vector Regressor with a Radial Basis Function (RBF) kernel, a GP using a RBF kernel with Automatic Relevance Determination (ARD), a GP using a Matérn 5/2 kernel with ARD. The MFGP model also trained on the gauge data is also shown for reference. The metrics include the average RMSE, the 5^th percentile RMSE (RMSE5), the 95^th percentile RMSE (RMSE95), and the R² score. The bolded values represent the best scores amongst the model benchmarks (not including the MFGP model).

Table B: As Table A, but for the 'data sparse' setup over Europe.

Table C: As Table B, but for the Upper Beas and Sutlej Basins.

Comment 3

“Finally, in the model comparison shown in Table 1-3, GP is only trained on ERA5 data at station locations. Why not use all ERA5 data in the study region to train and test the model? I would expect better model performance even for GPs.”

The model has information about ERA5 at all the training and test locations for the validation experiments. These locations fall within the ERA5 grid boxes, so there is theoretically little to no additional information to be gained by including neighbouring grid box values. To check this, we ran an experiment that was trained on all ERA5 data for the Beas and Sutlej Basins.

Results from this experiment are shown in Table D (below) and confirm that there is no added benefit in including this data. However, we have clarified this in the revised manuscript by adding additional text in Section 3, which states ‘The experiments were also conducted with all the ERA5 data for the study area (not shown), but showed no significant improvement over using the ERA5 data at the station locations only’. 

Note that we did not rerun these experiments over Europe as we would have needed to apply methodological approximations to overcome the memory and computational bottlenecks that come with this larger domain.

Table D: Comparison of MFGP performance using ERA5 for the whole study area (all ERA5) and using only ERA5 at the training and test site locations (limited ERA5) over Upper Beas and Sutlej Basins. The metrics include the average RMSE, the 5^th percentile RMSE (RMSE5), the 95th percentile RMSE (RMSE95), the R² score and the mean log loss (MLL). The bolded values highlight the best scores.

Comment 4

“I believe the authors use the Nash-Sutcliffe efficiency (NSE) in the manuscript, rather than R2, which should always be non-negative values.”

We confirm that we are using the coefficient of determination or R² score. This metric, defined and explained in Appendix B, and is given by:

R² = 1 - SS_res/ SS_tot = 1 - ∑_i(y_i - f_i)² / ∑_i(y_i - y_avg)²

where f_i is the i^th predicted value,  y_i is the i^th observed value, and y_avg is the mean of the observations. SS_res is therefore the sum of the squared residuals and SS_tot is the total sum of squares. A negative R² is possible and would indicate that the model is predicting worse than the precipitation mean. Although negative R² scores are unlikely in interpolation settings, they are possible when making predictions outside of the training distribution. To address the reviewer's concerns, we make this interpretation of negative R² clearer in the main body of the paper and Appendix B.