Review: Keck et al., 2023 - A landslide runout model for sediment transport, landscape evolution and hazard assessment applications

Keck et al. 2023 present a new cellular automation algorithm as extension to the Landlab framework for calculating runout distances of different types of mass movements. In addition, a phyton script is described to calibrate the model via a Monte Carlo sampling algorithm and statistical evaluation. Finally, several historical events were used to test the model performance. Cellular automata are not the latest approach in describing mass movements, but there are some benefits, such as high computational performance, provided that the methods are implemented efficiently. On the downside, cellular automata also have some shortcomings compared to fluid-dynamic models, which have found their application, particularly in the last two decades, to describe mass movements in both the scientific and engineering-geological fields - more about this later. Keck et al., 2023 identify three application scenarios in their introduction: (1) landscape evolution models, (2) landslide risk assessment, and (3) sediment budgets, although only case two is considered further in the article. The manuscript is quite long with over 1000 lines (including references, without supplement), extremely technical, since primarily the model and the algorithms were described, and for me difficult to follow over long passages.  After a careful evaluation of the article, it is not clear to me for which target audience the manuscript was written. Especially the first half reads like a hard-to-understand manual for tinkering with a cellular automaton.  As part of the Landlab framework, this code will find its users, so an article in Esurf would be justified – but …

From my point of view, however, the manuscript needs many fundamental changes and I recommend that the manuscript should be circulated among all authors several more times during the revision process before thinking about resubmission. The mathematical description is very buggy and must be rigorously reviewed by all authors, corrected, and distinctly simplified. It is required that the model is tested before applied to real world examples. Use cases that are addressed in the introduction must at least be explained, if not shown (landscape evolution, sediment budgets). Claims of better computational performance compared to fluid dynamic models need to be proven or weakened as an argument.

In its current form, I can't imagine the readership benefiting from this study. In the following, I explain what I consider to be the biggest problems. In addition to major issues, there is ambiguous wording in the text (especially where other studies are referenced (e.g. Lines 123 and following, 616 and following)) and spelling errors in the pseudo codes. Since I expect the manuscript to be fully revised as part of the internal revision cycle, I will omit the line-by-line comments.

Although very critical, I hope my review will be received constructively and will help the authors to improve this work.

Mathematical description, algorithms, variables, indices

It took me several attempts to work my way through sections 2.2 and following. I stopped reading several times in frustration because I got lost in all the variables and indices. Unfortunately, the textual description was not very helpful - and this although cellular automata are not a rocket science. In a first step I recommend the authors to make a table with ALL variables and indices (including a short description and units). I can't for the life of me imagine that anyone reading the article can keep track of all the variables that have been introduced.

In part, even the authors seem to have lost control of it. Here are just two examples: The variable \theta stands for local slope and appears for the first time in equation 12, but is explained much later after equation 14. \theta, however, appears in equation 25 in a completely different context. Here \theta stands "for a set of model parameters".  However, (local) slope already appears in equation 4, but with the variable S_{n}, which should also correspond to \theta in equation 12. However, this is not quite clear, because slope can be the topographic slope or the slope of topography + flow height of the landslide. I find it quite confusing that the variable n stands for "node index of debritons" - especially for programmers the n and m often stand for the number of rows and columns of a 2-dimensional array, or at least for the number of elements of a list. Anyway, in line 438 the variable n is then recycled and stands for the number of topographic changes from then on.

Some anomalies make me doubt that the derivation of the equations is sound. I'm only giving two examples here, but the list can be extended at will.

In equation 24, here the flow height is artificially limited - which is wild in itself, the variable h' is introduced as the limited flow height. The following line tells “Note that h' is used only in equations (8) and (11).”  How can I use a variable in (8) and (11) that is not defined until (24).  Beyond that, in equations (8) and (11) there is neither h nor h' but in line 261 after equation (12) it says that h is equivalent to qs^{I}_{n}. According to equation 2 this should be the total incoming flux. A material flow is not a flow height and the above derivation cannot be understood.

In equation 4, the slope of the cell S_{n} is compared with S_{c}*delta x, where S_{c} is a critical slope and delta x is the cell size. From my understanding slope is dimensionless [ ] but the cell size has the unit [m]. How can a dimensionless property be compared with a property that has the units [m]?

The indices gave me a big headache: A_{p,n|N_{a}} or qs_{Dnj,n} may sound logical to the authors who spent months on the paper - but the paper should be readable in one pass by people who were not involved in the study.  The indices have not only subscripts but also superscripts (which could easily be confused with exponents). If one of these indices then erroneously slips into the wrong position (compare lines 191 and 193), the matter becomes even more confusing. These issues make it extremely difficult to evaluate the model from the theoretical side. My recommendation would be to use the variables consistently throughout the manuscript, reduce the number of variables significantly, and simplify the indexes.

A simplification of the model description could be achieved by describing the theoretical basics precisely (both in equation and with text) but refraining from the description of the algorithms implemented in the code. This also applies to the calibration routine (Monte Carlo Sampling). Thus, one could get rid of a multiplicity of variables and the complicated indexing. The description of the algorithms with pseudo-codes is not absolutely necessary for the understanding of the model from my point of view and could possibly be transferred to the supplement. However, I think that scientists who are interested in the algorithms are more likely to read the Python code directly instead of pseudo-codes and long descriptions. 

Since the current manuscript follows the algorithms of the code in structure, flow routing, erosion, and deposition are treated in parallel. For the reader's understanding, it would be easier to untangle these strands and first describe how flow routing works and then describing what the rules are for deposition (without node lists und indices – they are not required). Erosion and entrainment could then be treated separately in their own section. I don't think you need to go into the individual nodes or node lists here - as it is handled in the code. Because even if the algorithms are correctly described with all indices (I can't evaluate this) and pseudo-codes would correctly represent the flow of the model, one cannot be sure as a user whether the methods were implemented correctly in the code. In the end, however, this is exactly what the majority of potential users need - a code that behaves as one would expect from the previously described theoretical basics.

Testing the model

In the best case, this model is used for the development of hazard zone plans, which control the development of settlements and infrastructure in alpine areas. Would I trust a hazard zone plan based on this model with the responsibility for a settlement with hundreds of inhabitants?  There is no margin for errors - I would have my doubts. Starting at line 449, we find the section "Model Validation", but here the new model is already used for real world examples. It is good to show a real-world application, but instead of modeling six landslides and also describing them in detail, one landslide for calibration and one next to it for testing the calibration would suffice. The work on the other landslides would not be lost, because they could still be published in a case study or something similar in an article on its own. The model does show partial good agreement with historical events after calibration, but in my view a new model needs a much more rigorous testing regime where flow routing, deposition, and erosion characteristics can be evaluated under controlled conditions.

To show how the model works with its rules I would recommend experiments with simple synthetic digital elevation models.

• A ramp with different slopes to see how a mass movement laterally diverges due to flow routing  – please compare with models applying a depth averaged version of the Navier Stokes equations.
• A ramp that turns into a plane with a kink to see what geometry the deposit forms. In addition to the plane, you could also take a ramp that is less steep than the critical angle for deposition and a ramp with counter rise.
• A ramp that turns into a steeper ramp (convex transition). Here it would be interesting to see how the flow depth changes. In nature (and in fluid dynamic models) the landslide mass accelerates at the transition to the steeper region and the flow depth decreases. I doubt that this can be represented in this model. However, in the model presented, the "erosion rate" is primarily controlled by flow depth.

Use Cases

In the introduction three different applications of the model are mentioned: (1) landscape evolution models, (2) landslide risk assessment, and (3) sediment budgets. That sounds very promising at first and would make the model the Leatherman of landslide research. However, it is not elaborated how far this model should work as part of a landscape evolution model. Over which time scales is it planned to operate? Is the model numerically efficient enough to be part of a LEM? The same applies to sediment budgets. How would the application of the model contribute to sediment budgeting?

Computational efficiency and parameterization complexity

In the manuscript, a number of claims are made regarding the computational performance and parameterization complexity of cellular automata compared to fluid dynamic models. Due to the simple structure of cellular automata, these models should be numerically more efficient than models from the field of computational fluid dynamics. From the description of the algorithms, it is difficult to estimate how efficient the code really is. Repeated sorting of large lists can be a performance brake, as can the handling of special cases (e.g. flow routing). The downslope motion occurs iteratively across cells, which is not fundamentally different from a forward-in-time calculation with time steps in fluid dynamic models. The examples presented in the manuscript are relatively small landslides and the chosen resolution of 10 m is quite coarse, so the number of cells of the presented models is small. Wouldn't it make sense to exploit the power of the model and calculate with a spatial resolution of 1 m, so that even smaller but quite relevant features of the topography are correctly represented? Over the last decade I performed countless landslide simulations on standard hardware with fluid dynamic models. The computational performance of fluid dynamic models is now superior (e.g. adaptive mesh refinement, parallelization, efficient equation solver) and easily allows >1000 model runs for an event, but at much higher resolution (1 m) and on a larger scale than presented in the study, where the modelled landslides tend to be small (< 1 million m³).  In summary - the claim that the presented code is more efficient than fluid dynamic models could be easily tested against one of the many freely available landslide models based on computational fluid dynamics. Results should be presented here as well, or the authors should weaken the claim of the better performance a bit.

While the claim of better computing power is theoretically understandable, I have to disagree about the claim of simpler parameterization. Since fluid dynamic models are mass AND momentum conserving, the parameterization is much easier than for rule-based cellular automata from my experience. In the presented model you even have to think about artificially limiting the flow height to avoid unreasonably high flow heights (Line 311, Eq. 24), or you have to derive the flow velocity via a series of empirical laws with numerous assumptions (steady state velocity) in order to calculate the erosion rate (Line 274, Eq. 16). In the minimal version of a landslide model assuming a Voellmy flow resistance law (without entrainment) it needs two parameters that have been calibrated for a wide range of mass movements and work well for debris flows, snow avalanche and even major landslides.

Closing remark

I would be grateful if the authors would take the suggestions and especially make the mathematical description readable for a person who is not involved in the study. In addition to solving the major problems, it would also be good if the experienced authors of this study in particular would invest a few evenings and would check line by line for consistency and question whether the formulations in the text are generally understandable. I don't see any insurmountable obstacles to fixing the problems but I'm afraid there is still a lot of work to be done to get the study ready for publication -good luck!

Keck et al., present a novel cellular automaton landslide runout algorithm which can be implemented in the Landlab modelling package. Keck et al., construct their model using rules primarily on sediment flux, and slope angle to determine runout extent, plus erosion & deposition along the path. They calibrate and validate the model using six interesting and varied case studies and a novel Markov Chain Monte Carlo algorithm to sample model parameter space. The paper would be of interest to the ESurf readership if substantial improvements are made, as the cellular automata model within the paper represents an important addition to the Landlab modelling package and would be of interest and use for varied geomorphic applications.

I present below more detailed general on suggested improvements to the manuscript and line specific comments.

General Comments:

The description of the model and resulting equations is very technical, and easy to get lost and confused in the indices and variables described. As it stands, it hard to follow and assess the validity and theoretical basis of all the equation components. The methodology would benefit from a higher-level description of the method with the base equations outlined, that are explained in a logical order and with plain wording. The introduction contains an overview of controls on landslide erosion and deposition, but it would be beneficial to explicitly justify the choice of model rules within this higher-level description of the method. This would allow the reader to much more quickly understand the important rules and variables that control modelled runout. While I understand why the notation and indexing was included in the formulas, as these are important details for how the model is coded and works, I believe this level of detail may be best included in a supplementary method that mimics the higher-level base equations outlined in the main paper. A table/glossary of indices and variables would be invaluable. A list of clear inputs to the models would also be advantageous. The figures throughout the text are well illustrated and clear and could be better utilized to explain concepts in text. Similar comments apply to the model calibration methodological section, where multiple indices are quickly introduced and again it’s easy to get lost.

Overall, the paper is long and could do with re-focusing on the important points that the articles want to make. Within the discussion section, the model validation section could be streamlined, and more time spent discussing the application of this model to landscape evolution models, sediment budget models and landslide hazard and risk analysis. Additionally, other cellular automata models are listed in the introduction – how does this model differ from those models, both in terms of rules and model performance?

The section on topographic controls is an interesting part of the paper but comes a bit out of nowhere and is not really set up in in the introduction or talked about in the rest of the paper apart from the conclusion. The relationships outlined in Figure 9 are weak and would benefit from a larger inventory dataset to derive the relationships between model performance and topographic indices. This may or may not be included in this paper, depending on the need to streamline and shorten the paper.

Line Comments:

Line 181: Describe how local deposition and erosion are calculated before Equation 1. Equation 1 would be better as one of the last equations.

Line 193: Typo in the notation of qso_n

Line 206: Can you initially explain why you chose Carretier et al., (2016) and Campforts et al., (2020) nonlinear, nonlocal deposition scheme?

Line 208: Would be good to define delta X for the non-mathematically minded, and also include a plain language description of what you are trying to achieve?

Line 214: What is meant by the single row of downstream nodes?

Line 221: What assumption (2 or 3) are you referring to?

Line 224: Bit lost in equation 6

Line 242: Define variable k here, not in three paragraphs time.

Line 273: What is the reasoning for using Julien and Paris (2010)?

Figure 4: If delta X is set by the resolution of grid cell, what is the impact of ground model resolution on runout model outputs

Line 342: What is meant by dividing the squared, cumulative, volumetric square error (SE) by the observed model volume? Instead of using simulated volume vs. observed volume? This needs to be explained more clearly.

Line 350: What do the variables of Δƞ_Mi and Δƞ_Oi exactly represent?

Line 431: Mapping landslide runout risk – should be something like ‘Mapping landslide runout hazard’, as you are not describing risk in this paragraph.

Line 563: What does Ē/ƪ represent? Would be good to clearly state this.

Line 636: Incomplete sentence.