This manuscript presents a novel approach for modeling gravel-bed river evolution, specifically the employment of sediment deposition process, within the context of Landscape Evolution Models (LEMs). The authors introduce the RiverBedDynamics component within the Landlab platform to improve the accuracy and applicability of LEMs for gravel-bed rivers. By integrating sediment transport dynamics, grain size evolution, and non-steady flow conditions, the model addresses key aspects of gravel-bed rivers in existing LEMs, particularly in simulating the complex interactions between river channels and surrounding landscapes.

The manuscript is well-structured, and the motivation for developing the new component is clearly articulated, supported by existing LEM frameworks and their limitations. The authors demonstrate how the RiverBedDynamics component can overcome these limitations by incorporating sediment transport dynamics, non-steady flow conditions, GSD evolution, and predicting bed surface changes across entire watersheds. Although, the manuscript can be accepted in its current form, some additional clarifications might be helpful for the readers to have a thorough understanding of the work. 

The authors presented a simple approach to employ deposition process based on the mass conservation at each cell. Some previous works have also employed the sediment deposition using different frameworks such as size-dependent sediment deposition probability in a Lagrangian framework in CIDRE (Carretier et al. 2023, https://gmd.copernicus.org/articles/16/6741/2023/), and Eularian framework in River.lab (Davy et al, 2017 https://doi.org/10.1002/2016JF004156), and  Minor et al. 2022 (https://doi.org/10.1029/2021JF006546). These previous works employing deposition process deserves a credit and need to be discussed how their sediment deposition framework complement or differs from this this work. The uniqueness of RiverBed Dynamics lies in its easy coupling with any flow solver, incorporation within the Landlab framework, and prediction of temporal changes in GSD (although GSD doesn’t change much during the simulation period, Figure 9d). In another work, Lei et al (https://doi.org/10.1029/2023WR035983) studied the change in GSD in response to different sediment supply. The authors can try to alter the sediment supply rates and supply GSD to study the sediment sorting (not necessary for this review as this is a model description paper, just a suggestion) or atleast discuss in the paper about the possible causes of minimal changes in the GSD.

Line 125: key aspects can be mentioned as the bracketed term.

Line 158-159: Rewrite the sentence

Figure 2: The fonts in Figure 2 are too small to read. The Figure with enlarged fonts in landscape format (rotating 90^O) might improve the readability.

Line 196-200: Can you please briefly mention what variables have been calculated at nodes and links specifically.

  Line 208-209: Why authors prefer central difference scheme as it calculates the velocity averaged over three links (x-1, x, x+1). Instead the velocity calculated between two links as . The velocity calculated between two links might be better representative of local channel hydraulics, (i.e., flow resistance). However, the form adopted might be better for numerical stability. Justify the choice in the manuscript.

Line 238-247: Again, please clarify what parameters being calculated at nodes and links.

Equation 15: , Authors can clarify the role of F_i in the numerator as transport of i^th size can be estimated as transport of i^th size divided by the transport of all size classes without F_i.

Line 335-336: The active layer depth may change with the strength of flow (discharge, and depth), especially the unsteady flow. The assumption of constant active layer might intercept the interaction between bedload and the substrate layer, and might not be the better representative of the actual process in the field. The authors should justify the limitations of the approach and what implication it might have on the bed evolution and grain size distribution.

Line 360-370: Does the two stratigraphic layers may have different GSD depending on the time evolution of bedload? How does the stratigraphic deposition improve the process representation? What difference does it make with a single deposition layer instead of a deposition layer of a constant depth of 1 m. Depending on the discussion, it seems all the stratigraphic layers might possess merely the same GSD. Please clarify.

Table 1: Option surface_water__velocity_prev_time_link was not explained or mentioned anywhere else in the paper.

Figure 9d: Its difficult to see the differences in GSD at different times. The authors can try to reduce the linewidth and reduce the frequency of plots and see if it improves the readability of the Figure (i.e. plotting only o, 8.1,16.2, and 93.9 days). Also Figures 9c and 9d shows minimal changes in GSD after the 60 days. Does this relate to the sediment sorting/armouring being developed after 60 days? The authors can add another subplots showing the vertical stratigraphic distribution of geometric mean size (at the same location x=1000m) and the GSD at different stratigraphic layers. Or Does it relates to continuous sediment supply with same GSD as that of bed. Also the authors can discuss the possible implications of altering the sediment supply rates and GSD on bed evolution.

How computational intensive is the application of RiverBedDynamics to real-world problems. Briefly mention about the computation times for different simulations in channel and synthetic watershed, and how feasible is it to apply it to a flood event in a river both spatially and temporally with computational requirements.

Line 672-677: The procedure of adjusting the water depth after the erosion or deposition is not very clear. How running the model for few internal cycles can give the better water surface elevations over the eroding or depositing cells.

Since many variables and their short representation have been used in the paper and its often confusing what variables are being referred to. It may be better, if the authors can provide the list of variables used and their full names at the end of manuscript.

In this manuscript, the authors present a new component, RiverBedDynamics, for the Landlab landscape evolution modelling framework. RiverBedDynamics facilitates modelling of transport-limited fluvial processes (i.e. deposition, erosion and transport of sediment), including grain-size dynamics. In combination with existing Landlab components, RiverBedDynamics can handle non-stready flow conditions. The component features several sediment-transport parameterizations, and results are benchmarked against analytical or earlier numerical solutions where possible. The authors include a brief example which nicely illustrates interactions between bed-elevation and the hydrogaph during a flood event. I did not try to run the code myself but it appears to be well archived and documented.

The authors make a strong case for the novelty of the component’s capabilities, which do represent a significant advance on existing landscape evolution models. The various elements are comprehensively described, and the examples do a nice job of illustrating them. For these reasons, I think the manuscript is appropriate for publication in GMD. I do have some comments/suggestions, which I outline below, first more generally, then line-by-line. Many of my comments relate to clarifying the authors’ motivations and intentions in developing the component – I hope this will help improve the accessibility and impact of the manuscript for as broad an audience as possible. The remaining comments are mostly related to (minor) technical aspects I found confusing or incomplete. With these changes, I think the manuscript would make an excellent contribution.

Specific comments

Anticipated spatial and temporal scales of use

I was unsure, reading the manuscript for the first time, on what spatial and temporal scales and/or resolutions the authors anticipate RiverBedDynamics being applied. Later in the manuscript, the authors do mention that their “mechanistic approach may be better suited for modeling relatively small-time scale processes”, that various adjustments/enhancements might be needed “to expand the component's applicability to longer timescales”, but also that the component is “flexible enough for short- and long-term simulations”. There is some inconsistency here, but my general impression is that the authors have aimed to design a component that can work across time scales, but is probably most suited to short time scales for now. I think this would be good to clarify much earlier in the manuscript, since it influences some of the assumptions made in constructing the model and also the types of problem that can be addressed using it. Not doing so lead to some confusion on my part which I will give some examples of below.

In the discussion of how to parameterize bed shear stress, emphasis is placed on a scenario where the spatial resolution is much finer than individual channels. Elsewhere, they mention potentially accounting for the effects of things like boulders in the future. Both of these examples point to a focus on very short, process-level spatial scales. Furthermore, the emphasis on accounting for non-steady water flow implies a focus on short timescales. The example simulation is indeed carried out on the timescale of individual flood events. However, the spatial resolution of 30x30 m is well above the process scale, and I imagine much larger than most of the channels in such a small domain. The rationale for this choice is not explained, which made we wonder if there are other limitations (e.g. computation time, numerical stability) that actually make running the model on shorter spatial scales impractical. So, while several lines of argument suggest a focus on short spatial and temporal scales, the example shown seems to depart from that in the spatial domain.

A related issue is that, if the focus is on short timescales of individual flood events, the wider motivation for the component becomes less clear to me, since, on these timescales, the drainage planform does not really evolve. Why, then, do these computations need to be carried out in a full LEM across two dimensions? (As opposed to, for example, the network based approach of Pfeiffer et al., 2020, which the authors refer to.) The authors refer at various points to the need to allow for interactions between different landscape processes, but, on these short timescales, I’m not sure what processes they are thinking about. As far as I can tell, they assume that sediment transport by flowing water applies across the entire landscape – if other processes (e.g. soil creep) can be incorporated with other components, it is not mentioned.

A lot of this confusion would be cleared up if the authors discuss early on that they are aiming to build a component that works across spatial and temporal scales, with a focus on shorter timescales to begin with (if that is indeed their aim). It would be good to also to clarify (a) the choice of spatial resolution in the example and (b) the need for a two-dimensional approach on short timescales.

Potential applications

The authors clearly explain some limitations of earlier landscape evolution models, and how their component is a more realistic treatment of process in gravel-bed rivers. But they give few (if any) examples of the kinds of problems they anticipate their component being used to address. The reader is therefore left wondering somewhat why the component is needed (aside from whether or not it is a “major step forward” in a technical sense). It would be great if the authors could give some examples of the kinds of problems that cannot be properly addressed with previous models, but can with this one. I think this would help engage the reader, increase the impact of the paper and inspire future work.

Introduction

I found the Introduction somewhat long-winded and unfocused. Especially at the beginning, there is quite a high level of detail about features of landscape evolution modelling that are not the main focus here. For example, the first paragraph contains quite a lot of detail about different approaches to flow routing, while the second discusses discretization, neither of which are central to the paper. Only in the third paragraph do we start to hear about the specific knowledge gaps that will be addressed by the new component. I suggest restructuring, and potentially shortening, the Introduction, to place more emphasis on the specific issues to be addressed and why there are important. This would also be a good place to clarify already the spatial and temporal scales considered and give some examples of potential applications.

Line-by-line comments

L72-74. The authors focus on correctly modelling bed shear stresses as the motivation for modelling the evolution of grain-size distributions. Another motivation might be understanding fluvial deposits (indeed, later on we learn that the component can keep track of the developing stratigraphy).

L108-110. Here the authors refer to the importance of modelling elevation changes across the entire catchment, not just within channels, and broader landscape interactions. Could some examples of these other processes and interactions be provided? The sediment-transport physics used in the component applies to flowing water, which might not apply outside channels – are other components needed to included other processes?

L113. “What’s needed” – colloquial.

L156. “Let’s examine” – colloquial. I am not necessarily opposed to this conversational style, but it is hardly used elsewhere in the manuscript, so it stands out here.

L160-167. This is the only place in the manuscript that snippets of code are given – in isolation, without more context as to the structure of Landlab components, I don’t think it’s particularly useful. Consider removing.

L182. “In our implementation…” I was a bit confused by this sentence. The section is about the component, so I initially took “our implementation” to refer to the component. But then the sentence seems contradictory, since the implementation both utilizes a specific flow-router and can be used with any flow router. Maybe with “our implementation” you are referring to the wider model you construct and use later in the examples, distinct from the individual component? Consider rephrasing for clarity.

L186. “Bed surface grain size properties variables” is a bit of a mouthful – consider either “properties” or “variables”.

L200. The term “unsteady friction slope” is new to me, could a definition be provided?

L222-236. As I mentioned above, I think the discussion of water depth vs. hydraulic radius for calculating shear stress has implicit assumptions about the scales on which the model will be run. It is good that both options are available, but the emphasis on water depth is justified by a scenario in which the spatial resolution is much finer than that of individual channels (which is not carried out in the later example, where the resolution is 30x30 m). The authors also mention that the use of water depth makes more sense for application across the entire landscape, where flow is not always channelized. This seems to suggest an assumption that significant sediment transport on hillslopes (outside channels) is driven by flowing surface run off, whereas we normally think about processes like soil creep in LEMs. If that is the assumption, that would be good to state/discuss explicitly. Can the other components be brought in to account for different/additional hillslope processes?

Figure 4. The size differences between the in/out arrows are not immediately obvious visually. An alternative, to make the figure more intuitive, might be to use 2D arrows where both the width and length scale with the sediment flux.

L312. Could the authors explain their rationale for using an explicit scheme here? In the flow routing, an implicit scheme is used. I assume this choice has implications for the stability of the solution (some later statements related to the example imply that too). It might be useful to give an approximate stability criterion, if one exists, to give readers an idea roughly what combinations of spatial and temporal resolution will be practical when running the model.

L320-321. I was surprised that boundary conditions are only needed at the outlet – I am used to solving Exner-type equations in one spatial dimension (i.e. longitudinal profiles), where we definitely do need boundary conditions at both ends of the domain. Is there an implicit assumption here, e.g. no sediment flux (zero gradient) over the edges of the watershed? If so, that might be perfectly reasonable, but it would be good to state it explicitly.

L321. The default boundary condition for the outlet is “zero gradient”. Does this imply zero sediment flux out of the watershed? That seems strange to me if so, and the fixed elevation option (base level) would be more familiar. Could a brief description of the physical meaning of the zero gradient option be included? Also, I don’t follow the subsequent description of the implementation of the zero gradient option. Please reconsider the phrasing there.

L331. I was initially unsure why the different sediment layers needed to be defined and tracked. It becomes clear later (e.g. L352-352), but only after going through quite a lot of details. Perhaps the start of this paragraph can be rephrased to make the motivation/necessity clear from the beginning.

L359. The stratigraphy tracking is interesting and potentially useful in geological applications – it could be mentioned earlier (e.g. in the Introduction) as part of the general motivation behind the component.

L410. The previous section was also number four.

L435-437. Can you comment more on the stability issue? Is there a predictable stability criterion (e.g. related to the explicit solution to the elevation equation)? Or does the instability come more from the coupling of the two equations, and need to be found by trial and error? In either case, this would be useful practical information for potential users.

L504-506. I am curious about the physical reason behind the (quite significant) differences in equilibrium slope for the various models. Could you briefly comment on that? In general, we don’t learn much about the differences between the parameterizations (other than that they deal with grain size in different ways). I appreciate it is not the purpose of the paper to review these models, but some basic information might be useful (also above where the models are introduced), to aid readers in choosing an appropriate parameterization for their purposes.

L514-518. Similarly, it is interesting to see the downstream part fining and then coarsening – could you briefly comment on the physical reason for that?

L555-556. I would expect the 30x30 m resolution to be much larger than most (all?) channels in a 6x6 km watershed. Is water depth or hydraulic radius used for the shear stress calculation? This might be a useful example to help readers think through which of the two options is appropriate given the resolution/scale of a given problem. Are there practical limitations that meant you chose this relatively coarse resolution?

L558-559. “both having uniform rainfall and the same total volume of water precipitated (24 mm) over all cells” – quite a complex phrase, do you just mean spatially uniform rainfall?

L578-579. It is cool to see the feedback between the bed-elevation changes and the hydrograph – this would be an example of the kind of problem that can be explored with the component, that you could mention already in the Introduction.

Figure 10. The inset hydrographs seem to have flipped vertical axes – is that deliberate?

L636-637. “Second, all our test cases involved channels without macro-roughness elements such as large boulders, vegetation, or any type of flow obstructions that can significantly alter the flow direction.” This kind of comment suggests to me that the authors expect the component to be run at very fine resolutions, which contrasts with the relatively coarse resolution of the example – this is the kind of thing that could be cleared up by explicitly stating somewhere the scales you are thinking about, and also explaining your choice of resolution in the example.

L664-665. “The coupled OverlandFlow-RiverBedDynamics approach in our LEM employs a decoupled method to solve for river bed evolution.” The authors refer to the flow routing and surface evolution as being “coupled” or in “continuous feedback” (L181) at various points in the manuscript, but here say that they are also “decoupled”. I understand what is meant (that over many time steps the two are coupled/in feedback, but in a single time step they are applied sequentially), but this wording might cause confusion. Consider revising (here and elsewhere in the manuscript).

L671-676. I don’t understand the correction being applied here. The authors write that “only the water surface elevation is affected after a change in bed surface elevation, not water depth or discharge”, but in the illustration, it is clear that different water depths are obtained at the end of the correction. Is the change in bed elevation re-calculated with the new water depths? If so, it is not mentioned. If not, I’m not sure what the purpose of the correction is in the first place.

L710-711. “Additionally, incorporating bank erosion and channel migration capabilities could improve predictions and make long-term simulations more realistic.” This seems quite a big leap to me, since the current framework uses a uniform gird and does not explicitly distinguish channels and their widths from the rest of the landscape. Perhaps the authors could expand on how this might be achieved, or give some references to studies that attempt a similar thing. Or consider removing.