I agree with the authors when they state that "translating the definition of PMP into a statistical model is interesting" (line 350).  They could also say "estimating PMP is a really hard problem".  To the authors' credit, they begin with the commonly accepted definition of PMP as an upper bound, and then construct a statistical model which fits this definition.  Fitting a model with a finite upper bound is challenging because precipitation data usually suggests that the distribution is unbounded, and further has a heavy tail.  It is probably not surprising that the authors ultimately find their approach to be unsuitable for implementing in practice, and conclude that the best statisical approach is to eschew the upper bound requirement and instead implement extreme value (EV) methods.  

Unfortunately, I think the manuscript's structure does not tell its story well.  Primarily I view the paper as an interesting way to discuss the challenges of PMP estimation, and talking about their particular model is one part of this larger story.  It strikes me that the take away message does not appear in the abstract or in the body until Section 6.  I think it would be better to move these messages up front.  The story I imagine is something like this:

1.  Statistically estimating PMP is hard because its definition assumes a bounded tail, but precipitation data suggests the tail is unbounded.  Because statistical estimation is hard, other methods like moisture maximization and Herschfield's scaling get used.  Uncertainty and climate change are hard to incorporate into these non-statistical methods and frequently-used moisture maximization approaches involve several subjective judgements.

2.  Starting with the ideas which underlie moisture maximization, we develop a sensible statistical model which assumes an upper bound.

3.  We perform simulation studies and use the method to fit PMP at two locations in Quebec, but find that estimates for the upper bound have unsuitable uncertainty.

4.  We conclude with a discussion and offer our suggestion for best practices.

I think all the pieces of this story are in the paper, but I do not think the current focus of the paper gets the essential message across very well.

I find the notation in the paper to be inconsistent.  In equation (1), Y_i denotes precipitation of storm i, but in equation (2) I believe Y_i has been replaced by P_i.  Equation (4) supposedly comes from Eq (1), but has quantities EP_i and EP_max, which are presumably PW_i and PW_max in equation (1)?  The ratio EP_i/EP_max is known/assumed to be less than 1, correct?  If so, please say this explicitly.  I believe equation (6) is used as the basis for the statistical model:  Y_i = EP_i/EP_max * r_i * PMP.  If I am following correctly, Y_i is random and observed. I think EP_i and r_i which underlie Y_i are random, but unobserved.  EP_max is a parameter but not known, and PMP is the parameter we wish to estimate.  So in the end, the authors propose a model for the observed precipitation Y_i, but use moisture maximization logic to include PMP as a parameter.  They choose a beta/Pearson 1 as their distribution to fit.  A cynical comment could be "the authors use a data-independent argument to conclude the data arise from a distribution, but which fits the data poorly."  I think the story to be told here is that if one begins with a supposition of an upper tail, and one tries to then fit a model based on that assumption, things are really hard.  

If I understand correctly, the authors propose a beta/Pearson 1 distribution and fit *all* of the nonzero rainfall data to it.  There is talk of thresholding on page 8, but it seems to be more tied to the discrete nature of the measurements rather than to thresholding for focusing on extremes.  An EV approach would pick a high threshold or take block maxima and fit an EV model, presumably a reverse-Weibull guaranteeing an upper bound.  Would such a method be better suited for estimating an upper bound than fitting a beta to the entire distribution?

The authors show qq plots for the EV models in Figure 6.  QQ plots for the beta fit are noticeably absent.  

l173.  Why is beta known to be greater than 1?

l304:  Figure??

The paper under consideration presents a statistical approach using the Pearson Type I distribution to estimate the upper bound of historical rainfall, incorporating uncertainty bounds. The authors aim to quantify uncertainty and address subjectivity inherent in the various stages of the World Meteorological Organization (WMO)-recommended Probable Maximum Precipitation (PMP) estimation methods. However, the WMO-recommended moisture maximization method focuses on maximizing highly efficient storms based on the physical mechanisms of those storms. From a statistical perspective, precipitation’s tendency to exhibit a heavy-tailed distribution poses challenges in defining an upper bound, and this study similarly encounters this issue. Here are my major concerns:

• Previous studies demonstrated that precipitation naturally often exhibits a heavy-tailed distribution (shape parameter greater than 0) that brings rare storms over global scale, the proposed method struggled to estimate upper bound of those places by majority of parameter estimation methods. In this study, the method of moments partially quantifies the range with limited data, but the range is unrealistically large, for example, the range estimated for St-Hubert station varies between 165 to 9006.
• A simulation study is conducted with distribution assumed convex density and found more than 40,000 sample size (In arid/semi-arid region that equivalent to more than 1000 years wet days) is required to stabilize the estimate. Given this, it is surprising that the authors did not attempt to expand the sample size for the two stations by incorporating numerical model ensemble precipitation products. Doing so could have supported their findings.
• This study compares their estimated upper limit with moisture maximization based PMP value. The PMP values using moisture maximization were found to be 282 mm for Montréal and 436 mm for St-Hubert, whereas the observed 24-hour maximum precipitation for these stations was 81.9 mm and 106.5 mm, respectively. Thus, the maximization ratio will be 3.44 and 4.09. The reason behind the exceptionally high maximization ratio may be due to the selection of storms and/or estimation of storm associated precipitable water. This study includes low magnitude storms (0.9 quantile might give more than 500 samples but previous studies mostly consider the highest 50 or less storms) and did not separate those storms that could lead to higher maximization ratio. Previous studies mostly limit the maximization ratio 2.0 (that only for orographic storms). Imposing a similar limit could provide some physically possible value around 200mm that aligns with 10,000-year return level (POT based) and PMP value would not much different within 26 km distance. Since the moisture maximization method provides an unrealistically high PMP value, comparing with this value to validate the method is questionable. It is recommended to use multiple study sites and consider those sites where maximization ratio lies below 2.0 and compare within those sites.
• National Academies of Sciences, Engineering, and Medicine (2024) recommends for risk-informed extreme value analysis methods that account for low exceedance probabilities and provide robust uncertainty and nonstationary quantification. It remains unclear how the proposed method offers advantages over or resolves issues better than these recommended approaches.
• The choice of Pearson Type I distribution over other distributions is missing.

Minor comment:

• Line 304: The placeholder “Figure??” needs Figure number. Additionally, there is inconsistent notation in Equation 4 compared to Equation 1. The term "EP" and “EP_max” should be clearly defined to maintain consistency and avoid confusion.

Interesting paper. However, I feel it may become more of use with additional investigations of refinements needed to the method that has been proposed. See my comments below.

l125: Is EPmax defined for a calendar day for nearby regions or defined as the maximum at that point location? Equation 1 referred to PWmax which is the atmospheric moisture. EPmax is I suspect the same but should be consistent or else defined. If non-seasonal, please refer to WMO guidelines for seasonal variations in PMP.

l140 - interesting. However, the sampled ri are non-iid, which complicates their use in defining the Beta distribution I think, Plus, there is an assumption that the sampled ri has an upper limit of 1. Given this limiting value will dictate/influence the PMP estimate, the uncertainty associated with this assumption is important to characterise.

Also, how can stationarity be assumed given there is a clear temporal trend in precipitable water time series. Would violation of the stationarity assumption distort the beta distribution parameters?

l375 - The authors recommendation makes sense. In addition to the issues they have mentioned, I also feel that the lack of independence (unless they parameters are being fitted using iid data above a threshold) and the presence of a trend are limiting factors. I wonder if using a nonstationary model and regional data can help overcome these limitations.