Bertoncini and Pomeroy quantify uncertainty in precipitation estimates using a network of in-situ precipitation gauges in the triple continental divide area of the Canadian Rockies, a region where precipitation can vary immensely across elevation bands and differently, depending on storm systems. Using the WMO guidelines for station density in mountainous areas, the authors transform and back-transform precipitation data (for normalized distribution), quantify a cumulative distribution function, and use a kriging and lapse rate approach to calculate and track precipitation standard deviation and coefficient of variation across space and time. In this way, the authors are able to determine areas where precipitation estimates are more and less uncertain and the areas where added in-situ observations would be most valuable in the future (e.g., relatively higher elevations). Much of the analysis, and thus manuscript text, includes a very clear description of the methodology used. I include no major changes to the workflow and thank the authors for their thorough depiction of the work in text and figures and for making the important and relevant connection to downstream hydrology. Limitations to the work, including the select reanalysis product, could be further discussed. As could a connection to other spatially distributed precipitation products. Otherwise, the small number of minor suggested changes I have made are with respect to clarifying language around some of the statistics (e.g., when uncertainty “rose” vs. “fell”) and domain description. The attached line-by-line comments should provide more clarity with these items, with the goal of better emphasizing the importance and value of this work, which I envision will serve as a frequent reference for many future projects.

Bertoncini and Pomeroy present and interesting approach for quantifying the uncertainty associated with sampling the spatiotemporal and elevational variability of precipitation in complex terrain, in addition to offering to the community a unique dataset for the Canadian Rockies.

General comments:

The paper is well written. It is relatively easy to understand the methods and follow the discussion, and the figures are well done. The literature review is interesting, but not exhaustive on the important topic of precipitation lapse-rate estimation methods. In particular, Dura et al. (2024, EGUsphere) proposes a few more relevant references. The dataset obtained is unique and helps to shed light on the important issue of network design in complex terrain. However, I must say that I am not totally convinced by some aspects of the methodology, partly because the authors do not evaluate the accuracy of the method through objective verification scores. This is important because the interpolation method proposed in the paper is new, even though it borrows from ideas already published in the literature to deal with elevational variability in precipitation intensity. In addition to the lack of verification, I also question some of the underlying hypotheses of the method itself. Overall, I believe that this paper has merit and should be published, but that important questions need to be answered first, and these will likely require additional experiments, as well as possible adjustments to the methodology.

In the next section, I propose specific comments for the introduction and methodology sections. I do not have much to say about the discussion part of the paper, although it might need to be revisited if changes are made to the methodology.

Concerning objective verification: verifying a gridded product that incorporates all available observations is obviously a challenge. I suggest that several test cases be chosen that correspond to various type of precipitation events. For these events, the interpolation procedure could be carried out using a subset of the available data, and the precipitation maps (with and without considering the altitudinal gradient in precipitation) could be compared to each other and to the precipitation observed at the remaining stations. Care should be taken to create a subset of data for verification that covers the distribution of elevations in the domain of interest. In addition to comparing kriging estimates to the observed precipitation, it is also necessary to compare the values of SD obtained from the kriging method to the distribution of the errors made when making predictions at unknown locations. In the NST transformed space, I would expect about two third of the estimations to fall within plus or minus one standard deviation from the kriging estimate. This should be checked. It is particularly important in this study because the discussion does not only focus on the estimate of precipitation, but also on its associated uncertainty.

Specific comments:

Line 60: "on the quantity variance". What do you mean by "quantity variance"? Not clear. Rephrase.

Line 62: "...semivariogram, which is the relationship between the variance in the observed quantity with the measured distances." This is not true. The semivariogram is defined as 1/2 E [ (Xr - Xs)^2 ] for random variables Xr and Xs at locations r and s. In its simplest expression, the semivariogram is only a function of the distance d(r,s), but this is not the only option. The variance of Xr and Xs need not even exist for the semivariogram to exist, and in all cases the semivariogram does not measure the variance of the observed quantities Xr and Xs (well, actually it does if the distance is larger than the range of the variogram). You could write something like "semivariogram, which is the relationship between the second moment of the differences between the observed quantity at two locations and the distance between these two locations".

Lines 64-69: physically-based residual kriging, a.k.a optimal interpolation, should also be mentioned here. See for example Pelli et al. (2022, SERRA). In the Canadian context, Brasnett (1999, JAMC) used this approach to combine a background field from a numerical model with observations, taking into account elevation differences explicitly through a variogram, to interpolate snow depth observations.

Lines 81-82: be more specific here with respect to how meteorological models and hydrological models take into account lapse rates. These two categories of models do things very differently. Discuss briefly advantages/disadvantages of how it’s currently done in these models. Add references.

Table 1: add references for the various instruments, either in the table or in the text.

Lines 159-160: I am not convinced that the use of ERA5-Land wind speed is acceptable for correcting precipitation accumulations for wind-induced undercatch. Have you compared ERA5-Land wind speed predictions with observations in this region? ERA5-Land has a resolution of 9 km but relies on an atmospheric model which has a resolution of 30 km. The orography field used in numerical models is generally smoother than the nominal resolution, so the effective resolution of ERA5 is even lower than 30 km. Vanella et al. (2022, J.Hydrol.) have compared both ERA5 and ERA5-land wind speed predictions to observations in different climates and topography in Italy and concluded that both products strongly underestimate wind speed (by 28% to 42% depending on the region). Is it the same in your region? How does this affect your precipitation estimates?

Lines 162-164: I find it surprising that you decided to use the BC and COOP data as is because no existing equations exist to do the bias correction. I assume that the BC gauges are impacted by wind? I read that these gauges use a pressure transducer to estimate precipitation. They might also be impacted by changes in atmospheric pressure. Wrt ruler-based snowfall, how is the density of fresh snow estimated in order to obtain a water equivalent? Does it / should it take into account wind speed? If you cannot correct the data, it might be a better idea to reject the data when the wind speed is above a threshold such that biased estimates are expected.

Line 193: the cubic-root transform is also used, see for example Lespinas et al. (2015) as well the more general Box-Cox transform, see van Hyfte et al. (2023, Tellus A).

Line 205: how is the semivariogram model chosen? With what frequency is each model chosen?

Lines 207-208: you mention that grid longitude and latitude is used to interpolate. When computing distances, does the method take into account the distortion caused by the fact that degrees of longitude are further apart in the south than in the north of the domain? Given the shape of your domain, this is important.

Line 208: back-transforming the data introduces a bias in the interpolation, because the kriging method is unbiased in the transformed space, but not in the original space. See van Hyfte et al. (2023) for details and formulas for taking this into account in the context of the Box-Cox transformation. Can you comment on the magnitude of that bias? Note that there is also a bias in the back-transformed standard deviation.

Line 216: you mention using 53 gauge pairs to estimate the lapse rate. I assume that you are recomputing the lapse rate on a daily basis? If so, does that number vary over time or are these 53 gauge pairs available for the whole time period covered by the product? Are these gauge pairs distributed relatively evenly over the domain? If not, do you think that this is problematic? Please provide a map of the location of these gauge pairs. Perhaps this information can be added to fig 1.

Line 226-227: Because you have interpolated the raw data to obtain the final daily horizontal precipitation field, you now need to rely on a reference elevation field interpolated from gauge elevations. I assume that you are recomputing this reference elevation field every day, since the network changes over time. Can you please confirm? Note that an alternative approach would be to bring all station observations to a reference elevation using equation (3) for each station, e.g. 1500m and generate a precipitation field valid at that altitude. Then, a daily lapsed precipitation could be generated on the SRTM elevation model by applying again equation (3) for each SRTM grid cell using a reference elevation of 1500m for each grid point. It would be interesting to compare the two approaches.

Lines 233-234: this sentence does not read well. I do not fully understand what is meant here. In particular, the use of “Therefore” at the start of the paragraph suggests that what is being said follows logically from the end of the previous paragraph. However, the link is not obvious. Please rephrase. Perhaps there is a sentence missing just before this one?

Lines 236-238: I believe that the term “coefficient of variation” and its abbreviation “CV” is an abuse of terminology. Indeed, if I understand well, the authors are computing the average of the SD field (over a year) and dividing by the average of the precipitation. However, the average daily value of SD does not correspond to the uncertainty associated with the mean annual precipitation. For example, if errors in daily precipitation estimates were independent, the standard deviation of the error on the annual average would correspond to the square root of the sum of the square of the SD values, divided by the square root of 365 days. Hence, if SD values were relatively constant over a year (they are obviously not), the standard deviation of the error on the annual average would be reduced by a factor of about twenty. The coefficient of variation of this error would be reduced by the same factor. Whatever the covariance structure of errors in daily precipitation estimates, the coefficient of variation would be obtained by manipulating the error variances (the square of SD), and not the standard deviation of the errors. However, the covariance between errors in daily precipitation estimates has not been modelled in this study – it is thus unknown. For this reason, it is an abuse of terminology to talk about the term “coefficient of variation”. What can be done? Well, the coefficient of variation can be computed for each day and each grid point, by dividing SD by the precipitation amount. Obviously, this is problematic for days without precipitation. It would still be possible to compute quantiles of the daily CV values and map one of them, for example a quantile halfway between the frequency of zeros and one. Another solution is simply to change the nomenclature and call the quantity reported in this paper something else than CV. This is what I suggest.