Reply to review by Marijn van der Meij. The replies are shown in italics below each comment.

The author presents a detailed and well-documented model named GOEMod that simulates the effect of stream incision on landscape evolution through groundwater and overland flow. The model simulates a set of hydrological and geomorphological processes that affect the groundwater level, stream incision and drainage density. He identified a process he calls groundwater capture, where streams with high incision rates draw the water table below neighbouring streams. With a sensitivity analysis the author identifies the soil transmissivity as main parameter controlling the drainage density.

The manuscript is well written and the model is very well documented. As the author mentions in his Introduction, groundwater flow is often simplified or not at all included in current landscape evolution models. This work shows a potentially important control of groundwater flow on landscape evolution. 

However, I think that the manuscript can improve on three points and would benefit greatly from additional discussion of the model and its results. The three points are the balance of the manuscript, the selected set of parameters and a comparison of model results with field evidence. I address these points in more detail in the general comments below. Next to that, I added several specific comments and technical corrections to the attached document. I look forward to seeing the revised manuscript.

With best regards,

Marijn van der Meij

1 Balance of the manuscript

In the title and abstract, the manuscript is presented as a study to identify the relation between groundwater flow and drainage density. For this purpose, the model GOEMod was developed. The manuscript itself however leans very much towards the model development and description of model behaviour and only shortly addresses the validity of the model and lacks a comparison with field evidence to support the model findings. I suggest to address the other comments below to improve the balance between model description and interpretation of the findings in the manuscript.

Reply: I agree and have added a section that discusses the comparison between the model and field evidence for changes in drainage density and its correlation with transmissivity. This discussion is based on the references that were mentioned in the 2nd paragraph of the introduction in the previous version of the manuscript. However I would like to reserve a more elaborate discussion for follow up work. Compilations of transmissivity data are very scarce which makes a comparison with drainage density difficult. However, I am also about to submit a manuscript with a new large-scale transmissivity dataset that may be very helpful there.

2 Selected parameters

The setting of the model is a humid climate zone and is roughly based on The Netherlands. The model builds further on stream density studies from the Netherlands (lines 40-41 & 50-52). Also, data from The Netherlands are used to parametrize the model. This data includes rainfall data, porosity, transmissivity and stream slope. However, there are several parameter choices in the model that deviate from this setting to simplify calculations or achieve a better model result:

    The transmissivity, porosity and slope parameters are adjusted (lines 284-289 & 398-400);

    The runtime of the model is set to 100.000 years;

    The elevation downstream of the model domain does not change over time (lines 225-226) and precipitation is also constant (lines 114-116).

Over the long simulation span of 100.000 years, downstream elevation cannot be assumed constant due to fluctuating sea water levels. Also, rainfall cannot be assumed stable over this time period. Fluctuations in these parameters will likely have a larger influence on stream incision than groundwater capture over these geological time scales.

This make me wonder if groundwater capture as resulting from the model is a realistic process in real landscapes, or whether it can only be simulated in modelling studies due to the unlikely conditions when under which it seems to occur. I would like to see some discussion in the manuscript on how realistic the parameter choices are and how likely it is that groundwater capture plays a role in the evolution of real landscapes.

Reply: The previous version manuscript unfortunately used a value of sediment transport coefficient that was much too low for transport-limited systems. An analysis of published sediment discharge datasets yielded a new more realistic and much higher sediment transport coefficient that results in a much faster incision and a much faster adjustment of stream networks. The stream network now reaches a dynamic steady-state over a timespan of ~2000 years. Because of this the revised manuscript now shows the results for model runs of 10,000 years as a standard, which is a timespan that sea level and climate can be expected to be relatively stable. The transmissivity and slope parameters in the revised manuscript now also closely follow the average values for the Netherlands reported by de Vries (1994) and others.

3 Field evidence

Following from the previous comment, I’m also missing a comparison of model results with field data to further support the concept of groundwater capture and the effect of transmissivity on drainage density. These points are briefly mentioned in the Introduction (lines 20-21 & 28-32), but are not discussed in the light of the model results. I would like to see such a comparison in the paper to further support the model results, for example in the form of a case study or as a meta-analysis of areas with different transmissivity and drainage densities.

Reply: Agreed. I have added a new section in the revised manuscript that compares the model results and published estimates of drainage density and transmissivity or permeability

Reply to additional review comments in the attached pdf:

Comment on Line 38: I'm curious about your ideas how to improve the hydrological modelling in landscape evolution models, for example using your approach. Could you elaborate on that in the discussion?

Reply: I have briefly commented on ways to improve the hydro(geo)logy in landscape models in the revised manuscript.

Comment L. 58: This focus requires a bit more introduction.

Reply: Agreed, this choice is now introduced in the last paragraph of the introduction.

Comment line 74: Is there a specific landscape where this initial relief is based on?

Reply: This admittedly is purely a guess. I have not been able to find data on typical relief for undissected landscapes (without streams).

Comment line 75: Can you switch sections 2.3 and 2.4? Section 2.3 already addresses the precipitation that is described in 2.4 

Reply: Agreed, i have switched the sections.

Comment line 76: Clever way to determine the partitioning of precipitation. Can you elaborate a bit on why you have to determine this per individual event?

Reply: The reason is that erosion by overland flow is also calculated per rainfall event. And therefore an calculation of overland flow volume per event is necessary.

Comment line 91: How is the ET determined in the model?

Reply: The ET is a fixed value as is now mentioned more clearly in the revised manuscript.

Comment Fig. 3: Can you change this to 1? A fraction of 1.05 doesn't make sense

Reply: Thanks for spotting this. The colorbar scale has been updated.

Comment Fig. 5: Can you clarify that the colours of the traingles refer to the base flow of the streams?

Reply: I have added a clarification to the figure caption.

Comment line 240: What is the effect of this and the other simplifications/assumptions on your model results? And is it possible to improve this description in (a later version of) the model?

Reply: It would perhaps be possible to relax this assumption. The resulting equation for a more realistic distribution of erosion is possible, but the integrand may not be solvable, and at the very least would reserve quite a bit of dedicated time and a lot of coffee. I would like to reserve a more thorough exploration of this and other simplifications for a future manuscript.

Comment line 275: If I understand correctly from the manuscript, the overland flow and groundwater recharge are calculated on an event basis, while hillslope erosion and stream incision are calculated on a yearly basis. Is this correct? If yes, could you elaborate how you bridge these time scales? If not, could you clarify how you dealt with the different time scales of the processes?

Reply: That is correct. Overland flow and recharge are calculated at an event basis. Recharge is summed over a year to yield an average recharge rate. However, overland flow and erosion are calculated on an event basis, and are then summed over a year to yield an average overland flow erosion rate. I have added this clarification to the revised manuscript.

Comment Table 2: How many values did you select for each parameter? And did you also vary multiple parameters in one run to study their interactions?

Reply: I have added the number of values to the table. The sensivitiy analysis used a relatively simple one at a time sensitivity analysis. It would be interesting to explore changes in multiple parameters and their interactions, but I would like to reserve this for future work

Comment Fig. 13: In the text you name this the drainage density instead of the stream density

Reply: Thanks for noticing, i have corrected this.

Comment line 370: Can you elaborate on how these simplifications affected the model results and how they might change with more complex descriptions in the model?

Reply: A brief assesment of the robustness of the conclusions is already contained in this section. A more elaborate discussion & more complex model physics are planned for follow up work.

Comment line 399: Aren't these parameters typical for the Netherlands, where you based your simulations on?

Reply: Not really, at least for the parts of the Netherlands that are above sea level, where tranmsissivity tends to be higher (approx. 1-5 10^-2 m2/s, see for instance Bogaart et al. 2003). And the average topographic gradient is somewhat higher too with 1/2500 m/m.

Comment line 404: Is GOEMod also able to simulate this process?

Reply: It is. The revised manuscript contains a new set of model simulations that explores the persistence of drainage networks and that includes model runs where base level or precipitation increases, which results in the reactivation of abandoned streams.

Reply: Many thanks for the thorough and very helpful review.

Reply to review by Stefan Hergarten. The replies are shown below the comments in italics.

Dear Elco Luijendijk;

first, I would like to mention that I enjoyed reading your paper about drainage density. Maybe I am even a bit biased because I really like this king of modeling -- playing with the equations, simplifying them in such a way that they become tractable, and then looking what we can get out of them. I noticed that mention the limitations in the sense of some oversimplification at several times. In this case, I do not worry so much about the simplicity of the involved processes. At some points, even more simplifications would have been possible, and a further reduction of the variety of processes could even improve our understanding. Anyway, this is my personal point of view and not a criticism on your work as a reviewer.

Nevertheless, I am still somewhat critical concerning the results and the depth of the discussion. Your model is quite specific. So I see the main merit not in providing a general tool, but in addressing the question for the principal dependencies of drainage density. Playing devil's advocate, I could suggest to assume a given across-slope profile (maybe not uncorrelated random numbers, rather something like a random walk) and apply your steady-state groundwater flow model with constant recharge. If the transmissivity increases, the gradients in hydraulic head decrease, which reduces the number of points where the groundwater level reaches the surface. This would exactly be a decrease in drainage density with increasing transmissivity. And then I would argue that this is already the main result of your paper and that you could omit the rest of the model.

Reply: That is true in principle, i.e., topography, recharge rate and transmissivity alone would determine drainage density. So in principle a dynamic model like the one presented in the manuscript would not be necessary to estimate transmissivity and one could use a static model like the model the reviewer proposes to test what drainage density would be for different sets of hydrogeological parameters. This is actually more or less what has been done in several studies cited in the introduction (Luo et al., 2010; Luo and Pederson, 2012; Bresciani et al., 2016), where existing drainage networks and topography have been used to constrain transmissivity. However, the drainage density that one would obtain from such a static model depends strongly on the topography and the degree of incision. All else being equal, a topography with a high relief would result in a much lower number of active streams than if one would analyze a topographic cross-section with a low relief. And relief is in itself of course affected by the number of streams, because this determines how much water each stream discharges and how much incision power these streams have. So there is a feedback here between topography, groundwater flow and drainage density that I tried to explore using the model presented in this manuscript. The main research question that i tried to answer was if a process can be identified that governs the relation between transmissivity and drainage density.

While this is a exaggerated, my feeling is indeed that you do not get the maximum out of your model. After following the theoretical part, I was a bit disappointed by the results and discussion sections.

There is one central point where I cannot assess how powerful the model is. You perform simulations over a given time span. Some of these (highest transmissivity) already arrive at one single stream over the entire slope (so 0.1 streams per km at 10 km domain size). I would like to know whether the model is able to predict the existence of more than one stream over very long time (so in some kind of steady state) or whether it always arrives at a single stream if we just wait long enough. If the latter was the case, this would reduce the merit of the model seriously.

Reply: Thats a very good point. To test this I included a model run in the revised version of the manuscript that show the evolution of drainage density over a longer timescale of 1 million years. The results show that these runs do reach a steady-state with more than one stream. The reason is that larger streams that incise the fastest at the start at some stage lose their advantage because the higher incision rate results in a lower slope of the stream in the out of plane direction. This lowers the incision and allows smaller streams to catch up. So this is a negative feedback that provides a limit on drainage density. I have added a more elaborate discussion of this effect to the revised manuscript. This process is obviously very dependent on the evolution of the regional baselevel. In case where the baselevel does not limit stream slope and incision, the system does evolve to one stream. However, the model results in the revised manuscript use a new more realistic value of sediment transport coefficient that results in much faster incision and a steady-state with > 1 streams within ~2000 years.

In the following, I discuss a few more specific points about the results, in particular about the parameter study.

(1) Figure 13(a) shows some kind of power-law dependency of the drainage density on both the transmissivity and and the initial slope. It looks as if the exponent was -1, so the drainage density seems to be proportional to 1/transmissivity and to 1/slope. This would be the chance to leave the purely qualitative level of the discussion and to get a least some quantitative results. I would urge you to investigate these dependencies and to try to explain these directly from the model equations.

Reply: Thanks for the suggestion. I would however like to reserve a more quantitative analysis of the results for a follow up manuscript. In particular the steady-state drainage density may be predicted by analytical equations that i am currently developing, but that will take some time to finalize.

(2) The consideration of erosion with a fixed base level seems to be a bit unfortunate to me. Practically, you mix the buildup of across-slope topography by incision with the decay of slope gradient. As a rough estimate: 5000 m downstream length, 5e-3 initial slope is 25 m difference in topography. Incision is some meters, so not so strong that the decay of overall slope is already dominant, but in turn not negligible. It would be better to separate both effects.

Reply: Thanks for the remark. I have included a base level change in the model runs in the revised manuscript. For the area in the Netherlands that this model is loosely based the rate of baselevel change is very low (20 mm/kyr). However, the sensivity analyses in the revised manuscript explore the effect of higher baselevel decrease rates, which do have an important effect drainage density.

(3) You consider a variation of the slope exponent n. In the context of these models of fluvial incision and sediment transport, it makes no sense to change the exponent of the power-law relation and to keep the factor of proportionality constant. Beyond this, the ratio of the discharge exponent m and the slope exponent n is constrained much better than the value of n itself, so changing n alone really makes no sense. So you should keep m and n constant and consider the sediment transport kf as a variable parameter in your sensitivity analysis.

Reply: Thanks for the comment. I have repeated the sensitivity analysis with kf as a variable parameter instead of n and have added these results to the revised manuscript. I have calculated the variability of the kf using compiled sediment discharge data by Lammers and Bledsoe (2018, https://doi.org/10.1002/esp.4237). In the process I also discovered that the value used in the previous version of the manuscript was much too low for transport-limited streams, and have now used a more realistic and higher value. As a result streams incise faster and the stream network reaches a steady-state in within ~2000 years in the new set of model runs used in the revised manuscript.

(4) Just as an option for future work: periodic boundary conditions may be simpler here.

In summary, I feel that you developed a very nice model, but earning the fruits of your work is somewhat weak. Nevertheless, I would not argue against publishing it if you were running out of time. In this case, point (3) would be the point where I would definitely request a revision.

I hope you find my suggestions helpful.

Best regards,

Stefan

Reply: Many thanks for the very helpful review