This paper explores various approaches to improve the prediction of ice supersaturated regions that is provided by numerical weather prediction models. This piece of work is crucial in informing ongoing efforts and trials in attempt to mitigate aviation’s contrail climate forcing. The paper is well-structured, employs suitable methodology, and presents reasonable results. It aligns with the scope of Atmospheric Chemistry and Physics. I recommend it for publication pending consideration of the following methodological questions and suggestions for improvements:

1. [Section 2.1] The description of the MOZAIC/IAGOS dataset might create confusion for the reader. For example, MOZAIC was transferred to the European infrastructure IAGOS in 2011, and the authors stated that only data between 2000 and 2009 was used. If this is the case, why not just name the in-situ measurements as "MOZAIC", instead of "MOZAIC/IAGOS"? In addition, what is the rationale for not including more in-situ measurements from the recent IAGOS dataset (i.e., between 2009 – present day)?
2. [Section 2.1] What is the temporal resolution (between flight waypoints) that is provided by the MOZAIC/IAGOS in-situ measurements? Have the authors considered and attempted to minimize the potential for autocorrelation between waypoints?
3. [Section 2.2] It is unclear what is the specific product of the ERA5 reanalysis data that was used in this manuscript. Is it the ERA5 high resolution realization (HRES) reanalysis? In addition, what is the spatiotemporal resolution of the ERA5 data that was downloaded and used in this analysis? How does the spatiotemporal resolution of the meteorological data influence the results presented in this paper?
4. [Section 2.2] Given that the spatiotemporal resolution of the meteorological data is an important factor that can influence the quality of ISSR and RHi estimates, have the authors considered re-running the analysis using “model level” data in addition to “pressure levels”?
5. [Section 2.2] It is worth including some details on the specific parameterization used to estimate the saturation pressure over ice (p_ice), which is required to estimate the RHi. Does the use of different parameterizations, i.e., Sonntag (1994) or Murphy & Koop (2005) lead to differences in the presented results?
6. [Section 2.2] A recent open-source repository found that the interpolation method across the vertical level (i.e., linear interpolation, log-log interpolation, or cubic spline interpolation) could lead to differences in the RHi estimates from the ERA5 humidity fields (https://py.contrails.org/notebooks/specific-humidity-interpolation.html). For example, given the non-linear lapse rate of the specific humidity, a linear interpolation across the vertical level may lead to overestimation of the specific humidity. Therefore, the authors should consider exploring the impact of the interpolation methodology on their presented results.
7. [General comment] The authors correctly highlighted that contrail mitigation is most effective when long-lived persistent (warming) contrails are avoided. Since contrails forming near the RHi threshold (close to 100%) tend to be shorter lived relative to those formed at higher RHi’s, one potential strategy is to focus on regions where ice supersaturations are notably higher than threshold conditions. For instance, could there be an improvement in the ETS scores if the authors consider raising the threshold for the “predicted contrail formation” to be, say, RHi_ERA5 > 110%, instead of 100%?
8. [Structure] Section 3 of the manuscript includes both the methodology and results, thereby making the section unnecessarily long. The authors should consider separating the methodology and results into two different sections to improve the readability of the manuscript.

This paper examines the key challenge of accurately predicting ice supersaturated regions (ISSRs) at cruising altitudes using reanalysis data and in situ measurements as ground truth. This paper is a continuation and culmination of some of the authors' efforts and previous research to use dynamical proxies to improve the ability to predict relative humidity with respect to ice. It aims to provide an end-to-end statistical framework for better estimation of ISSR formation, focusing on the different statistical design options and their respective performance. The conclusion, although disappointing, is a significant result for the contrail research community.

The paper is well organised, uses appropriate statistical methods and provides plausible and important results. We recommend its publication as it is, with only a minor technical correction:

[Section 3.1] "This means that the likelihood ratio must be greater than 3.85". I assume here that 3.85=0.5/0.13. But if I'm not mistaken, to make ISS more likely than no ISS, we need the posterior odds ratio to be greater than 1.0, which translates into a likelihood ratio threshold of 1.0/0.13=7.69. Since ln(7.69)=2.0, this would be consistent with the threshold on the logit

Here are also some points to consider for further scientific discussion or future research, particularly with regard to the statistical methodology (Part 3):

1. [Section 2.1] As mentioned by the other reviewer, it would be interesting to understand how the spatio-temporal resolution of the meteorological data affects the statistical results developed in the Part 3.
2. [Section 3.1] It could be interesting to use the same validation framework (split train/test datasets) to estimate the ETS using the individual Bayesian estimators directly as binary classifiers. Section 3.1 could be developed in a very similar way to 3.2, as methods based on Kullback-Leibler distance can be used as model/variable selectors (although less frequently than mutual information). Distance between variables is not needed if you are using forward additive selection or backward elimination.
3. [Section 3.1] It might be interesting to use slightly more complex Bayesian estimators such as naive Bayes, which tend to work surprisingly well even in a situation where the variables are correlated, provided it is improved by optimising the threshold applied to the posterior output. Accounting for variable correlation with a simplified DAG and simple conditional probability parameterisations is also an option.
4. [Section 3.2.2] The models chosen in this section could be automatically determined by forward additive selection or backward elimination using the mutual information as variable selector.
5. [Section 3.2.4] Ridge or lasso regularisation could be used to better control the risk of overfitting due to correlated variables.
6. [Section 3.2.4] Although it is highly unlikely that it would have changed the results at all, modern tree-based machine learning techniques (random forests/LightGBM) could be used as they are designed to deal with non-linearities, overfitting/correlated variable problems and implicit variable selection at the same time.
7. [Section 3.2.4] Optimising the output probability threshold might improve the results slightly (but certainly not significantly).
8. [Section 4.2.3] As mentioned by the authors, it is not surprising that adjustment techniques like [Teoh et al., 2022] are not really needed here, as the non-linear nature of the classifier has learned to infer the implicit bias non-parametrically.
9. [General] As mentioned by the other reviewer, it might be interesting to reformulate the problem as a multi-class problem, either using situations of even higher supersaturation or using the ontology developed in [Wilhelm et al., 2020] (no persistent contrails, persistent contrails, most warming contrails). A direct regression framework between reanalysis and in-situ relative humidity is also an option.