The authors present a study in which they couple SWAT+ with  gwflow model and examine the parameters, sensitivity, and uncertainties for four study areas in the United States. The areas differ not only in terms of climatic conditions, but also in terms of hydro(geo)logic processes such as water withdrawal, GW-SW interactions, etc.

They demonstrate how a systematic workflow using relatively recently developed tools such as iES can be used within the PEST framework to take full advantage of PEST and its utilities in calibrating groundwater and surface water models and in analyzing the uncertainties in the predictions made by these models before and after calibration.

Although I believe that more studies of this type are needed, and that uncertainties and sensitivities are still too infrequently determined systematically, I do have comments that should be included.

In general, I find that the current manuscript does not discuss enough the results in a broader context (a discussion section would be extremely useful). For example, it could be discussed why bed_thickness is an important parameter in almost all study areas. Here a stronger link to existing literature should be made. Also, a comparison/discussion to integrated models could be made (such as, HGS, Parflow, etc.) and the assumptions of the coupling SWAT+&gwflow for the groundwater-surface water interaction could be discussed.

Model simplifications can be extremely important (even if difficult to quantify) and should be discussed more. Otherwise, it seems that with the presented approach all uncertainties are known very precisely, which is not the case.

Simplifications and limitations in the calibration should also be discussed more intensively. Calibration against GW annual average values are certainly not optimal. Especially since it has often been shown that a calibration against groundwater levels alone is often not sufficient (see for example Schilling et al 2019).

Also, I am not sure how the calibration has dealt with parameters that are correlated. Is this corrected by the ensemble approach or do the dependencies still exist? Some more information would be useful . Also I don't quite understand how the flux from adjacent aquifer was determined. Is this flux also calibrated? From Table 4 it does not seem so but since I am not sure of all the abbreviations, it might be good to at least include in the SI the full name with a small description of the parameters from Table 4.

Other minor comments.

Line 60ff: How is this different from the Null-Space Monte Carlo approach (e.g. Alberti et al.2018, Herckenrath et al., 2011, Moeck et al., 2020) One-two sentences would certainly be useful.

Line 95ff: The calibration/validation selection is relatively classic (which is ok) but perhaps more information could be gleaned from the data when "extremes" are considered. This is certainly something that could be addressed in a discussion.

Line 191ff: Monthly flow data is used, but the time series can be decomposed to get different components such as monthly or event quantities; base and fast flow, flow duration statistics, etc. to use in a multi-component objective function. If done properly, the processing behind each component should distill information of a different type from the measurement data set. Each of these components should therefore inform different parameters or groups of parameters. This can be done for surface water but also for groundwater data (Heudorfer et al., 2019).

Figure 2: How are the data be used between the different weather stations (simple interpolation)?

Figure 14: Even if the general trend is clear and also described how the calculation takes place, I have my difficulties with the values on the x and y axes. For example, the values of mu vary from 0 to 1250 or 30000 depending on the study area. Does this say anything about the significance/sensitivity? Some more context would be extremely useful in the text.

Figure 16: Seasonality in some parameters should be discussed in the text and also why there is a trend in bed_thick for Winnebago and bed_k for Cache.

Figure 17: Even with the posterior ensemble some observed streamflow data cannot reproduced. Does that mean the model is too simple for some events?

Table 7: why is stream seepage negative? Moreover, the recharge values seems to be very small compared to a visual (and only simple) comparison to the Reitz et al. (2017) data. It does not mean that the Reitz data are better (or the method) but please check if the data you simulated are in a plausible range.

***********************************

Schilling, Oliver S., Peter G. Cook, and Philip Brunner. "Beyond classical observations in hydrogeology: The advantages of including exchange flux, temperature, tracer concentration, residence time, and soil moisture observations in groundwater model calibration." Reviews of Geophysics 57.1 (2019): 146-182.

Moeck, C., Molson, J., & Schirmer, M. (2020). Pathline density distributions in a null‐space Monte Carlo approach to assess groundwater pathways. Groundwater, 58(2), 189-207.

Alberti, L., L. Colombo, and G. Formentin. 2018. Null-space Monte Carlo particle tracking to assess groundwater PCE (Tetrachloroethene) diffuse pollution in north-eastern Milan functional urban area. Science of the Total Environment 621: 326– 339.

Herckenrath, D., C.D. Langevin, and J. Doherty. 2011. Predictive uncertainty analysis of a saltwater intrusion model using null-space Monte Carlo. Water Resources Research 47.

Heudorfer, B., Haaf, E., Stahl, K., & Barthel, R. (2019). Index‐based characterization and quantification of groundwater dynamics. Water Resources Research, 55(7), 5575-5592.

Reitz, M., Sanford, W. E., Senay, G. B., & Cazenas, J. (2017). Annual estimates of recharge, quick‐flow runoff, and evapotranspiration for the contiguous US using empirical regression equations. JAWRA Journal of the American Water Resources Association, 53(4), 961-983.

1. Please make sure the abbreviation is defined at the very beginning, GLM, NHD, SWAT, etc.
2. For table 1, what is the reason of using 1000 m instead of 500 m resolution for Arkansas Headwaters?
3. Figure 1, please add unit of DEM?
4. Figure 7, it looks like the predicted values are smaller than observed value in most cases, what is the reason? It needs more discussion to evaluate the model performance.
5. Figure 9, why there is only 2 points in the mean absolute error plot?
6. Could you add more detail how you calibrate the model. Until how many iterations you consider good calibration?
7. More studies with similar hydrology parameters should be tested to see if the method could have consistent performance.