General comment

Comment #1:

The paper presents a study that developed a Surface water-Groundwater model in an Alpine outwash plain. The model's calibration and validation relied on a sophisticated data acquisition system and external information from nearby meteorological stations and previous studies that aimed to identify the bedrock. The study also revealed the seasonal impact on groundwater flow through the use of transient modelling.

The authors then applied their results to estimate the potential impact of outwash cascade effects on the formation of future outwash plains in the current glacier position.

The methodology is interesting, and the results are quite promising. Therefore, I think that the paper has the potential to deserve publication. However, a few points, listed below, should be addressed before possible publication.

Answer #1:

We thank the referee for his/her careful review of our paper and his/her detailed, constructive comments. In the following comments, we will address the main referee’s comments.

Comment #2:

In line 165, the authors utilized the unconfined aquifer formulation for groundwater modeling. Given the shallow depth of the aquifer body, it's possible that the vadose zone could significantly impact your results. Have you assessed whether neglecting the vadose zone is applicable to your domain? It would be valuable to explore this possibility further to ensure that the assumptions align with the specific conditions of your study area.

Answer #2:

In this study, we indeed decided to neglect the vadose zone. Due to the very coarse texture of the outwash plain sediments, it is expected that the capillary fringe is thin and that the water content in the unsaturated zone above the water table decreases rapidly. In our case, the soil texture is composed of more than 90% sand and usually less than 1 to 2% clay. This leads to a soil texture dominated by large pores where the capillary effect is limited. This statement has been tested by building a very simple HYDRUS-1D model which solves the Richard’s equation and allows accurate modelling of the interface between a saturated and unsaturated zone. It shows that the water content above the saturated zone drops rapidly to values close to 0.1, in the first 50 cm. Similarly, the unsaturated hydraulic conductivity drops sharply and also leads to a very dry soil surface, where actual evapotranspiration is strongly limited. This will be made clearer in the revised version.

Comment #3:

Additionally, during the winter period in the upstream portion of the domain, most of layer 1 is unsaturated. It may be worthwhile for the authors to consider solving a Richard-based model for the first layer and using a simpler formulation for the other layers. Given the ample data available for the studied domain, have the authors considered using UFZ MODFLOW package or other models that allow for this type of modeling, such as ParFlow (Maxwell et al., 2005 https://doi.org/10.1175/jhm422.1) or OpenGeoSys (Kolditz et al., 2012  https://doi.org/10.1007/s12665-012-1546-x)?

Answer #3:

We did not consider to implement any unsaturated flow package based on the assumption stated above. In addition, even though the processes in the superficial unsaturated layer may not be fully well represented, the focus of this study is on the aquifer-scale groundwater flow and exchanges with the stream. Due to the automatic calibration based on an extensive database of observations, we believe that the seasonal aquifer behavior is correctly modelled and the use of a more complex (unsaturated) model would not change the modelled behavior and would not provide additional information regarding the research questions.

We believe that unsaturated processes would only be required for processes such as groundwater - plant interactions such as root water uptake or evaportranspiration, which can be safely neglected in our case.

Comment #4:

Did the authors investigate whether water evaporation, which may be more significant during summer days when solar irradiation is stronger, was negligible? Additionally, it is not clear to me if the authors evaluated the rain runoff during raining events. It would be valuable to understand how rainfall was accounted for in the study and whether it was incorporated into the groundwater model.

Answer #4:

We indeed did not clearly state this in the paper and it should be added (will be done in the revised version). As discussed above, evaporation from the top sediments is highly limited by the coarse nature of the sediments and the limited capillary effect. As such, we neglected this process.

We also did not specify that we also neglected the direct recharge effect of rain falling on top of the aquifer. This simplification was made as we could never observe any groundwater response after a rain event in any groundwater wells. The groundwater heads only responded to rain events indirectly due to the increase in stream level which then propagated in the aquifer as a pressure wave diffusion process (which is correctly represented in the current model).

The reason for this lack of response is indeed likely due to the effect of the dry unsaturated zone, which delays rainwater infiltration in the unsaturated zone. Since we did not model this unsaturated layer, we could not model unsaturated infiltration and therefore neglected surface rain input

Finally, again, we believe that those processes do not impact the aquifer-scale groundwater processes analyzed in this study.

We will add in the methodology comments of the revised document a paragraph to better detail the input boundary conditions. We will also add in the discussion a more detailed discussion based on three previous comments of the reviewer.

Detailed comments

Comment #5:

Line 172: could the author report more detail on the coupling method between surface and subsurface model that the authors had employed among the possibilities offered by MODFLOW? This information could be possibly helpful in providing a deeper understanding of the results obtained.

Answer #5:

We used the most state-of-the-art packages offered in MODFLOW to represent surface-groundwater interactions. The streamflow (SFR) package is the only package which allows to simulate the dynamic stream stages at every node based on a channel topography and  on Manning’s equation. This thus leads to a dynamic estimation of the pressure gradient between the surface and subsurface and thus to a realistic estimation of the water flux.

One key challenge with this package is the definition of the stream channels, which we attempted to represent in a relatively simple way but with enough details to provide a realistic representation of the processes. We will try to provide some more details about this package.

Comment #6:

Could the author provide more detail on how the two objective functions described between lines 195 and 201 are used for calibration? In particular, the authors defined a single objective function that utilizes a weighted sum of the multi-objective function. How was the weight of the OBJ function selected? It would be helpful to provide a clear explanation of this process to ensure the reader understands how the calibration was performed.

 Answer #6:

Concerning the single objective function, we first weighted both objective functions (SWE and snow cover) equally, that is their total initial weighted sum of residuals is equal. In this way, both the objective functions are optimized to a similar degree.. Based on the initial calibration results, we then decided to increase the weight of the snow cover objective function to improve the representation of the snow cover fraction, which tended to decrease too quickly in the model during the late season. This iterative adjustment of the weighting resulted in a total weight of 1 and 2.4 for the SWE and snow cover objective functions. This final weighting allowed to  match well the snow cover fraction over the whole season, while maintaining a root mean square error of about 100 mm for the estimation of the SWE. The final weights correspond to a heuristic choice but that we can (and will) justify in the revised version.

We will provide the main equations and results of this model in an appendix of the revised manuscript. 

Comment #7:

is it possible to include section number 4 in Figure 2?

Answer #7:

Yes, we will try to better illustrate hillslope recharge on figure 2 in the revised document.

Comment #8:

In line 381, you mention that exfiltration is also correlated with stream water infiltration, citing Figure 8c. However, it is not clear from this figure how these two variables are related. It seems that exfiltration has a similar hysteresis effect to upstream discharge, as you note. To gain a better understanding of the relationship between these two variables, it would be helpful to include a graph comparing upstream discharge and exfiltration, if possible.

Answer #8:

We agree that this statement is not clearly illustrated. We will provide an additional figure and a clearer discussion of this processes in the revised document.

Comment #9:

In line 411, could the authors provide more explanation about the selection of the porosity value?

Answer #9:

Yes, we forgot to specify this in the methodology and will be added in the methodology. Porosity was estimated by measuring the saturated water content using a Decagon 5TM at five locations in the upper sediment layer. Additionally, the sediment texture was analyzed by laser granulometry, resulting to a large dominance of the sand fraction (>90%). Such coarse sediment texture usually have a low porosity of 0.25 to 0.3.

General comments

Comment 1:          

Müller et al present a MODFLOW-based numerical model of a proglacial aquifer system in a high elevation Swiss catchment. A wide variety of field observations are used to establish the structure of the model, and to subsequently parameterise and calibrate it. Indeed, a multi-variate objective function is used. The model is capable of reproducing the available observations rather well. A particularly novel aspect is the identification of overdeepenings beneath the present glacier where proglacial outwash aquifers could form in future. The calibrated model is transferred to these possible future locations of outwash plain formation in order to estimate the potential first-order influence on outwash plains on groundwater storage and release in the future and at larger scales. I find the paper competently written and well organised, and commend the authors for their commitment to Open Science (through the provision of a Jupyter notebook, PEST control files etc.). In summary, I am confident that the work will represent a useful contribution to the literature.

Answer 1:               

We thank the referee for his careful review of our paper and his/her detailed, constructive comments. In the following comments, we address the main referee’s comments.

Comment 2:          

Some of the description of the numerical model structure, as well as the forcing data applied under the future scenario (which seem to be missing).

Answer 2:               

We thank the reviewer for his/her careful inputs in the manuscript. As also highlighted by the first reviewer, some description and justification for the model boundary conditions can be improved.

We discussed in the answer to the first reviewer’s comments the choices of not including an unsaturated layer as well as neglecting direct rain input.

The reviewer highlighted some concerns in the annexed pdf which we address hereafter.

The reviewer questioned the simplifications made regarding the stream channel structure. This is indeed a challenging point as MODFLOW does not allow a very flexible definition of those channels. We however used the latest package development to account for a realistic definition of the river cross-sections which was essential to model the surface-groundwater interactions. We assessed the impact of using only the major static channels and this is discussed in lines 450 to 460 of the current manuscript. Although this method is not ideal, we believe that our model structure and calibration procedure to mitigate at least partially such simplifications and provides a sufficient representation of the aquifer-scale interactions.

The reviewer also highlighted some unclear details about the lateral recharge. In this study, we mostly included surface water flow from the hillslope resulting from snow melt and rain. Those inputs were defined as point recharge as observed in the field. No surface flow was observed in the south side of the hillslope, likely due to less solar input and less melt. A diffuse subsurface recharge from the hillslope was not included due to a lack of data regarding this process. We however assume that diffuse recharge is likely not a dominant process in summer  due to the coarse and steep nature of the hillslope and crystalline bedrock as discussed in previous work (https://hess.copernicus.org/articles/26/6029/2022/hess-26-6029-2022.html). In winter lateral subsurface recharge may play a more important role when melt is assumed to be very limited. This process is however difficult to measure in the field. From our work in this catchment in general(see also https://hess.copernicus.org/articles/26/6029/2022/hess-26-6029-2022.html)), we have highlighted a winter baseflow of about 0.3 mm/day, potentially due to bedrock exfiltration. The hillslope side of the outwash plain have an area of about 1.5x10⁶ m²; if we assume a diffuse recharge of 0.3 mm/day, this represents about 5 L/sec of recharge over the whole outwash plain. This flux is at least 10 times lower than the stream infiltration estimated for winter. This short discussion will be further detailed in the revised manuscript.

The reviewer also pointed out some unclear details regarding future scenarios. In the case of the future scenario, we chose here a quite extreme case, where no lateral recharge occurs (no recharge from the hillslopes). This choice was made to highlight the effect of cascading outwash plains only and to highlight their specific hydrological functioning in case of an extreme drought case.

More generally, in this study, we used a rather unconventional approach where we did not model the entire catchment functioning, but rather focused on a specific hydro-geomorphological structure within it and attempted to include all major processes that influence the groundwater storage of this structure. As highlighted by the reviewers, the use of a more complex modelling framework such as HydroGeoSphere, could have provided a more integrated understanding of the processes. We however lack much information about the subsurface geological structure and subglacial processes in other parts of the catchment and believe that the development of a more complex model may be very challenging and may introduce other sources of uncertainties which may more difficult to interpret.

With the use of this hybrid model (highly detailed 3D aquifer modelling coupled with a rather simple lateral recharge routine), we have more control and understanding of the key processes which occur in the outwash plain aquifer specifically.  Since our goal was to focus on the outwash plain and not the whole catchment, we believe this approach is well suited to our research questions. This will be made much clearer in the revised version.

Comment 3:          

In addition, I am unsure whether the 3 day model einitialization period is sufficient (perhaps it is, since the aquifer system seems to form a kind of bedrock-dammed “bucket” in which groundwater levels near the surface are maintained). Several additional (generally minor) comments are provided in the attached PDF.

Answer 3:               

The 3 day initialization period seems indeed rather short but was verified before running the full calibration. The initial groundwater head was initialized 1m below the ground, which is not far from the measured groundwater stage. Due to the very conductive nature of the sediments, the aquifer rapidly reaches an equilibrium with the stream and 3 days (even 2) was found to be sufficient to obtain such equilibrium.

Comment 4:          

More generally, in agreement with Reviewer 1, I also think it would have been interesting to apply a “fully-integrated” hydrological model code such as ParFlow, ATS, or HGS to the problem. Some studies applying these models at catchment scale in mountain settings have emerged (https://www.sciencedirect.com/science/article/pii/S0022169417301002; https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020WR029390); multi-scale work using nested integrated models could perhaps be proposed as a future research direction for such research, since this would enable many strong a prior assumptions made regarding relevant processes, channel locations, recharge locations etc. to be relaxed, whilst retaining detail in the aquifer of interest.

Answer 4:               

We thank the reviewer for this broader comment. As discussed before, we believe our approach also has some benefits compared to a more complex fully integrated model, as it allows us to focus more on one structure, with more control over the specific processes occurring in this aquifer.

We also agree that future work proposing a more integrated modeling of the whole catchment would be interesting but would require more detailed geological data, especially regarding the nature of the bedrock fractures and groundwater flow within such aquifer. We finally believe that using different methods with different levels of complexity, as well as detailed analysis of field data of various sources are all required to provide a fully sound understanding of alpine catchments, as each method is limited by the choices of model structure and assumptions.

We will add a paragraph in the discussion/conclusion part of the revised document to highlight our choices and the potential use of fully-integrated models.

Annexed pdf

We thank the reviewer for his/her further detailed inline comments. All remaining comments will be addressed in the revised version of the paper.