# Generation of realistic synthetic catchments to explore fine continental surface processes

Raphaël Bunel, Nicolas Lecoq, Yoann Copard, Eric Guérin, Marco van de Wiel, Nicolas Massei

► To cite this version:

Raphaël Bunel, Nicolas Lecoq, Yoann Copard, Eric Guérin, Marco van de Wiel, et al. Generation of realistic synthetic catchments to explore fine continental surface processes. Earth Surface Processes and Landforms, 2021, 46 (3), pp.593 - 610. 10.1002/esp.5048. hal-03558806

HAL Id: hal-03558806

https://hal.science/hal-03558806v1

Submitted on 2 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Generation of realistic synthetic catchments to explore fine continental surface processes

Bunel, R, Lecoq, N, Copard, Y, Guérin, E, Van de Wiel, M & Massei, N

Author post-print (accepted) deposited by Coventry University's Repository

**Original citation & hyperlink:**

'Generation of realistic synthetic catchments to explore fine continental surface processes', Earth Surface Processes and Landforms, vol. 46, no. 3, pp. 593-610.

https://dx.doi.org/10.1002/esp.5048

DOI 10.1002/esp.5048

ISSN 0197-9337

ESSN 1096-9837

Publisher: Wiley]

This is the peer reviewed version of the following article 'Generation of realistic synthetic catchments to explore fine continental surface processes', Earth Surface Processes and Landforms, vol. 46, no. 3, pp. 593-610, which has been published in final form at https://doi.org/10.1002/esp.5048. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving.

Copyright © and Moral Rights are retained by the author(s) and/ or other copyright owners. A copy can be downloaded for personal non-commercial research or study, without prior permission or charge. This item cannot be reproduced or quoted extensively from without first obtaining permission in writing from the copyright holder(s). The content must not be changed in any way or sold commercially in any format or medium without the formal permission of the copyright holders.

This document is the author's post-print version, incorporating any revisions agreed during the peer-review process. Some differences between the published version and this version may remain and you are advised to consult the published version if you wish to cite from it.
1 **Generation of realistic synthetic catchments to explore fine**

2 **continental surface processes.**

3 **Raphaël Bunel¹*, Nicolas Lecoq¹, Yoann Copard¹, Eric Guérin², Van de Wiel**

4 **Marco³,**

5 **Nicolas Massei¹**

6 ¹University of Normandy, UNIROUEN, UNICAEN, CNRS, M2C, 76000 Rouen,

7 France

8 ²University of Lyon, LIRIS, CNRS, UMR5205, 69000 Lyon, France

9 ³Centre for Agroecology, Water and Resilience, Coventry University, Coventry,

10 United Kingdom

11 * Corresponding author. Tel.: +33 235146807.

12 E-mail address: raphael.bunel1@univ-rouen.fr (R. Bunel).

13
14 Abstract

15 Understanding, analysing, and predicting the erosion mechanisms and sedimentary
16 flows produced by catchments play a key role in environmental conservation and
17 restoration management and policies. Numerical case-testing studies are generally
18 undertaken to analyse the sensitivity of flood and soil erosion processes to the physical
19 characteristics of catchments. Most analyses are conducted on simple virtual
catchments with physical characteristics that, unlike real catchments, are perfectly
20 controlled. Virtual catchments generally correspond to V-shaped valley catchments.
21
22 However, although these catchments are suitable for methodical analysis of the
23 results, they do not provide a realistic representation of the spatial structures of the
24 landscape and field conditions. They can, therefore, lead to potential modelling errors
25 and can make it difficult to extend or generalize their results. Our proposed method
26 bridges the gap between real and traditional virtual catchments by creating realistic
27 virtual catchments with perfectly controllable physical characteristics. Our approach
28 represents a real alternative to traditional test case procedures and provides a new
29 framework for geomorphological and hydrological communities. It combines a field
30 procedural generation approach, geographic information system processing
31 procedures, and the CAESAR-Lisflood landscape evolution model. We illustrate how
32 each of these components acts in the process of generating virtual catchments. Five
33 physical parameters were adjusted and tested for each virtual catchment: drainage
34 density, hypsometric integral, mean slope of the main channel, granulometry, and land
35 use. One of our virtual catchments is compared with a real catchment and a virtual
catchment produced by a standard method. This comparison indicates that our
36 approach can produce more realistic virtual catchments than those produced by more
37 traditional methods, while a high degree of controllability is maintained. This new
38
method of generating virtual catchments therefore offers significant research potential
to identify the impacts of the physical characteristics of catchments on hydro-
sedimentary dynamics and responses.

42

Keywords: realistic virtual catchments, controllable physical characteristics,
procedural generation, GIS processing, CAESAR-Lisflood.

45
46

# 1. Introduction

47

48

49

50

51

52

53

Understanding, analysing, and predicting the erosion mechanisms and sedimentary flows produced by catchments play a key role in environmental conservation and restoration management and policies. The hydro-sedimentary dynamics of catchments are controlled by natural climatic and geological forcings, which lead to the long-term structuring of landscape forms (i.e. geomorphology), as well as by anthropogenic forcings, which modify catchment characteristics over much shorter time scales (e.g. land use, Schumm and Parker, 1973; Nicholas and Quine, 2007).

54

55

56

57

58

59

60

61

62

63

64

65

Several studies have been conducted to understand the influence of these forcings on the evolution of catchments and river systems (e.g. Macklin and Lewin, 1993; Brewer and Lewin, 1998). Most of them have been based on the acquisition of extensive field data. However, studies based on these data sets suffer from several limitations related to the spatial and temporal scale of the considered processes. First, field data, acquired on an ad hoc basis, do not exhaustively cover the great diversity of environments in terms of geomorphology, geology, pedology, and land use. It is therefore impossible to use these data to accurately analyse the sensitivity of a system to any of these parameters (Coulthard and Van De Wiel, 2012). Second, there is discontinuity between the temporal scale of the measurements and the temporal scale of the change: landscape evolution typically occurs on time scales too long to be understood based on field measurements only.

66

67

68

69

70

Sedimentary records represent alternatives for analysing the evolution of hydro-sedimentary dynamics and the forcings responsible for them. Nevertheless, these archives, which are generally incomplete because of erosion or non-sedimentation of particles, provide only a partial description of the past hydro-sedimentary dynamics of a catchment (Straub et al., 2020). Their low resolution and entanglement of forcings
make it difficult to quantify the influence of each forcing and the effects of interactions (Notebaert et al., 2011; Hancock and Coulthard, 2012). Therefore, many uncertainties remain regarding the hydro-sedimentary dynamics of these catchments.

Numerical modelling can be an alternative solution to understand this dynamic by allowing the sensitivity of a system to a single parameter to be analysed over very long time scales. This can be done by adjusting only one parameter at a time and keeping all other parameters constant. This approach thus makes it possible to analyse separately how environmental conditions (e.g. land use) and internal processes (e.g. flows and sediment transport) influence the morphological evolution of catchments over the long term. Monitoring this evolution on large spatial and temporal scales requires the use of models with low computation time requirements to limit the duration of simulations. Landscape evolution models (LEMs) have been designed to operate at low computational costs by reducing the complexity of the processes involved in soil flow and erosion. These LEMs therefore allow simulations to be carried out over long periods of time and at the scale of entire catchments, unlike reductionist models that describe the physical processes in more detail (Bishop, 2007). The application of LEMs to virtual or case-test catchments is generally carried out to better understand the impacts of the physical characteristics of the catchments on their hydro-sedimentary dynamics (Baartman et al., 2013). This analysis would not be possible with real catchments because of their high complexity and heterogeneity. The use of virtual catchments with selected characteristics makes it possible to isolate the influence of each parameter under a range of environmental conditions and to avoid interference between parameters. However, although these catchments make methodical analysis of the results possible, they do not provide a realistic representation of the spatial structures of the landscape and field conditions. They are generally square or
rectangular in shape and consist of a single valley formed by two symmetrical inclined planes. These catchments can therefore lead to potential modelling errors (e.g. production of unusual flows and sediment flows) and can make it difficult to generalise the results.

In this study, we propose a method to generate catchments with realistic and perfectly controllable morphometric characteristics to provide the opportunity to test, in subsequent work, the sensitivity of hydrological and sedimentary processes to the physical characteristics of the catchments. This type of empirical approach based on the use of synthetic cases is indeed fundamental to conducting such sensitivity studies, as existing databases are never sufficiently exhaustive in terms of basin shapes and hydro-sedimentary measures. No sensitivity analysis is conducted for this study as it focuses strictly on the development and evaluation of a method to produce realistic synthetic catchments. The proposed methodology creates virtual catchments by combining the terrain generation approach of Génevaux et al. (2013), geographic information system (GIS) processing procedures, and the CAESAR-Lisflood LEM (Coulthard et al., 2013). The latter was selected following a comparative study of different LEMs, as detailed in Section 2b and 2c below. Aside from size and shape of the catchment, the proposed method allows five key parameters to be adjusted independently (i.e. changing one of the parameters does not change the value of the others): drainage density, hypsometric integral, average slope of the main channel, land use and grain size. The first three are defined using a modified version of the Génevaux et al. (2013) algorithm and additional GIS procedures, whereas the remaining two are adjusted within CAESAR-Lisflood. This combination of approaches allows to produce a fully defined virtual catchment for each of the five physical parameters considered. It is therefore able to guarantee a systematic catchment
generation procedure, ensuring a realistic representation of field conditions, spatial
landscape structures, and modelled flows (Penny et al., 2013; Baartman et al., 2018).
This new method of generating virtual catchments offers great research potential to
identify, in future LEM studies, the impacts of the physical characteristics of
catchments on flow dynamics and sediment flows, and therefore on their
geomorphological evolution.

## 2. Review of methods used to create virtual catchments

Since the advent of computer modelling, LEMs have been used to answer a wide range
of geomorphological questions, often with the aim of gaining an exploratory
understanding of the processes that lead to landscape development (Ahnert, 1977;
Willgoose et al., 1991; Tucker and Slingerland, 1994). These models simulate spatially
and temporally the processes of erosion and sedimentation and thus the dynamic
adjustment of landscape shapes. Currently, the LEMs with the most active user
communities are SIBERIA (Willgoose et al., 1991), CHILD (Tucker et al., 2001),
LAPSUS (Schoorl et al., 2002), and CAESAR-Lisflood (Coulthard et al., 2013). These
LEMs have become laboratories of choice (Tucker and Hancock, 2010) where
hypotheses related to the influence of morphometry, land use, and grain size on
landscape dynamics can be tested. In this section, we review the methods used to
define morphometry (Section 2a), grain size (Section 2b) and land use (Section 2c) for
the virtual catchments used in these LEM analyses. It is important to note that section
2a refers to existing approaches to produce virtual catchments and is independent of
the LEM under consideration. In contrast, Sections 2b and 2c focus on the specific
approaches adopted by the different LEMs listed above to consider the influence of
grain size and land-use in the simulation of water flow and sediment transport. A
comparative analysis of each LEM's specific capabilities is carried out with the aim of
ultimately specifying the most satisfactory LEM to represent land use and grain size in
virtual catchments.

a. Review of existing methods for creating virtual catchments

The morphological complexity of the produced virtual catchments used in LEM
analyses has been addressed in previous studies using different parameters such as
basin relief, soil roughness (Baartman et al., 2013), mean slope of the catchment,
slope profile, drainage length, and catchment area (Zhang and Han, 2017).

Virtual catchments generally are square or rectangular in shape. They are derived from
mathematical expressions producing a valley formed by two symmetrical inclined
planes (Schoorl et al., 2000; Coulthard and Van De Wiel, 2007, 2010; Engida, 2010;
Baartman et al., 2013; Wang, 2013). The valley may include a floodplain, alluvial
terraces, or terraced crops (Baartman et al., 2013) (Figure 1A). Different slope profiles
can be chosen: one with a wide valley bottom or one with a restricted valley bottom
and an extended plateau (Beaujouan et al., 2001) (Figure 1B).

Virtual catchments can also be produced by modifying the digital elevation model
(DEM) of a real-world reference catchment. For example, catchment characteristics
such as average slope, drainage length and catchment shape can be adjusted from a
real-world DEM by changing the elevations and the axial extents of the catchment
through the use of GIS packages (e.g. Zhang et al., 2017; Zhang and Han, 2017;
Coulthard and Van De Wiel, 2017). However, this approach generates results that are
only modified copies of the input DEM. With this method, for example, it is not possible
to influence the shape of the hydrographic network.

No case-testing studies have attempted to examine the influence of the shape of the
river system, although this is a parameter that strongly influences the response of
catchments to precipitation (Carlston, 1963). A manually sketched hydrographic network can eventually be imposed in inclined planes (Thommeret et al., 2010). However, this operation can lead to results that, even after DEM conditioning, are not hydrologically and geomorphologically correct because the density and structure of this arbitrary network does not reflect the topography and microtopography of the catchment.

One alternative solution is to use a LEM to generate a "plausible" topography that would be influenced by the chosen rates of tectonic uplift and fluvial incision (e.g. DeLong et al., 2007; Hobley et al., 2017). In this case, the LEM is used to create an initial catchment topography that can be used for subsequent sensitivity analysis in other LEM simulations. In this approach, fluvial incision is modelled using the stream power law, which expresses the variation in erosion rate as a function of an erosion coefficient K, drainage area A, channel slope S, and constants m and n (equation 1) (Whipple and Tucker, 1999):

$$E = KA^mS^n. \quad (1)$$

The values of the coefficients K, m, and n play a major role in the evolution of the landscape as they control the relationship between the slope of the landscape, the drainage area, and erosion rates, as well as the prevalence of concentrated (fluvial) processes over diffuse (slope) processes. However, the parameters K, m, and n remain unknown constants, and their impacts on the produced relief still raise many questions (Croissant and Braun, 2014). They are based on questionable assumptions about the physics of erosion processes, which significantly limit the use of this model (Beer and Turowski, 2015; Lague et al., 2014). For example, a linear relationship between the erosion rate and the slope of the channel is generally assumed (n = 1), although several studies have demonstrated the existence of incision thresholds (e.g.
195 Harel et al., 2016; Lague et al., 2014). In addition, the application of this law requires
196 very fine spatial and temporal discretisation, which imposes significant computational
197 costs for large-scale landform generation (Tucker and Hancock, 2010).

198 Finally, Ebert et al. (1998) have defined a fractal approach that enables the automatic
199 generation of landscape by combinations of noise functions. The creation of
200 landscapes is very fast with this method: the algorithms are relatively simple and
201 require few parameters adjustment compared to approaches based on physical
202 processes. However, these approaches generate large-scale landscapes that seem
203 unrealistic (geologically fresh looking) and there is a significant lack of control over the
204 positioning of landscapes features, such as rivers and crest lines. Specific methods
205 have been designed to include river networks into the fractal generation procedure
206 (Kelley et al., 1988). However, due to their stochastic nature, these algorithms produce
207 river networks and catchments that are not hydrologically and geomorphologically
208 consistent.

209 The difficulty of efficiently generating large catchment areas that meet
210 geomorphological requirements and the lack of controllability are therefore important
211 problems for most current methods.

212 b. Review of LEM methods for adjusting soil granulometry in virtual
213 catchments

214 Although most LEMs integrate erosion and sedimentation processes into their
215 operations, this is often done in a fairly basic way (Coulthard et al., 2007) (Table 1).

216 For example, SIBERIA and LAPSUS are based on a simplified form of Einstein's
217 equation (1950) that assumes an increase in sediment transport with flow and slope.

218 Within these models, this equation takes the form (equation 2):
219

$$q = \beta_1 Q^{m_1} S^{n_1} \quad (2)$$

where $q$ is the sediment transport rate per unit width (m³ s⁻¹ m⁻¹); $Q$ is the flow rate per unit width (m³ s⁻¹ m⁻¹); $S$ is the slope (m m⁻¹); and $\beta_1$, $m_1$, and $n_1$ are constants (Hancock and Willgoose, 2002). The value of $q$ calculated from equation (2) is used to deduce the mass balance of each cell. Based on this calculation, the elevation of each cell is adjusted, and this adjustment reflects erosion and sediment deposition. However, the sediment transport formulation in SIBERIA and LAPSUS considers only one grain size fraction (Willgoose, 2005). This means that at present, these models are not able to simulate certain sedimentological and morphological properties of river beds related to the selective sorting of sediments. For example, although they are able to simulate the morphology of alluvial cones (Veldkamp et al., 2017), they struggle to reproduce some of the sedimentological features associated with these deposits (e.g. fining of particles from upstream to downstream).

Unlike SIBERIA and LAPSUS, CHILD employs a complete formulation of the Einstein transport equation using sediment density ($\rho_s$), water density ($\rho$), gravitational acceleration ($g$), grain size fraction ($D$), and dimensionless flow by compression ($\varphi$) based on shear stress and flow velocity (equation 3):

236

$$q_i = \phi \sqrt{\frac{(\rho_s - \rho)gD_i^3}{\rho}}. \quad (3)$$

Furthermore, CHILD does not treat the channel bed as a single grain size fraction (Tucker et al., 2001). Indeed, in CHILD, selective sediment transport is possible with the use of a two-fraction approach (sand and gravel) based on Wilcock (1998). However, this two-fraction system is not sufficient to simulate phenomena such as the formation of armour plating at the bottom of the channel, for which selective erosion phenomena are important (Wilcock, 2001).

237

238

239

240

241

242
243 The most detailed representation of sediment granulometry is performed by CAESAR-
244 Lisflood, which simulates the transport of nine granulometric fractions using one of two
245 equations: the Einstein equation (1950) and the Wilcock and Crowe equation (2003)
246 (Van De Wiel et al., 2007). These formulas have different fields of application as they
247 are derived from experiments carried out on different materials. Einstein's equation
248 (1950) was developed from channels with beds mainly composed of grain sizes
249 between 0.785 and 28.65 mm. The Wilcock and Crowe (2003) formula was developed
250 from channel experiments involving sediment mixtures with grain sizes ranging from
251 0.5 to 64 mm.

252 The nine particle size fractions are included within a three-dimensional layer system,
253 thus enabling simulation of the development of stratigraphic units and depositional
254 sequences (Coulthard and Van De Wiel, 2007; Van De Wiel et al., 2007). This
255 configuration also allows the representation of phenomena for which selective erosion
256 processes are important (e.g. development of an armour layer at the bottom of the
257 channel in a limited transport model). CAESAR-Lisflood is also able to consider regolith
258 erodibility in presence of negligible soil cover and particles fracturing. This is based on
259 an experimental soil development model developed by Vanwalleghem et al. (2013).
260 This more complete representation of erosion and selective particle transport
261 processes is a definite advantage of CAESAR-Lisflood over other LEMs (Temme et
262 al., 2017).

263 c. Review of LEM methods for land use adjustment in virtual catchments

264 Vegetation plays a major role in the evolution of the landscape, influencing the
265 hydrological processes that take place in the soil and the degree of cohesion of the
266 soil. From a hydrological point of view, vegetation can, among other things, absorb the
267 kinetic energy of raindrops and reduce surface runoff by promoting infiltration (Brandt
and Thornes, 1987; Rey, 2003). Vegetation also reduces soil moisture through evapotranspiration, which increases the soil's resistance to shear (Matsushi and Matsukura, 2006). From a mechanical point of view, the plant root system improves cohesion of the soil, increasing its shear strength (Preston and Crozier, 1999).

CHILD, SIBERIA, LAPSUS, and CAESAR-Lisflood have been used in a series of studies to assess the effects of vegetation on the topography of a former mining quarry (Evans et al., 1998), on drainage density (Collins et al., 2004; Istanbulluoglu and Bras, 2005), on channel length profile (Collin and Bras, 2010), and on ravine formation (Bastola et al., 2018). In these studies, depending on the LEM used, the effects of vegetation on hydrological processes and soil cohesion were taken into account in different ways.

Limited detachment conditions related to the presence of vegetation can be reproduced in SIBERIA, CHILD, and LAPSUS using the erodibility factor $β_1$ of the sediment transport equation (2). This parameter generates high erosion rates when it takes on high values (Schoorl et al., 2002; Schoorl et al., 2014). LAPSUS has an additional parameter, the sedimentation factor, to reproduce the effect of vegetation on carrying capacity (Schoorl et al., 2002; Schoorl et al., 2014). Low values of this parameter indicate longer-distance transport of sediments, and high values generate large volumes of deposits (provided that the transport capacity has been exceeded).

In CAESAR-Lisflood, the increased resistance of vegetated surfaces can be achieved through a linear plant growth model, although it is mainly a riparian vegetation model.

In summary, the influence of vegetation on soil mechanical properties is generally better taken into account in LAPSUS and less well represented in CAESAR-Lisflood.

However, only LAPSUS and CAESAR-Lisflood, and to a lesser extent CHILD, are able to reproduce the effects of vegetation on flow hydraulics. SIBERIA only focuses on the
293 erodibility of the regolith. The CHILD and CAESAR-Lisflood LEMs use the Manning
294 coefficient in the calculation of flow velocity and shear stress (Istanbulluoglu and Bras,
295 2005). The adjustment of this coefficient can be done spatially by considering a zonal
296 distribution of vegetation. In CAESAR-Lisflood, the influence of vegetation on soil water
297 properties can also be achieved by adjusting a parameter acting on the storage and
298 release of water inside the soil (detailed in section 3f). In LAPSUS, simulation of
299 vegetation effects on soil hydrology is achieved using two parameters that can vary for
300 each land use category: the infiltration capacity (Debolini et al., 2015) and runoff
301 coefficient (Lesschen et al., 2009).

302 CHILD, LAPSUS, and CAESAR-Lisflood can address the effects related to the spatio-
303 temporal evolution of vegetation as the values of their parameters can vary both in time
304 and space. For example, unlike SIBERIA, they can model the growth and senescence
305 of vegetation in different places. The spatio-temporal representation of vegetation is a
306 strong point of LAPSUS, CHILD, and CAESAR-Lisflood as it is known to have a great
307 influence on sediment transport and system connectivity (Foerster et al., 2014;
308 Coulthard and Van De Wiel, 2017; Mishra et al., 2019). However, the annual time step
309 of LAPSUS complicates its application for the analysis of the (very non-linear) effects
310 of vegetation on flow and erosion processes (Baartman et al., 2013; Hancock et al.,
311 2015). A finer time step model such as CHILD or CAESAR-Lisflood is, therefore, better
312 suited to analyse the influence of vegetation on landscape evolution (Temme et al.,
313 2017).

314 Combined, the considerations from Sections 2b and 2c show that the CAESAR-
315 Lisflood model allows the most detailed representation of erosion and deposition
316 phenomena, with complete formulation of the Einstein equation and multiple
317 granulometric fractions. It also is the only LEM that provides a detailed representation
of the influence of vegetation on catchment hydrology in response to a short rain event (Table 1). For these reasons, the CAESAR-Lisflood model is considered to be the most appropriate model to represent the influence of grain size and plant cover in the virtual catchments.

### 3. Generation of realistic virtual catchments with perfectly controllable physical characteristics.

#### a. Methodological approach overview

There are many parameters for characterising the morphometry of a landscape (Hengl and Reuter, 2009), but the morphological complexity of a landscape is mainly governed by the shape of the hydrographic network and the complexity of the relief (Cavalli and Marchi, 2008; Baartman et al., 2013). The method proposed in this study produces virtual catchment areas through five steps (Figure 2). First, a baseline virtual catchment DEM with a desired shape and size is created, and the location of the catchment outlet is specified. In this same step, which follows a modified algorithm of Génevaux et al. (2013), the drainage density is also adjusted. This algorithm for generating the virtual catchment with the desired drainage density is described in section 3b.i. In a second step and third step, the hypsometric integral and the mean slope of the main channel of the virtual catchment are adjusted (3b.ii et 3b.iii). Finally, land use and grain size of the soils of each virtual catchment are then adjusted during their implementation in the CAESAR-Lisflood LEM (3c.i et 3c.ii). Aside from size and shape of the catchment, there thus are five physical parameters that control the morphology of the virtual catchment: drainage density, average slope of the main channel, hypsometric integral, land use, and soil granulometry. The choice of these five parameters is further justified in sections 3b and 3c. An important point is that each parameter can be adjusted independently of the other four, which, in the context of a LEM sensitivity analysis,
343 allows precise examination of the influence of each parameter on the outputs of any
344 subsequent analyses.

345 b. <u>Generation of virtual catchment DEMs</u>

346 i. Drainage density adjustment

347 The baseline for the virtual terrains is generated based on a modification of the method
348 by Génevaux et al. (2013). The catchment size and catchment shape (e.g. elongated,
349 rounded, stunted, axi-symmetrical or curved catchments), and the relief can take any
350 form desired by the user and are controlled by input parameters and two user-sketched
351 slope maps, i.e. a terrain slope map and a river slope map (figure 3). Génevaux et al.'s
352 (2013) algorithm was originally designed to generate several river systems within an
353 island by randomly drawing outlet points around the island. For the current study it was
354 modified to be able to impose a single outlet and thus force all river segments to drain
355 to this outlet.

356 Virtual catchments produced by the proposed method can be any size. Large basins
357 are likely to be driven by marked internal processes ("autogenic processes") related to
358 the storage and transient remobilisation of sediments. This reinforces the inertia and
359 non-linearity of the entire system towards external forcing(s) (Jerolmack and Paola,
360 2010; Van De Wiel and Coulthard, 2010). Thus, in the context of a LEM sensitivity
361 analysis, using large catchments should make it easier to highlight dynamic differences
362 in the morphometry of catchments. To ensure the efficiency of LEM calculation, the
363 modelling of large spatial scales requires the use of grids with relatively low resolution
364 (~30 m or more), as high-resolution data result in a large increase in computing time.
365 Indeed, according to Perron (2011), an increase of the resolution by a factor *n* results
366 in an increase in the calculation time of about *n*³. A low-resolution DEM inevitably omits
certain complexities of the river system, such as the presence of sedimentary bars and islets in the channel or bank collapses (Van De Wiel et al., 2007). However, the main characteristics of the hydrographic network remain visible at low spatial resolution. In addition, it has been shown that LEMs can reliably predict sediment transport even at low grid resolutions (Hancock et al., 2010). Considering these aspects, we herein illustrate our method using DEMs with a surface area of 35 km² and a resolution of 30 m. The method itself, however, can generate virtual catchments of any size and any resolution.

In addition to the key geometric properties of the catchment, Génevaux et al.'s (2013) algorithm also permits to specify and adjust key geometric characteristics of the hydrographic network, namely the drainage pattern (e.g. parallel or trellip) and the drainage density. Shape characteristics of the hydrographic network are commonly used by hydrologists and geomorphologists to study the hydrological cycle, the exchanges between the land surface and the atmosphere, and the response of the land surface to environmental changes (Tarboton et al., 1989; Rodriguez-Iturbe et al., 1992). Many parameters exist to characterise the complexity of the hydrographic network, and many of them are correlated (e.g. drainage density, texture of drainage, frequency of drainage, frequency of first-order sections, and fractal dimension), although some, such as the confluence ratio and the length ratio, do not follow this correlation scheme. However, these parameters were not considered as their range of values remains poorly constrained. It is essential to constrain the parameters of the virtual basins in a range of the natural variation of the used parameters to avoid modelling errors (e.g. aberration of simulated flows). Therefore, the drainage density, a geomorphological measure that characterises the proximity of the channels to each other, was considered herein to define the shape of the hydrographic network.
392 Drainage density is generally calculated as the length of river channels divided by catchment area [L/L²] (Horton, 1945).
393

394 The river slope map (Figure 3b) adjusts the distribution and the structure of the river
395 network (e.g. dendritic, parallel or trellis) and drives the position of plateaus, hillslopes
396 and valleys: rivers and valleys will develop preferentially in areas where specified river
397 slopes are low. Therefore, this map makes it possible to intuitively define how the
398 network will develop without having to sketch out the details of river locations. On the
399 other hand, the terrain slope map (Figure 3a) controls the slope of valleys and the
400 importance of floodplain development: low values will result in gentle cross slopes and
401 significant floodplain development. To create symmetrical dendritic channel networks,
402 the two slope maps have simple symmetrical gradients: the river slope map should be
403 defined to favour the placement of high-altitude areas at the periphery of the domain
404 (high river slope values) and the spreading of the river network from the tip of the
405 domain to areas with low river slope values (Figure 3b). To produce an accentuation
406 of stream incision from the upstream to the downstream part of the catchment, the
407 terrain slope map would need slightly greater values towards the centre thalweg
408 (Figure 3a). The production of virtual catchments with a more orderly drainage network
409 (e.g. trellis, rectangular, parallel, contorted or annular network) would involve the
410 production of more constraining slope maps. In these cases, the user should sketch
411 the river slope map so that the low slope values are shown at the desired locations for
412 the main valleys. For example, for an annular pattern, the river slope map should
413 represent alternating rings of low and high slope values. Similarly, it is quite possible
414 to place the main channel on one side of the catchment so that short, steep tributaries
415 are formed on one side and longer, low-sloped rivers on the other.
416 The algorithm generates from the river slope map a hydrographic network in the form
417 of 3D polylines set. This geometric graph develops from an outlet and colonizes the
418 upstream available space using a probabilistic approach (Génevaux et al., 2013).
419 Network growth is based on a weighted function combining the elevation of candidate
420 network nodes (calculated from the river slope map) and their priority defined from the
421 Horton-Strahler stream classification. The river system expansion process also
422 considers several production rules acting on the minimum allowed spacing between
423 river sections and on the relative frequency of Horton-Strahler streams in the final
424 network. The algorithm then uses the geometry and topology of the generated stream
425 network to divide the domain into Voronoi polygons, each polygon consisting of ridge
426 points, water inflow points and an outflow. The elevation of the ridge points is
427 calculated from the elevation of the neighbouring cells centre and from the terrain slope
428 map which will influence the gradient of the hillslopes. The flow discharge at the outlet,
429 calculated from the catchment area using a power law, is used to classify river sections
430 according to the Rosgen (1994) approach and thus refine the flow paths. Finally, the
431 virtual catchment DEM is generated by combining the river network graph, the Voronoi
432 diagram, flow discharge and Rosgen types assigned to each river node. The
433 roughness of the landscape is also adjusted by a correlated random noise which varies
434 according to the geomorphological context. The construction of the DEM is based on
435 a procedural model combining terrain primitives and allows the creation of hills,
436 mountains and valleys as well as the detailed geometry of streams. This vector model,
437 which allows the elevation to be obtained at any point in the catchment, is finally
438 discretized and converted into a raster data model.

439 The adjustment of the drainage density is done during the creation of the DEM.
440 However, the drainage density is not specified directly as in input parameter in the
algorithm, but rather follows from a parameter defining the minimum allowed spacing between river sections. Hence, the drainage density value is experimentally regulated by a trial and error procedure requiring the subsequent extraction of the DEM hydrographic network from the generated virtual catchment area. Several methods exist to extract the hydrographic networks of DEMs. For low-resolution DEMs (> 10 m), the detection of channel heads based on their geometric characteristics (e.g. plan-view curvature or long profile) is not reliable (Passalacqua et al., 2010; Grieve et al., 2016; Clubb et al., 2014), but the use of a simple contributory area threshold is sufficient. This threshold represents the minimum contributing area necessary to initiate the hydrographic network. The threshold value strongly affects the density of the extracted network. This value cannot be determined from the analysis of the drained slope-area relationship as it has been shown that the slope reversal observed on this type of graph represents the transition from convex to concave topography, not the transition from diffuse to concentrated erosion (Montgomery and Foufoula-Georgiou (1993), Perera and Willgoose (1998) and Hancock (2005). The method representing the cumulative distribution of drainage areas (CAD) is capable of identifying this transition zone and consequently the channel heads (McNamara et al., 2006). The CAD is the probability *P* of crossing a drainage area value at a given point in the catchment area. The probability of crossing *P* is obtained by dividing the number of DEM cells with a drainage area *a* greater than or equal to a given area (*a*\*) by the total number of cells (*N*<sub>T</sub>) (equation 4):

$$P (a \geq a^*) = \frac{\sum_{i=1}^{n} N_i}{N_T} \qquad (4)$$

$$N_i = 0 \text{ if } a_i < a^*, \quad N_i = 1 \text{ if } a_i \geq a^*$$

where $a_i$ is the drainage area of the *i*th pixel.
The CAD has a shape that is common to many catchments. This distribution can be divided into three regions that express different flow regimes (Perera and Willgoose, 1998) (Figure 4D).

Region 1, associated with low drainage areas, represents the part of the catchment dominated by diffuse flows (slopes). Region 2 corresponds to the part of the catchment affected by concentrated flows. This region is described by a straight line in a logarithmic scale and therefore follows a power law. Region 3 consists of the areas near the outlet of the catchment. The discontinuities in this region come from the junctions of tributaries with the main channel. In this region, the probability of crossing decreases rapidly with the increase of the drainage area.

Figure 4 shows the CADs of three virtual catchments that differ only in their drainage density. The catchment with the most branched drainage system has a probability of crossing in Region 1 that decreases slowly with increasing drainage area (Figure 4C). This reflects the fact that a high percentage of the catchment area has a low drainage area. However, the location of the transition between regions 1 and 2 is an indicator of the flow regime and, consequently, of erosion processes in the catchment area. The border between CAD regions 1 and 2 is located to the left, especially because diffuse flow transport is important (Perera and Willgoose, 1998). Log-log linearity (region 2 – power law) is observed for very small drainage areas for the catchment with high drainage density because for this catchment, flows are mainly concentrated. Region 1 is therefore small for this catchment with a slope break that is difficult to identify (Figure 4C).

The river systems can then be derived from the flow accumulation thresholds identified on the CADs (boundary between regions 1 and 2). The very different values of the flow
489 accumulation thresholds from one basin to another obviously result in very different
490 drainage networks in terms of density (Figure 5).

491 The drainage densities of the three catchments generated with our method are
492 respectively low (2.19 km/km², Figure 5A), moderate (5.51 km/km², Figure 5B), and
493 very high (8.35 km/km², Figure 5C) according to Smith's classification (1950).

494 The three hydrographic networks produced all extend to the limit of the catchment
495 although in reality the development of the networks is done both in terms of density
496 and channel length. Thus, the proposed method does not allow considering an areal
497 growth of the network which supposes an increase in channel length with the increase
498 in drainage density.

499 ii. Adjusting the hypsometric integral

500 The hypsometric integral is the second parameter chosen because it is the only
501 parameter that integrates the distribution of altitudes within a catchment area. The
502 hypsometric curve of a catchment represents the relative proportion of the surface area
503 of the catchment above a given altitude, and is widely used in geomorphology to
504 determine the erosive stages of catchments (Strahler, 1952; Schumm, 1956). This
505 area-altitude curve is suitable to analyse the relief distribution in a catchment and to
506 perform comparisons between different catchments. The shape of the hypsometric
507 curve reflects the geomorphological development stage of the catchment. Convex
508 hypsometric curves are characteristic of catchments with young relief and high erosion
509 potential, whereas concave curves indicate a stage of more mature relief and very low
510 erosion potential (Strahler, 1952). The hypsometric integral corresponds to the area
511 under the hypsometric curve (Strahler, 1952; Schumm, 1956) and is also used to
analyse the maturity of the catchment. For example, values close to 1 are associated
with convex curves, and reflect immaturity of the relief.

Harlin (1978) developed an approach that considers the hypsometric curve as a
cumulative probability distribution. It consists of representing the hypsometric curve by
a third-order polynomial function, allowed by the relatively simple shape of the
hypsometric curve and the low number of inflections (equation 5) (Harlin, 1978),

$$
f(x) = \frac{h}{H} = a_0 + a_1 \frac{a}{A} + a_2 \left(\frac{a}{A}\right)^2 + a_3 \left(\frac{a}{A}\right)^3 \quad (5)
$$

519

520 where *h* is the altitude, *H* is the total elevation difference, *a* is the surface area above
521 altitude *h*, and *A* is the total surface area.

522 The hypsometric integral can be calculated by integrating *h/H*. The resolution of this
523 integral consists of using linearity to separate the integral into a sum of integrals
524 (equation 6) (Harlin, 1978):

524

$$
\mathrm{HI} = \int_{0}^{1} a_{0} + a_{1}x + a_{2}x^{2} + a_{3}x^{3} f(x) dx = \sum_{k=0}^{3} \frac{a_{k}}{k+1}. \quad (6)
$$

525

526 The proposed approach to adjust the hypsometric integral is based on a baseline DEM
527 produced in the previous step with the modified version of the Génevaux et al. (2013)
528 algorithm. It contains three main steps:

528

529 (i) The attribute table of a DEM includes altitude and surface values, and it is therefore
530 possible, by a simple operation, to calculate the relative surface area values (a/A).

531 (ii) Once the relative area values are obtained, the relative altitude values (h/H) can be
532 calculated. To do this, it is sufficient to describe the hypsometric curve in the form of
533 the third-order polynomial law given in the equation 6. The adjustment of the

532

533
coefficients $a_0, a_1, a_2$, and $a_3$ of the polynomial function then allows it to act on the
relative altitude values ($h/H$) and thus on the shape of the hypsometric curve and its
integral.

(iii) It is then possible to spatially calculate the altitude data by applying the equations
obtained in the previous step. These equations use as inputs the spatial values of
relative surfaces.

To illustrate the method, we defined three hypsometric curves with very different
shapes and hypsometric integral values (Figure 6D).

The first is strongly concave and has a very low hypsometric integral (HI = 0.25). The
second has a convex shape with a hypsometric integral equal to 0.65. The third has a
stronger convex shape with a hypsometric integral equal to 0.75. The equations
describing these three hypsometric curves can then be used to generate three
catchments showing very different hypsometry (Figures 6A, 6B and 6C). The low
altitudes of catchment area A are very widely represented (low hypsometric integral)
and are representative of a basin for which erosion processes would have been intense
(Strahler, 1952). The altitudes of catchment area B are on average higher because
higher altitudes occupy a larger area. This trend increases for catchment C with high
altitudes very widely represented and low altitudes restricted to the most downstream
areas. Therefore, catchments A and B are characteristic of catchments with young
relief and high erosion potential (Strahler, 1952).

iii. Adjustment of the mean slope of the main channel

There are many parameters to describe the relief and its elevation. Most of them are
correlated, such as the total relief, main channel slope, and average slope of the
terrain. The adjustment of the slope of the main channel was preferred because this
morphometric parameter is one of the first to be considered for the analysis of erosive power of a river (Whipple and Tucker, 1999). The river slope map provided in the first step to create the initial DEM does not allow for easy and accurate adjustment of the main channel mean slope, as its primary purpose is to indicate where the drainage network will be established.

The mean slope of the main channel (MCS) is the length of the longest flow path divided by the difference in elevation between its two ends (Asquith and Slade 1997).

The average slope of the main channel is adjusted by redefining the altitude range of the catchment based on the original topography produced during the adjustment step of the hypsometric integral. This resizing is done without modifying the distribution of values within the DEM (use of a linear mathematical function) so as not to affect hypsometry or previously defined hydrographic properties (equation 7):

$$z = \frac{(x - x_{min}) \times (z_{max} - z_{min})}{(x_{max} - x_{min})} + z_{min} \quad (7)$$

where z represents the new altitude of the catchment at a given point. This value is between a minimum altitude $z_{min}$ and a maximum altitude $z_{max}$; x is the initial altitude of the catchment between $x_{min}$ and $x_{max}$.

Our method allows the automatic adjustment of $z_{max}$ from an expected MCS slope gradient. This adjustment is performed via the following equation (equation 8):

$$z_{max} = \frac{(-x_{min} \times \left(\frac{x_{max} - P_{max}}{x_{max} - x_{min}}\right) + ((Pd \times L) + x_{min}))}{1 - \left(\frac{x_{max} - P_{max}}{x_{max} - x_{min}}\right)} \quad (8)$$

where $P_{max}$ represents the initial maximum altitude of the main channel (m) and L indicates the length of the main channel. The "longest flow path" function available in
579 ArcHydro can be used to estimate the length of the channel. The altitudes at each end
580 of the flow path can then be obtained using the "Extract Multi Values to Points" function.

581 It is important to emphasise that the resizing of altitudes is based on a linear
582 mathematical function which allows to fully preserve the original terrain slope patterns
583 and ratios. This means that virtual catchments with high longitudinal slopes will
584 necessarily have significant transverse slopes. This results in a greater incision in the
585 drainage system and limited development of the floodplain. However, shallow
586 longitudinal slopes will generate significant floodplains that are conducive to the lateral
587 development of watercourses. This correlation between longitudinal and transverse
588 slopes is generally found in river and hydrological system classification schemes
589 (Rosgen, 1994; Winter, 2001). Indeed, Winter (2001) defined two main types of valleys
590 according to their geographical locations: highland and lowland valleys. The former
591 tends to have strong longitudinal gradients with steep slopes and almost non-existent
592 flood plains. In contrast, lowland valleys have relatively low longitudinal gradients with
593 gentle transverse slopes and significant development of floodplains.

594 To illustrate our method, three virtual catchments were produced that differ only in the
595 mean slope of the main channel (Figure 7).

596 The catchment areas are represented with the same scale based on the strongest
597 relief to make the altitude differences visible. The location of the main channel and its
598 ends are also shown on each map. The average slope of the main channel of
599 catchment area A is 0.005 m m⁻¹ which is representative of low-lying areas (Mesquita,
600 2008). That of the second catchment area is two orders of magnitude higher at 0.1 m
601 m⁻¹. Finally, the last catchment area has a main channel slope equal to 0.2 m m⁻¹,
602 which is a value generally observed in mountainous regions (Shen et al., 2016).
603

### c. Adjusting grain size and landuse in the LEM CAESAR-Lisflood

604

605

606

607

608

609

610

611

Finally, as mentioned above, land use and granulometry have many effects on surface processes and thus on landscape evolution (Istanbulluoglu and Bras, 2005; Veldkamp et al., 2017). For example, deforestation or reforestation of surfaces modify sedimentary movements and geomorphological connectivity across the basin (Liébault et al., 2005). Soil granulometry also has a significant impact on sediment production by promoting or not promoting the development of armour layers at the bottom of the channel (Hancock, 2009). This justifies the importance of these two parameters within the proposed method.

612

#### i. Land use adjustment

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

Land use adjustment of the virtual catchments is performed using the parameter *m* and the Manning coefficient of the CAESAR-Lisflood model. The Manning coefficient is used to calculate the depth and velocity of flows and is defined based on land use in the catchment area (see section 2.c). Its value varies between 0.02 (urban areas) and 0.2 (forests) (Lewis et al., 2013). The parameter *m* is based on a modified version of the semi-distributed hydrological model TOPMODEL (Beven and Kirkby, 1979). It influences the characteristics of the modelled flood hydrograph by determining how quickly a catchment responds to rainfall input (Welsh et al., 2009). Large values of *m* increase the water retention capacity of the soil, leading to reduced flood peaks and a slow decline in the hydrograph's recession curve. These values therefore represent well-vegetated catchments (Welsh et al., 2009). Conversely, low *m* values increase soil transmissivity, induce hydrological reactivity of the system to rain, and thus represent catchments with sparse vegetation. A value of 0.02 is commonly used to represent dense forest cover, and a value of 0.005 is specified to represent poorly vegetated areas such as grasslands. Intermediate values represent a full range of
vegetation cover (Beven et al., 1984). An improvement of the hydrological component of the CAESAR-Lisflood model was undertaken by Coulthard and Van De Wiel (2017) to input the spatial variations of the parameter *m* and the Manning coefficient. Our test-case generation method exploits this function by adjusting the land use of a catchment according to its geomorphology. Moreover, numerous studies show close correspondences between landscape forms and land use in various regions of the world (Mesquita, 2008; Segundo et al., 2017; Elhag and Boteva, 2017; Rudi et al., 2018). We propose three land use scenarios showing different location of their forest cover depending on whether it is at the valley bottom, on the slopes, or on the plateaus: (i) in the first scenario, forest cover is concentrated at the bottom of the valley, and the surface is bare for the rest of the catchment area; (ii) in the second scenario, forest is located on the slopes with a total absence of vegetation on the remaining surfaces; and (iii) the last scenario corresponds to the case where the forest is located on the plateaus without any vegetation at the bottom of the valley or on the slopes. The three geomorphological contexts of plateau, valley, and slope were extracted from the DEM using the semi-automatic classification method proposed by Jasiewicz and Stepinski (2013). This method classifies the numerical model into eight geomorphological classes (excluding the "pit" and "flat" categories eliminated during the post-treatment of the DEM; see section 4). This process is based on the concept of geomorphons. Geomorphons are elementary structures that allow different geomorphological contexts to be expressed, and are obtained by comparing the altitude of a given cell with the altitudes of the eight surrounding cells and calculating the zenith and nadir angles for the eight elevation profiles. The eight categories produced by the Jasiewicz and Stepinski method (2013) were then reclassified to generate maps with three major classes: valley, plateau, and slope (Figure 8).
In the example in Figure 8, the plateaus, defined by the “peak” and “ridge” classes, represent only 13% of the area because of their incision by many valleys. The slopes, simple ("slope"), concave ("hollow"), or convex ("spur") from the shoulder ("shoulder") to the foothills ("footslope"), represent a total of 65% of the surface area of the catchment. The valley bottoms are relatively wide and represent 22% of the surface area of the virtual catchment.

## ii. Adjustment of soil granulometry

It is possible to define different soil granulometry scenarios using CAESAR-Lisflood. This LEM allows the representation of sediment erosion, transport, and deposition through nine granulometric classes, included in a superimposed layer system (Van De Wiel et al., 2007). The calculation of the sediment flow $q_i$ for each cell and each grain size fraction $D_i$ is done either with Einstein's formula (1950) (equation 3) or with Wilcock and Crowe's formula (2003) (equation 9):

$$q_i = \frac{F_i U_i^3 W_i^*}{(s - 1)g} \quad (9)$$

where $F_i$ is the fractional volume of the particle size class $i$ in the active layer, $U_i^*$ is the shear rate, $W_i^*$ is a function linking the fractional flow to the total flow, and $s$ is the ratio of sediment density to water density (Van de Wiel et al., 2007). The sediment flow can be restricted by the availability of a given grain size fraction in the topmost surface layer (limited detachment model). The calculation of shear stress ($\tau$) is essential for the application of Einstein's (1950) and Wilcock and Crowe's (2003) formulas. The value of $\tau$ is determined from the density of the water ($\rho$), a drag coefficient ($C_d$), and the velocity of the flows ($u$), integrated based on the runoff height (equation 10):

$$\tau = \rho C_d u^2. \quad (10)$$
Eroded material is transported either as bedload or in suspension (although only the finest fraction can be transported in suspension) and is distributed to neighbouring cells according to the local slope of the ground in the first case and the speed of flow in the second case (Van De Wiel et al., 2007). At the end of each iteration, all the material transported as bedload is deposited in the receiving cell. In addition, the deposition of suspended sediments ($V_{dep}$) is calculated from their falling velocity $v_f$ and their concentration $\kappa$ (equation 11) (Van de Wiel et al., 2007):

$$V_{dep} = \kappa v_f D x^2 dt. \quad (11)$$

The soil granulometry of the catchments is initially homogeneous and is spatially distributed as selective erosion processes begin. The removal of fines from the channel bed during the first iterations generally results in a period of high sediment production (Coulthard et al., 2002).

Figure 9 shows three soil granulometry scenarios that can be implemented for virtual catchments. Scenario 1 produces soils composed of 80% silt and 20% sand (sandy bed river). The granulometry of scenario 2 is represented by depletion of fine particles and enrichment of coarse particles (60% silt; 25% sand, and 15% pebbles – gravel bed river). Scenario 3 follows the same trend with lower and higher proportions of fine and coarse particles respectively. The sediments are composed of 40% silt, 25% sand, 20% pebbles, and 15% cobbles (cobble bed river). In each scenario, particles smaller and larger than 0.63 mm are respectively transported in suspension and as bedload.

## 4. Post-treatment of DEMs in virtual catchments

Before hydrological modelling can be applied to a DEM, it must be adjusted to be hydrologically correct. This post-treatment consists of removing peaks and valleys to allow water to flow freely to the outlet of the catchment.
701

a. Altitude smoothing using the ANUDEM algorithm

702

703

704

705

706

The local soil roughness, controlled by a function for generating coherent spatially auto-correlated noise (Perlin noise), is not a directly and easily adjustable element in the method of Génevaux et al. (2013). It seems exaggerated because random variations in altitude are observed in some areas, leading to the formation of depressions and peaks (Figure 10A).

707

708

709

710

711

712

713

714

715

716

717

718

Areas of incorrect altitude must be automatically detected and corrected, but all this must be done without degrading slope continuity, overall morphology, shape and drainage structure. That would not be possible with a conventional pit removal procedure which would result in the creation of numerous flat areas. There is an interpolation technique specifically designed to create hydrologically correct land surfaces while preserving both the ridge lines and the hydrographic network. This is the ANUDEM (Australian National University Digital Elevation Model) adaptive interpolation and smoothing method developed by Hutchinson (1989). The underlying interpolation algorithm of the ANUDEM programme is described in Hutchinson et al. (2013). The DEM is obtained with an iterative interpolation procedure of the spline type that respects the sudden changes in relief observed at the bottom of valleys. The interpolation procedure follows a multi-grid approach to optimise the execution rate.

719

720

721

722

723

ANUDEM imposes a global drainage network using an algorithm that can consider different input entities (e.g. cliffs, lakes, and drainage networks). The parameter of the ANUDEM method relevant to reducing DEM noise is the discrete error parameter. This parameter controls the smoothing of the surface: the higher the value, the greater the filtering of peaks and valleys, whereas the lower the value, the greater the roughness.

In this study, the default value of this parameter (1) produces satisfactory results.

724
Figure 10 illustrates the smoothing effect, showing that most depressions and peaks are removed with a decrease in roughness while perfectly preserving landforms.

### b. Correction of depressions and flat areas in DEMs

DEMs contain depressions or flat areas that can be artefacts. These elements cause discontinuity and local entrapment of surface runoff during hydrological modelling of catchments, which may not reflect the effect of overall morphology (Petroselli, 2012).

Flat areas are relatively rare in reality, they occur totally randomly and we have no control over their frequency in the catchment and the effects they might have on water and sediment discharges (especially bedload). As a consequence, several algorithms for correcting depressions and flat areas have been published: (i) The Filling method (hereinafter Fill) (Jenson and Domingue, 1988) available in many GIS software packages such as ArcGIS and GRASS, (ii) the TOPAZ method (Martz and Garbrecht, 1999) available in the RICHDEM python tool suite (Barnes et al., 2014)

and (iii) The PEM4PIT method, distributed as an extension of the ArcGIS software (Grimaldi et al., 2007).

A methodological comparison of these three methods was carried out by Fernandez et al. (2016) to evaluate their effectiveness in correcting flat areas present in DEMs. In this study, the results were analysed in relation to the extent of the area affected by the correction, changes in elevation and slope, and the distribution of flow velocities.

The results show the disadvantages of the standard Fill correction method. For advanced methods (TOPAZ and PEM4PIT), the correction processes are strongly influenced by the relief, sizes of the depressions, and their distribution (Fernandez et al., 2016). There are no methods that work optimally for all contexts. However, Fernandez et al. (2016) highlighted the greater versatility of the TOPAZ tool, which is
effective in regions with very irregular terrain and in regions with low slopes. We have
therefore chosen this tool for the correction of our virtual catchments DEM.

The TOPAZ tool was therefore applied to fill the depressions using an approach where
the propagation of flows is limited to the cardinal directions (east, north, west, and
south; i.e. a D4 regime) as seen in most LEMs. Indeed, with a traditional D8 method
(flows in the cardinal and diagonal directions), depressions would remain and the
continuity of flows in the CAESAR-Lisflood LEM would not be as well-ensured.

Although these treatment procedures improve the quality of the DEM, they are likely
to cause a slight change in drainage density and hypsometric integral. However, these
changes remain minor and do not substantially affect the properties of the virtual
catchment.

## 5. Results

Figure 11 compares the Austreberthe catchment in Normandy, France (Figure 11C)
with a virtual catchment produced by our method (Figure 11B).

In order to make the comparison between the two catchments more explicit, the
drainage density, hypsometric integral and mean slope of the main channel were set
at values commonly observed in the Normandy regional context (Laignel, 2003;
Lequien, 2006; Mesquita, 2008; El Janyani, 2014). The virtual catchment generated by
our method is visually similar to the morphology of the real catchment. The two river
systems are very similar in organisation: they are built around continuous ridge lines
sharing several sub-catchments. Both have a dendritic structure for which a
hierarchical distribution of the main channel and its tributaries is observed. The
distribution of the valley, plateau and hillslope areas follow a coherent pattern
consistent with that observed for the real catchment. Our method generates
773 catchments whose relief appears to be shaped by erosion and weathering. These
774 characteristics and this type of organisation have not been found in previous virtual
775 catchments (e.g. V-shaped valley catchments, Figure 11A).

776 However, the terrain generated by our method still lacks continuity between valley and
777 plateau areas, which explains the prominence of ridge lines, whether at the head of
778 the catchment area or further downstream (Génevaux et al., 2013). It would be
779 possible to prevent the creation of these huge ridges by limiting the maximum local
780 slope to a threshold value. The narrowing of the valley section commonly observed
781 with the decrease in drainage area is not or only slightly visible for our catchment (Zink
782 et al., 2012). On the other hand, the hydrographic network of the virtual catchment
783 appears somewhat less branched and consists of numerous tributaries. Its bifurcation
784 ratio (3.90) is higher than that of the Austreberthe catchment (3.44), indicating that low-
785 order streams are more frequent in the river system. In addition, the stream length ratio
786 is higher for the virtual catchment (2.30 versus 1.84), suggesting that high-order
787 streams are proportionately longer. However, the values of drainage densities,
788 hypsometric integrals, and average slopes of the main channels of the two catchment
789 basins remain very close. Indeed, according to Smith's (1950) classification, both
790 catchments have a low drainage density (~2km/km²). Their relatively high hypsometric
791 integrals (~0.7) reflect convex hypsometric curves. Thus, high altitudes are
792 predominantly represented and low altitudes are restricted to downstream areas.
793 Finally, main channel slopes are both relatively low (~0.005m/m) potentially resulting
794 in long water and sediment transit times.

795 On the other hand, land use was defined in CAESAR-Lisflood based on the
796 geomorphology of the virtual catchment (Figure 8) since there is a strong landform-
797 land use relationship in Normandy: indeed, most of the vegetation is concentrated on
the hillslopes (Mesquita, 2008). A homogeneous grain size representative of the soils of the region (80% silts and 20% fine sands) was also defined in CAESAR-Lisflood to finalize the creation of the virtual catchment.

## 6. Discussion

This comparison confirms that the proposed method can successfully generate coherent and plausible virtual catchments on a large scale, while offering a high degree of control (e.g. the density of the hydrographic network is controlled by the user through the adjustment of a single parameter). We can control the structure of the terrain and the hydrographic network through simple and intuitive parameters. We are able to define the location of plateau and valley areas, the location of the catchment outlet, the location and trajectory of channels through the definition of intuitive slope maps.

The strength of our method is to allow the separate adjustment of the drainage density, hypsometric integral, mean slope of the main channel, land use and soil granulometry in order to analyse precisely their influence on the evolution of the landscape. All the virtual catchments produced by our method are therefore designed assuming independence between the parameters considered. However, in reality, correlations exist between each of these parameters but also with other parameters. Indeed, drainage density is influenced by factors such as soil permeability (influenced by soil granulometry), climate, vegetation, length and average slope of the slopes, etc. In general, low drainage density values are observed in areas based on highly permeable materials (generally coarse grading) with low vegetation cover and low relief (Devlin, 2015). While high drainage density values indicate areas with impermeable soils with sparse vegetation, mountainous relief (Nautiyal, 1994; Lazzaro et al., 2015) and fine grading (Devlin, 2015). Strahler (1952) also found that the hypsometric integral is inversely correlated with drainage density and channel slope. Separate adjustment of
the different parameters can therefore lead us to produce unrealistic catchments.
Thus, in addition to ensuring that each parameter is constrained within a range of
values encountered in reality (Smith, 1950; Strahler et al., 1952; Shen et al., 2016),
the user must verify the plausibility of a combination of parameters. In this sense, a
bibliographical work is a preliminary step in the design of virtual catchments with the
proposed method.

It should be also noted that the imposed drainage density slightly affects the shape of
the catchment initially defined by the two slope maps. Indeed, the algorithm relies on
the geometry of the river network to generate the Voronoi diagram and thus to calculate
the drained areas. Catchments with low drainage density tend to deviate further from
the shape initially defined since their Voronoi diagram have a coarser structuring. It is
therefore impossible with the proposed approach to generate a set of catchments that
have exactly the same shape, and that really only differ in density. However, the
differences remain insignificant since the overall shape of the catchment area is
completely preserved.

7. Conclusion

This work proposes a new method for generating virtual catchments to simulate hydro-
sediment transfer and landscape evolution, leading to a high level of control over the
generation process. Our method consists of different components for the individual
adjustment of five physical characteristics: drainage density, hypsometry, average
slope of the main drain, granulometry, and land use. The generation of the virtual
catchments is given by coupling a modified version of the field generation method of
Génevaux et al. (2013) with GIS processing procedures and the CAESAR-Lisflood
LEM. As with any catchment generation method, parameter adjustment requires a
certain level of knowledge to find a consistent and plausible set of values. The values
presented in this paper show that our method can be used to produce realistic virtual catchments efficiently, simply, and with a high degree of control, unlike standard approaches that produce unrealistic V-tilt type catchments. The difference shown between Figure 11A and Figure 11B highlights the significant impact of the proposed approach and shows that the new method provides catchments with more hydrologically correct drainage networks than previous methods. Hence, our method seems to be effective in producing virtual catchments that can be used in the future to precisely examine the influences of different physical parameters of catchments on hydro-sediment dynamics. This illustrates the diversity of issues that can be addressed later through this work, such as the respective roles of landscape evolution and climate variability on the evolution over time of the hydro-sedimentary response and its non-linear nature (feedback).

## Acknowledgments

This work was financed by the Seine Normandy Water Agency (AESN) and also obtained financial support from the OZCAR Research Infrastructure (French Critical Zone Observatories, Research and Applications) through the Réseau de Bassins Versants (RBV).

We thank Sara J. Mason, M.Sc., from Edanz Group (https://en-author-services.edanzgroup.com/) for editing a draft of this manuscript. We also are grateful for the reviewers' constructive comments, which improved the manuscript.

## Data Availability Statement
871 The data sets used and/or analyzed during the current study are available from the

872 corresponding author on reasonable request.
874 Ahnert, F., 1977. Some comments on the quantitative formulation of geomorphological processes in a theoretical model. Earth Surf. Process., 2, 191-201.

876 Asquith, W.H., and Slade, R.M., Jr., 1997. Regional equations for estimation of peak-streamflow frequency for natural basins in Texas: U.S., Geological Survey Water-Resources Investigations Report, 96-4307, 68p.

879 Baartman, J.E.M., Masselink, R., Keesstra, S.D., and Temme, J.A.M., 2013. Linking k landscape morphological complexity and sediment connectivity, Earth Surf. Process. 881 Landf., 38, 1457-1471.

882 Baartman, J.E.M., Temme, J.A.M., and Saco, P.M., 2018. The effect of landform variation on vegetation patterning and related sediment dynamics, Earth Surf. Process. 884 Landf., 43, 2121-2135.

885 Barnes, R., Lehman, C., MullA, D., 2014. Priority-flood: an optimal depression-filling and watershed-labeling algorithm for digital elevation models, Comput. Geosci., 62, 17-127.

888 Bastola, S., Dialynas, Y.G., Bras, R.L., Noto, L.V., and Istanbulluoglu, E., 2018. The 889 The role of vegetation on gully erosion stabilization at a severely degraded landscape: 890 A case study from Calhoun Experimental Critical Zone Observatory, Geomorphology, 891 308, 25-39.

892 Beaujouan, V., Durand, P., Ruiz, L., 2001. Modelling the effect of the spatial distribution 893 of agricultural practices on nitrogen fluxes in rural catchments, Ecological Modelling, 894 137, 93-105.

895 Beer, A.R., and Turowski, J.M., 2015. Bedload transport controls intra-event bedrock 896 erosion, Earth Surf. Dynam. Discuss., 3, 1-30.

897 Beven, K.J., and Kirkby, M.J., 1979. A physically based, variable contributing area 898 model of basin hydrology, Hydrological Sciences Bulletin 24, 43-69.

899 Beven, K.J., Kirby, M.J., Schofield, N. and Tagg, A.F., 1984. Testing a physically- 900 based flood forecasting model (TOPMODEL) for three U.K. catchments, J. Hydrol., 69, 901 119-143.

902 Bishop, P., 2007. Long-term landscape evolution: linking tectonics and surface 903 processes, Earth Surf. Process. and Landf., 32, 329-365.

904 Brewer, P.A., and Lewin, J., 1998. Planform cyclicity in an unstable reach: complex 905 fluvial response to environmental change, Earth Surf. Process. Landf., 23, 989-1008.

906 Brandt, J., Thornes, J., 1987. Erosional energetics, Energetics of Physical 907 Environment, 51-87.

908 Carlston, C. W., 1963. Drainage Density and Streamflow, U.S. Geol. Surv. Prof. Pap. 909 422-C.
910 Cavalli, M., and Marchi, L., 2008. Characterisation of the surface morphology of an
911 alpine alluvial fan using airborne LiDAR, Nat. Hazards Earth Syst. Sci., 8, 323-333.

912 Clubb, F.J., Mudd, S.M., Milodowski, D.T., Hurst, M.D., and Slater, L.J., 2014.
913 Objective extraction of channel heads from highresolution topographic data, Water
914 Resour. Res., 50, 4283-4304.

915 Collins, D.B.G., Bras, R.L., and Tucker, G.E., 2004. Modeling the effects of vegetation-
916 erosion coupling on landscape evolution, Journal of geophysical research, 109, 1-11.

917 Collins, D.B.G., and Bras, R.L., 2010. Climatic and ecological controls of equilibrium
918 drainage density, relief, and channel concavity in dry lands, Water resources research,
919 46, 1-18.

920 Coulthard, T.J., Macklin, M.G., and Kirkby, M.J., 2002. A cellular model of Holocene
921 upland river basin and alluvial fan evolution, Earth Surf. Process. Landf., 27, 269-288.

922 Coulthard, T.J., Hicks, D.M., and Van De Wiel, M.J., 2007. Cellular modelling of river
923 catchments and reaches: advantages, limitations and prospects, Geomorphology, 90,
924 192-207.

925 Coulthard, T.J., Neal, J.C., Bates, P.D., Ramirez, J., de Almeida, G.A.M., and
926 Hancock, G.R., 2013. Integrating the LISFLOOD-FP 2D hydrodynamic model with the
927 CAESAR model: implications for modelling landscape evolution, Earth Surf. Process.
928 Landf., 38, 1897-1906.

929 Coulthard, T.J., and Van De Wiel, M.J., 2007. Quantifying fluvial non linearity and
930 finding self-organized criticality? Insights from simulations of river basin evolution,
931 Geomorphology, 91, 216-235.

932 Coulthard, T.J., and Van De Wiel, M.J., 2012. Modelling river history and evolution,
933 Phil. Trans. R. Soc. A, 370, 2123-2142.

934 Coulthard, T.J., and Van De Wiel, M.J., 2017. Modelling long term basin scale
935 sediment connectivity, driven by spatial land use changes, Geomorphology, 277, 265-
936 281.

937 Croissant, T., and Braun, J., 2014. Constraining the stream power law: a novel
938 approach combining a landscape evolution model and an inversion method, Earth Surf.
939 Dynam., 2, 155-166.

940 Debolini, M., Schoorl, J.M., Temme, A., Galli, M., and Bonari, E., 2015. Changes in
941 agricultural land use affecting future soil redistribution patterns: a case study in
942 southern Tuscany (Italy), Land degradation & development, 26, 574-586.

943 DeLong, S.B., Pelletier, J.D., and Arnold, L., 2007. Bedrock landscape development
944 modeling: Calibration using field study, geochronology, and digital elevation model
945 analysis, Geological Society of America, 119, 157-173.

946 Devlin, J.F., 2015. HydrogeoSieveXL: An Excel-based tool to estimate hydraulic
947 conductivity from grain-size analysis, Hydrogeology Journal, 23, 837-844.
948 Ebert, D., Musgrave, K., Peachey, D., Perlin, K., and Worley, S., 1988. Texturing and Modeling: A Procedural Approach. Academic Press Professional.

950 Einstein, H.A., 1950. The bed-load function for sediment transport on open channel flows. Tech.Bull., vol. 1026. USDA, Soil Conservation Service, p. 71.

952 Elhag, M., Boteva, S., 2017. Understanding of the geomorphological elements in discrimination of typical mediterranean land cover types, IOP Conf. Ser.: Earth Environ. Sci, 95.

955 El Janyani, S., 2013. Incidence des bétoires et de la karstogénèse des plateaux crayeux de la Haute-Normandie sur le fonctionnement hydrologique de l'équifère de la craie : Modélisation hydrogéologique des influences climatiques à différentes échelles spatio-temporelles. PhD thesis, University of Rouen.

959 Evans, K.G., Willgoose, G.R., Saynor, M.J., House, T., 1998. Effect of vegetation and surface amelioration on simulated landform evolution of the post-mining landscape at ERA Ranger Mine, Northern Territory, Supervising Scientist Report 134, Supervising Scientist, Canberra.

963 Engida, A.N., 2010. Hydrological and suspended sediment modeling in the lake Tana Basin, Ethiopia, PhD Thesis, Université Joseph-Fourier – Grenoble I.

965 Fernandez, A., Adamowski, J., and Petroselli, A., 2016. Analysis of the behavior of three digital elevation model correction methods on critical natural scenarios, Journal of Hydrology: Regional Studies, 8, 304–315.

968 Foerster, S., Wilczok, C., Brosinsky, A., Segl, K., 2014. Assessment of sediment connectivity from vegetation cover and topography using remotely sensed data in a dryland catchment in the Spanish Pyrenees, Journal of Soils and Sediments, 14, 1982-2000.

972 Génevaux, J.-D., Galin, E., Guérin, E., Peytavie, A., and Benes, B., 2013. Terrain generation using procedural models based on hydrology, ACM Trans. Graph., 32.

974 Grieve, S.W.D., Mudd, S.M., Milodowski, D.T., Clubb, F.J., and Furbish, D.J., 2016. How does grid-resolution modulate the topographic expression of geomorphic processes? Earth Surf. Dynam., 4, 627–653.

977 Grimaldi, S., Nardi, F., Di Benedetto, F., Instanbulluoglu, E., Bras, R.L., 2007. A physically based method for removing pits in digital elevation models, Adv. Water Resour., 30, 2115–2158.

980 Hancock, G.R., and Willgoose, G.R., 2002. The use of a landscape simulator in the validation of the SIBERIA landscape evolution model: transient landforms, Earth Surf. Process. Landf., 27, 1321-1334.

983 Hancock, G.R., 2009. A catchment scale assessment of increased rainfall and storm intensity on erosion and sediment transport for Northern Australia, Geoderma, 152, 350-360.
986 Hancock, G.R., Lowry, J.B.C., Coulthard, T.J., Evans, K.G., and Moliere, D.R., 2010. A catchment scale evaluation of the SIBERIA and CAESAR landscape evolution models, Earth Surf. Process. Landf., 35, 863-875.

989 Hancock, G.R., and Coulthard, T.J., 2012. Channel movement and erosion response to rainfall variability in southeast Australia, Hydrol. Process., 26, 663-673.

991 Hancock, G.R., Lowry, J.B.C., and Coulthard, T.J., 2015. Catchment reconstruction – erosional stability at millennial time scales using landscape evolution models, Geomorphology, 231, 15-27.

994 Harel, M.-A., Mudd, S.M., and Attal, M., 2016. Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184-196.

997 Harlin, J.M., 1978. Statistical moments of the hypsometric curve and its density function, Mathematical Geology, 10, 59-72.

999 Hengl, T., Reuter, H., 2009. Geomorphometry: Concepts, Software, Applications, Elsevier, Developments in Soil Science, 33, 796p.

1001 Hobley, D.E.J., Adams, J.M., Nudurupati, S.S., Hutton, E.W.H., Gasparini, N.M., Istanbulluoglu, E., and Tucker, G.E., 2017. Creative computing with Landlab: an open-source toolkit for building, coupling, and exploring two-dimensional numerical models of Earth-surface dynamics, Earth Surf. Dynam., 5, 21-46.

1005 Horton, R.E., 1945. Erosional development of streams and their drainage basins: hydrophysical approach to quantitative morphology, Bulletin of the Geological Society of America, 56, 275-370.

1008 Hutchinson, M.F., 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits, Journal of Hydrology, 106, 211-232.

1010 Hutchinson, M.F., Stein, J.L., Gallant, J.C., and Dowling, T.I., 2013. New Methods for Incorporating and Analysing Drainage Structure in Digital Elevation Models, Geomorphometry, Nanjing, 2013.

1013 Istanbulluoglu, E., and Bras, R.L., 2005. Vegetation-modulated landscape evolution: Effects of vegetation on landscape processes, drainage density, and topography, Journal of Geophysical research, 110, 1-19.

1016 Jasiewicz, J., and Stepinski, T.F., 2013. Geomorphons — a pattern recognition approach to classification and mapping of landforms, Geomorphology, 182, 147-156.

1018 Jenson, S.K., and Domingue, J.O., 1988. Extracting topographic structure from digital elevation models. Photogramm, Engng Remote Sensing, 54, 1593-1600.

1020 Jerolmack, D.J., and Paola, C., 2010. Shredding of environmental signals by sediment transport, Geophysical Research Letters, 37, 1-5.

1022 Kelley, A., Malin, M., and Nielson, G., 1988. Terrain simulation using a model of stream erosion, Computer Graphics, 22, 263-268.
1024 Lague, D., 2014. The stream power river incision model: evidence, theory and beyond, Earth Surf. Process. Landf. 39, 38-61.

1026 Laignel, B., 2003. Caractérisation et dynamique érosive de systèmes géomorphologiques continentaux sur substrat crayeux. Exemple de l'Ouest du Bassin de Paris dans le contexte nord-ouest européen. HDR Thesis, University of Rouen, 138p.

1030 Lequien, A., 2006. Dynamique erosive actuelle et transferts fluviatiles dans l'Ouest du bassin de Paris : exemple des bassins versants littoraux en Haute-Normandie : le Dun, la Ganzeville et l'Yères. PhD thesis University of Rouen. 160p.

1033 Lazzaro, M.D., Zarlenga, A., Volpi, E., 2015. Hydrological effects of within-catchment heterogeneity of drainage density, Advances in Water Resources, 76, 157-167.

1035 Lesschen, J.P., Schoorl, J.M., and Cammeraat, L.H., 2009. Modelling runoff and erosion for a semi-arid catchment using a multi-scale approach based on hydrological connectivity, Geomorphology, 109, 174-183.

1038 Lewis, M., Bates, P., Horsburgh, K., Neal, J., and Schumann, G., 2013. A storm surge inundation model of the northern Bay of Bengal using publicly available data, Q. J. R. Meteorol. Soc., 139, 358-369.

1041 Liébault, F., Gomez, B., Page, M., Marden, M., Peacock, D., Richard, D., Trotter, C.M., 2005. Land-use change, sediment production and channel response in upland regions. River Res. Appl., 21, 739-756.

1044 Lin, Z., and Oguchi, T., 2011. Analyses of watershed longitudinal/transverse profiles and stream-net structure using high-resolution DEMs, Geomorphometry, Redlands, CA, 2011.

1047 Macklin, M.G., and Lewin, J., 1993. Holocene river alluviation in Britain, Z. Geomorph. Suppl., 88, 109-122.

1049 Martz, L.W., and Garbrecht, J., 1999. An outlet breaching algorithm for the treatment of closed depressions in a raster DEM, Comput. Geosci., 25, 835-844.

1051 Matsushi, Y., and Matsukura, Y., 2006. Cohesion of unsaturated residual soils as a function of volumetric water content, Bull. Eng. Geol. Env., 65, 449-455.

1053 McNamara, J.P., Ziegler, A.D., Wood, S.H., and Vogler, J.B., 2006. Channel head locations with respect to geomorphologic thresholds derived from a digital elevation model: A case study in northern Thailand, Forest Ecology and Management, 224, 147-156.

1057 Mesquita, J., 2008. Facteurs de contrôle multi-échelles des modalités de réponse hydrologique de bassin versant en substrat carbonaté. Comparaison Haute-Normandie (FR) et Région d'Austin Texas (US). PhD Thesis, University of Caen, 264p.

1060 Mishra, K., Sinha, R., Jain, V., Nepal, S., Uddin, K., 2019. Towards the assessment of sediment connectivity in a large Himalayan river basin, Science of the Total Environment, 661, 251-265.
1063 Nautiyal, M.D., 1994. Morphometric Analysis of A Drainage Basin Using Aerial Photographs: A Case Study of Khairkuli Basin, District Dehradun, U.P, Journal of the Indian Society of Remote Sensing, 22, 4.

1066 Nicholas, A.P., and Quinne, T., 2007. Modelling alluvial landform change in the absence of external environmental forcing, Geology, 35, 527-530.

1068 Notebaert, B., Verstraeten, G., Ward, P., Renssen, H., and Van Rompaey, A., 2011. Modeling the sensitivity of sediment and water runoff dynamics to Holocene climate and land use changes at the catchment scale, Geomorphology, 126, 18-31.

1071 Passalacqua, P., Do Trung, T., Foufoula-Georgiou, E., Sapiro, G., and Dietrich, W. E., 2010. A geometric framework for channel network extraction from lidar: Nonlinear diffusion and geodesic paths, J. Geophys. Res.-Earth, 115, F01002.

1074 Penny, G.G., Daniels, K.E., and Thompson, S.E., 2013. Local properties of patterned vegetation: quantifying endogenous and exogenous effects, Phil. Trans. R. Soc. A, 371, 1-29.

1077 Perera, H., and Willgoose, G., 1998. A physical explanation of the cumulative area distribution curve, Water Resources Research, 34, 1335-1343.

1079 Perron, J.T., Numerical methods for nonlinear hillslope transport laws, Journal of Geophysical Research, 116, 1-13.

1081 Petroselli, A., 2012. LIDAR Data and Hydrological Applications at the Basin Scale, GIScience & Remote Sensing, 49, 139-162.

1083 Preston, N.J., and Crozier, M.J., 1999. Resistance to shallow land slide through root-derived cohesion in East coast hill country soils, North Island, New Zealand, Earth Surf. Process. Landf., 24, 605-676.

1086 Rey, F., 2003. Influence of vegetation distribution on sediment yield in forested marly gullies, Catena, 50, 549-562.

1088 Rosgen, D.L., 1994. A classification of natural rivers. Catena, 22, 169-199.

1089 Rudi, G., Bailly, J.S., Vinatier, F., 2018. Using geomorphological variables to predict the spatial distribution of plants species in agricultural drainage networks, PLoS ONE, 13, 1-20.

1092 Schoorl, J.M., Sonneveld, M.P.W., and Veldkamp, A., 2000. Three-dimensional landscape process modelling: the effect of dem resolution, Earth Surf. Process. Landf., 25, 1025-1034.

1095 Schoorl, J.M., Veldkamp, A., and Bouma, J., 2002. Modeling water and soil redistribution in a dynamic landscape context. Soil Science Society of America Journal, 66, 1610-1619.

1098 Schoorl, J.M., Temme, A.J.A.M., and Veldkamp, T., 2014. Modelling centennial sediment waves in an eroding landscape - catchment complexity, Earth Surf. Process. Landf., 39, 1526-1537.
1101 Schumm, S.A., 1956. Evolution of Drainage Systems and Slopes in Badlands at Perth
1102 Amboy, New Jersey, Geological Society of America Bulletin, 67, 597-646.

1103 Schumm, S.A., and Parker, R.S., 1973. Implications of complex response of drainage
1104 systems for Quaternary alluvial stratigraphy, Nature Physical Science, 243, 99-100.

1105 Segundo M. I. G., Bocco, G., Velasquez, A., Gajewski, K., 2017. On the relationship
1106 between landforms and land use in tropical dry developing countries. A GIS and
1107 multivariate statistical approach, Investigaciones Geograficas, 93.

1108 Shen, X., Anagnostou, E.N., Mei, Y., Hong, Y., 2016. Data Descriptor; A global
1109 distributed basin morphometric dataset, Scientific data, 4:160124.

1110 Smith, K.G., 1950. Standards for grading texture of erosional topography, Am. J. Sci.,
1111 248, 655-668.

1112 Strahler, A., 1952. Dynamic Basis of Geomorphology, Geological Society of America
1113 Bulletin, 63, 923-938.

1114 Straub K.M., Duller R.A., Foreman B.Z., Hajek E.A., 2020. Buffered, Incomplete, and
1115 Shredded: The Challenges of Reading an Imperfect Stratigraphic Record, Journal of
1116 Geophysical Research: Earth Surface, 125(3), 1-44.

1117 Temme, A.J.A.M., Armitage, J., Attal, M., van Gorp, W., Coulthard, T.J., and Schoorl,
1118 J.M., 2017. Developing, choosing and using landscape evolution models to inform
1119 field-based landscape reconstruction studies, Earth Surf. Process. Landf., 42, 2167-
1120 2183.

1121 Thommeret, N., Bailly, J.S., and Puech, C., 2010. Extraction of thalweg networks from
1122 DTMs: application to Badlands, Hydrol. Earth. Syst. Sci., 14, 1527-1536.

1123 Tucker, G.E., and Slingerland, R.L., 1994. Erosional dynamics, flexural isostasy, and
1124 long-lived escarpments: A numerical modelling study, Journal of Geophysical
1125 Research, 99, 12229-12243.

1126 Tucker, G.E., Lancaster, S.T., Gasparini, N.M., Bras, R.L., and Rybarczyk, S.M., 2001.
1127 An object-oriented framework for distributed hydrologic and geomorphic modelling
1128 using triangular irregular networks, Comput. Geosci., 27, 959-973.

1129 Tucker, G.E., and Hancock, G.R., 2010. Modelling landscape evolution, Earth Surf.
1130 Process. Landf., 35, 28-50.

1131 Van Der Wall, J.A., and Ssegane, H., 2013. Do polynomials adequately describe the
1132 hypsometry of monadnock phase watersheds? Journal of the American Water
1133 Resources Association, 49, 1485-1495.

1134 Van De Wiel, M.J., Coulthard, T.J., Macklin, M.G., and Lewin, J., 2007. Embedding
1135 reach-scale fluvial dynamics within the CAESAR cellular automaton landscape
1136 evolution model, Geomorphology, 90, 283-301.

1137 Van De Wiel, M.J., and Coulthard, T.J., 2010. Self-organized criticality in river basins:
1138 Challenging sedimentary records of environmental change, Geological Society of
1139 America, 38, 87-90.
1140 Vanwalleghem, T., Stockmann, U., Minasny, B., and McBratney, A. B., 2013. A
1141 quantitative model for integrating landscape evolution and soil formation, J. Geophys.
1142 Res.-Earth, 118, 331-347.

1143 Veldkamp, A., Baartman, J.E.M., Coulthard, T.J., Maddy, D., Schoorl, J.M., J.E.A.,
1144 Storms, A.J.A.M., Temme, van Balen, R., van De Wiel, M.J., van Gorp, W., Viveen,
1145 W., Westaway, R., and Whittaker, A.C., 2017. Two decades of numerical modelling to
1146 understand long term fluvial archives: Advances and future perspectives, Quaternary
1147 Science Reviews, 166, 177-187.

1148 Wang, Y., 2013. Simulating complex hydro-geomorphic changes in lake-catchment
1149 systems, PhD Thesis, University of Southampton.

1150 Welsh, K.E., Dearing, J.A., Chiverrell, R.C., and Coulthard, T.J., 2009. Testing a
1151 cellular modelling approach to simulating late-Holocene sediment and water transfer
1152 from catchment to lake in the French Alps since 1826, Holocene, 19, 785-798.

1153 Wilcock, P.R., 1998. Two-fraction model of initial sediment motion in gravel bed rivers,
1154 Science, 280, 410-412.

1155 Wilcock, P.R., 2001. Toward a practical method for estimating sediment transport rates
1156 in gravel-bed rivers, Earth Surf. Process. Landf., 26, 1395-1408.

1157 Wilcock, P.R., and Crowe, J.C., 2003. Surface-based transport model for mixed-size
1158 sediment, Journal of Hydraulic Engineering, 129, 120-128.

1159 Willgoose, G., Bras, R.L., and Rodriguez-Iturbe, I., 1991. Results from a new model of
1160 river basin evolution, Earth Surf. Process. Landf., 16, 237-254

1161 Willgoose, G., 2005. Mathematical Modeling of Whole Landscape Evolution, Annual
1162 Review of Earth and Planetary Sciences, 33, 443-459.

1163 Winter, T. C., 2001. The Concept of Hydrologic Landscapes. Journal of the American
1164 Water Resources Association, 37, 335-349.

1165 Wipple, K.X., and Tucker, G.E., 1999. Dynamics of the stream-power river incision
1166 model: Implications for height limits of mountain ranges, landscape response
1167 timescales, and research needs, Journal of Geophysical research, 104, 17661-17674.

1168 Zhang, J., and Han, D., 2017. Catchment Morphing (CM): a novel approach for runoff
1169 modeling in ungauged catchments, Water Resources Research, 53, 10899-10907.

1170 Zhang, J., Dawei, H., Song, Y., and Dai, Q., 2017. Exploration of virtual catchments
1171 approach for runoff predictions of ungauged catchments, Hydrol. Earth Syst. Sci.
1172 Discuss.

1173 Zink, J.M., Jennings, G.D., and Price, A., 2012. Morphology characteristics of
1174 Southern Appalachian Wilderness Streams. Journal of the American Water Resources
1175 Association, 48, 762-773.

1176

1177
1178 Figure 1. Examples of virtual catchments used in numerical case-testing studies. (A)

1179 Baartman et al (2013) and (B) Beaujouan et al (2001).
<table><thead><tr><th></th><th>SIBERIA</th><th>CHILD</th><th>LAPSUS</th><th>CAESAR-Lisflood</th></tr></thead><tbody><tr><td>Sediment transport equation</td><td>Simplified Einstein (1950)</td><td>Einstein (1950)</td><td>Simplified Einstein (1950)</td><td>Einstein (1950) or Wilcock and Crowe (2003)</td></tr><tr><td>Number of grain sizes</td><td>One</td><td>Two (sand and gravel)</td><td>One</td><td>Nine</td></tr><tr><td>Vegetation influence on soil mechanical properties?</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Vegetation influence on soil hydrological properties?</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Spatial and temporal evolution of land use?</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Temporal resolution</td><td>Annual</td><td>Dynamic; can be as low as seconds</td><td>Annual</td><td>Dynamic; can be as low as seconds</td></tr></tbody></table>

1180

<caption>1181 Table 1: Technical comparison of selected attributes of four LEMs</caption>
1182 **Figure 2.** Overview of the proposed virtual catchments generation method. (A)

1183 Generation of topography using a modified version of the Génevaux et al (2013)

1184 method and GIS processing procedures. (B) Production of different land use and

1185 granulometry scenarios in CAESAR-Lisflood LEM
1186 **Figure 3.** Constraint slope maps used to generate all of the virtual catchments

1187 presented in this paper. (A) Terrain slope map (B) River slope map.
1188 **Figure 4.** Cumulative area distribution associated to the three virtual catchments

1189 shown in figure 5.(A) Weakly, (B) moderately and (C) strongly divided hydrographic

1190 network. (D) Cumulative area distribution of the Middle Creak catchment in Australia

1191 (after Perera and Willgoose, 1998). The vertical lines in each graph indicate regions

1192 with different flow regimes (see text for further explanation).
1193 **Figure 5.** Maps of three virtual catchments produced with our method by using

1194 drainage densities (Dd) equal to (A) 2.19km/km², (B) 5.51km/km² and (C) 8.35km/km².

1195 Cumulative area distributions associated to these catchments are respectively shown

1196 in figure 4.
1197 **Figure 6.** DEMs derived from figure 5A's DEM by adjusting integral hypsometric to (A)

1198 0.25, (B) 0.65 and (C) 0.75 values. The hypsometric curves of the three virtual

1199 catchments are displayed in graph D.

1200
1201 **Figure 7.** DEMs of three virtual catchments derived from figure 6B's DEM by adjusting

1202 mean main channel slopes to (A) 0.005, (B) 0.1 and (C) 0.2 values. Note that the

1203 altitude values are represented with the same range of values for all three maps.
1204 **Figure 8.** Geomorphological map obtained for a virtual catchment example by

1205 reclassifying the height geomorphological contexts generated by the Jasiewicz and

1206 Stepinski method (2013) into three contexts: plateau, hillslope and valley.
1207 **Figure 9.** Examples of soil granulometry scenarios that can be implemented in LEM

1208 CAESAR-Lisflood. Setup 1, 2 and 3 present different proportions of fine and coarse

1209 particles.
1210 **Figure 10.** DEM noise reduction by applying the ANUDEM interpolation and smoothing

1211 method. (A) Initial DEM (B) DEM produced by ANUDEM re-interpolation/smoothing.
1212 **Figure 11.** Topography comparison of (A) commonly produced virtual catchment (Van

1213 De Wiel & Coulthard, 2010) (B) virtual catchment produced by our method (C) and a

1214 real catchment (Austreberthe, Normandy, France)