This paper introduces a novel machine learning algorithm for the simulation of 2D surface flooding. The algorithm takes inspiration from scientific computing techniques to develop the machine learning architecture, to arrive at a setup that provides a time-dynamic and spatially distributed simulation, and a speedup in the order of factor 100 compared to a numerical simulation model. To my knowledge, this is the first application of this kind of methodology for flood problems. A number of questions therefore remain unsolved, and the considered case examples do not yet reflect all complexities that we encounter in the real world. On the other hand, this work is likely to inspire a variety of derivative works in hydrology, where the potential impact of these techniques so far is not really recognized. Because this work is a frontrunner, I have a number of suggestions that are mainly aimed at making the work accessible to a broader audience. I think these can be adressed in a minor revision, and I recommend the paper for publication, subject to these modifications.

Main comments:

1. Limitations - and important feature of numerical methods is that they (mostly) preserve e.g. mass and momentum. As far as I can see, this is not the case for the algorithm proposed here, and it is not straightforward to see how this can be implemented in the encoder/decoder architecture. This should be clearly mentioned as a limitation.

2. A major challenge with this type of model is actually implementation. For example, efficient data pipelines for custom graph operators are not a straightforward problem. I would strongly suggest a section or appendix summarizing the main computational challenges and how you suggest to adress them.

3. Along the same lines, I appreciate that the authors have tried to keep the Methdology description generic. However, this also makes it very hard to read in some places. I think you could greatly help the readers with an Appendix that includes a detailed variant of Figure 2, where you include equations 5 to 10, and where hyperparameters (G, p, etc.) reflect that actual values used in your implementation.

4. Regarding the hyperparameters, I strongly miss a table that summarizes the final set of hyperparameters. These are now spread out through the paper. In addition, the complexity of the different MLPs used throughout the methodology is currently entirely unclear.

5. The benchmark models are introduced in a few lines of text, and not easily understood. Also here, I suggest including an appendix that details these.

Detailed comments:

line 106: a is already used as symbol for area, use v for velocity?

Figure 2:

    -top part of the figure

        -mention what are the static inputs and the outputs in your work in the figure

        -include a recursive arrow from output 1 back to dynamic inputs, to make it clear that the model prediction is recycled

    -bottom part of the figure

        -include the following captions above MLP, GNN and MLP: "process individual nodes", "process neighborhood of each node", "process individual nodes"

        -consider referencing the equation numbers inside the figure, to make it clear which illustration relates to what step

        -include j in one of the orange circles

        -edge feature encodings should receive and output arrow from the MLPs as well. As a whole, maybe this figure would be easier to understand if you distinguish h_si, h_di (why not h_dy?) and eps_ij separately (three sets of arrows)

    -caption

        -h_si and h_di are not clear from the caption (and explained somewhere much later in the main text)

line 138: explain that input features and hyperparameters will be explained in section 3.2

line 145: Ut−p:t are the dynamic node features --> Ut−p:t are the dynamic node features (hydraulic variables) for time steps t-p to t

line 148: define I_epsilon

line 154: explain that G is a hyperparameter

Equation 7: At this point it would be really good to know how terrain differences come into play, and you don't explain it until Section 3.2. I think a short explanation is needed here, because it actually hinders the understanding of the methodology

line 184: is the activation function only applied on the final graph layer?

line 191: Do the neural networks in the graph layers include bias terms? You refer to sparsity in multiple places in the paper, but you never explain it and why it is relevanmt (a large area of the image is 0 and does not need to be processed?)

Eq. 12: it is confusing that u is not the same vector as when introducing the SWE. I have no good suggestion for improving this though.

line 212: define unit normal vector (probably already needed in the context of Fig. 2)

Section 3.3: Consider merging this with 4.2. Searching for the actual parameters used in training is yet another unnecessary complication

Algorithm 1: Nice!

line 266: I don't understand why H is now fixed (previously, you introduced the curriculum). Is this the maximum of H considered?

Section 4.3: I think it would be interesting to see some illustrations of how e.g. MAE evolves over simulation time. This would help us understand better if the method is stable

line 298: I'm not sure what propagation rule you refer to. As mentioned above, these model variations need to be documented better. 

Figure 5: Show units also on the legends of the difference plots

Figure 7+8: From these figures, speedup in the order of factor 100 seems more realistic to me (not 600 as mentioned somewhere in the paper)