Dear Referee,

We would like to thank you for your comments and review aimed at improving our paper to make it more interesting to a wider audience. We will modify our paper following your recommendations. We agree that the results need to be presented in a broader, less local-context way and we will include modifications to position our paper within the context of cold and humid climate research on groundwater recharge (GWR). To reflect this, the title will be modified to “Simulation of long-term spatiotemporal variations in regional-scale groundwater recharge: Contributions of a water budget approach in cold and humid climates”. The introduction will be modified to position the study within the recent research on GWR in cold and humid climates, such as Grinevskiy et al., 2021; Kløve et al., 2017; Nygren et al., 2020; and Pozdniakov et al. 2020 (complete references at the end of our response). As well, the sub-sections 5.1 and 5.2 will be modified to compare the regional-scale results to GWR estimates obtained in likewise cold and humid climates such as in Scandinavia and northern Russia. It will include a comparison of the GWR trends within the warming climate observed in these regions since the middle of the 20^th century. Thus, it will be possible to discuss further the contributions and limitations of water budget approach in these climates. Hereafter we describe the main modifications that will be made to the document based on your comments. The general comments are addressed first (labelled C#), followed by the minor remarks (labelled mC#).

C1: Lines 172-175. The Authors mention three basic hypothesis. I have serious concerns on the third one: the watershed response time is shorter than one month, thus compensating for the absence of water routing. This is a very strong hypothesis, as it allows neglecting transient dynamics of the aquifer. Before proceeding with any computation, Author should demonstrate and convince the reader that this a reliable assumption.

A1: We agree with the Referee that this hypothesis needed to be justified. However, monthly time steps are frequently used in groundwater flow modelling and are considered a reasonable timescale for groundwater flow processes in cold and humid climates (Dripps and Bradbury, 2007; Guay et al., 2013). To avoid misleading the reader, L172-177 will be rephrase: “The HydroBudget model is calibrated on superficial watersheds, based on the hypotheses that: 1) surface watersheds match hydrogeological watersheds, 2) the rivers drain unconfined aquifers. Under these conditions, for any given watershed, potential GWR should be similar to river baseflow at the outlet, and the sum of runoff and potential GWR should be equal to the total flow at the outlet. Dripps and Bradbury (2007), and Guay et al. (2013) have shown that results of GWR water budget models (SWB and HELP, respectively), compared better to groundwater flow models (GFLOW and CATHY, respectively) at the monthly time step than at a daily time step. Thus, runoff, AET, and potential GWR were computed with HB with a daily time step, and the results were aggregated on monthly time steps. Calibration and validation of the model and analysis of the results were all performed on the monthly aggregated results.”

C2: The model accounts for 8 parameters to be calibrated. It seems to me that some of them are set of parameters. For example the Runoff factor (as explained by the Authors “Partitioning between runoff computed with the RCN method and infiltration into the soil reservoir”) do depend also on the land cover or not? in the first case, you are calibrating each runoff factor for each soil. Am I wrong? Maybe I do not understand correctly.

A2: We thank the Referee for this comment, and we agree that there is a need for clarification. The runoff computed with the RCN indirectly depends on the land cover since simple land cover classes are used to compute RCN values that are later used to compute runoff for each grid cell. However, a single value of the runoff factor parameter (f_runoff) is uses for the entire watershed (one value for all grid cells). As detailed in Dubois et al. (2021; eq(1) p°9) : R = f_runoff x R_RCN with R the runoff for each grid cell and R_RCN the runoff computed for each grid cell. To ensure that this is clear, we will change the description of f_runoff  in Table 2 for “Correction factor for runoff computed with the RCN method for the partitioning between runoff and infiltration into the soil reservoir”. As well, L150 will be rephrased as follows: “It is based on simplified process representations and is driven by eight parameters that need to be calibrated. These parameters are held constant for all the grid cells and through time (Table 2).”

C3: The coupling of time steps (daily time step for soil modelling and monthly time step for GWR) is not clear to me. Please, give more details on that.

A3: We thank the Referee for pointing out the missing information concerning the aggregation on a monthly time-step of the outputs computed with a daily time-step. The model uses daily precipitation and temperature data to compute daily runoff, evapotranspiration, and potential GWR. These daily results are aggregated at a monthly time-step to perform 1) calibration and validation and 2) results analysis. We integrated a clarification of this point in the modified L172-177 presented in A1.

C4: The targets for calibration are both the total surface flow and the baseflow. However, the first is observed, whereas the second is estimated (through the Lyne and Hollick filter). In my opinion this is a weak point of the whole calibration procedure: generally speaking I do not like to calibrate a model using output of another model. Even if they match each others, what can I say on the reliability of both? The Authors should at least convince the reader on the reliability of the baseflow estimates.

A4: We thank the Referee for this comment and agree that it is a sensitive point in our research. In eastern Canada and in other similar geo-climatic environments, the topography-driven water table make streams and rivers discharge areas of groundwater systems. In these conditions, baseflow estimated with recursive filters have long been used as proxies of groundwater recharge for model calibration. Although it is an approximate proxy, the lack of local data/studies that would demonstrate its relevance and limits encourage authors to continue using river baseflows in GWR and other groundwater-related studies (Rivard et al., 2009). In section 5.4 on method limitations, paragraph L472-483 acknowledges the use of baseflows as an imperfect proxy for calibration. This paragraph also underlines the usefulness of using a baseflow estimate based on recursive filters for simulation in an area where mainly long-term river flow measurements are available. Nevertheless, to avoid any misleading, L178-180 will be rephrased as follows: “The daily flow rates at the outlets were those from the 51 gauging stations (Fig. 1). Baseflows were estimated from the flow rate time series calculated following the proposed approach of Ladson et al. (2013) based on the Lyne and Hollick (1979) recursive filter, using a stochastic calibration and 30 filter passes. As aquifers generally discharge into superficial water bodies in southern Quebec, baseflows estimated with recursive filters are considered a satisfactory, although approximate, proxy of GWR dynamic.”

C5: The previous one is a very important point, also because the objective function (eq. 1) is a linear combination of two different metrics, one referring to the total flows, the second one only to the baseflows. In my opinion, dependence of the calibration on the weights adopted to define the objective function should be explored much more in details. The explanation you gave for your choice (lines 203-205) is too simplistic.

A5: We thank the Referee for pointing this out and we agree that our explanation about the weights on the objective functions was too simplistic and may have been misleading. The caRamel algorithm performs the multi-objective optimization using the objective functions separately, here KGE_qtot and KGE_qbase, and plots the evolution of the Pareto front (model runs function of the objective functions). Therefore, the weights used to compute the KGE_mean are only used to extract the best compromise from the ensemble of model runs computed during the optimization and do not impact much the choice of the best compromise. To illustrate this and provide a convincing argument, supplementary table A1 will be completed to values of the objective function based on other combinations of KGE_qtot and KGE_qbase (see below). L203-205 will be rephrased as follows: “The weights attributed to each objective function in KGE_mean were arbitrarily chosen to select the calibrated parameter set that maximizes the quality of simulated baseflows, considered as a proxy for GWR (KGE_qbase), without losing the benefits of the multi-objective optimization. Supplementary Table A1 shows that the objective functions and the calibrated parameters were moderately sensitive the weights.”

*F_T is held constant at the nominal value of 16.4 d

C6: I found the sensitivity analysis performed on the W6 group of gauging stations very interesting and potentially the core business of the paper (which otherwise it is only a modelling study, important for Quebec, but not interesting for the international reader of HESS). Why the Authors decided to carry out the sensitivity analysis only on one group of gauging stations. In my opinion, it could be much more interesting to perform the same analysis to several groups of gauging stations to explore the possible differences among watershed and the dependences on climate forcing, soil, land cover etc In other words, has the ranking proposed for W6 a general validity? why?

A6: This is an excellent suggestion. We are currently carrying out the sensitivity analysis for the rest of the gauging station groups. As soon as the runs will be over and the results are interpreted, we will post here a supplementary answer with the update. We thank the Referee for this helpful comment.

C7: Presentations of the results are sometimes not easily readable. For example, in section 4.3 one can find several information already shown in figure 6 and table 4). I think that the description of the results (sections 4.2, 4.3, 4.4) can be much simplified. On the contrary, results on the temporal evolution (analysis of trends) deserves for much more attention and a dedicated picture to presents results.

A7: We agree with the Referee that results presentation would benefit from being simplified. To make the manuscript clearer and more concise, we will withdraw fig. 6 since the presented information are redundant with table 4. As suggested, we will develop the trend analysis since the length of the time series would allow to follow the GWR evolution since the 1960s. We will include our proposition of the new sections in our future answer, along with the sensitivity analysis.

C8: The entire discussion on the temporal patterns of groundwater recharge relates the time variation of meteorological forcing to discharge, neglecting possible changes in the land cover. I do not know if such an assumption is correct, however it should be explicitly stated

A8: This is a very pertinent comment. We have worked under the assumption that no major land cover change occurred during the simulation period. However, we did not mention this in the text and we thank the Referee for pointing this out. The following sentence will be added at the end of the paragraph L172-177: “As well, the calibration and validation of the model were performed under the hypothesis that no major land use change occurred during the simulation period (i.e. land use was held constant).”

C9: Results of figure 6 are surprising: proportion between runoff, AET, GWR does not depend on the watershed (differences among watersheds are in the order of 1-2%). Figure 5(b) shows different patterns: for example over W2 GWR/P spans between 0.1 and 0.15 for most of the grids; on the contrary, over W4, it seems that GWR/P is >0.3 for half of the cells. Maybe I missed something, but it does not seem to me that Fig 5 and Fig 6 are coherent.

A9: We thank Referee for this comment, and we have double-checked our results. Considering the spatial variations in precipitation within the study area (952 mm/yr in W1, 1 039 mm/yr in W2, and 1 123 mm/yr in W4; table 1; fig. 2), variations are found between absolute values of potential GWR (109 mm/yr in W1, 119 mm/yr in W2, and 147 mm/yr in W4; table 4). The ratios of average potential GWR/average precipitation for these three watersheds are 0.11₄, 0.11₅, and 0.13₁, respectively. When computed from the spatially-distributed GWR/P ratio presented in fig. 5b, the averaged ratios for the same three watersheds are 0.11₄, 0.11₃, and 0.13₁, respectively. Results are equal and we think that the difference in precipitation explains that the two figures are coherent although it seems that the GWR rates would be higher in W4.

Minor comments:

mC1: Line 102-109 All the information given here are summarized in Table 1. This lines appear to me not useful.

L97-108 will be rephrased as follows: “The study area is located in the province of Quebec (humid continental climate), between the St. Lawrence River and the Canada-USA border, and between the Quebec-Ontario border and Quebec City (35 800 km²) (Fig. 1). It is comprised of the watersheds of eight main tributaries of the St. Lawrence River (numbered 1 to 8 from west to east) (Table 1). Watersheds W1 (Châteaugay River), W2 (Richelieu River), and W4 (Saint-François River) are partially located in the USA (42%, 83%, and 15% of their total areas respectively). Topography is flat with low elevation areas close to the St. Lawrence River and higher elevations in the Appalachian mountain range, associated with steeper slopes. Land cover includes agriculture, forest, wetlands, urban uses, and surface water (Fig. 2a). Agriculture dominates in the watersheds located in the St. Lawrence Platform, while forest occupies most of the Appalachian watersheds.”

mC2: Line 112. In my opinion it would be better to add also the mean bias as metrics of the goodness of interpolation, to avoid systematic under or overestimation.

L110-112 will be rephrased as follows : “The high density of measurements during this period generated minimal error on the interpolated data, with a root mean square error (RMSE) of 3 mm/d for precipitation (mean bias of +0.1 mm), 2.5°C for minimal temperature (mean bias of -0.5°C), and 1.5°C for maximal temperature (mean bias of +0.1°C; Bergeron, 2016).”

mC3: Figure 1. In my opinion the map in the middle is not useful: I suggest to eliminate it, retaining the location map and the watersheds map.

Fig. 1 will be modified as presented in the supplementary PDF file.

mC4: Figure 4. As GWR is your main output, I would present several graphs as the panel (b) for several stations (one as an example is not useful)

Fig. 4 will be as presented in the supplementary PDF file.

The associated paragraph L247-267 will be rephrased as follows: “The calibrated HB model was used to simulate potential GWR for the entire study area on a 500 m x 500 m grid for the 1961-2017 period. Examples of the resulting monthly vertical inflows, baseflows, and potential GWR are illustrated in Fig. 4 for the downgradient stations in W3, W7, and W8. The simulated potential GWR compared favorably with baseflows estimated using the Lyne and Hollick (1979) digital filter. Maximum values were reached simultaneously in April, during the spring month(s) of maximum VI. A second baseflow and GWR peak was observed and simulated in November-December of most years. Lowest values were reached in July to September (high AET rates) and February (minimum VI). Similar matching results (timing and amplitude) were obtained for the river flow (not presented here). Simulated AET was null in the winter until the spring thaw, after which it quickly reached its highest value (> 100 mm/month) in July and decreased at the end of August, reaching null values again in November. Comparable results were obtained for the other gauging stations.”

References:

Grinevskiy, S. O., Pozdniakov, S. P., and Dedulina, E. A.: Regional-Scale Model Analysis of Climate Changes Impact on the Water Budget of the Critical Zone and Groundwater Recharge in the European Part of Russia, 13, 428, https://doi.org/10.3390/w13040428, 2021.

Kløve, B., Kvitsand, H. M. L., Pitkänen, T., Gunnarsdottir, M. J., Gaut, S., Gardarsson, S. M., Rossi, P. M., and Miettinen, I.: Overview of groundwater sources and water-supply systems, and associated microbial pollution, in Finland, Norway and Iceland, Hydrogeol J, 25, 1033–1044, https://doi.org/10.1007/s10040-017-1552-x, 2017.

Nygren, M., Giese, M., Kløve, B., Haaf, E., Rossi, P. M., and Barthel, R.: Changes in seasonality of groundwater level fluctuations in a temperate-cold climate transition zone, Journal of Hydrology X, 8, 100062, https://doi.org/10.1016/j.hydroa.2020.100062, 2020.

Pozdniakov, S. P., Vasilevsky, P. Yu., Grinevskiy, S. O., Lekhov, V. A., Sizov, N. E., and Wang, P.: Variability in spatial–temporal recharge under the observed and projected climate: A site-specific simulation in the black soil region of Russia, Journal of Hydrology, 590, 125247, https:/

Dear Referee,

As a complement to our previous response, we would like to share here the results of the sensitivity analysis performed on the eight group of gauging stations. We have used the method initially presented for W6 on all the watersheds in the study area. We obtained the following ranking (Figure 1 will be added to our revised manuscript), with the relative sensitivity of simulated river flows to parameters changes presented on the left panel (a) and the sensitivity of simulated potential GWR on the right panel (b).

L220-230 will be rephrased as follows: “The model sensitivity to its parameters for the eight groups of gauging stations was obtained with 60 repetitions of the design (540 model runs). The relative sensitivity of the model to some parameters varied markedly between the groups of gauging stations for the simulated river flow (Figure 1a) but appears to be more constant for the simulated potential GWR (Figure 1b). River flow simulation was mostly sensitive to the snow-related parameters (T_M and C_M), except for the western watersheds where the f_runoff was more important. The simulated flow rates were less sensitive to the other parameters. The simulated potential GWR was most sensitive to f_runoff and least sensitive to snowmelt parameters (T_M and C_M) for all the watersheds. The ranking from the second to the fifth highest sensitivity of potential GWR varied from a group of gauging stations to another but was relatively similar. Although the model clearly showed limited sensitivity to the soil freezing time (F_T) for the two simulated variables, they seemed slightly more important for the eastern watersheds W7 and W8. Overall, the potential GWR was more sensitive to parameter variations than river flow since all the μ* for the river flow were lower by a factor 2 to 10 than for the potential GWR (values not presented here).”

As well, L445-450 will be modified as follows: “The impact of long and cold winters was included in HB through the widely used degree-days method that represents snowpack evolution (Massmann, 2019), and through the representation of freezing soil conditions with a threshold temperature and a duration of the threshold temperature to freeze the soil (TT_F and F_T respectively). The sensitivity analysis shows that the simulated potential GWR is sensitive to TTF, while both flow rates and potential GWR have limited sensitivity to FT (Figure 1). The colder watersheds seemed more sensitive to these parameters while the simulation of river flow in the warmer watersheds were less sensitive to the snow-related parameters. This result underlines the importance of including soil freezing in GWR modeling for cold regions. Specific studies on winter recharge would allow to deepen the processes involved and to propose more elaborate representations of this phenomenon.”

Dear Referee,

We would like to thank you for your feedbacks and propositions of improvement for our manuscript. We will take into account all your comments in our modifications. Hereafter we describe the main modifications that will be made to the document based on your comments. The general comments are addressed first (labelled C#), followed by the minor remarks (labelled mC#).

C1: In response to reviewer #1, comment 5, you provide a helpful table of how calibration results change with different weights in the objective function. You say that the calibration is “moderately sensitive to the weights” – but what implication does this have for the results? Do these results all fall within the sensitivity analyses? The key question here is does changing the objective function change the GWR estimates and trends? Would your interpretations and conclusions change?

A1: We would like to thank the Referee for this helpful comment. To demonstrate that these weights have only a “moderate” impact on the results, we are planning to further amend the supplementary Table A1 presented in our previous answer to include the interannual potential GWR. As well, we computed the trends on the potential GWR (Mann-Kendall test, significant if the p-value < 0.05) and since they were systematically positive but not significant, we did not include them in the table. The L203-205 will be rephrased as follows: “The weights attributed to each objective function in KGE_mean were arbitrarily chosen to select the calibrated parameter set that maximizes the quality of simulated baseflows, considered as a proxy for GWR (KGE_qbase), without losing the benefits of the multi-objective optimization. Supplementary Table A1 shows that the objective functions, the calibrated parameters, and the interannual potential GWR values were only moderately sensitive to the weights, with a variation of potential GWR between 176 mm/yr and 194 mm/yr when for the three tested weights (40%-60%, 50%-50% and 60%-40% for KGE_qtot and KGE_qbase respectively). No significant trends in annual potential GWR (Mann-Kendall test, significant if the p-value < 0.05) were found within these 75 model runs (all positive but not significant).”

Table A. 1: Best compromises, interannual potential GWR, and  optimized parameters obtained with the best compromises for the gauging stations of W6 and ranges obtained with the 25 best compromises (before regionalization) using different calibration weights and different baseflow separation methods.

C2: Further to #1, and in line with reviewer #1’s comment 4, what effect does the selection of baseflow separation method have on the results? As you point out in the discussion and in the response to reviewer #1 there can be significant variance in baseflow estimates between different methods. If you selected a different method, would the interpretation change? I acknowledge that these analyses may be time consuming, but I think they speak to the robustness of the approach. These baseflow estimates come with such high uncertainty that I feel they should almost be treated as another parameter – how would your results change if baseflow varied?

A2: We understand the Referee’s concern about that the choice of the baseflow separation method and its impact on the simulation results. We chose to use the Lyne and Hollick filter (1979) with the proposition of Ladson et al. (2013 – cited in our manuscript) with a stochastic calibration and 30 passes of the filter. When compared to other options, this filter appeared to provide the “most likely” smoothed groundwater signal at a daily time step (Figure 1 – the figure will not be presented in our manuscript). When aggregated to monthly time steps, the higher baseflow peaks of the Eckhardt (2005) and Chapman (1991) filters translated in higher baseflows during the high flow periods (mainly during spring) (Figure 2 – the figure will not be presented in our manuscript). Nevertheless, we compared the effect of the choice of the baseflow separation method on the calibration and the GWR estimates and trends for the gauging stations from W6 in the same manner as we compared the effect of the calibration weights to the results. We performed an automatic calibration with the Echkardt and Chapman filters for the baseflow estimates and extracted the 25 best compromises to compare the results. The comparison is reported in the supplementary Table A1 presented above. As the absolute value of baseflow changes with the filter method, the interannual GWR and calibration parameters vary in consequence (Figure 3 - the figure will not be presented in our manuscript). Similarly to the simulations from the model calibrated with the Lyne and Hollick filter, all trends  in simulated GWR based on the Eckhardt and Chapman filters are positive but not significant (Mann-Kendall test, significant if the p-value < 0.05). The objective functions remain similar, between 0.7 and 0.8, except for KGE_qbase in validation which is slightly lower for Eckhardt and Chapman than for Lyne and Hollick.

To include these new results in the manuscript, the paragraph L178-182 will be rephrased into: “The daily flow rates at the outlets were those from the 51 gauging stations (Fig. 1). Baseflows were estimated from the flow rate time series following the proposition of Ladson et al. (2013) of a standard approach of the Lyne and Hollick filter (Lyne and Hollick, 1979), using a stochastic calibration and 30 filter passes. To compare the effects of the baseflow filters on GWR simulation, baseflows were also estimated with the Eckhardt filter (2005) and the Chapman filter (1991) for the gauging stations of W6. Total flows and baseflows were divided by the area of the given watershed to provide flow values in mm/yr, and thus facilitate the comparison of calibration results between watersheds of very different sizes.”

L472-483 will also be rephrased into: “The simulated transient potential GWR was calibrated with baseflows computed using regressive filters from river flow rate time series. Although Partington et al. (2012) showed that the association between baseflows and GWR is not always satisfactory and different baseflow separation methods can lead to differences in volumes and timing (Gonzales et al., 2009; Zhang et al., 2017), baseflows are generally considered to be an acceptable proxy for GWR in cold and humid climates similar to that of Canada (Chemingui et al., 2015; Rivard et al., 2009). They are widely used for the calibration of GWR simulation (Batelaan and De Smedt, 2007; Croteau et al., 2010; Dripps and Bradbury, 2007; Gagné et al., 2018; Rivard et al., 2013). Using baseflow estimated based on the Lyne and Hollick (1979), Eckhardt (2005), and Chapman (1991) filters influenced the proportion of baseflow in total river flow, and consequently led to variations in simulated potential GWR (Supplementary Table A1). However, the quality of the simulations were comparable due to parameter variations that remained in the possible range of the parameters (Table 2) and trends in the simulated potential GWR remained positive but not significant with the three methods. The structure of water budget models makes them compute GWR as the residual of the water budget, thus propagating the computational error from the other terms of the water budget, namely runoff, interception, evapotranspiration, subsurface runoff (depending on the model complexity) onto GWR estimation (Crosbie et al., 2015; Scanlon et al., 2002). Using baseflows as calibration data for GWR simulation limits the modeling error of GWR rates, but makes them dependant on the baseflow separation method.”

C3: I really appreciate how this work can take trends in baseflow (estimated from measurements at sparse gage stations) and use that to estimate the spatial distribution of potential groundwater recharge across a watershed. However, I think the manuscript would benefit from a more explicit discussion of the uncertainties that propagate through the workflow. Again, this kind of analysis and discussion really helps to make the results more robust and applicable.

A3: We thank the Referee for pointing out that an explicit discussion about uncertainty was missing. To quantify the uncertainty, potential GWR was simulated over the study area with the 100 best sets of regionalized parameters and the uncertainty was estimated with the standard deviation between the 100 model runs for each month. As well, to give an idea of the error comparison between the input data and the simulated variables, the mean bias will be compared (the first Referee suggested to add the mean bias of the input data).

This explanation about the uncertainty computation will be added after L205: “The 100 best compromises of each group of gauging stations were used to produce the 100 best regionalized parameter sets and the HB model was run with these parameters, estimating uncertainty from the standard deviation.”

As well, these sentences will be added after L245: “The mean bias on simulated river flow and potential GWR, computed over the entire period of observation data, and averaged by group of gauging stations, was comprised between -9 mm/month and 5 mm/month (Table 3). The uncertainty computed with the 100 best regionalized parameter sets was ≤ 10 mm/yr for the three simulated variables (Table 4).”

Finally, this sentence will be inserted after L440: “Furthermore, the uncertainty analysis showed that the calibration method allowed reaching acceptable uncertainty levels in the simulated variables and that the mean bias in the simulated river flow and potential GWR was similar to that of the mean bias in the input data.”

Table 3 will be modified as follows:

Table 3: Objective functions for the simulated outputs, for calibration and validation periods, and mean bias over the entire period of measurement

*The presented values are for the stations of which the watershed are completely located in Quebec

Table 4 will be modified as follows to include uncertainty in runoff, AET and potential GWR:

Table 4: Simulated runoff, actual evapotranspiration (AET), and potential groundwater recharge (GWR) and uncertainty for the study area between 1961 and 2017, in mm/year and in percentage for winter, spring, summer and fall, and for the eight watersheds (W1 to W8).

*Part of the watershed is located in the USA - the presented values are for the Quebec part only

Winter: December, January, February; spring: March, April, May; summer: June, July, August; fall: September, October, November

C4: You compare results from this work to results from previous studies – a lot of them spatially-distributed numerical models. Can you provide some idea of how the spatial distribution of GWR compares between them, as opposed to just the mean/variance/ranges? Are you picking up the same high/low recharge patterns? The same temporal distributions? This is hinted in section 5.1 but I Think a more robust comparison would improve the manuscript.

A4: We understand the Referee’s point of view and agree that such analysis would be interesting. However, following comments from the first Referee to make our paper more interesting for a wider audience and less focused on the study area, we will develop the comparison with other research studies on GWR in cold and humid climates. Therefore, we do not wish to elaborate further on the comparison with Quebec-based studies.

Minor comments:

mC1: Line 60-61: It is bizarre that the reference that “acknowledges the lack of representiveness of the daily results….” Is from 2007, yet the reference for the model they are acknowledging is from 2010.

Indeed, the paper from Dripps and Bradburry (2007 – full reference in our manuscript) was the reference for the SWB model. Therefore, the citation of Westenbroek et al., (2010) is not necessary and will be removed from the L61.

mC2: Line 80-81: I think this is needed in a lot of places, not just southern QC!

L80-81 will be rephrased as follows : “There is clearly a need for reliable regional-scale estimates of GWR that can be updated relatively easily using widely available data in cold and humid climates where groundwater is often the main source of drinking water in rural areas (Groupe Agéco, 2019; Kløve et al., 2017).”

mC3: Line 304-305: So it is associated more with the precip trends than with the soil type?

Clayey areas are mainly located in the St. Lawrence Lowlands, the flattest of the study area, that also happen to receive the less rainfall (mainly < 1 000 mm/yr) in the study area. The combination of the three (precipitation distribution, soil type, and flat topography) explains that besides having the smallest potential GWR rates of the study area, clayey areas are also associated with the smallest runoff rates and AET rates. Therefore, we do not think that the text need to be modified.

mC4: Section 4.4; 1^st para: This is a rough paragraph to read. I think a table or graph would be more appropriate.

The first paragraph of section 4.4 L310-321 will be rephrased as follows: “HydroBudget simulated the temporal evolution of the water budget in the study area since the 1960s, thus producing an exceptionally long simulated time series of runoff, AET, and potential GWR for the area (Fig. 8). The effect of interannual variability in precipitation was clear in the simulated runoff with low runoff rates (< 350 mm/yr) produced during the driest year and high runoff rates (> 550 mm/yr) produced during the wettest year (Fig. 8a, b). The simulated AET varied less than runoff, mainly between 450 mm/yr and 560 mm/yr (Fig. 8c). Variations of potential GWR were relatively high, comprised between 90 mm/yr and 200 mm/yr, and seemed more influenced by precipitation than by temperature variations (Fig. 8c).”

mC5: Line 343-344: “somewhat higher” is almost double.

L343-346 will be rephrased as follows: “This value is much higher than that obtained using HB in the same area (109 mm/yr), but the resulting preferential recharge areas located close to the Canada-USA border are similar with both approaches (i.e., 70 mm/yr to 250 mm/yr in Chemingui et al. (2015) and 70 mm/yr to 345 280 mm/yr with HB).

mC6: Line 381-382: awkward sentence structure. Maybe provide range in main sentence text (between 89 and 198) and then average in brackets?

L381-382 will be rephrased into: “The annual variability in potential GWR is closely linked to that of total precipitation, with maximum variation between 89 mm/yr (1968) and 198 mm/yr (1983).”

mC7: Line 387: It is hard to support that this is the novel contribution when you follow it up with the fact that it is supported by a lot of previous literature. I would consider rephrasing this to highlight the specific contribution beyond what the existing literature provides.

L387-390 will be rephrased into: “This spatiotemporal link between GWR and precipitation and temperature patterns is coherent with that reported in other studies, both in similar and different climate and geological environments (Abdollahi et al., 2017; Ashaolu et al., 2020; Chemingui et al., 2015; Fu et al., 2019; Hayashi and Farrow, 2014; Hu et al., 2019). However, one of the contributions of this work was to show that temperature influenced GWR rates at the watershed scale and at the regional scale, both temporally (seasonality) and spatially (increase of GWR from warmer to colder climates).”

mC8: Line 399: This is a regional contribution. While significant to those interested in this region, authors should try to focus on the broader applicability and contributions.

The paragraph L398-407 will be rephrased into: “The increasing trends in simulated winter potential GWR are most likely related to that of winter temperature and precipitation between 1961 and 2017. A contribution of this work is to show that the absence of decreasing trends in potential GWR despite the overall statistically significant increases in temperature and AET indicates that the increase in precipitation (vertical inflows in Table A. 2) is large enough to compensate for the increases in AET. These observations were also true when comparing two decades. Similarly, Levison et al. (2016) and Rivard et al. (2009) did not found significant trend in baseflow during the second half of the 20^th century. Inversely, Nygren et al. (2020) found that groundwater levels significantly decreased in Sweden and Finland between 1980-1989 and 2001-2010 as a result of shifting from cold to temperate climate, inducing a change of the GWR dynamic from snowmelt-dominated to rain-dominated. In comparison, the temperature increase observed in southern Quebec over the 1961-2017 period was probably not large enough to produce such a change. An original contribution of this work was to show that the simulation of GWR with a water budget model allows associating changes in the climatic conditions to long-term trends in the simulated variables and their impacts on the regional hydrological dynamic.”

mC9: Line 439: I generally try to avoid using unquantifiable terms like “good” – how else can you describe this that can be backed up by the results?

The term “good” L439 will be changed for “satisfying”.

mC10: Line 490: ‘work’ instead of ‘word’

The term “word” L490 will be changed for “work”.

mC11: Line 515: Could this have been interpreted from trends in baseflow? Or does this work provide additional info that goes beyond that provided by baseflow results alone?

We do think that this work provides more information on trends in GWR that an analysis of trends in baseflow for several reasons: 1) GWR is simulated continuously over the study area for the 1961-2017 period while baseflows estimated from the river flow time series are only available for periods of time (gaps in the time series and abandonment of gauging stations). Rivard et al (2009 – reference in our manuscript) did not find any consistent trend in Eastern Canada based on baseflow analysis. As well, computing trends on the simulated runoff and AET allows understanding better how the trends in the climate data (increase of precipitation and warming temperature) translated into the regional hydrological dynamic (increase of runoff and AET, but not GWR). 2) River flow measurements during winter are highly uncertain due to ice cover and ice flowing, therefore analysing trends from the raw data only during these periods could be misleading. This point has been added in our proposition of modification of paragraph L398-407.

Supplementary references:

Chapman, T. G.: Comment on “Evaluation of automated techniques for base flow and recession analyses” by R. J. Nathan and T. A. McMahon, 27, 1783–1784, https://doi.org/10.1029/91WR01007, 1991.

Eckhardt, K.: How to construct recursive digital filters for baseflow separation, 19, 507–515, https://doi.org/10.1002/hyp.5675, 2005.

Groupe Agéco: Recherche participative d’alternatives durables pour la gestion de l’eau en milieu agricole dans un contexte de changement climatique (RADEAU1) (participative research for sustainable options in water management in agricole region and within a climate change context), Ministère de l’agriculture, des pêcheries et de l’alimentation, Fonds Vert, Quebec City (Canada), 2019.

Kløve, B., Kvitsand, H. M. L., Pitkänen, T., Gunnarsdottir, M. J., Gaut, S., Gardarsson, S. M., Rossi, P. M., and Miettinen, I.: Overview of groundwater sources and water-supply systems, and associated microbial pollution, in Finland, Norway and Iceland, Hydrogeol J, 25, 1033–1044, https://doi.org/10.1007/s10040-017-1552-x, 2017.

Nygren, M., Giese, M., Kløve, B., Haaf, E., Rossi, P. M., and Barthel, R.: Changes in seasonality of groundwater level fluctuations in a temperate-cold climate transition zone, Journal of Hydrology X, 8, 100062, https://doi.org/10.1016/j.hydroa.2020.100062, 2020.