Dear Marvin Höge,

We appreciate your suggestions on our manuscript, which we indeed consider all valuable in enhancing its quality. In response to your comments, we aim to address your discussion points as following.

In the introduction, it is mentioned that some consider machine learning methods as “black boxes”. Explainable artificial intelligence (xAI) tools or similarly called methods are supposed to help here. Therefore, a brief coverage of advances in this field – ideally in the field of groundwater research if available - would be worth mentioning.

We acknowledge the importance of addressing the implementation of Explainable AI techniques in our research, and we intend to integrate a discussion on this topic into our manuscript. This addition will contribute to a more thorough understanding of the interpretability and transparency of our modeling approach.

The main point of concern, however, is as follows: Overall, the performance of the employed CNN model, as presented, e.g., in Figure 6, is not fully convincing. Even the well-performing models - according to NSE and R2 - show mainly a sinusoidal pattern with only slight variations – yet, it is these variations that would be interesting to be modelled. Otherwise, a sine function-based model with a mean trend might often be sufficient and provide the same goodness-of-fit values. Therefore, one can assume that the used model architecture (together with only precipitation and temperature as inputs) is not complex enough to capture more of the dynamics. The subsequent analysis of performance with respect to geospatial and time series features therefore appears to be weaker than it could be. It appears to be difficult to deduce relations between features and model performance if the model does not perform convincingly in the first place. All correlations reported are rather weak with the strongest anti-correlation being -0.62.

Regarding your main concern about the model's performances, we hypothesize that the monthly resolution used in this study may contribute to the sinusoidal pattern seen on the well-performing model example. The seasonality evidenced on the monthly temperature and precipitation used as inputs can certainly affect the model behaviour. It is worth noting that the CNN model has previously demonstrated effectiveness in predicting weekly groundwater level time series with high overall accuracy, as shown by Wunsch et al. (2022) "Deep learning shows declining groundwater levels in Germany until 2100 due to climate change". Therefore, we are confident that our approach works well for higher temporal resolutions and we propose a shift to a weekly temporal resolution for modeling the time series. This adjustment is intended to mitigate the potential seasonality introduced by the monthly resolution. By comparing the model's performance at weekly resolution, we seek to address this issue and establish a more robust foundation for examining geospatial and time series correlations also in form of possibly higher correlations.

It is important to emphasize that while adopting a different model may yield improvements in overall performance, we have consistently applied the same model to all monitoring stations. Consequently, local or more specific variations in performance metrics across stations are primarily attributed to external factors rather than model selection. Certain models may better adapt to specific locations while displaying lower performance when applied to the entire dataset. In this study, our primary objective is to analyze these relative performance differences in the presence of external influences. Therefore, even if an alternative model were to enhance overall performance, our analysis will retain its current form, focusing on the examination of these external influences.

Along these lines, the used geospatial features are all interesting but the timeseries features are too many and some are hardly tangible. A thorough explanation of their meaning, range of values, etc. would be beneficial. If time series features shall be part of the analysis, I recommend to focus on only a few of the ones provided. Yet, in this case, the question remains: What is exactly gained from relating these time series features to model performance? For instance, a time series can be rated as complex while an adequate model could still be able to predict it.

The selection of time series features represents the outcome of a pre-evaluation aimed at elucidating their physical significance in the context of groundwater level time series. We agree with your suggestion that a more extensive interpretation of these features should be incorporated into our manuscript. We will augment our discussion with a detailed explanation of each feature, emphasizing their relevance to the physical aspects of groundwater level dynamics.

Regarding the detailed comments addressed on the attached PDF, we will go through them and incorporate them when in agreement.

In conclusion, we are committed to addressing your comments as comprehensively as possible.

Dear Jonathan Frame,

Thank you very much for your comments on the manuscript; we value your incomes and would like to incorporate the suggestions when suitable and in accordance with the current objective of the paper.

“this paper would greatly benefit from additional description of the training procedure and evaluation when dealing with gap-filled or processed data. Another benefit would be further evaluation of the sensitivity of model performance to the gap filling and data processing measures.”

We will certainly elaborate more on the data pre-processing, by including a more detailed description of the data exploration and gap-filling methods used before applying the CNN model. Regarding the training procedure and sensitivity of model performance to the gap filling, we believe that it can be helpful as further and future research but by only using time series with good data quality in terms of data gap lengths and frequency we seek to avoid major influence of data imputation approaches. Therefore, we think that a sensitivity analysis is a bit out of the scope of this study and could be focussed on in a follow-up analysis.

Line 44: These claim should be cited: “In terms of accuracy and calculation speed, the CNN models outperform the LSTM. NARX models performed, on average, better than CNN.”

We agree that this needs to be cited.

Line 47: “Most studies have successfully applied these techniques for GWL forecasting using only meteorological variables as inputs.” You might be interested in this paper: Gholizadeh et al., “Long short-term memory models to quantify long-term evolution of streamflow discharge and groundwater depth in Alabama” Science of The Total Environment Volume 901, 25 November 2023, 165884, where the did in fact include site geospatial characteristics to make predictions of wells that were held out from the training set (ungauged).

Gholizadeh et al. 2023 applied an LSTM model including static inputs that refer to the aquifers' hydrogeology as an attempt to model ungauged locations. As the model does not include groundwater levels as inputs, the authors attribute the satisfactory model performance to input features such as hydraulic conductivity, soil depth, soil porosity, and maximum water content. These findings can contribute to the central discourse of the paper.

Line 116: Can you please provide the total number (and percentage) of gap filled values referred to here: “To provide the CNN model with continuous time series, we performed a data imputation process through a Multiple Linear Regression (if enough dynamically similar wells based on the Euclidean Distance”. Can you also explain if these values were removed, or should be removed, from the loss during training, and also removed from the evaluation?

Thank you for raising this point. We will include precise numbers regarding the data gaps and imputation. From the 505 groundwater level time series, 241 (48%) are complete, 254 (50%) have gaps of 2 consecutive values and 10 (2%) have gaps of 3 consecutive values. Overall, the time series have less than 5% gap-filled values. We did not remove them from the training phase since the number of filled values is not considerably high.

Line 120: Similarly can you provide the total number (and percent) of data points that were modified as outliers described here: “We removed these anomalies by finding the highest slope in the cumulative sum”? Is this a standard approach? I don’t think this description is satisfactory. I see from your code that you identify these based on “initial point where the values increase by 0.5 of the standard deviation“. This is an important point that should be explained in the paper, as well as the decision to use this processing method.

Only 28 wells were identified to have jumps/steps on the temporal record. The cumulative sum is commonly used to detect changes in the mean or variance along the time series and is not referred to as outliers. Here, we intended to detect jumps/steps on the observed values that can hinder the model training or might introduce confusion to the model due to potential changes in the dynamic of the ground water levels. The optimal fraction of standard deviation was determined through trial and error by visually inspecting the detections and selecting the value that best adjust to most of the jumps. We will definitely include these explanations in the revised manuscript. Thank you very much for pointing this out.

Line 135: Is it really necessary to give the equations for r-squared and NSE, as you don’t provide the equations for MSE or BIAS? There is also more unfamiliar calculations made in Tables 2 and 3 with no equations provided, and also the main CNN model is not described with equations. I guess I would suggest just removing equations 1 and 2, avoiding an asymmetry in descriptions, rather than adding equations for all the rest of the calculations.

We agree to remove the equations to make the manuscript more consistent regarding the inclusion of equations.

Line 176: “Occasionally, in poorly performing models, the pattern of the GWL observations has been generally learned but with a strong Bias.” This is a little concerning, and I think it would be work describing in more detail. Similar to your NSE/r-squared cutoffs above, can you provide a quantification of these problematic BIAS wells, something like in how many wells does the prediction not intersect the observation?  What causes this BIAS, is it an unusually high section in the training period? I wonder if there is anything in the data preprocessing that plays into this issue

This comment relates a lot to the main concern, raised in the first reviewer comment (RC1). We will address this issue more in detail through re-running the model on a weekly resolution, which we expect will improve model performance.