SPGM: A Scalable PaleoGeomorphology Model

Numerical models of landscape evolution are playing an increasingly important role in providing an improved understanding of geomorphic transport processes shaping Earth’s surface topography. Improving theoretical underpinnings coupled with increasing computational capacity has led to the development of several open source codes written in low-level languages. However, adapting these codes to new functionality or introducing greater flexibility often requires significant recoding. Here we present a multi-process, scalable, numerical model of geomorphological evolution, built with a modular structure and geared toward seamless extensibility. We implement recent algorithmic advances that reduce the computational cost of flow routing – a problem that typically scales quadratically with the number of unknowns – to linear in time while allowing for parallel implementations of geomorphic transport processes. Our scalability tests demonstrate that such parallelizations can achieve an order of magnitude speedup on a typical desktop computer, making large-scale simulations more tractable.


Introduction
The interplay of a variety of physical processes acting over a range of space and time scales is responsible for the evolution of Earth's surface topography. These physical processes broadly fall in four categories: (i) river incision that leads to advective transport of material over distances of up to several thousands of km, (ii) diffusive processes such as soil creep and landslides that operate over comparatively shorter length scales, (iii) large-scale surface deformation arising from tectonic forces in the crust and lithosphere, and (iv) transient long-wavelength but small-amplitude show significant speedups. Braun and Willett [9] presented a parallelization approach based on open-multi-processing (OpenMP), producing near-linear scalability on a shared-memory machine. However, both MPI-and GPU-based parallelization strategies involve increased development and debugging times and are likely to compromise ease of code adaptability. We therefore adopt the parallelization approach described in Braun and Willett [9] for its simplicity and ease of implementation.

Software framework
SPGM has been implemented using C++ [14] and the basic abstractions utilized are shown in Fig. 1. The Config class is responsible for parsing input configuration files and for providing collaborating classes access to named parameters and parametergroups. Each of the physical processes implemented in SPGMwhich we interchangeably refer to as 'modules', e.g. Precipitation -inherit from the abstract Process class. The Process class contains a reference to the Config class and its sub-classes (blue rectangle in Fig. 1) implement the Execute() function, which contains the core numerical algorithm pertaining to a given physical process. The ModelBuilder class instantiates a number of 'Processes' based on an input configuration file and populates an instance of the Model class with them. The SurfaceTopology class manages the underlying geometry of the computational mesh and recomputes drainage networks at the beginning of each time step, as shown in the program flowchart in Fig. 2.
In the following sections we describe the parsing of configuration files where general parameter values, as well as timeseries data for specific parameters can be specified. We also briefly describe the algorithms used in mesh generation, flow-routing and output generation.

Configuration files and parameter time series
We use simple, plain-text configuration files to list a number of required parameters, followed by several groups of parameters corresponding to physical processes modeled. Parameter values such as rate of precipitation or uplift can be specified in the configuration file as either: (i) a scalar value, (ii) a time-series in a text file containing a sequence of (time, value) pairs, or (iii) a timeseries in a text file containing a sequence of (time, file) pairs-where each file lists nodal values at that time, thus allowing space-time variability (see Supplementary Data). Linear interpolation is used in the second and third options for model times that fall in between times at which parameter values are provided.

The computational mesh
The initial topography (xyz triplets) and boundary conditions (BC) can be specified in a simple four-columned text file. Only arbitrary Dirichlet boundary conditions are currently supported and nodes with a fixed elevation are marked by '1' in the fourth (BC) column-free nodes are marked by '0'. Output from a given model can also be used as the initial topography of another model (see Appendix C.2 for more details).
We compute a Delaunay triangulation and the corresponding Voronoi tessellation of the initial point-data on the x-y plane using the divide and conquer method [15], following the algorithm outlined in Lischinski [16]. Numerical solutions obtained in the various modules, described in latter sections, directly update the height-field on this computational mesh. Each point in the triangulation is connected to a set of neighbors -also known as 'natural neighbors' -and its associated Voronoi polygon marks the drainage area associated with the point. For highly irregular distributions of initial point-data, iterative Laplacian smoothing can be applied as: where, at each iteration, φ is the smoothing factor, N is the number of neighboring vertices to node i,x j is the position of the jth neighboring node andx i is the new position for node i.

Flow routing
Computing upstream drainage area is computationally expensive -typically scales quadratically with the number of unknowns [17] -and efficient algorithms are necessary as upstream drainage area needs to be recomputed at every time-step in response to changes in elevation through time that affect flow directions. We implement a variant of the recursive algorithm in Freeman [17] that has been adapted in Braun and Willett [9] to compute the order, S, in which nodes need to be traversed such that the contributing drainage area for downstream nodes can be incrementally computed-thereby reducing the original computational cost to linear in time.
Computation of S starts from outlet nodes, specified as BCs in the initial topography file, whose elevations are fixed through time. As an important byproduct of this procedure, each catchment is labeled with a unique catchment tag. Node indices in S are grouped into catchments and node indices for each catchment are listed upstream from their corresponding outlet nodes such that each node of confluence appears before all other nodes that contribute flow to it. Consequently, nodes are traversed according to S in reverse order in order to compute discharge at each nodesee Braun and Willett [9] for a detailed account of the procedure.

Output generation
Model outputs are written out in standard VTK format to facilitate the usage of advanced visualization capabilities provided by popular visualization packages such as Paraview [18]. At a given time-step, the mesh, drainage networks and related scalar fields including 'h', 'bc', 'cid' and 'dh' that represent nodal heights, boundary conditions, catchment tags and the difference in nodal heights compared to that in the initial topography, respectively, can be output as VTK files.

Modules
In this section we describe the various modules that each represent individual processes responsible for topographic evolutionsee Fig. 1. We describe modules that account for: (i) precipitation that leads to surface runoff, (ii) uplift, (iii) river incision processes that arise from surface runoff, and (iv) diffusive hill-slope processes. River channels can, however, transition between two different modes, namely, detachment-limited and transport-limited states. In a detachment-limited state the rate of incision is limited by low erodability, whereas inefficient transport of sediments in a transport-limited state limit the evolution of river channels. In other words, a detachment-limited scenario is purely erosive and deposition of sediments only occurs in a transport-limited scenario. However, generally both these modes are necessary to account for the diversity of river channel morphology [19]. We therefore implement separate modules to represent both modes of river incision.

Precipitation
This module simply assigns a parcel of water, w i , to the ith node and is derived as: where a i , λ and ∆t are the area of the Voronoi polygon corresponding to the ith node, rate of precipitation that can vary in space and time and the time-step in years, respectively.

Uplift
This module incorporates uplift (tectonic and dynamic being indistinguishable) based on uplift rate, U, which can be specified as a constant or as a time-series-see Section 2.1. The latter option allows for variable uplift rates in both space and time.

Detachment-limited fluvial transport
In a detachment-limited scenario, fluvial incision can be described by a power law function [20]: where z, t, K f and A are elevation, time, erodability constant and contributing drainage area, respectively. The value of K f depends on a number of factors, including lithology and climate, while its dimension depends on dimensionless constants m and n (that both have small positive values in which their ratio, m/n, is generally ≈0.5 [20]. The solution to this equation is obtained from the implicit scheme described in Braun and Willett [9] that ensures numerical stability even when large time-steps are used-although, larger time-steps may result in lower accuracy (see Fig. 4d in Braun and Willett [9]. The solution scheme traverses the nodes according to S, thus following the donor-link to incrementally progress upstream, starting from outlet nodes (specified as BCs), which have a fixed elevation. When n = 1 in Eq. (3), it can be solved explicitly to compute the evolving elevation of a node. Otherwise, when n ̸ = 1 in Eq. (3), an iterative Newton-Raphson scheme [21] is usedsee Braun and Willett [9] for a detailed derivation of the approach. The solution scheme is parallelized using OpenMP over the total number of catchments identified, since nodes that belong to a given catchment can be processed in isolation from nodes in other catchments-see Section 5 for scalability results.
We test this module with a simple model, where the initial topography of a 100 km × 100 km region is set to a uniform elevation of 500 m. Boundary conditions are set on nodes along y = 0, with their elevation set to 0 and an uplift function is imposed from time t = 0 as: The model is computed with m = 0.5, n = 1, K f = 6 × 10 −4 yr −1 , ∆t = 100 yrs and U 0 = 5 × 10 −3 m yr −1 (see Supplementary Data). These parameters were chosen such that topographic evolution (Fig. 3A) reaches a steady-state at ≈3 × 10 4 years (Fig. 3B) and visually resembles results shown in Fig. 4 of Braun and Willett [9].

Transport-limited fluvial transport
We traverse nodes according to S in reverse order (see Section 2.3) to incrementally compute volumetric discharge, D, through each node. We implement the mathematical model in Kooi and Beaumont [22], where local sediment flux, Q , along river networks can potentially be at a disequilibrium with carrying capacity, Q e , which is defined as: where m and n are dimensionless constants that vary around 1see Tucker and Slingerland [23] for a detailed derivation of Q e .  River networks evolve towards equilibrium at a rate proportional to the disequilibrium present [22]. Sediments are deposited when Q > Q e and changes in elevation are given by: Material becomes entrained when Q < Q e and changes in elevation are given by: where l i and l s are the local channel length and the erosion lengthscale, respectively. Local sediment flux, Q , thus evolves towards an equilibrium with the carrying capacity, eroding and depositing material in the process. We assume l s to be larger for bedrock, compared to that for previously deposited alluvial sediments-see Appendix A for more details. We also record sediment accumulation throughout model evolution so that appropriate values of l s are used for fluvial entrainment at a given location, depending on the presence of previously deposited sediments. We traverse nodes according to S in reverse order and solve Eqs. (6) and (7) explicitly, based on simple continuity of mass. This solution scheme is equally amenable to parallelization as that in Section 3.3 and scalability results are discussed in Section 5.
We test this module with a simple model, where the initial topography in a 1000 km × 1000 km region is derived based on the complementary error function, with a small amount Perlin noise [24] added to induce the development of realistic river networks. Boundary conditions are set on nodes along y = 0, with their elevation set to 0. The model is computed with parameters chosen arbitrarily to demonstrate a simple case of topographic evolution (m = 1, n = 0.5, K f = 5 × 10 −2 , ∆t = 1000 yrs and λ = 0.2 m yr −1 ); these parameters result in a realistic landscape and sediment distribution after 20 × 10 6 years (Fig. 4A), with a reasonable evolution of mean topography over this period (Fig. 4B).
Note that since sediment transport and deposition in a detachment-limited scenario depends on the surface runoff through each node, this module must appear in the parameter file (see Supplementary Data) after the precipitation module.

Short-range hill-slope processes
Short-range material transport due to processes e.g. soil creep, landslides, etc. is approximated by anisotropic diffusion: where K d is a space-time varying diffusivity field. Diffusivities for both bedrock and sediments are specified in parameter files  (see Appendix A for more details). We assume an absence of any sediments at the start of the model-consequently, as material gets entrained and sediments are routed throughout model evolution, the distribution of diffusivities vary in space-time, as described in Section 3.4. We cast Eq. (8) as a simple finite element problem [25], based on P1 elements on the triangulated mesh. We use an implicit Euler time integration scheme and solve the resulting sparse system of equations using the conjugate gradient iterative solver implemented in the Eigen library [26]. A detailed verification procedure for accessing the numerical implementation is given in Appendix B.
To further demonstrate the effects of this module, we impose it on the simple model in Section 3.4-see Appendix A for descriptions of additional parameters pertaining to this module. Model results are shown in Fig. 5, where the spatial distribution of diffusivities (Fig. 5B) reflect sediments deposited at the foothills.

Illustrative examples
In Appendix C we describe two models-the first shows the effect of transient dynamic topography on the evolution and reorganization of drainage networks; the second shows the influence on the evolution of drainage networks when river systems transition from a detachment-limited to a transport-limited state. These models are based on the simple model described in Section 3.5 and serve to demonstrate how the various modules can be brought together to explore a broad range of geomorphological scenarios. Annotated parameter-files and initial topographies for models used to produce Figs. C.9 and C.10 can be found online in Supplementary Data.

Discussion
In order to demonstrate the parallel scalability of the fluvial transport modules, we have computed models with a similar initial topography as that in Appendix C.1, but with a range of grid densities. Fig. 6A and 6B show strong scaling results for the detachmentlimited and transport-limited modules, respectively. The scaling tests were computed on a desktop machine with two Intel Xeon E5-2650 processors, each with eight cores (2.6 GHz clock-speed). The parallelized solution schemes in both the detachment-limited and the transport-limited modules show near-linear scalability for the densest mesh (1600 × 1600), while scalability gradually drops off with increasingly smaller mesh sizes in both. Our results suggest that an order of magnitude speedup can be achieved on typical desktop machines-thus making continental scale models with tens of millions of nodes more computationally feasible.
Studying continental scale geomorphological evolution over tens of millions of years, particularly in order to better understand the influences of long wavelength but small amplitude dynamic topography, has been computationally prohibitive. Moreover, since paleotopography is poorly constrained over geological timescales, models of geomorphological evolution generally start with synthetic initial topography conditions-thus, the ability to Erosion length-scales need to be specified for both bedrock and sediments in the parameter file-see text for more details. b Diffusivities need to be specified for both bedrock and sediments in the parameter file-see text for more details. c The dimensions of K f depend on the formulation used. carry out systematic model parameter space explorations is critical for examining model sensitivities and uncertainties under different plausible scenarios. The implementation of efficient, parallel algorithms to model fluvial transport in SPGM makes model parameter space exploration more feasible and has the potential to provide new insights into the influence of dynamic topography on landscape evolution over geological time.

Conclusions
In this paper, we present the implementation of a parallel, multi-process [27] numerical model, where physical processes that contribute to mass redistribution are integrated separately to record geomorphological evolution. The modular structure of the code, with clearly defined interfaces between mesh generation, implementation of optimized numerical algorithms and output generation is particularly geared toward ease of adaptability for studying a wide range of geomorphological scenarios. We present simple examples of geomorphological evolution that demonstrate the capabilities of the code while providing guidelines for setting up more complex models for exploring the influences of forcing functions that can vary in space-time.

Acknowledgments
MG was supported by Statoil ASA, NSF EAR -1358646 and EAR -1247022. R.D.M was supported by Statoil ASA through ARC IH130200012. Table A.1 lists physical parameters associated with each module, along with their corresponding names in parameter files. Additionally, annotated parameter-files and initial topographies for models used to produce Figs. 3-5 can be found online in Supplementary Data.

B.1. Verification of nonlinear diffusion algorithm
The Method of Manufactured Solutions (MMS) [28,29] offers a convenient way to verify the accuracy of implementations of numerical algorithms. A major advantage afforded by MMS is the ability to verify codes with nonlinearities for which analytical solutions are not readily available. In the generally applicable version of the method, one simply includes in the code a source term, S, that originates from the choice of a non-trivial, but analytical solution, which also defines boundary conditions-see Roache [29] for a detailed review of the procedure. In other words, one starts with an analytical solution and derives S by passing the solution through the governing equation, in order to satisfy the solution.
For our purposes, we choose Eq. (B.1) as the analytical solution to Eq. (8) and apply Eq. (B.2) as the space-time varying diffusivity field to derive S (Eq. (B.3)), using the symbolic mathematics package, SymPy [30]. We add S to the right hand side of Eq. shows the rate of convergence as a function of grid spacing.

C.1. Effects of plume-induced dynamic topography
We combine the model in Section 3.5 with a transient uplift function that represents the motion of a mantle plume beneath  continental lithosphere (e.g. Braun et al. [31]. Laboratory models of mantle plumes suggest that initiation of plumes involve a phase of rapid uplift, followed by a more gradual decline [32]. Here we approximate the domal rate of uplift (representing dynamic topography from an upwelling mantle plume) by a Gaussian function, centered at the plume axis. The center of the plume moves along a near-diagonal path across the domain and the associated rate of uplift initially increases, followed by a gradual waning of both strength and radial extent of influence. Fig. C.9 shows various model attributes after a 20 Myr period of evolution and animation S1 demonstrates the significant influence that plume-related dynamic topography can exert on the evolution and reorganization of drainage networks that in turn influence patterns of erosion and sediment deposition through time. Preexisting topography evolves at a slower rate compared to regions influenced by the plume, suggesting that our results are qualitatively in agreement with those presented in Braun et al. [31], where it was shown that erosion of dynamic topography increases linearly with its wavelength. It is also important to note that drainage reorganization occurs predominantly in downstream channel networks where channel gradients are gentler, whereas upstream channel networks remain mostly intact over the last 10 Myr period (animation S1). This is likely to have implications for drainage reversals in river systems, induced by transient plume-related dynamic topography.

C.2. Transition from detachment-limited to transport-limited scenario
In order to demonstrate the contrasting effects of detachmentlimited and transport-limited modes of fluvial erosion, we take the initial topography and dynamic uplift history from the model described in Appendix C.1 and parameterize two stages of drainage network evolution. During the first 10 Myr period we impose a detachment-limited state, where all eroded material is removed from the computational domain. Drainage networks then switch to a transport-limited state, potentially due to a change in climatic conditions and eroded sediments are deposited downstream over the subsequent 10 Myr period. This is achieved by running the model using two separate configuration files, each spanning a 10 Myr period of model evolution-the final topography from the first stage is used as an input to the next. Fig. C.10 shows instances of topography at 5 Myr intervals and animation S2 demonstrates the sharp contrast in the morphology of drainage networks over the two separate stages of evolution. During the first 10 Myr period, the plateau front is rapidly eroded away, whereas over the next 10 Myr period, incisive drainage networks reach deep into the interior of the plateau.
While river systems are likely to transition between states more gradually, this example serves to demonstrate that such transitions can be modeled in SPGM and that they can be made more gradual by cascading several intermediate stages, which would feature a varying K f to effect a smoother transition.