general comments:

The topic of the manuscript is of high relevance for the management of large aquifer systems. The presented approach provides a very suitable methodology to develop a groundwater model with predictive potential for a complex heterogeneous large-scale groundwater system. Such a model supports the understanding of the system dynamics and enables to identify parameters impacting diverse system responses. For the first time stochastic calibration and informed global sensitivity analysis are used to calibrate the groundwater model of the main aquifer system in the Po River watershed. The potential of the presented methodology is well concluded.

The manuscript is well structured and well readable. The title clearly reflects the contents of the paper. The abstract provides a concise and complete summary. The scientific methods and assumptions are valid and clearly outlined. 

The authors give proper credit to related work and clearly indicate their own contribution. The number and quality of references appropriate.

I recommend the publication only after the revisions described below.

specific comments:

In the following I would like to recommend several revisions in order to improve the manuscript.

The lateral extent and the base of the groundwater system should be described more clearly (l. 180/181). This should cover a more detailed the description of the interfaces with the sub-basins (l. 212/213) particularly the vertical distribution inflow boundary condition (l. 219). Furthermore, the basic geologic concept behind the vertical discretization is missing.

In lines 265/266 a reference to the formulas where the targets of the calibration kc and rq might be added. 

It is not clearly described whether the proportion of geomaterials, Fig. S1, is an a priori information or the result of calibration. As this is an important information anyway Fig. S1 should be included in the manuscript and not part of the supplementary material. In order to support the descriptions in l. 351 (and similar descriptions), it would be helpful to have a figure with the distribution of the geomaterials available. 

The concept behind the combination of the information in Fig. 4 is not really clear. Why are Fig. 4a and 4b combined with 4c and 4d? The information in Fig 4a might not really be important for the purpose of the manuscript. Fig 4b might be improved if the difference between simulate and observed heads are provided instead of the observed heads only. In Fig 4d the dark grey and dark red areas could not easily be distinguished as described in l. 359.

The concept behind the combination of the information in Fig. 5 is not really clear. The description of Fig 5a, l. 372-380, is too coarse. The definition of the macro areas is not clearly motivated. There are more areas related to the several macro areas as described. In Fig 5b it does not become clear what the colour distribution in the represents. The comparison described in l. 367-369 does not really become clear from Fig. 5c.

Fig. 6 might be reorganized as different information, v and h, is combined. It does not become clear why different cross section are used in Fig 6 and Fig. 5.

In order to emphasis the importance of the 3d approach it might be useful to describe the Morris indices for all model layers, especially the lower ones. If helpful an additional figure might be provided which might be added as supplementary material. 

Fig. S2 might be included in the manuscript eventually in combination with Fig. 8 as this is an important result of the study.

A reference for the Penman-Monteith model should be added (l. 157) if this not covered by the reference ‘Allen et al. (1998)’. 

technical corrections:

The sizes of the following figures should be increased. The names of the rivers should be readable in Fig. 1. The cross sections in Fig. 5b are not clearly visible. The sediment types in Fig. 5c are not clearly visible and a corresponding legend is missing. Details in the graphs in Fig. 6 are only hardly visible.

The formula '𝑄_𝑠 = 𝑟_𝑞𝑅′_𝑠𝑆_𝑠' should treated as separate equation (l. 214). The further equations should be renumbered then. 

Within Fig 8 the number ‘6.5e5’ is printed.

This paper presents a pioneering framework designed to model large-scale groundwater systems using a comprehensive three-dimensional approach. This approach not only accounts for the dynamics of river-aquifer interaction but also integrates three-dimensional spatial distributions of geo-materials. Notably, the framework excels in achieving a delicate equilibrium between the requisite simplification essential for large-scale groundwater models and the anticipated outcomes. Of particular significance is the utilization of a multi-objective optimization approach for model calibration and sensitivity analysis.

From a general standpoint, I have two major concerns before recommending the paper publication:

This study, is build upon the substantial groundwork laid by Manzoni et al. (2023), a fact that merits explicit acknowledgment in both the introduction and the abstract. Indeed, in my opinion, the contributions made in Manzoni et al. (2023) represent a cornerstone of the framework presented herein, serving as more than just a dataset but rather as a fundamental component upon which this work is built. 

The aspects presented in lines 59-61 may not inherently appear novel 'per se'. The authors should thus better emphasize the unique contributions and novelties of their work to distinguish it from existing literature.

Other specific comments:

In light of Manzoni et al. (2023), the statement highlighting the limitations of data-driven models due to constraints in the quantity and quality of available training data, particularly in large-scale scenarios where data accessibility across the entire domain might be limited, warrants revision. Specifically, it should be contextualized to reflect insights gleaned from Manzoni et al. (2023) and potentially revised to underscore the advancements or strategies proposed in this study to address such challenges.

The clarity of Equation “𝑄𝑠 = 𝑟𝑞𝑅𝑠𝑆𝑠” is  ambiguous. Assuming my interpretation is correct, where the subscript 𝑠 represents the 𝑠-th subbasin, further clarification from the author is needed regarding whether the efficiency term is constant across time or space. Additionally, a more detailed discussion on the assumption of the correction coefficient, 𝑟_𝑞, being constant for all subbasins would enhance understanding. This could include explanations of the rationale behind this assumption, any potential implications or limitations, and how it aligns with existing literature or empirical evidence.

The description of the domestic water flux associated with groundwater resource utilization lacks the necessary details for reproducibility and applicability in other areas. To enhance the transparency and replicability of the work, it is suggested that the authors provide additional information regarding the data utilized, including its sources and any assumptions or hypotheses underlying its selection or interpretation. This could involve detailing the methodology for acquiring the total volumetric flow rate associated with groundwater extractions of drinking water within a municipality, as well as specifying the criteria used to define the surface area covered by the municipality. By providing this supplementary information, the study's findings can be better contextualized and applied to different geographical regions or scenarios.

The section pertaining to multi-objective calibration appears to lack clarity and would benefit from a more precise formulation of the mathematical framework. It is recommended to clearly define the two objective functions, delineate the independent variables, and specify any constraints involved in the optimization process. Providing explicit details about the mathematical formulation will enhance understanding and facilitate the replication of the calibration methodology.

The authors should provide further clarification regarding the statement "We note that CCDE does not require defining a single weighted multi-objective function, as otherwise required by the standard DE." It would be beneficial to elaborate on the distinction between the adopted algorithm (CCDE) and the standard NSGA-II algorithm. For example, the NSGA-II algorithm does not inherently require a weighted multi-objective function either. Algorithms such as NSGA-II typically operate based on the concept of dominance and the Pareto front without necessitating explicit weight assignment during the optimization process. Therefore, expanding on the characteristics of CCDE and how it differs from conventional multi-objective optimization algorithms like NSGA-II would provide valuable insights. Furthermore, it appears that the original paper by Storn and Price (1997) primarily focuses on single-objective optimization, similar to the work by Trunfio in 2015. Given that the current study deals with non-contrasting objective functions, it may explain the absence of weight utilization. However, providing additional details on the rationale behind the choice of CCDE and its suitability for handling the specific characteristics of the optimization problem in this context would enhance understanding.

It is suggested that the authors include a figure illustrating the framework chain, depicting the data inputs required as well as the expected outcomes. This visual representation will aid in comprehending the workflow of the methodology and provide a clear overview of the research process. Additionally, the figure can serve as a useful reference for readers to understand how various components of the framework interact and contribute to the overall analysis.