eSCAPE: Regional to Global Scale Landscape Evolution Model v2.0

eSCAPE is a Python-based landscape evolution model that simulates over geological time (1) the dynamic of the landscape, (2) the transport of sediment from source to sink, and (3) continental and marine sedimentary basins formation under different climatic and tectonic conditions. eSCAPE is open-source, cross-platform, distributed under the GPLv3 license and available on GitHub (escape-model.github.io). Simulated processes rely on a simplified mathematical representation of landscape processes the stream power and creep laws to compute Earth’s surface evolution by rivers and hillslope transport. 5 The main difference with previous models is in the underlying numerical formulation of the mathematical equations. The approach is based on a series of implicit iterative algorithms defined in matrix form to calculate both drainage area from multiple flow directions and erosion/deposition processes. eSCAPE relies on PETSc parallel library to solve these matrix systems. Along with the description of the algorithms, examples are provided and illustrate the model current capabilities and limitations. eSCAPE is the first landscape evolution model able to simulate processes at global scale and is primarily designed 10 to address problems on large unstructured grids (several millions of nodes).


Introduction
Since the 80s, many software have been designed to estimate long-term catchment dynamic, drainage evolution as well as sedimentary basins formation in response to various mechanisms such as tectonic or climatic forcing (Braun and Sambridge, 15 1997;Coulthard et al., 2002;Davy and Lague, 2009;Simoes et al., 2010;Salles, 2016;Grieve et al., 2016b;Hobley et al., 2017). These models rely on a set of mathematical and physical expressions that simulate sediment erosion, transport and deposition and can reproduce the first order complexity of Earth's surface geomorphological evolution (Tucker and Hancock, 2010;Shobe et al., 2017).
In most of these models, climatic and tectonic conditions are imposed and often consist in rather simple forcing such as 20 uniform spatial precipitation and vertical displacements (uplift or subsidence) far from reflecting the complexity of the natural system. In addition such approaches are unable to properly explore potential feedback mechanisms between each of the Earth components. In fact, only a handful of these models are able to account more completely for the dynamics of the lithosphere and mantle, the role of sedimentation and provide a more quantitative representation of climate relative to its interactions with topography (such as orographic rain) (Beaumont et al., 1992;Salles et al., 2011;Thieulot et al., 2014;Yang et al., 2015;Salles et al., 2017;Beucher et al., 2019). When made possible, it is often realised through the coupling of specialised numerical models involving the expertise of geodynamicists, geophysicists, Earth surface and atmospheric scientists.
Several geodynamic numerical models constrained by both geological and geophysical observations have been developed and some of these global models (Zhong et al., 2000;Moresi et al., 2003;Heister et al., 2017) have shown how mantle convection drives the motion of tectonic plates and dictates the long-term evolution of the Earth. Similarly progresses in the understanding 30 of past, present, and future climates have been made by the development of mathematical models of the general circulation of a planetary atmosphere or ocean that simulate climate at an increasing level of detail (Dutkiewicz et al., 2016;Brown et al., 2018).
Yet, we are still missing a tool to evaluate global scale evolution of Earth surface and its interaction with the atmosphere, the hydrosphere, the tectonic and mantle dynamics. Such a tool will certainly provide new insights and help to better characterise 35 many aspects of the Earth system ranging from the role of atmospheric circulation on physical denudation, from the influence of erosion and deposition of sediments on mantle convection, from the location and abundance of natural resources to the evolution of life.
The model presented in this paper is a first step toward the development of a parallel global scale landscape evolution model. It allows to couple the Earth's surface with global climatic perturbations and geodynamic forces acting within the Earth's 40 interior. Landscapes and sedimentary basins evolution in eSCAPE are driven by a series of standard stream power incision and diffusion laws (Howard et al., 1994;Tucker and Slingerland, 1997;Chen et al., 2014) designed to address problems from regional to global scales and over geological time (10 5 -10 9 years). Due to the inherent assumptions made in the set of equations used, eSCAPE is not intended to estimate the evolution of individual fluvial channels but to quantify large scale and long term evolution of Earth's surface regions (Salles et al., 2017;Armitage, 2019). 45 First, this paper presents the implicit, iterative approaches that are used to solve the multiple flow direction water routing and the erosion deposition processes (section 2). Then in section 3, I provide a list of all the parameters required to run the eSCAPE model and I discuss the input and output formats. In addition, three examples based on both generic and global scale experiments are described in detail and showcase the code main capabilities. Finally in section 4, I analyse the scalability of eSCAPE and discuss some of the limitations and future implementations that are necessary to improve the performance of the 50 code on parallel architectures.
2 Modelled processes and algorithms eSCAPE (Salles, 2018) is a parallel landscape evolution model, built to simulate landscapes and basins dynamic at various space and time scales over unstructured grids. The model accounts for river incision using stream power law, hillslope processes and sediment transport in land and marine environments. It can be forced with spatially and temporally varying tectonics 55 (horizontal and vertical displacements) and climatic forces (temporal and spatial precipitation changes and sea-level fluctuations). eSCAPE is primarily written in Python with some functions in Fortran and takes advantage of PETSc solvers (Balay and most approaches (Wallis et al., 2009;Wallace et al., 2010;Tarboton, 2013;Bellugi et al., 2011;Braun and Willett, 2013) are based on an initial ordering process followed by efficient priority-queue implementations with some variants such as the sub-basin acyclic graph partitioning method in Salles and Hardiman (2016) or the breadth-first traversal approaches proposed by Barnes (2019). Except for the approach proposed by Barnes (2019), the previous methods scale well as long as the number of processors used is modest but quickly deteriorates as inter-processors communication cost increases.

70
In addition, when using the aforementioned implementation strategies, several problems might arise in (1) load balancing, when catchments size greatly changes in the simulated domain or (2) handling very high resolutions where multiple processes are needed for a single catchment. In addition, most of these methods assume a single flow direction (SFD - Fig. 1a). This assumption makes the emergent flow network highly sensitive to the underlying mesh geometry and most dendritic shape of obtained stream networks is often an artefact of the surface triangulation. To reduce this effect, authors have proposed to 75 consider not only the steepest downhill direction but also to represent other directions appropriately weighted by slope (multiple flow direction -MFD). Using MFD algorithms prevent locking of erosion pathways along a single direction and help to route flow over flat regions into multiple branches (Tucker and Hancock, 2010). Yet, graph traversal approaches cannot be easily modified to incorporate MFD algorithms as catchments are no longer strictly isolated in low slope areas and flow pathways often diverge (Fig. 1b).

80
To overcome these limitations, Richardson et al. (2014) proposed to use linear solvers. The approach consists in writing the FA calculation as a sparse matrix system of linear equations (Eddins, 2007;Schwanghart and Kuhn, 2010). It can take full advantage of purpose-built, efficient linear algebra routines including those provided by parallel libraries such as PETSc (Balay et al., 2012). eSCAPE computes the flow discharge (m 3 /y) from FA and the net precipitation rate P using the parallel implicit drainage area (IDA) method proposed by Richardson et al. (2014) but adapted to unstructured grids (Fig. 1).

85
The flow discharge at node i (q i ) is determined as follows: where b i is the local volume of water Ω i P i where Ω i is the voronoi area and P i the local precipitation value available for runoff during a given time step. N d is the number of donors with a donor defined as a node that drains into i (as an example the donor of vertex 5 in the SFD sketch in Fig. 1a is 1). To find the donors of each node, the method consists in finding their receivers first. Then, the receivers of each donor is saved into a receiver matrix, noting that the nodes, which are local minima, are their own receivers. Finally the transpose of the matrix is used to get the donor matrix. When Eq. 1 is applied to all nodes and considering the MFD case presented in Fig. 1a, the following relations are obtained: q 4 = b 4 + q 1 w 1,4 + q 2 w 2,4 q 5 = b 5 + q 1 w 1,5 + q 4 w 4,5 q 6 = b 6 + q 4 w 4,6 + q 5 w 5,6 + q 7 w 7,6 q 7 = b 7 + q 10 w 10,7 q 8 = b 8 + q 3 w 3,8 + q 4 w 4,8 + q 6 w 6,8 + q 7 w 7,8 + q 10 w 10,8 q 9 = b 9 + q 3 w 3,9 + q 8 w 8,9 + q 10 w 10,9 (2) and sum to one for each node: n w m,n = 1 In eSCAPE, the number of flow direction paths is user-defined and can vary from 1 (i.e. SFD) up to 12 (i.e. MFD) depending of the grid neighbourhood complexity. The weights are calculated based on the number of downslope neighbours and are proportional to the slope (Quinn et al., 1991;Tarboton, 1997;Richardson et al., 2014).

100
In matrix form the system defined in Eq. 2 is equivalent to Wq=b or: the vector q corresponds to the unknown flow discharge (volume of water flowing on a given node per year) and the elements of W left blank are zeros.
As explained in Richardson et al. (2014), the above system is implicit as the flow discharge for a given vertex depends on its neighbours unknown flow discharge. The matrix W is sparse and is composed of diagonal terms set to unity (identity matrix) and off-diagonal terms corresponding to at most the immediate neighbours of each vertex (typically lower than 6 in constrained Delaunay triangulation).
In eSCAPE, this matrix is built in parallel using compressed sparse row matrix functionality available from SciPy (Jones et al., 2001). Once the matrix has been constructed, PETSc library is used to solve matrices and vectors across the decom-110 posed domain (Balay et al., 2012). The performance of the IDA algorithm is strongly dependent on the choice of solver and preconditioner. In eSCAPE, the solution for q is obtained using the Richardson solver (Richardson, 1910) with block Jacobi preconditioning (bjacobi). This choice was made based on the convergence results from Richardson et al. (2014) but can be changed if better solver and preconditioner combinations are found. Iterative methods allow for an initial guess to be provided.
When this initial guess is close to the solution, the number of iterations required for convergence dramatically decreases. I take 115 advantage of this option in eSCAPE by using the flow discharge solution from the previous time step as an initial guess. This allows to decrease the number of iterations of the IDA solver as discharge often exhibits small change between successive time intervals.

Erosion and sediment transport
River incision, associated sediment transport and subsequent deposition are critical elements of landscape evolution models.

120
Commonly these are defined based on either a transport-limited (Willgoose et al., 1991) or a detachment-limited (Howard et al., 1994) approach. On one hand, the transport-limited hypothesis assumes that rivers may be able to transport sediment up to a concentration threshold (often referred to as the stream transport capacity) linked to discharge, slope, sediment size, and channel form (channel depth/width ratio) and that an infinite supply of sediment is available for transport. On the other hand, the detachment-limited hypothesis supposes that erosion is not limited by a transport capacity but instead by the ability 125 of rivers to remove material from the bed. Even though validations of each hypothesis have been conducted based on field studies calibration (Snyder et al., 2003;Tomkin et al., 2003;van der Beek and Bishop, 2003;Valla et al., 2010;Hobley et al., 2011) there are many evidences suggesting that both transport and detachment limited behaviours take place simultaneously in natural systems and models accounting for transition between the two have been proposed in the past (Beaumont et al., 1992;Braun and Sambridge, 1997;Coulthard et al., 2002;Davy and Lague, 2009;Hodge and Hoey, 2012;Salles and Duclaux, 2015;130 Carretier et al., 2016;Turowski and Hodge, 2017;Lague, 2010;Shobe et al., 2017;Hobley et al., 2017;Salles et al., 2018). For simplicity, the approach proposed in this paper is similar to the initial version of eSCAPE (v1.0.0 -Salles (2018)) and is based on a standard form of the stream power law assuming detachment-limited only behaviour. In the future, a better representation of erosion and sediment transport could be added such as the SPACE approach proposed by Shobe et al. (2017).
As mentioned above and following Howard et al. (1994), I consider that sediment erosion rate is expressed using a stream 135 power formulation function of river discharge and slope. The volumetric entrainment flux of sediment per unit bed area E is of the following form: where K is the sediment erodibility parameter, Q is the water discharge, S is the river slope. In eSCAPE, I incorporate local precipitation-dependent effects on erodibility (Murphy et al., 2016) and use the flow discharge defined in previous section 140 Q = P A to represent rainfall gradients effect on discharge. A is the flow accumulation and P the upstream annual precipitation rate. m and n are scaling exponents. In our model, K is user defined and the coefficients m and n are set to 0.5 and 1 respectively (Tucker and Hancock, 2010). E is in m/y and therefore the erodibility dimension is (m·y) −0.5 .
The entrainment rate of sediment (E) is approached by an implicit time integration and consists in formulating the stream power component in Eq. 5 in the following way: where λ i,rcv is the length of the edges connecting the considered vertex to its receiver. Rearranging the above equation gives: ( with the coefficient K f,i|rcv = K √ Q i ∆t/λ i,rcv . In matrix form the system defined in Eq. 6 is equivalent to: Γη t+∆t = η t .
Using the case presented in Fig. 1a, the matrix system based on the receivers distribution is defined as: This system is implicit and the matrix is sparse. The SciPy compressed sparse row matrix functionality (Jones et al., 2001) is used to build Γ on local domains and the Eq. 7 is then solved using Richardson solver with block Jacobi preconditioning 155 (bjacobi) using an initial guess for the solution set to vertices elevation.
Once the entrainment rates have been obtained, the sediment flux moving out at every node Q out s equals the flux of sediment flowing in plus the local erosion rate. Q out s takes the following form: Ω is the voronoi area of the considered vertex and F f is the fraction of fine sediment that remains in suspension. The solution 160 of the above equation requires the calculation of the incoming sediment volume from upstream nodes Q in s . At node i, Eq. 8 is equivalent to: where e i = (1 − F f )E i Ω i and N d the number of donors. Assuming that river sediment concentration is distributed in a similar way as the water discharge we can write a similar set of equalities as the ones in Eq. 2. Then a matrix system as proposed for 165 the FA (Eq. 4) can be obtained. The new system is then solved using the PETSc solver and preconditioner previously defined.

Priority-flood depression filling
In most landscape evolution models, internally-draining regions (e.g., depressions and pits) are usually filled before the calcu- Priority-flood algorithms consist in finding the minimum elevation a cell needs to be raised to (e.g., spill elevation of a cell) to prevent downstream ascending path to occur. They rely on priority queue data structure used to efficiently find the lowest spill 180 elevation in a grid. Depending on the chosen method, priority queue implementation approaches affect the time complexity of the algorithm (Barnes et al., 2014). In eSCAPE, the priority-flood + variant of the algorithm proposed in Barnes et al. (2014) is implemented. It provides a solution to remove automatically flat surfaces and it produces surfaces for which each cell has a defined gradient from which flow directions can be determined.
In eSCAPE, this part of the algorithm is not parallelised and is performed on the master processor. It starts from the grid 185 border vertices and processes vertices that are in their immediate neighbourhoods one by one in the ascending order of their spill elevations (Barnes et al., 2014). The initialisation step consists in pushing all the edges nodes onto a priority queue. The priority queue rearranges these nodes so that the ones with the lowest elevations in the queue are always processed first.
To track nodes that have already been processed by the algorithm a Boolean array is used in which edge nodes (that are by definition at the correct elevation) are marked as solved. The next step consists in removing (i.e. popping) from the priority 190 queue the first element (i.e. the lowest node). This node n is guarantee to have a non-ascending drainage path to the border of the domain. All non-processed neighbours (based on the Boolean array) from the popped node are then added to the priority queue. In the case where a neighbour k is at a lower elevation than n its elevation is raised to the elevation of n plus before  Figure 2. Illustration of the two cases that may arise depending on the volume of sediment entering an internally drained depression (panel a). The red line shows the limit of the depression at the minimal spillover elevation. b) The volume of sediment (V in s ) is lower than the depression volume Vpit. In this case all sediments are deposited and no additional calculation is required. c) If V in s ≥ Vpit, the depression is filled up to depression filling elevation (priority-flood + ), the flow calculation needs to be recalculated and the excess sediment flux (Q ex s ) is transported to downstream nodes. being pushed to the queue. Once k has been added to the queue, it is marked as resolved in the Boolean array. In this basic implementation of the priority-flood algorithm, the process continues until the priority queue is empty (Barnes et al., 2014). 195

Depression filling and marine sedimentation
The filling algorithm presented above is used to calculate the volume of each depression at any time step. Once the volumes of these depressions are obtained, their subsequent filling is dependent of the sediment fluxes calculation defined in section 2.2 (Fig. 2a). In cases where the incoming sediment volume is lower than the depression volume (Fig. 2b), all sediments are deposited and the elevation at node i in the depression is increased by a thickness δ i such that: where η f i is the filling elevation of node i obtained with the priority-flood + algorithm and the ratio Υ is set to V in s /V pit . If the cumulative sediment volume transported by the rivers draining in a specific depression is above the volume of the depression (V in s ≥ V pit - Fig. 2c) the elevation of each node i is increased to its filling elevation (η f i ) and the excess sediment volume is allocated to the spillover node (Fig. 2c). The updated elevation field is then used to compute the flow accumulation following the approach presented in section 2.1. The sediment fluxes are initially set to zero except on the spillover nodes and using Eq. 9 the excess sediments are transported downstream until all sediments have been deposited in depressions, have entered the marine environment, or have moved out of the simulation domain.
In the marine realm, sedimentation computation follows a different approach to the one described above. First, the flow accumulation is computed using the filled elevation in both the aerial and marine domains and a maximum volumetric deposition 210 rate ζ i is calculated based on the depth of each marine node: with η sl the sea-level position. Using similar solver and preconditioner as the ones proposed for the flow discharge calculation, we solve implicitly a matrix system equivalent to the one in Eq. 4 with the same weight (W ) and a vector b equals to q s,i − ζ i .
From the solution, only positive sedimentation rates are initially kept and the sedimentation thicknesses for these nodes are set 215 to ζ∆t. Then remaining sediment fluxes on adjacent vertices are found by computing the sum of ζ and obtained sedimentation rates and by considering again only positive values.

Hillslope processes and marine top sediment layer diffusion
Hillslope processes are known to strongly influence catchment morphology and drainage density and several formulations of hillslope transport laws have been proposed (Culling, 1963;Tucker and Bras, 1998;Perron and Hamon, 2012;Howard et al., 220 1994;Fernandes and Dietrich, 1997;Roering et al., 1999Roering et al., , 2001. Most of these formulations are based on a mass conservation equation and with some exceptions such as CLICHE model (Bovy et al., 2016), these models assume that a layer of soil available for transport is always present (i.e. precluding case of bare exposed bedrock) and that dissolution and mass transport in solution can be neglected (Perron and Hamon, 2012).
Under such assumptions and via the Exner's law, the mass conservation equation widely applied in landscape modelling is of 225 the form (Dietrich et al., 2003;Tucker and Hancock, 2010): where q ds is the volumetric soil flux of transportable sediment per unit width of the land surface. In its simplest form, q ds obeys the Culling model (Culling, 1963) and hypothesises a proportional relationship to local hillslope gradient (i.e. q ds = −D∇η, also referred to as the creep diffusion equation): in which D is the diffusion coefficient that encapsulates a variety of processes operating on the superficial soil layer. As an example, D may vary as a function of substrate, lithology, soil depth, climate and biological activity (Tucker et al., 2001;Tucker and Hancock, 2010). The creep law is found in many models such as GOLEM (Tucker and Slingerland, 2017), CHILD (Tucker and Slingerland, 1997), LANDLAB (Hobley et al., 2017) or Badlands (Salles and Hardiman, 2016;Salles et al., 2018), and 235 in Willgoose et al. (1991), Fernandes and Dietrich (1997), Tucker and Slingerland (1997), Simpson and Schlunegger (2003).
In eSCAPE, hillslope processes rely on this approximation even though field evidences suggest that the creep approximation (Eq. 12) is only rarely appropriate (Roering et al., 1999;Tucker and Bradley, 2010;Foufoula-Georgiou et al., 2010;DiBiase et al., 2010;Larsen and Montgomery, 2012;Grieve et al., 2016a). In the future, a possible improvement could be based on the nonlinear hillslope transport equation incorporating a critical slope to model hillslope soil flux (Roering et al., 1999(Roering et al., , 2001.

240
For a discrete element, considering a node i the implicit finite volume representation of Eq. 12 is: N is the number of neighbours surrounding node i, Ω i is the voronoi area, λ i,j is the length of the edge connecting the considered nodes and χ i,j is the length of voronoi face shared by nodes i and j. Applied to the entire domain, the equation above can be rewritten as a matrix system Qη t+∆t = η t where Q is sparse. The matrix terms only depend on the diffusion 245 coefficient D, the grid parameters and voronoi variables (χ i,j , λ i,j , Ω i ). In eSCAPE, these parameters remain fixed during a model run and therefore Q needs to be created once at initialisation. At each iteration, hillslope induced changes in elevation η are then obtained in a similar way as for the solution of the other systems using PETSc Richardson solver and block Jacobi preconditioning.
In addition to hillslope processes, a second type of diffusion is available in eSCAPE and consists in distributing freshly de-250 posited marine sediments in deeper regions. This process is the only one treated explicitly and in this case the length of the diffusion time step ∆t m must be less than a CFL factor to ensure numerical stability: where D m is the diffusion coefficient for the newly deposited marine sediments. Even with a reasonable small time step, the Eq. 14 can produce incorrect results. Following Bovy et al. (2016), the following set of inequalities are also added: where q out ms,ij is the flux of sediment from the marine top layer leaving node i towards the downstream neighbours j, h i is the depth of the marine top layer, η m is the elevation associated to the highest downslope neighbour of i and α is a factor lower than 1. These inequalities are always satisfied if positive outgoing fluxes are scaled by a factor β given by:

260
In eSCAPE, marine diffusion of freshly deposited sediment is performed explicitly using the CFL condition described in Eq.
14 and the restriction proposed in Eq. 16.

Usability and applications
In this section, I present the main files used to run eSCAPE and to visualise the generated outputs. I then illustrate the capability of the code using a series of 3 examples presenting two generic models and one global scale experiment.

Input parameters and visualisation
eSCAPE uses YAML syntax for its input file. YAML structure is done through indentation (one or more spaces) and sequence items are denoted by a dash. When parsing the input file, the code is searching for some specific keywords defined in tables 1, 2 and 3. Some parameters are optional and only need to be set when specific forces (table 2) or physical processes (table 3) are applied to a particular experiment.

270
All the input parameters that are defined in external files like the initial surface, different precipitation or displacement maps are read from VTK files. These input files are defined on an irregular triangular grid (TIN). Examples on how to produce these files are provided in the eSCAPE demo repository on Github and Docker. The only exception is the sea-level file which is a two-columns CSV file containing on the first column the time in years and ordered in ascending order and on the second one the relative position of the sea-level in metres (curve in table 2).

275
The domain and time keywords (table 1)    This surface is exposed to an uniform precipitation regime of 1 m/y and is uplifted linearly from its fixed western side to the eastern one that experiences an uplift of 5 mm/y (Fig. 3a). The proposed setting is similar to the one in Braun and Willett (2013) and the value of the bedrock erodibility parameter K is set to 2 × 10 −4 in order to reach steady-state during the simulated 10 5 300 years. Under such conditions, the model is purely erosional and therefore neither the aerial and marine sedimentation nor the depression filling algorithm are considered. In addition hillslope processes are also turned off, meaning that this example only relies on the implicit parallel flow discharge and erosion equations defined in sections 2.1 and 2.2.
Three cases are presented after 10 5 years for different time steps ∆t varying from 10 4 to 10 3 and 10 2 years in Figure 3a, 3b and 3c respectively, implying that the number of steps is 10, 100 and 1000. In both cases the implicit schemas converge 305 for the chosen solver and preconditioner (i.e. Richardson with block Jacobi). The solutions for the mean landscape elevation ( Fig. 3d) show that the landscape reaches steady state in all cases and overall the final elevations are in good agreement with a maximum elevation of 482 ± 3 m and a number of catchments n c almost identical between models (87 ≤ n c ≤ 94). Yet as the time step increases the differences between models increase over time. By the end of the simulation, the mean elevation difference between the case with ∆t equals to 10 2 y and the one at 10 3 y is around 2.5% whereas the difference with a ∆t 310 of 10 4 y is above 30% (Fig. 3d).  by specifying the number of directions ( Fig. 1a and flowdir parameter in table 1) appropriately weighted by slope that rivers could potentially take when moving downhill.
lowest regions (at 0 m elevation) located on the edges of the domain and increasing to 1000 metres towards the center. The triangulated circular grid of 50 km radius is built with a resolution of approximately 200 m. The three experiments with varying water routing directions are ran for 100,000 years with a ∆t of 1000 years under a 1 m/y uniform precipitation. In addition to stream incision (bedrock erodibility K set to 2 × 10 −5 ), hillslope processes are also accounted for using a diffusion coefficient D of 10 −2 m 2 /y.

335
After 20,000 years, the dendritic flow accumulation pattern observed on the surface for the SFD case (Fig. 4b) is analogue to many natural forms of drainage systems but is actually a numerical artefact and depends on the random locations of the nodes in the surface triangulation. By increasing the number of possible downstream directions, this sensitivity to the mesh discretisation is significantly reduced (as illustrated in Fig. 4b where a second direction is added). In addition, routing flow to more than one destination node allows a better representation of channels pathways divergence into multiple branches over flat 340 regions (Tucker and Hancock, 2010).
Landscape evolution models tend to be highly dependent on grid resolution and this dependency is mostly related to the approach used to route water down the surface (Schoorl et al., 2000;Pelletier, 2004;Armitage, 2019). As discussed by Armitage (2019), enabling node-to-node MFD algorithm decreases the dependence of landscape features (e.g. valley spacing, branching of stream network, sediment flux etc) to grid resolution. As shown in Figure 4b and c, SFD algorithm leads to increase 345 branching of valleys whereas the MFD approach, by promoting wider flow distribution, produces smoother topography where local carving of the landscape is reduced. Armitage (2019) also showed that when using models that operate at scale larger than river width resolution, node-to-node MFD algorithm creates landscape features that are not resolution dependent and that evolve closer to the ones observed in nature. Therefore it is recommended to use more than one downhill direction (flowdir) in eSCAPE when looking at global and continental scale landscape evolution or for cases where multiple resolutions are 350 considered within a given mesh.

Global scale simulation
The last example showcases a global scale experiment with eSCAPE. The simulation looks at the evolution of the Earth 200,000 years into the future starting with present day elevation and precipitation maps. The model is ran forward in time without changing the initial forcing conditions and is primarily used to highlight the capabilities of eSCAPE and does not 355 represent any particular geological situation.
The elevation is obtained from the ETOPO1 1 arc-minute global relief model of Earth's surface that integrates land topography and ocean bathymetry (Amante and Eakins, 2009). For the rainfall, I summed all the WorldClim gridded climate monthly dataset to obtained a global yearly rainfall map with a spatial resolution of about 1 km 2 (Fick and Hijmans, 2017). From these dataset, I then built the initial surface and climate meshes at 16 km resolution consisting in approximately 3 millions 360 points ( Fig. 5a and 6a). The model inputs are temporally uniform, but any other climatic scenarios could have been chosen for illustration as well as tectonic conditions (both vertical: uplift and subsidence and horizontal: advection displacements). In addition to these grids, the following parameters are chosen: flowdir is set to 5 (similar to the MFD flow routing approach), the time step ∆t equals 500 years, bedrock erodibility K is 5 × 10 −5 , the diffusion coefficient D equals 10 −1 , the fraction of sediment in suspension Ff is 0.3 (F f parameter in section 2.2 and defined in table 3) and the marine fresh sediment diffusion 365 D m is set to 5 × 10 5 (see section 2.5 and sedimentK parameter in table 3). The simulation took 2 hours to run on a cluster using 32 processors.
From these set of input parameters, eSCAPE can predict the global evolution of topography (Figures 5 and 6) and quantifies the associated volume and spatial distribution of sediments trapped in continental plains or transported into the marine realm.
By recording eroded sediment transport over each drainage basin, eSCAPE provides estimation of sedimentary mass fluxes 370 carried by major rivers into the ocean (section 2.2). As shown in Figures 5 and 6, the predicted locations of largest basin outlets match quite well with observations and many of the biggest simulated deltaic systems are related to sediment transported by some of the world's largest rivers (Milliman and Syvitski, 1992;Syvitski et al., 2003). The model can also be used to evaluate the evolution of drainage systems, the stability of continental flow directions, the exhumation history of major mountain ranges, the timing and geometry of sedimentary bodies formation (e.g. deltas or intra-continental deposits) as well as basins  (Saad, 2003) and I decided to use the block Jacobi preconditioner as it is one of the simplest methods and produces in combination with the Richardson solver good scalability (Richardson et al., 2014). Yet other combinations such as the Richardson solver with the Euclid preconditioner (i.e. HYPRE package, Falgout et al. (2012)) might exhibit better scalability in some cases.

395
The analysis of the profiling work realised in Fig. 7a suggests that for purely erosive models (similar to the ones presented in section 3.2.1) most of the computational time is spent solving the Richardson iterative method (PETSc solver KSP). From the graph on the right hand side of the panel (Fig. 7a), one can deduce that performance improvements are obtained when the problem size increases. However the scaling performance decreases when reaching 32/64 processors depending on the problem size. This does not agree with some of the conclusions from Richardson et al. (2014), where the scaling of the implicit drainage 400 accumulation algorithm continues even for large numbers of processors (>192). To improve performance, I will be exploring two directions. First the problem might be related to the chosen solver and preconditioner combination and I will run new tests using the HYPRE package as discussed above. Secondly, the poor performance for larger processors might also be linked to issues related to either Python PETSc wrapper or installation problems and incompatibilities between some of the compilers and packages that I used. In the future, different software libraries and compilers versions from GNU and Intel will be tested 405 and might help to improve the performance for increasing number of processors.
For experiments accounting for marine deposition and pit filling, a similar trend is found when comparing performance against processors number (right hand side in Fig. 7b). However the results from the profiling (left hand side in Fig. 7b) suggest that more than half of the computation time is now spent on non-PETSc work with the biggest proportion related to the pit filling function. In eSCAPE, the priority-flood algorithm is performed in serial (see section 2.3). This is the major limitation 410 of the code as shown by the time spent in broadcasting the information from the master node to the other processors (MPI Bcast in Fig. 7b). To take advantage of parallel architectures, several authors (Wallis et al., 2009;Tesfa et al., 2011;Yıldırım et al., 2015;Zhou et al., 2017) have proposed partitioning implementation of depression filling algorithms. However, most of these methods require frequent interprocess communication and synchronisation. This becomes even more problematic in the case of eSCAPE where the depression filling algorithm needs to be performed at every time step (Barnes, 2019). Barnes

415
(2017) presented an alternative to the aforementioned parallelisation methods that limits the number of communications. Yet this approach is not fully satisfactory as it only fills the depressions up to the spilling elevation but does not provide a way of implementing efficiently the variant of the algorithm proposed in (Barnes et al., 2014). Finding a strategy to perform a parallel version of the variant of the priority-flood algorithm that could provide flow directions on flat surface will likely improve greatly the performance of eSCAPE.

Conclusions
In this paper, I describe eSCAPE, an open-source, Python-based software designed to simulate sediment transport, landscape dynamics and sedimentary basins evolution under the influence of climate, sea level and tectonics. In its current form, eSCAPE relies on the stream power and creep laws to simulate the physical processes acting on the Earth's surface. The main difference with other landscape evolution models relies on the formulation used to solve the system of equations. The approach 425 builds upon the Implicit Drainage Area calculation from Richardson et al. (2014) and consists in a series of implicit iterative algorithms for calculating multiple flow direction and erosion deposition that written in matrix form. As a result, the obtained systems can be solved with widely available parallel linear solver packages such as PETSc.
Performance analysis shows good parallel scaling for small number of processors (under 64 processors as shown in section 4) but some work is required for larger numbers. The code profiling suggests that the main issue is in the inter-processes commu- the future, a parallel approach allowing depression filling and flow direction computation over flat regions will be critical to improve the overall performance of the code.
Examples are provided in the paper and available through the Docker container. They illustrate the extent of temporal and spatial scales that can be addressed using eSCAPE. As such this code is highly versatile and useful for geological applications 435 related to source to sink problems at regional, continental and for the first time global scale. It is already possible to use eSCAPE to simulate global geological evolution of the Earth's landscape at about 1 km resolution providing accurate estimates of quantities such as large-scale erosion rates, sediment yields and sedimentary basins formation. In the future, the code will be coupled with atmospheric and geodynamic models to bridge the gap between local and global scales predictions of Earth past and future evolutions.