A new efficient calculation procedure was developed for estimating the input parameters for groundwater model simulation to explore the locations of contamination sources. The authors demonstrated the step-by-step performance evaluation of the procedure through hypothetical field conditions to near real-world situations. Although the manuscript contains several points of uncertainty or illegibility, as noted below, the overall structure makes it possible to understand the authors’ great work as described above.

・ In Subchapter 4.3, the authors described Case 3 as a real-world groundwater problem; however, insufficient information was provided for the regional groundwater flow system. At a minimum, “aquifer structure (number of layers),” “number of grid cells in the depth direction,” “water input to the aquifer (i.e., presence or absence of local precipitation),” and “reasonable information to define the regions as four; for example, geological information except the northern and southern boundary areas as mentioned" would make it easier to agree with the process by which the model was discretized from the actual groundwater flow field to the structure necessary to achieve the objectives.

・The title in Chapter 3 appears irregular with the other chapter titles because it includes the specific names of algorithms (SA-CSA-TS). Abstracting the title, for example, to “Overview of the Developed Algorithms, ”would make the overall context of the article easier to understand.

・The manuscript contained many figures, which reduced readability. It would be more readable to consolidate Figures 14, 15, and 17, which are similar; a multi-panel figure with branch numbers 14a, 14b, and 14c can be used. This approach could also be applied to pairs of figures that are “of the same type but depict different experimental cases,” such as Figures 12 and 13 or Figures 10 and 11. If possible, the authors should consider categorizing the figures into those that are strongly relevant to the derivation of conclusions and those that are not; the latter should be moved to the supplement.

・Limitations of this approach at the present time must be addressed. The estimation of three-dimensional flow and dispersion of contamination, which is one of the major purposes of groundwater flow modeling, is not captured by this approach.

・ In this manuscript, a model with a simplified depth layer is constructed. If the layer (as grids) is increased for depth, would the computation time remain within a practical range? I believe that the importance of this question is related to the applicability of the developed method beyond the estimation of contamination locations under the simplified aquifer system (mentioned in the last line of conclusion).

This paper proposed a new synergistic method called SA-CSA-TS based on (1) multi-chain (2) surrogate-assisted (3) hybrid optimization framework integrating the Cooperative Search Algorithm (CSA) and Tabu Search (TS) for Groundwater Contamination Source Identification (GCSI) problem. Three cases with synthetic and practical data were carried out, demonstrating that the proposed SA-CSA-TS method can consistently identify both the values of hydrogeological parameters and the locations and time-varying release fluxes of contaminant sources. The novel contributions of this paper are significant enough to be published on HESS. Only a few minor revisions are required.

1. Section 2.2.2, Page 5. Why was CSA chosen instead of other heuristic optimization algorithms, such as the well known Genetic Algorithm, Simulated Annealing, etc.? The reasons for choosing CSA need to be explained. After reading [Feng et.al., 2021], I got that CSA is a new optimization algorithm with outstanding performance, but it’s not well known in the hydrology community. Please provide more information about CSA.
2. Section 2.3, Page 6 - 7. This paper established a surrogate-assisted optimization method for GCSI problem, which need to construct multiple surrogate models for each well. The computational burden of surrogate model might be too heavy, if the surrogate model itself is too expensive. Consequently, the surrogate model selected in this method should be cheap and effective. The reason of selecting RBF has been sufficiently demonstrated in the discussion part in section 5.2, line 385-390, page 17. But the introduction of the compared surrogate models in section 2.3 is somehow insufficient and misleading.

• The relationship between the Gaussian Processes Regression (GPR) and Kriging interpolation should be explained. The Gaussian Processes Regression originated from machine learning, and Kriging interpolation originated from geostatistics. The two are essentially the same method in mathematics, and the main difference between them is that they have different expressions of the same concept, because they grew up in difference community. It is necessary to introduce what kind of covariance function is used in the Gaussian Processes Regression here and how the hyper-parameters are set. If the mathematical equations are the same with Kriging, only one should be kept.
• Radial Basis Functioncan be used as a covariance function in Gaussian Processes Regression (GPR), and a kernel in Support Vector Regression (SVR). Please list out the setup details of GPR and SVR (e.g. the parameters assigned in scikit-learn function call. Please double check the details if the default parameters were used.), in order to demonstrated what were actually compared in this work.
• GPR(merge with Kriging), SVR and RBF methods has many hyper-parameters, which may have significant influence to the final fitting performance. The hyper-parameters in GPR (mainly in the covariance function) can be manually specified, or automatically optimized, as elaborated in [Rasmussen and Williams, 2006]. The hyper-parameters of SVR (use RBF kernel) C, gamma, epsilon, are usually optimized with a crude grid-search, or manually specified according to expert knowledge. Please double check the hyper-parameters used in the comparison, are they automatically optimized? Or had been specified with default values? The fitting RMSE varies greatly, possibly because the default hyper-parameters are used, which do not appropriate for the problem in this paper. A good practice is automatically optimize, or manually tune the hyper-parameters after initial sampling and first construction of surrogate models.

1. Section 3, page 8. This section is also methodology section, and this is actually the novel contribution of this paper. Some content in this section overlaps with section 2. It is better to move this section as section 2, in order to highlight the research significance. Move the introduction of surrogate models, CSA and TS, as sub-sections within this section.
2. Table 3, page 16. Add RMSE values in this table.
3. These two references are the same book? Please double check.

Rasmussen, C.E., Williams, C.K.I., 2006a. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, USA.

Rasmussen, C.E., Williams, C.K.I., 2006b. Gaussian processes for machine learning, Adaptive computation and machine learning. MIT Press, Cambridge, Mass.