The manuscript presents a much-needed methodology for cross-validation of spatial data. In my opinion, the strongest point is the use of the W statistic to identify the best CV split. However, there are a few points which I feel should be addressed in the discussion.

1. The proposed methodology using clustering algorithms seems valid, but how can we know if it provides the best possible result? An algorithm that optimizes the W statistic directly as a function of the CV fold indices would be more desirable, instead of relying on the clustering algorithm´s internal metric as a proxy. As a suggestion for future work, I recommend using a genetic algorithm to assign CV indices to the data points directly.
2. The W statistic explained 60% of the variability in map accuracy, but would this be consistent across different datasets? At least one more case study would be needed to verify this.

Minor comment:

Line 90: cross “the”.

The study proposes a novel cross-validation method for spatial data that aims to deliver more representative measurements of spatial map accuracy than commonly-used methods. This is a relevant concern for GMD readers with the rise in use of machine learning methods for geoscientific modelling. Issues with model evaluation in the spatial setting have been identified in a number of recent studies. The paper is well-written and contributes a practical solution for a common issue.

In my opinion, the most exciting/innovative idea in this work is the concept of defining the evaluation method based on the desired data for which he model is intended to return predictions. This would require researchers to more carefully define the purpose of their models before and during the model creation process, which should be common practice. In reality, this is often not done, or done in a ‘standard’ way which doesn’t accurately reflect the intended use of the model.

The method presented in this paper is a very practical solution to this, where the desired target dataset is an input of the evaluation algorithm and therefore researchers are required to clearly consider and define it. I think this is a significant contribution to model development methodology and should be more clearly emphasised in the manuscript. The possibilities, benefits and disadvantages of this concept could also be discussed - for example, when models are used in production, the prediction area is a moving target; would that require continual re-evaluation?

The paper suggests that kNNDM is, essentially, a computationally-cheaper alternative to the previously-published method by the authors, NNDM LOO. In the article, the only limitation of leave-one-out CV methods described is that of computational time. However, to my knowledge, even if computation is not considered, LOO CV methods may not be the optimum method due to higher variation in the resulting models (due to the bias-variance tradeoff). Could this explain why kNNDM 10-fold seems to perform better in the case of strong clusters (Figure 5)? For me, this would be more convincing than the computation speedup comparison, which is relatively trivial given that LOO CV is the most extreme version of k-fold CV.

Following on from this, it seems likely that the value of k would impact the results. Use of 10 folds is very common; is there theoretical justification for this? It would be useful to see some comparisons of the results with multiple values of k.

In Figure 1, it is shown that the W statistic will also be larger if training points are regularly distributed, as well as when clustered. Does this mean that the null hypothesis might be rejected for regularly distributed datapoints? Does this explain why NNDM LOO performed better for regularly distributed data (Figure 5)?

Finally, I would recommend testing the method on at least one additional dataset, as the results presumably depend on the spatial autocorrelation present in the dataset used.

Minor comment: I assume the hyperparameters of the models are not tuned as it is not mentioned, but this could be stated explicitly.