Ruetenik et al. introduce and apply a method to find optimal parameter values in Landscapes Evolution Models using a global 10Be compilation. They use a Montecarlo approach to run multiple one-time-step simulations to find the values of m, n, K, D minimizing the differences between simulated erosion rates and apparent erosion rates measured by CRN. They finally group the results by lithological and climatic domains to investigate the behaviour of optimal parameters function of these variables. I believe it offers a novel way to calibrate model parameters and even highlight interesting patterns about the behaviour of the stream power law and hillslope diffusion parameters correlated with precipitation patterns. I suggest this work is suitable for publication in esurf and will benefit the geomorphological community, however, I have two concerns about the method itself and the analysis of the results as well as other minor comments that I suggest should be addressed prior to publication. The manuscript is relatively short, which is not an issue, but could benefit from additional information and sensitivity analysis to clarify the method and its potential/limitation. Please find bellow my different comments. 

 Boris Gailleton 

 Main concerns:

One of the main advantages of the method according to the authors is the consideration of all the pixels in the investigated catchments, as opposed to only considering fluvial processes like Harel et al., 2016 or Marder and Gallen (2022). However, this also raises a number of questions as 10Be data is sampled at the river mouth. First, rivers are by nature only representing a minority of pixels in the catchments and are therefore under-represented statistically-speaking. Secondly, the connectivity with the sampling site will be significantly different between river processes and hillslope processes: rivers “instantly” transmit the sediments to the sampling sites while hillslopes may be disconnected. The models proposed in this study seem to assume that any change in the current topography through time is instantly transmitted to the outlet of the river, which may be true on a landscape simulated only with the stream power law and linear hillslope diffusion, but more difficultly on a real landscape. Finally, the authors present a method running LEMs on one snapshot of the (current) topography, while cosmogenic nuclides integrate larger time scales. Stochastic events (e.g., monsoon, landsliding) and/or complex sediment dynamics (e.g., recycling, residence time, intermediate traps) induce potentially significant variability in the signal (see Dingle et al., 2018 for an extreme example: https://doi.org/10.5194/esurf-6-611-2018, 2018). I acknowledge these are challenging points to include in the actual simulations - especially because the process laws of the models do not explicitly account for the sediment flux (they only express vertical topographic changes) - but I suggest it is not sufficiently discussed or addressed in the manuscript in its current form.

My other main concern is about the results, expressed in Figure 2, 3, 4, 5 and 6. It is not entirely clear whether the method picks a set of parameters, runs this same set for all the 3618 selected catchments, and calculates the global performance of this set of parameters; or if the Montecarlo process is done for each basin separately and then averaged. In other words, it is not entirely clear what these points exactly represents. In any cases, I suggest it would be relevant to show and discuss the breakdown of the global datasets: one could expect an average best parameter to be misrepresentative of the whole variability: multiple clusters could eventually emerge. This is also applicable to the Lithology and Climate comparisons: grouping metamorphic rocks into a unique category can be misleading as one would expect a very different rock strength via K (or even erosion processes via n) in a schist domain vs a gneiss domain – while both metamorphic. On the opposite, showing that the best-fit parameters are relatively constant in different contexts would be an important finding.

Other comments:

- There is a missing subsection title in the method part, when describing fluvial and hillslope processes.

- I find the use of the term “concavity” for the non-linear hillslope diffusion equation a bit confusing as it commonly refers to river processes – as stated in the manuscript. I would recommend using “hillslope concavity” or another term to avoid any confusion.

- The “non-linear hillslope diffusion” term can also refer to Roering et al. (1999) type of laws – that include a critical slope component in the equation. While p indeed makes the equation non-linear, it would be relevant to explicitly state the difference with the more common (in my experience at least) alternatives like Roering et al., 1999.

- The ratio between m/n is set to 0.45. While the average is in accordance with most of the literature, this value can vary quite significantly (this is acknowledged in the manuscript, and for example Wobus et al, 2006 suggest a 0.3-0.6 range using theoretical prediction, Harel et al., 2016 0.51 +/- 0.14 based on geomorphometrics and Gailleton et al., 2021 details significant variations at global scale). I am not suggesting to rerun the full analysis for ranges of m/n, but I would recommend to at least test the sensitivity of the prediction to this ratio in at least one site or one subset of the study. I believe the study would be more robust and if the method is resilient to this ratio, it would represent a notable advantage over others. 

- To echo a comment about smoothing the DEM from Reviewer #1: the location of the river network is of prime importance as it define the transition between advective and diffusive processes. However, the choice of methods and their implications is not very detailed for (i) the preprocessing of the dem to remove local minima (or smooth the data) and (ii) the flow topology (D8). These can significantly affect the location of the rivers and the proportion of river pixels, especially in the case of lower-relief areas. Another point that would deserve a mention is the flow direction. Steepest descent algorithm reduces the river as 1D series of pixels with the same area no matter the drainage area. Because the process equations used to calculate E_predicted depends on the location of the river I believe this is an important point to raise – for example mentioning alternative like D_inf or multiple flow directions which would change the ratio rivers/hillslope. Again, I am not suggesting to rerun all the analysis with a different topology, but more to address this point in the method or discussion part.

- While some parameters benefit from the Montecarlo “brute-force” approach, others can be estimated through other means and optimize the process – especially because this approach is applied on real topography rather than simulated ones. For example, it is possible to estimate ranges of possible m/n and Ac through geomorphometry (e.g., Mudd et al., 2018, Wobus et al., 2006). It is also possible to directly determine the location of the river network (and therefore a more natural delimitation between hillslopes and fluvial processes) and bypass the need of Ac (e.g., Clubb et al. 2014). 

  - The parameter p is both used for the precipitations in equation 8 and the (hillslope) concavity exponent. I suggest to at least use a capital letter P for precipitations or another notation to avoid confusion. 

  - It is unclear in the method section whether the model considers spatial variations in K_lith, D and p (precipitation). If it does not, it can be an important point of discussion: catchments showing significant lithologic contrasts can really obscure the reading of geomorphometrics (e.g., Gailleton et al., 2021b - https://doi.org/10.1029/2020JF005970). 

 Minor comments:

l.40: I agree. I would also add a technical limitation to the parameter values: some are inter-connected - for example a value of K is only valid in a given context, and if n, m or other associated laws are modified a new K needs to be constrained. 

l. 43: “vital” may be a bit of an exaggeration in my opinion, maybe “natural” or “crucial”? 

l.107: “discharge” more than “flood” variability? 

l.109: Harel et al. 2016 also reports significant variability in this value, also suggested by Gailleton et al. 2021 from geomorphometrics. 

l. 112: As long as K remains constant (can be good to remind that). 

l.139: While in line with the median in our global compilation, we also highlighted how variable this ratio could be (or at least its expression in the shape of the topography – which might not always reflect their translation into m/n when the system is significantly transient). 

l.200: It is not clear what K* represents at this point of the manuscript. Please clarify.

l.232: missing “s” in variation.

l.246: refer to equation.

l.293: I disagree with this sentence: Harel et al., 2016 calculate ksn from M_chi using the method detailed in Mudd et al., 2014. M_chi is the gradient of a statistically segmented Chi-elevation River profile and is equal to ksn if A_0 = 1, otherwise proportional to it to a factor A_0. M_chi is then directly relatable to K and n within the stream power referential via the same relationship” M_chi = (E/(K * A0^m))^1/n. 

l.295-300: while Harel et al., 2016 indeed only consider the fluvial section of the landscape, they do integrate spatial variation of ksn within a basin and calculate it for each river pixel. I suggest rephrasing/reworking the whole paragraph comparing the results with Harel et al., 2016.  

l.301: Braun and Willett, 2013 barely mention linear diffusion, I suggest to use more relevant reference(s).