Comment on esurf-2021-35

The manuscript presents the development of a landscape evolution model (LEM) based on analytic solutions to the stream power incision model (SPIM), a simplified representation of vertical fluvial incision. Despite its simplicity, the SPIM is widely used in landscape evolution models for its demonstrated capabilities in reconstructing fluvial relief in response to changing boundary conditions at a relatively large spatial scale (i.e., Venditti et al. 2020). Model simplicity lent itself to explicit analytic solutions. Recently, there has been a growing interest in the behavior of the analytical solutions as part of studies that explored the expected dynamics of river long profiles as a function of changing geomorphic parameters and boundary conditions (i.e., Royden and Perron, 2013) and inversion applications (i.e., Goren et al. 2014a, Fox et al. 2014, Rudge et al. 2015, and many more).

The manuscript presents the development of a landscape evolution model (LEM) based on analytic solutions to the stream power incision model (SPIM), a simplified representation of vertical fluvial incision. Despite its simplicity, the SPIM is widely used in landscape evolution models for its demonstrated capabilities in reconstructing fluvial relief in response to changing boundary conditions at a relatively large spatial scale (i.e., Venditti et al. 2020). Model simplicity lent itself to explicit analytic solutions. Recently, there has been a growing interest in the behavior of the analytical solutions as part of studies that explored the expected dynamics of river long profiles as a function of changing geomorphic parameters and boundary conditions (i.e., Royden and Perron, 2013) and inversion applications (i.e., Goren et al. 2014a, Fox et al. 2014, Rudge et al. 2015, and many more).
The manuscript is well-written and clear. The motivation behind the research is well justified: developing an infinitely accurate 2D forward model for the evolution of fluvial relief in response to changing boundary conditions (in the current study, only changing uplift rate is demonstrated).
However, model implementation and its demonstrated usefulness remain somewhat immature and could benefit from additional exploration. Furthermore, the presentation of the model operation suffers from some biases. Here I specifically refer to (1) The physical meaning of the greedy algorithm that searches for acceptable topologies; (2) The application of "heterogeneous but constant uplift rate ð(ð)"; and (3) The claim that the model also solves for hillslope dynamics.
For any given topology of a fluvial drainage network and associated boundary conditions, the analytic solutions of the SPIM predict the 2D relief of all the nodes in the model for all times. There is no need for numerical iterations. The donor-receiver relation is needed only to calculate the response time, tau(x), numerically (including for time-invariant and space variable erodibility and precipitation). Once the response time is known, the solution for each node can be derived independent of the other nodes (see, for example, solutions in eq. 6, Goren et al. 2014a, andeq. 13, Goren 2016, for the case where K and P vary in time). Therefore, for any given topology, Saleve produces a graphical representation of the analytic solution.
What the analytic solutions cannot predict is the drainage network topology. Here is where the model could become more significant. Saleve addresses it by updating the topology using a greedy algorithm that attempts to minimize chi (or potentially, chi*, Willett et al. 2014, eq. 5) gradients across divides (in the manuscript, presented as elevation gradients), by updating the topology following the steepest descent in consecutive "update steps". The execution of each of these steps and the number of steps needed to achieve a stable topology have no physical time associated with them. While this issue is acknowledged in the manuscript, I'm not sure it is sufficiently discussed. Updating the topology toward the steepest descent commonly mimics drainage reorganization. However, whether each node should drain in the direction of its steepest descent neighbor regardless of the model spatial resolution is questionable. I.e., for sufficiently large grid spacing, this is probably not a very good representation. Additionally, there is no process related to this reorganization and, as stated in the manuscript, no time scale associated with it.
Many LEMs use the steepest descent criterion to update the network topology (not all though, the DAC LEM (Goren et al., 2014b), for example, uses a physical criterion to decide if reorganization should take place or the older topology should be maintained) without discussing its meaning and related time scale. This omission is particularly apparent in Saleve because topological convergence based on the steepest descent algorithm while aspiring to minimize delta z across divides is the main numerical operation of the model.
Moreover, an interesting case study, which is not attempted in the current manuscript (the author might want to consider including it), is when U varies in space and time. A convergence issue might arise in such cases: Each topology is associated with different tau(x) for each node x and a different max(tau(x)). Each topology, therefore, samples a different time range of the uplift rate history (from the present to max(tau) in the past), both generally, for the whole landscape, and specifically for each node, following eq. 7. Is there a way to guarantee topological convergence for such cases? Couldn't there be scenarios where the topology jumps between different configurations that sample different U histories without converging?
Section 4 starts with a declaration that the model input would be "potentially heterogeneous but constant uplift rate ð(ð)". In practice, in this section, the model is run with a uniform and constant U. The outcomes do not differ from those of section 5 with a time-variable uplift rate. The reason is that section 4 implicitly assumes that before U = 10 mm/yr, U was 0. This means that there is a temporal change (and no spatial change) in U, that generates a knickpoint, much like those in section 5. In fact, Figure 4c shows 2 knickpoints: One that corresponds to the transition from U = 0 to U = 10 mm/yr and the second from U = 10 to U = 20. The present form of section 4 is therefore redundant.
Hillslopes are suggested to be represented by the SPIM with m=0, giving rise to a constant critical slope. To maintain a constant slope over the grid, with information propagation based on the response time, the erodibility of hillslopes, K_3, should be infinite. This is a reasonable approximation, although a physically weird concept, for critical angle hillslopes, but there is no need for a numerical solver to represent it. These are just constant slope lines that adjust instantaneously to changes in elevation at their base (l_2).
To summarize, the concept of an infinitely accurate LEM, based on analytic solutions is appealing. However, the implementation presented in the current contribution raises some doubts regarding the model's usefulness for the landscape evolution community.

Line comments
Page 7, line 4. The concept of "thresholding the response time so that for every node ð(ð, ð¡) = min(ð(ð), ð¡)." is not clear. If it meant to address the case when the response time is > max time for which there is information about the uplift, then this is like assuming U = 0 before information is available. This is a specific choice. Why assume U = 0 and not any other value?