This research article employs a probabilistic approach to test two hypotheses regarding the geological evolution of caves in an alpine karst catchment, thus in my opinion, the study fits well within the scope of the HESS journal. In general, the article is well written, and its objective is clear. The methodology consists of techniques already proposed and tested in previous studies, however, the application of such techniques in the context of geological evolution assessment is quite novel. In my opinion, the paper can be accepted for publication provided that the present comments/suggestions are correctly addressed:

• In the document the authors state that the simulation of the networks for a resolution of 50x50 meters can be performed in less than 2 minutes. In this sense, what is the limiting factor for not using a map with a higher resolution?
• In Figure 3, how is the expected network defined?
• In Figure 3, is it possible to draw the expected network with different plotting settings to enhance its visibility? As it is, the expected network can be confused with the simulated conduits.
• Could you provide an interpretation of the variability of the simulated conduits? Is it an indicator of the uncertainty of the expected conduit? Could you include a metric to measure such variability?
• In line 243 the authors mention that there is an additional scenario that is not explored in this study. Is there any justification for not analyzing it?
• In line 125 the authors wrote that in this study you are considering that conduits form preferentially in the direction of the maximum downward hydraulic gradient. Could you please add some references to studies where this assumption is also employed? Alternatively, could you please add some words explaining the reasoning behind this assumption?
• In Section 6 (Findings) you do some references to Figure 4. However, they seem to refer to elements in Figure 5.
• There is a missing closing bracket in line 112.

NB: I have received an email notification, but I did not read the referee1’s comments to avoid any influence in my judgment. Some repetitions could thus appear in the following (or I could have a completely different opinion, I don’t know).

The present work is devoted to the study of two hypotheses on the Hölloch cave karstogenesis process with a geometric modelling approach. The Hölloch cave is part of the large Gottesacker karst system and contains 2 parts from a karstogenesis point of view: an active one and an inactive one. The authors recall two explanations proposed by speleologists, which are not exhaustive and/or excluding, for the co-existence of the 2 parts. Then, they translate these hypotheses into input constraints (i.e., inlet/outlet pairings) for a stochastic modelling method of karstic conduit paths. Then, observing the results, they comment on both hypotheses.

The topic addressed in this paper is, I think, completely relevant for HESS. The English is good and the reading easy. However, I have unfortunately two major concerns that prevent me to recommend its publication: 1) the lightness of the contribution as compared to the recent papers already published by the same authors on the same topic and the same case study; 2) what I would qualify as an over-exploitation (or even mis-exploitation (?)) of the modelling approach. More details are further provided.

The submitted paper is a follow-up of two papers recently published in Hydrogeology journal in 2021 and 2022: Fandel, C., Ferré, T., Chen, Z., Renard, P., Goldscheider, N., 2020. A model ensemble generator to explore structural uncertainty in karst systems with unmapped conduits. Hydrogeol. J. 229–248. https://doi.org/10.1007/s10040-020-02227-6 and Fandel, C., Miville, F., Ferré, T., Goldscheider, N., Renard, P., 2022. The stochastic simulation of karst conduit network structure using anisotropic fast marching, and its application to a geologically complex alpine karst system. Hydrogeol. J. 927–946. https://doi.org/10.1007/s10040-022-02464-x . To perform properly this review, I have also carefully read both these papers to clearly understand the added value of this new one.

In this new submission, the concepts, tools, and data are identical to those presented in these two previous papers. The 2021 paper already presented the SKS approach (which was initially published by Borghi et al., 2012 and used in an inverse approach in 2016) to generate ensemble of models. The main contribution of this first paper was the application to the Gottesacker case study. The innovation of the previous 2022 paper was to use an anisotropic fast marching algorithm (provided by an already existing external library: skfmm) instead of the initial version that used an isotropic fast marching . This change was necessary to consider the hydraulic gradient when using a 2D method. A deepest comment of this paper is out of the scope of this review and will not be provided. But this previous paper also proposed a very short sensitivity analysis, and, in particular, it explored the effect of various pairing of sources and sinks on the resulting networks for the Gottesacker karst (fig. 18 of 2022’s paper). Note that the modelling results are compared to a “reference” which is not the real karstic system, but a conceptual representation proposed by Chen et al. (2018) allowing to reproduce the hydrological observations. In the present work, the two hypotheses for the formation of the inactive karst conduits are translated in the modelling approach as changes in the position and pairing of inlets/outlets. They are described in 22 lines page 9 + the figure 4. All the remaining of the methodology (page 1-8) is a rewording of the previous papers. As well as a large part of the discussion and of the “messages” of the paper. As the inlets-outlets pairing effect was already discussed in the previous 2022 paper, the example of a pairing guided by a speleogenetical history would have been more pertinent as 15 additional lines in the previous paper than as a “new” 16 pages paper, avoiding the large number of redundancies.

My second main concern is about the interest of using the model to answer the scientific question “should the inactive conduit having been formed during the glaciation or by a late opening of the QS spring?”. Indeed, the inactive conduit is oriented north-south. The modelling algorithm searches a path between an inlet and an outlet, with a secondary influence of fractures (randomly generated) and a pseudo-hydraulic gradient. If you do not put an inlet / outlet on the south of the paleo-spring, there is no reason (theoretically and numerically) to generate a path aligned with the inactive conduit: from top view, the path “south -> QO spring” is opposite to the input (= topographic) gradient. No need for the model to see that. Thus, as QO is the only outlet in hypothesis 1, only an inlet located approximately in its south could explain such a path. As the test performed by the authors does not propose this solution, it is obvious that they are not going to generate a consistent solution. In addition to this first point, in the hypothesis 2, some paths are generated accurately, but it is because the authors let the spring QA exist: thus, some paths are generated between QO and QA which is a main direction aligned with the inactive part. Again, this is consistent with the topographic gradient and the direction QO-QA, no need for the model to guess this hypothesis is consistent with the observations.

But, and this is really important, what about the direction of fluxes, never discussed (a limitation due to the fact they simplify in 2d)? Indeed, in hypothesis 2 these paths are obviously going from QO (higher) to QA (lower): in that case, QO is not the spring of the inactive part but an “entrance”. What about the field data? Are the inactive conduits indeed sloping towards QA or in the opposite direction? Were there indices of QO being an inlet and not an outlet? It is not said so, in the text, where QO is presented as an outlet. If the conduits are sloping towards QA, then QO cannot be their previous outlet, and the tests performed for hypothesis 1 were, from start, bound to fail (again, no need to perform any computation to conclude that). If, oppositely, the inactive conduits are sloping towards QO, then their “dead-end” should be the place of a past inlet. And paths will be generated between them and QO, independently of the presence or not of the paleo-glacier (thus independently of both hypotheses). Most of all, it will imply that what the authors consider as a proof for benefit of hypothesis 2 is wrong as hypothesis 2 generates conduits sloping in the opposite direction.

Still on this second aspect: by definition, stochastic simulations are stochastic, not determinist. Thus, they do not aim to find the “true location” of a conduit but to propose equiprobable paths given a limited knowledge. In the present approach, the paths are found by fast marching. The cost is a combination between aerial distance (in 2D here) and topographical gradient (because it is what is used here), and fractures. Only the fractures change from one simulation to another, thus, it just allows to create a kind of a glow around the two main elements of the cost function: gradient and distance between the two points of the considered pair. Thus, as soon as a hypothesis is consistent with these two elements, no further ensemble simulation is required. The stochastic approach does not seem to provide any particular interest in that precise case.

To finish with the modelling approach used here, the 2D approach simplification seems not appropriate as it completely ignores the staggering of the karstic conduits (illustrated in theory by the authors in their figure 1). The Borghi et al. algorithm was however 3D and other works by other authors proposed 3D approaches, why not having worked with them? In this precise case, the topographic gradient is used to mimic the hydraulic gradient, but what does the 2D maps represent? Do the authors consider that the conduits are vertical until they reach the 2D map level? In that case, should the paths computed by the model be seen as conduits developed along the piezometric surface? If yes, which one: today’s surface or a paleo-water table? The paper is clear: the authors use only something which they associate to today’s surface, but it is not consistent with their own initial illustration in figure1. If there is a scenario of water-level dropdown like the theory in figure1, this level is moving through time, and thus the “2D equivalent level” (and map) is changing… This is not at all considered by the model here. Another critical point related to this : if what we see in 2D maps are the conduits developed inside the limestone level after the conduits have vertically crossed the highest part to reach the water-table, then on this virtual surface, everything should be limestones => the geological units presented on this 2D maps have no meaning, more precisely, the sandstone and overlying units should be ignored: they have been crossed by the conduits vertically and on the virtual surface we consider, which is below, we are always inside the limestones (except when the limestones are completely eroded). Saying the same differently: the 2D units presented on the maps seem to be what is seen at the surface, while the network develops at depth, inside the Schrattenkalk. This problem was already there in the previous papers.

These are the main reasons for my recommendation. Below are some additional comments, provided following the order of the paper.

“On the flow” comments:

The figure 1 is misleading: it is presenting the principles of a two-phase karstification due to a drop in the base level, resulting in a two-level network. But the used approach in the paper is 2D and can not consider different levels of karstification as elevation is not taken into account in a map view. It is indirectly considered with the anisotropic fast marching approach, where only a global gradient, here parallel to the topography, is introduced to try mimicking the gravity effect. Using such a figure is not appropriate if it does not correspond to the spirit of what is feasible with the approach and if it does not correspond to what is tested/investigated in the following.   

Section 2, section 4 and section 5.2 almost say differently the same things: the hypotheses should not be repeated and split in various incomplete parts but regrouped in a single complete section.

Line 53: “this study presents a model-based approach to reconstructing the geologic processes driving cave formation”. I disagree, the modelling approach builds networks by considering the influence of user-defined geologic influencing factors. It does not reconstruct a geologic process, as no physical rules/equation are used. 

Line 57: why only two hypotheses? Why only the QS being younger? What about the “ages” of inlets? Inlets could also change depending on the erosion of overlying deposits…

Figure 2:

• use different colours for active vs. inactive conduits.
• Currently and considering the conceptual schema of Chen et al., QO is linked to QS: on the figure it is only a stream: what happens at depth? Do we have information about the connections between the QO point on map and QS point?
• The consistency between the 2D map and schematic cross-section is not clear : contrary to the map, on the cross section, the marls do not outcrop in the valley, as well as the sandstone on the Hoher Ifen, the marls at the same point, then again the sandstone at Gottesacker, etc… Looking at the geological map provided in the 2021 paper, the 2D map used here is also not consistent : limestones should outcrop almost everywhere. The geology is not complex, and the rasterizing effect does not justify all these differences. As this seems to be used in the cost function, this is questioning.
• On the cross section, N1, N6, N11 and N16 should be indicated to help the reading.

Holloch cave: show developed cross-sections of the conduits if possible (vertical organisation?)

Line 107: hypothesis 2 not clear: if the connection is from QO to QA, what should today explain the dis-activation as QA still exist?

Line 127: in 2022 paper, the hydraulic gradient was said to coincide with the elevation of the bottom surface of the limestone for the synthetic case. Here, it is said it is the topography, with the justification “it is simpler to calculate” (line 123): how do you explain such differences? The cross-section in figure 2 is not so evident with different level of erosion for limestone (the Hoher Ifen for example and its surroundings). “Simpler to calculate” is not relevant if the consequences are large (see remarks in above comments). 

Line 136: “The conduits leading to each outlet can be simulated in separate iterations, to represent springs of different ages”: in SKS there are not 1 iteration but several ones as you simulate one path between 2 points, and then a second one which is influenced by the first one (second iteration). The authors explain it well in their first paper, thus here this sentence is unclear: what was the real message to pass?

Section 5 and 5.1: a resume of previous papers (but see point 1 of the main comments).

Line 155: referring to Fandel et al 2021 for a detailed description of geology is strange (and partly unfair) as the geology was there taken from the previous works of Goldscheider et al., Chen et al., etc.

Line 157: the model is said “sliced using the topography surface”: why that? the karstic system does not develop at the topographic surface (see also detailed comment above).

Line 158: cells of 50x50m seem quite big for the level of detailed searched here… (the difference in source elevation is largely lower which can induce large border effect in the method, see previous comment on the 2D approximation).

Line 161-171: all already said in previous papers.

Line 163: “expected paths”. These paths are not the karstic conduits in itself (the superposition of the Holloch cave map show this clearly) but a equivalent conceptual model allowing to reproduce the hydraulic response. Why do the authors try to fit it instead of the conduit real paths if they have them? Also, if we have not the real paths, as the concepts guiding the conduit modelling approach is not the same than trying to find an equivalent hydraulic model, why trying to compare both model results? They do not have the same purpose.

Line 169: “These results support placing confidence in the ability of pyKasso-generated ensembles to simulate”:  I find it an “over-interpretation”. The end of the sentence: “particularly when the inlet/outlet assignments are fixed” confirms what I said above about the importance of inlet/outlets pairing issue.

Figure 3: Already there in Fandel et al 2022 (said by the author in the legend): why spending so much time on already presented works and results?

Line 192: what proves that the other inlets N1, N6 and N11 exist at the time of the glacier? Is there any field element for that?

Line 194: “The existing inlets remained the same”: why, what support this assumption?

Findings: See the critics in my top comments.

Discussion: apart for the large number of repeated points from previous papers:

• Line 251-257: globally break open doors.
• Line 261: “This is likely a limitation of the model’s ability and our simplified assumptions to predict exact conduit locations rather than general orientations and connections” As I said in the above comments, the goal of stochastic approach is precisely NOT to predict the exact conduit location. If you want an exact prediction, you have to use a deterministic approach. I am surprised by this sentence.

Section 7.1: Not very informative, it is globally a rewording of what was previously said.