Transmissivity and groundwater flow exert a strong influence on drainage density

Abstract. The extent to which groundwater flow affects drainage density and erosion has long been debated, but is still uncertain. Here, I present a new hybrid analytical and numerical model that simulates groundwater flow, overland flow, hillslope erosion and stream incision. The model is used to explore the relation between groundwater flow and the incision and persis5 tence of streams for a set of parameters that represent average humid climate conditions. The results show that transmissivity and groundwater flow exert a strong control on drainage density. High transmissivity results in low drainage density and high incision rates and vice versa, with drainage density varying roughly linearly with transmissivity. The model evolves by a process that is defined here as groundwater capture, whereby streams with a higher rate of incision draw the watertable below neighbouring streams, which subsequently run dry and stop incising. This process is less efficient in models with low trans10 missivity due to the association of low transmissivity and high watertable gradients. A comparison of different parameters shows that drainage density is the most sensitive to transmissivity, followed by parameters that govern initial slope and stream erosion. These results imply that permeability and transmissivity exert a strong control on drainage density, stream incision and landscape evolution and that models of landscape evolution may need to explicitly include groundwater flow.

because of the focus on humid regions where infiltration excess overland flow is of minor importance (Dunne, 1978;Bogaart et al., 2003). Groundwater flow and saturation overland flow contribute to steady-state and transient streamflow, respectively.

60
Both components of streamflow lead to erosion and incision of the stream. In addition, the areas outside streams erode by hillslope diffusion, which is a simplified representation of processes such as soil creep (Culling, 1960(Culling, , 1963.
The model starts with a rectangular model domain with a constant slope in one direction. The rectangular model domain contains a single cross-section that is oriented perpendicular to the slope and that is used to solve the groundwater flow, overland flow and the hillslope diffusion equations. The initial topography in the direction of the cross-section is randomly perturbed. 65 The topography evolves over time as a result of stream incision and hillslope diffusion. The model simulates groundwater flow, overland flow and hillslope diffusion in the 2D cross-section. All streams are assumed to run perpendicular to the cross-section and are perfectly straight and parallel. Streamflow and stream incision are calculated by multiplying the water and sediment flux in the 2D cross-section with the contributing area perpendicular to the cross-section. The resulting model is 2.5D in the sense that it simulates water fluxes and incision of a series of perfectly parallel streams that develop along an inclined topography. 70 The workflow and equations for each component of the model are discussed in detail in the following sections.

Initial topography
The model starts with a random initial topography, which is calculated as using a series of 200 linear segments with random placement and random perturbation of the elevation at the start and end points of the segments. For the model simulations shown in this study, the average initial elevation is 0 m and the initial relief is 0.5 m. The subdivision of precipitation between evapotranspiration, overland flow and groundwater flow in the model is calculated for individual precipitation evens. For each precipitation event groundwater recharge is assumed to equal the available storage in the unsaturated zone, i.e., all groundwater stored in the unsaturated zone is assumed to eventually percolate to the groundwater table. For each point in the model domain the available storage in the unsaturated zone is calculated using the depth of the 80 watertable and the specific yield of the subsurface as: where s is storage (m), S y is specific yield (dimensionless), z is the elevation of the land surface (m) and h is the elevation of the watertable (m). Groundwater recharge for a single precipitation event is calculated as: where P d is precipitation depth per event (m), R i is groundwater recharge depth per event (m), The time-averaged potential recharge rate R p (m s −1 ) is calculated as the sum of the individual recharge events as: where f i is the frequency of precipitation event i ( a −1 ), and ∆t r is the duration of the reference time period, which is one year (s). The actual recharge rate is calculate by subtracting evapotranspiration: where ET is the evapotranspiration rate (m s −1 ). Note that for simplicity the evapotranspiration rate is assumed to be constant and independent of the depth of the watertable.
Saturation overland flow is calculated as the amount of precipitation that exceeds the available storage (s) in the unsaturated zone. For each node in the model domain and for each precipitation event the saturation overland flow depth is calculated as: where Q si is the saturation overland flow depth per precipitation event (m).

Precipitation
Precipitation events are quantified using rainfall intensity statistics for the Netherlands (Beersma et al., 2019), using a precipitationfrequency curve shown in Figure 2. The rainfall-intensity curves follow a generalized extreme value distribution with the 100 parameters given by the following equations (Beersma et al., 2019),: where η is the location parameter, γ is the dispersion parameter and κ is the shape parameter of the distribution and D is the 105 duration of each rainfall event (s). For the model experiments shown in this study a rainfall duration (D) of 3 hours is used.
The precipitation depth for a single precipitation event is calculated as: where P d is the rainfall depth per event (m) and T is repetition time (a), which is the reciprocal of precipitation frequency f (a −1 ). The model simulates overland flow and groundwater recharge for an average year. The total number of precipitation 110 events in a single year is found by adding up a series of precipitation events until the sum of the individual events matches a desired volume for the total annual precipitation (P t ): The precipitation events per year are calculated starting with a frequency (f 1 ) of 1 a −1 . Subsequently higher frequency (and lower magnitude) events are added progressively until the desired amount of total precipitation per year is reached. The 115 precipitation statistics are based on an average humid climate such as the Netherlands, with a total precipitation (P t ) of 0.75 m a −1 .

Groundwater flow
The model is based on the assumption that groundwater flow can be considered to be in steady-state on the timescales of stream and hillslope erosion processes. This was judged to be reasonable because the groundwater flow system reacts much faster than the relatively slow rates of erosion. Given these assumptions, groundwater flow and the position of the watertable can be calculated using analytical solutions of steady-state groundwater flow.
First, the out of plane component of groundwater flow, i.e., groundwater flow parallel to the direction of the nearest stream, is calculated using Darcy's equation and assuming the out of plane hydraulic gradient is equal to the (out-of-plane) slope of the nearest stream: where T is transmissivity (m 2 s −1 ), and S is stream slope (m m −1 ). Out of plane groundwater flow can be significant for cases with high transmissivity or stream slope, as shown in Fig. 3.
The remaining in-plane groundwater flow (i.e., towards the nearest stream) is calculated using the value of recharge calculated in equation 3 and subtracting the out-of-plane discharge: where R e is the effective in-plane recharge m s −1 , and L u is the length of the upstream contributing area m.
In-plane groundwater flow is calculated using the Dupuit-Forchheimer equation, which describes depth-integrated, steadystate groundwater flow between two groundwater discharge points (Forchheimer, 1886;Bresciani et al., 2016): groundwater seepage takes place. The equation assumes that the lateral differences in hydraulic head (h) are much smaller than the thickness of the aquifer (Bresciani et al., 2016).

140
For points at the edge of the model domain that are only bound by a discharge point on one side, the equation reduces to: Where L b is the distance between the discharge point and the downstream model boundary (m). The average in-plane groundwater recharge rate R e was calculated as the average effective recharge rate for all the nodes in between two seepage nodes. The seepage nodes represent points where the watertable is at the surface and where groundwater discharge occurs. The 145 position of the seepage nodes is not known in advance, but is calculated using the following iterative procedure: 1. First one seepage node is picked at the lowest elevation in the model domain and the watertable is calculated using equations 13 and 14. In most cases the calculated watertable is still above the land surface in a large part of the model domain after the first iteration.
2. Subsequently a new seepage node is added at the node with the lowest elevation in the part of the model domain where 150 the modelled watertable exceeds the land surface.
3. The watertable is recalculated using this new additional seepage node.
4. The last two steps are repeated until the modelled watertable is below or at the land surface (i.e., h ≤ z) everywhere.
An example of the calculated watertable and seepage locations following the procedure is shown in Fig. 4.

155
The water flow in each stream consists of two components, steady baseflow supplied by groundwater discharge, and transient flow that consists of overland flow. The calculation of both components is described in the following two sections:

Baseflow
The baseflow in each stream node is calculated in two steps. First, streams nodes are found by finding the node with the lowest elevation for each series of neighbouring seepage nodes in the model domain. Note that the term seepage nodes is used here 160 to denote nodes where groundwater discharge occurs. The two-dimensional (in-plane) value of groundwater flow toward each stream node (q b ) is calculated for each stream by finding the nodes contributing groundwater to each stream and by summing the product of the recharge rate at each node R and the width of each node (∆x). The contributing area is found using the taking the watertable h as calculated using equations 13 and 14 and finding the two nodes on either side of each stream where the hydraulic gradient changes direction. The three dimensional (out of plane) value of baseflow was calculated by multiplying 165 in-plane baseflow (q b ) with an upstream length of each stream: where Q b is baseflow. An example of the calculate value of baseflow is shown in Fig. 5.

Saturation overland flow
The volume of water that is contributing to overland flow is calculated per precipitation event as: Where V 0 is the volume of water to be discharged in a stream channel (m 3 ), x 1 and x 2 are the position of the topographic divide on either side of the channel (m) and Q si is the rate of overland flow for each node as calculated using equation 5. An example of the resulting distribution of precipitation excess and overland flow is shown in Fig. 6.

Waterlevel in streams 175
The waterlevel in streams as a result of baseflow and overland flow is calculated using the Gauckler-Manning equation for stream discharge. The Gauckler-Manning equation for stream discharge is (Gauckler, 1867;Manning, 1891): where v is the mean flow velocity in a stream channel (m s −1 ), K n is an empirical coefficient (m 1/3 s −1 ), which is defined as 1/n where n is Manning's roughness coefficient (s m −1/3 ), R is hydraulic radius (m) and S is the slope of the water surface, 180 which is taken as equal to the slope of the channel bed (m m −1 ). With the common assumption that channels are much wider than deep, R ≈ h c and the equation can be simplified to: where h c is the water height in the channel (m). The discharge is equal to the product of the flow velocity and the crosssectional area. To simplify the equations the cross-section for each stream is assumed to be triangular. The linear in-plane slope 185 of the channel bed (i.e., perpendicular to the flow direction in the channel) is denoted as S t (m m −1 ). The cross sectional area A of the channel equals S t h 2 . Using this the discharge equation can be written as: where Q w is the discharge in the channel (m 3 s −1 ) and C 2 = K n S 1/2 S t . Rewriting this equation yields and expression for the waterlevel in streams as a function of discharge: 2.6.4 Transient stream discharge The discharge generated by overland flow operates on short timescales of hours to days. Modelling this process directly would require short timesteps that would make the model prohibitive computational expensive. Instead here I derive new equations for the discharge over time and the total discharge in a stream following a single precipitation event. To keep the solution mathematically tractable the assumption is made that each precipitation event generates a volume of overland flow that is added instantaneously to the stream channel. The volume is subsequently discharged over time.
The continuity equation for discharge of a stream channel is given by: where V w is the water volume in the channel (m 3 ). We define the volume to be discharged as the product of a cross-sectional 200 area that changes over time and a fixed stream length L u (m): assuming that the shape of the channel is a triangle and combination with the stream discharge equation (eq. 19) yields: Integration of this equation with boundary condition h c = h 0 at t = 0 yields: The initial heigh of the waterlevel in the channel (h 0 ) is related to the initial volume of water in the channel V 0 : Adding eq. 25 into eq. 24 and replacing the constants yields an expression for the decrease of waterlevel in a channel over time in response to drainage of an initial volume V 0 : This equation was validated by comparison with a numerical solution of for the waterlevel over time (Fig. 7). The numerical solution was calculated using the discharge equation (eq. 19) to calculate discharge over time. At the first timestep the initial  Table 1, a precipitation depth of 0.0365 m, which is the theoretical maximum 1-day precipitation depth for a return time of 1 year for the Netherlands, and a contributing area for overland flow with a length of 100 m.
overland flow volume V 0 is added. Subsequently discharge is calculated and subtracted from the initial volume V 0 . The heigh of the waterlevel at each timestep was calculated using equation 25. This process is repeated until the waterlevel is less than 1 215 mm.

Stream incision
The sediment flux in a stream channel at carrying capacity is given by (Tucker and Bras, 1998;Tucker and Hancock, 2010): is water discharge (m −3 s −1 ), S is channel slope (dimensionless) and m is the discharge exponent (dimensionless) and n is the slope exponent (dimensionless). The erosion by streams is subdivided into two components, the erosion by baseflow driven by groundwater discharge and the additional transient erosion by saturation overland flow during and directly after precipitation events. The equations for these two erosion components are discussed in the following sections.
Note that the incision of streams results in an adjustment of the stream slope (S). The initial value of stream slope is fixed.

225
After each timestep a new value of stream slope is calculated using the new elevation of the base of the stream, and the elevation of the downstream edge of the model domain, which is assumed not to change over time: where z s is the elevation fo the stream in the modelled cross-section (m), z b is the elevation of the downstream model boundary (m), and L d is the distance to the downstream model boundary (m). The elevation of the downstream boundary does 230 not change over time and is the product of the initial slope (S 0 ) and L d .

Stream incision by baseflow
The sediment flux in the steam channel carried by baseflow is calculated using equation 27 and using the value of baseflow calculated using equation 15 as the value for water discharge Q w . The incision of the stream is calculated by dividing the sediment flux by the channel width (w). In addition, the erosion is divided over the upstream length of the stream by assuming 235 that erosion increases linearly from zero at the start of the stream to the maximum value at the end of the stream. This yields the following equation for erosion of the stream channel: where z is the elevation of the base of the channel (m), t is time (sec), φ is porosity (dimensionless). Note that the linear increase of erosion by baseflow along the stream is a simplification that makes the problem more tractable mathematically. In Channel width was calculated as (Lacey, 1930): where k w and ω are empirical parameters. We use values of k w = 3.65 and ω = 0.5 (van den Berg, 1995).

Stream incision by overland flow
The erosion by streams caused by the discharge of overland flow is calculated by combining the equation for stream discharge due to overland flow over time (eq. 26) with the sediment discharge equation (eq. 27). The combination of these two equations yields:

250
Integrating this equation from t = 0 to t = ∞ yields an expression for the total volume of sediment eroded from the channel after a single precipitation and discharge event: where V s is the total volume eroded in a stream by overland flow following a single precipitation event, a, b and c are constants that are defined as: The eroded volume results in incision of the stream. Stream incision by overland flow is calculated by distributing the eroded volume (V s ) evenly over the width of the stream channel (w). In addition, erosion by overland flow is assumed to 260 increase linearly from the start of the channel to the position of the modelled cross-section. This is a similar simplifcation as the linear increase of erosion by baseflow discussed in the previous section, and was likewise also used to keep the problem mathematically tractable. This yields the following expression for the incision of the stream due to overland flow: where ∆z s is stream incision following a single precipitation event (m).

Hillslope diffusion
Erosion of the parts of the model domain outside of the streams follows the hillslope diffusion equation (Culling, 1960): Where K d is the hillslope diffusion coefficient (m 2 s −1 ). This equation was solved numerically using a standard implicit finite difference approach using a matrix solver implemented in the Python module Numpy (Harris et al., 2020).

Iterative solution
The solution of the equations for groundwater flow, overland flow, streamflow and erosion follow an iterative scheme that is detailed in Fig. 8. After setting up an initial topography, the initial watertable is calculated using the procedure described in  m 2 s −1 , which is lower than values of 0.012 to 0.03 m 2 s −1 reported by de Vries (1994). The relatively low value of transmissivity was chosen because higher values led to a very low drainage densities, which made them less suitable as a starting point for sensitivity analyses. For initial slope a value of 5 × 10 −3 was chosen, which is higher than the average stream gradient in the Netherlands of 5 × 10 −4 m m −1 . Lower values resulted in models where channel initiation was so slow that the modelled drainage density reflected the initial random topography, as will be further discussed in section 4.1. 290

Model sensitivity analyses
To  (Revil, 2002;Gleeson et al., 2014;El-husseiny, 2020). The variation of the slope exponent was based on a global compilation by Harel et al. (2016). 3 Results Figure 9 shows the result of an example model run using the base case parameters values as shown in Table 1. The figure shows 305 the evolution of the land surface from an initial randomly generated surface over a timespan of 100,000 years. The system starts out with a very high drainage density, with 65 streams. Subsequently the model evolves towards a system dominated by fewer streams over time, with the final number of 13 active streams after a runtime of 100,000 years.

Groundwater capture
The decrease in number of streams is caused by a process that is referred to here as groundwater capture, by which faster eroding streams draw the watertable below neighbouring streams and reduce the baseflow and saturation overland flow of these 310 streams until they fall dry. The process of groundwater catchment capture is illustrated in more detail in Figure 10, which shows the land surface, watertable and stream fluxes for three timeslices before, during and after a groundwater catchment capture event. After 73 kyr (Fig. 10a) there are still 15 streams that are active. However, the second stream from the left hand side has smaller catchment area and a lower baseflow, and therefore incises slower than its two neighbouring streams. After 79 kyr (Fig.   10b) the neighbouring streams have incised further and have drawn the watertable below the base of the stream. The second 315 stream does generate overland flow, but has ceased to generate baseflow. At 86 kyr the watertable has been drawn down further and the watertable is too deep for the second stream to generate saturation overland flow, i.e., all the precipitation that falls near the stream can be accommodated in the unsaturated zone above the watertable. As a results the stream has become inactive.
The stream channel has partly been filled by hillslope diffusion, because the sediment flux caused by hillslope diffusion is no longer compensated by the removal of sediment by streamflow.

320
The evolution of the stream network over time is summarized in Figure 11. Initially each small depression is potentially an active stream channel, which results in a large number of active streams. However as the incision of streams increases the number of active streams drops rapidly in the first 15 kyr, followed by a slower reduction in the number of active streams for the remainder of the model run. The continuing incision and the reduction of active streams by groundwater capture also means that the area susceptible to saturation overland flow decreases over time, because saturation overland flow takes place 325 near active stream channels where the watertable is located close to the surface. Apart from a short phase at the start of the model run streamflow and stream erosion are predominantly generated by groundwater flow, which is named baseflow once it enters the stream channel (Fig. 11b, c).

Sensitivity of drainage density to transmissivity
Groundwater capture is dependent on the transmissivity of the subsurface. This is illustrated by three model runs with different 330 values of transmissivity (Fig. 12). The model run with the lowest transmissivity shows the higher drainage density (Fig. 12a) and a model run with a high value of transmissivity shows a low value of drainage density ( (Fig. 12c). The reason is that a low transmissivity results in higher watertable gradients, which makes it much more difficult for streams to draw the watertable below neighbouring streams and to capture their groundwater discharge. In contrast, high transmissivity results in a relatively flat watertable, which means that small differences in incision can already lead to a watertable that is drawn below the base 335 of streams. There is a small positive correlation between transmissivity and stream incision rate. The lower number of streams at high values of transmissivity means that these streams receive more water and therefore have a higher erosion power.
However, at high transmissivity the groundwater flow direction also changes and the proportion of out-of-plane groundwater flow increases (Fig. 3). This results in a lower in-plane groundwater flow toward streams and lower incision rates at high transmissivities. These two mechanisms compete with each other. At transmissivities exceeding 10 −2.0 m 2 s −1 the increase in 340 out-of-plane groundwater flow becomes dominant and results in a decrease in incision rate (Fig. 13).

Sensitivity of drainage density to hydrological and erosion parameters
Comparison of drainage density and incision to a number of parameters that govern erosion and streamflow show that drainage density is the most sensitive to transmissivity, followed by stream slope and slope exponent (Fig. 13). Drainage density is positively correlated with slope exponent, hillslope diffusion coefficient and precipitation, and is negatively correlated with 345 transmissivity and initial stream slope.
The strong negative correlation of drainage density and initial downstream slope is illustrated in Fig. 14 and can be explained by two effects. First, higher stream slope means that individual streams have a higher incision power, which means that it is easier to draw the watertable below adjacent streams and capture their groundwater discharge. Second, a higher initial downstream slope means that more of the total groundwater recharge is directed downstream in the out-of-plane direction, 350 and there is less in-plane groundwater flow (Fig 3). This results in flatter watertables and an increase in groundwater capture.
The slope exponent results in similar magnitude of changes in drainage density and stream incision (Fig. 13), with low slope exponents resulting in high incision rates and few streams, and vice versa.
The effect of changes in hillslope diffusion rate on the development of streams is illustrated in Fig. 15. As many previous studies have highlighted (Tucker and Bras, 1998;Perron et al., 2008), hillslope diffusion competes with stream incision. High 355 rates of diffusion result in high sediment fluxes toward streams. It is difficult for small streams to compete with these high rates of sediment delivery, and as a result only larger streams survive. In contrast, low rates of hillslope diffusion allow more streams to persist.
Interestingly, precipitation and specific yield have a limited effect on drainage density (Fig. 13). Precipitation affects groundwater recharge. Given that groundwater flow depends on the ratio of recharge over transmissivity (see eq. 13), in principle 360 recharge has an effect of equal magnitude as transmissivity, but opposite sign. However, the variability of transmissivity of the subsurface is several orders of magnitude, whereas the variability of groundwater recharge is much lower, approximately one order of magnitude in humid regions (Moeck et al., 2020). This results in an overall much more important effect of transmissivity on drainage density.
Specific yield only has a subtle effect in that lower values of specific yield increase the volume of saturation overland flow 365 and make it more difficult for streams to run completely dry, even if they are disconnected from the watertable. However, given the subordinate importance of overland flow in generating streamflow and erosion in these model experiments, this effect is very modest and does not change the modelled drainage density or incision rate after a model runtime of 100,000 years ( Fig.   13).

Limitations of the model code
The model code presented here was intentionally kept as simple as possible to keep the solution mathematically tractable and the computational effort manageable. The main limitations are that the model is 2.5D and only represents a system with perfectly parallel and straight streams, groundwater flow is steady-state and the contribution of transient groundwater discharge or subsurface stormflow to streamflow is neglected, and the treatment of erosion by overland flow is highly simplified, with 375 an instant addition of all overland flow to the nearest active stream channel, and a highly simplified distribution of the eroded volume over the length of the channel.
In spite of these limitations the conclusions of the importance of groundwater flow for drainage density and stream incision are arguably relatively robust. The model results show that groundwater capture is a somewhat inevitable consequence of coupling groundwater flow, streamflow and erosion equations. Regardless of the exact equations and numerical implementation 380 used, the watertable will always crop out at perennial streams, and the incision of these streams will draw the watertable down.
Smaller streams losing their groundwater discharge and eventually falling dry is therefore a logical consequence of differences in incision rate between streams. These conclusions are also supported by previous model studies that found that groundwater flow had a strong effect on erosion (de Vries, 1994;Huang and Niemann, 2006;Barkwith et al., 2015;Zhang et al., 2016) and by observations that found a correlation between permeability, transmissivity, or lithology and drainage density at a number of 385 locations (Carlston, 1963;Luo and Stepinski, 2008;Bloomfield et al., 2011;Luo et al., 2016).
However, one thing that the model does not represent well is the initiation of stream channels. Due to the 2.5D nature of the model the initial topography is represented only in the modelled cross-section, and all streams are effectively parallel. This means that initially, small streams develop in most small depressions. The presence of many small streams divides the water to be discharged over many streams, which each have a very low incision power but are nonetheless considered active streams 390 in the model. In reality, many of these streams would not initiate, because the initial topography would not consist of a set of small depressions that are lined up and connected along an inclined slope. Instead, most of these depressions would be isolated. Therefore the number of active streams at low runtimes should be much smaller than the high number in the model runs reported here (see for example Fig. 11a). In reality, the stream network would develop much more gradually by headward erosion, as shown by natural experiments (Schumm, 1956;Wells et al., 1985;Kashiwaya, 1987;Talling and Sowter, 1999).

395
For the final model results presented in this study the unrealistic initiation of streams is less of a problem, because the stream networks adjusts relatively rapidly in the first few 100s to 1000s of years as stream incision progresses far enough for most of the initial streams to fall dry and for the initial topography to be irrelevant. Therefore for the final model runtime of a 100,000 years the unrealistic initiation of stream networks does not affect the results. The exception is however, formed by model runs with low transmissivity (< 10 −4 m 2 s −1 ) or low intial slope (< 10 −4 m m −1 ), for which the ability to capture groundwater 400 discharge and the incision power of the streams is too low to escape the initial phase with a dense stream network.

Groundwater capture
To my knowledge the process of groundwater capture that is illustrated here has not been proposed before. de Vries (1976) and de Vries (1994) discussed the opposite process, in which new streams establish themselves when the watertable reaches the land surface following an increase in precipitation and groundwater recharge. Twidale (2004) suggested that stream incision 405 may induce a positive feedback by increasing groundwater flow and runoff to streams that incise faster, but did not provide any further details.
The process of groundwater capture is expected to be important wherever streamflow is dominated either by groundwater discharge or saturation overland flow. In systems where infiltration-excess flow is dominant, which includes most semi-arid and arid regions (Kidron, 2021), groundwater capture could still occur and would determine which streams discharge groundwater 410 and are active perennially, but streams that are detached from the groundwater table would still be able to incise. In addition, groundwater capture may be less important for systems dominated by subsurface stormflow in perched aquifers that are detached from the watertable. The importance of this flow mechanism is however uncertain and under debate (Freeze, 1972;Chifflard et al., 2019).
Note that the difference between the groundwater capture process and the more well studied process of stream capture 415 (Whipple et al., 2017) is that groundwater capture does not require a physical connection between two streams, and may therefore not leave a clear trace in the landscape, as opposed to river capture. The capture of the groundwater discharge of one stream by another results in the first stream falling dry, but the hillslope between these streams remains intact, as can be seen in Fig. 10.

420
A new model is presented that couples equations that govern groundwater flow, overland flow, streamflow and erosion and represents humid regions where infiltration-excess overland flow is negligible. The coupling of these equations reveals a strong dependence of drainage density on groundwater flow and transmissivity. The dependence of drainage density on groundwater flow is controlled by a process identified here as groundwater capture, by which faster incising streams draw the watertable below neighbouring streams, which fall dry and lose the ability to incise any further. This process is more efficient when the 425 watertable is relatively flat, which can be either due to low groundwater recharge or high transmissivity, with transmissivity often being the dominant control due to its high variability. Sensitivity analysis show that the importance of transmissivity for drainage density is roughly equal to parameters that control the incision power of stream channels. The importance of transmissivity exceeds other hydrological and erosion parameters such as precipitation and hillslope diffusion coefficient. The results imply that groundwater and transmissivity may be a dominant control on drainage density in humid regions. Landscape 430 evolution models may need to include groundwater flow and the watertable. At the same time a closer integration of landscape and groundwater data and models may help improve knowledge of the evolution of the earth's surface.
Code availability. The model code presented in this manuscript is named the groundwater flow, overland flow and erosion model, GOEMod, and has been published at Zenodo (Luijendijk, 2021). The code is also available at GitHub (https://github.com/ElcoLuijendijk/goemod).