The direction of landscape erosion
 ^{1}BlueMarbleSoft, Kyushu, Japan
 ^{2}Department of Computer Science and Technology, University of Cambridge, Cambridge, UK
 ^{1}BlueMarbleSoft, Kyushu, Japan
 ^{2}Department of Computer Science and Technology, University of Cambridge, Cambridge, UK
Correspondence: Colin P. Stark (cstarkjp@gmail.com)
Hide author detailsCorrespondence: Colin P. Stark (cstarkjp@gmail.com)
The rate of erosion of a landscape depends largely on local gradient and material fluxes. Since both quantities are functions of the shape of the catchment surface, this dependence constitutes a mathematical straitjacket, in the sense that – subject to simplifying assumptions about the erosion process, and absent variations in external forcing and erodibility – the rate of change of surface geometry is solely a function of surface geometry. Here we demonstrate how to use this geometric selfconstraint to convert a gradientdependent erosion model into its equivalent Hamiltonian, and explore the implications of having a Hamiltonian description of the erosion process. To achieve this conversion, we recognize that the rate of erosion defines the velocity of surface motion in its orthogonal direction, and we express this rate in its reciprocal form as the surfacenormal slowness. By rewriting surface tilt in terms of normal slowness components and deploying a substitution developed in geometric mechanics, we extract what is known as the fundamental metric function of the model phase space; its square is the Hamiltonian. Such a Hamiltonian provides several new ways to solve for the evolution of an erosion surface: here we use it to derive Hamilton's raytracing equations, which describe both the velocity of a surface point and the rate of change of the surfacenormal slowness at that point. In this context, gradientdependent erosion involves two distinct directions: (i) the surfacenormal direction, which points subvertically downwards, and (ii) the erosion ray direction, which points upstream at a generally small angle to horizontal with a sign controlled by the scaling of erosion with slope. If the model erosion rate scales faster than linearly with gradient, the rays point obliquely upwards, but if erosion scales sublinearly with gradient, the rays point obliquely downwards. This dependence of erosional anisotropy on gradient scaling explains why, as previous studies have shown, model knickpoints behave in two distinct ways depending on the gradient exponent. Analysis of the Hamiltonian shows that the erosion rays carry boundarycondition information upstream, and that they are geodesics, meaning that surface evolution takes the path of least erosion time. Correspondingly, the time it takes for external changes to propagate into and change a landscape is set by the velocity of these rays. The Hamiltonian also reveals that gradientdependent erosion surfaces have a critical tilt, given by a simple function of the gradient scaling exponent, at which raypropagation behaviour changes. Channel profiles generated from the nondimensionalized Hamiltonian have a shape entirely determined by the scaling exponents and by a dimensionless erosion rate expressed as the surface tilt at the downstream boundary.
When geomorphologists describe the evolution of a landform, a direction of erosion is often invoked: for example, we speak of a bank cutting laterally, a cliff retreating, a knickpoint eroding upstream, or a river channel incising down into bedrock. Generally, such statements are taken at face value, and the erosion direction in each case is understood from context, e.g. erosion in a bedrock channel is broadly considered to take place subvertically downwards, hewing closely to gravity, except at knickpoints where it occurs subhorizontally upstream and along the channel walls where it acts subhorizontally and roughly orthogonal to streamflow. At the same time, we recognize that the geomorphic processes driving or mediating erosion are associated with particular directions relative to the geometry of the surface, which presumably has consequences for the direction in which that surface erodes: weathering acts roughly normal to an exposed surface, mechanical abrasion involves obliquely streamwise impacts that can be resolved into normal and tangential components (as can frictional wear by sliding ice or debris), and so on. There are obviously many directions involved in driving the evolution of a landscape, so what can we say about the direction of motion of the erosion surface itself? Our goal here is to answer this question using some concepts and tools of differential geometry and classical mechanics.
1.1 Tracking points on an erosion surface
Tracking the motion of a solid object is easy if the surface of the object is not eroded during motion: all that is needed is to tag the surface with markers and monitor their displacements. This is not possible for a surface undergoing erosion because all such markers are destroyed by the erosion process itself. We can nevertheless describe, in a mathematical sense, how points on an erosion surface move – if we know something about the process of erosion. The purpose of this section is to preview how this task can be achieved and to provide some conceptual context. The ideas outlined here are developed in full in the main body of the paper.
A moving erosion surface has only one intrinsic direction available at each surface point: the local normal to the surface. Describing motion in any other way entails the supply of extra information through the choice of an additional direction as a reference. Since gravity acts downwards, the usual choice is to assign vertical as the reference axis and to express erosion rate as a vertical velocity. On the other hand, for problems such as sea cliff retreat or riverbank erosion it can be more convenient to pick horizontal as the reference. Whatever the choice, basic trigonometry makes it easy to transform an erosion function between any of these geometries (but with a complication; see Sect. 3.1).
The minimal approach therefore avoids supplying a reference direction and treats surface erosion as acting intrinsically in the local normal direction. In light of this, we may be tempted to infer that points on an erosion surface move in the normal direction: in general, however, they do not.
To see why, let us examine a surface evolving by some unknown mechanism. Let us assume for simplicity that the surface is an always smooth 1D line in 2D x–z space (Fig. 1). Mark the surface at time T_{a} and again at a very small time interval later ${T}_{b}={T}_{a}+\mathrm{\Delta}T$. Each surface can be considered as a set of points: T_{a}={a} and T_{b}={b}, where a and b are 2D vectors.
In the absence of an equation of motion, we are free to pair each point a∈T_{a} with any otherwise unpaired point b∈T_{b} (Fig. 1, “free” inset). We could enforce a strict order to the pairings, but we would still have a very large number of choices. What matters is that the motion of surface points from T_{a} to T_{b} is defined by our choice of mappings, and for the moment this choice is arbitrary.
Now, there are two ways to assess the rate of surface motion: one familiar, the other much less so. The familiar quantity is the velocity vector, which we get by measuring the distance between, and direction defined by, each pair of points $\mathit{v}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}:=\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\left(\mathit{b}\mathit{a}\right)/\mathrm{\Delta}T$. The unfamiliar quantity is the normalslowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ (Sect. 3.1), which get by measuring the time ΔT it takes for the surface to move a given distance in its normal (intrinsic) direction. We visualize $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ as a series of small planes emanating from a, parallel to the local tangent to T_{a} and approaching T_{b}. The term “slowness” is used because its units are reciprocal speed; it could also be called the “pace” of erosion, meaning the time needed to erode a reference distance.
This brings us to our key premise: when we specify an erosion function, we are explicitly defining the behaviour of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ but only indirectly obtain the behaviour of v. That is because an erosion rate function measures the time it takes for the surface to move a given distance and not the travel time for points on the surface. If the process of erosion is isotropic, this subtle distinction is moot; if, however, the erosion rate depends on gradient, the distinction is fundamentally important (Fig. 1, “isotropic” and “anisotropic” insets).
We can understand why if we realize that by quantifying the elapsed time between successive erosion surfaces, the erosion rate function actually defines a metric, i.e. a tool for measuring the “length” of the covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$. If the erosion rate depends only on location, meaning that it is independent of surface tilt and thus isotropic, the corresponding metric is Euclidean, which makes $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ and v point in the same direction and leads each point on T_{a} to pair with its nearest neighbour (in a Euclidean sense of the term) on T_{b}. This is the most intuitive way of linking points on one surface to another, but is not correct for erosion in general. That is because if the erosion rate is also a function of gradient, the resulting metric will be anisotropic and nonEuclidean, $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ and v will point in different directions, and the way the metric measures the shortest distance between successive erosion surfaces will no longer be a simple use of Pythagorean geometry. Metrics of this kind – that depend on position and orientation – are called Finsler metrics. They constitute a way to measure travel time between two points when resistance to motion varies with direction in a nontrivial way. Physical analogues include measuring travel time when walking over hills or navigating a boat in a wind. In special cases they may reduce to, or at least incorporate, a Riemannian form.
Transformation of the erosion equation into a metric function takes a few steps. The first is to reparameterize the directional parts of the erosion equation using the components of the slowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}=[{p}_{x},\phantom{\rule{0.125em}{0ex}}{p}_{z}]$ while leaving any spatial dependence untouched. For example, if the erosion function depends explicitly on surface gradient tan β, where β is the angle of the surfacenormal relative to vertical, we can use the substitution $\mathrm{tan}\mathit{\beta}=\left{p}_{x}/{p}_{z}\right$. The normal erosion speed is replaced with the reciprocal magnitude of the slowness covector ${\mathit{\xi}}^{\u27c2}=\mathrm{1}/p=\mathrm{1}/\sqrt{{p}_{x}^{\mathrm{2}}+{p}_{z}^{\mathrm{2}}}$.
If this reparameterization is possible, we get an equation that can be rearranged into the form ${\mathcal{F}}_{*}(\mathit{a},\stackrel{\mathrm{\u0303}}{\mathit{p}})=\mathrm{1}$. This ℱ_{*} is a fundamental metric function, which measures the shortest time interval for the surface to erode a unit distance in a given direction. Among several special properties exhibited by this function, the crucial one is its order1 Euler homogeneity in $\stackrel{\mathrm{\u0303}}{\mathit{p}}$, which means that ${\mathcal{F}}_{*}(\mathit{a},\mathit{\lambda}\stackrel{\mathrm{\u0303}}{\mathit{p}})=\mathit{\lambda}{\mathcal{F}}_{*}(\mathit{a},\stackrel{\mathrm{\u0303}}{\mathit{p}})$.
Squaring and scaling the metric function defines a quadratic Hamiltonian $\mathcal{H}(\mathit{a},\stackrel{\mathrm{\u0303}}{\mathit{p}})\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}:=\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathcal{F}}_{*}^{\mathrm{2}}/\mathrm{2}=\mathrm{1}/\mathrm{2}$, which is the key result of this study. This “geomorphic” Hamiltonian provides us with equations of surface motion in the form of Hamilton's equations, which allow ray tracing and other methods to be used to solve for landscape evolution. It tells us not just that surface points move according to a Hamiltonian flow but also that they follow geodesic paths, i.e. paths of shortest erosion time.
Point velocities, and therefore point pairings, {a,b}, are given by one half of Hamilton's equations: differentiating the Hamiltonian by each of the erosion slowness covector components p_{x} and p_{z} in turn, we get a vector expressing the change of point position with time: $\mathit{v}=\partial \mathcal{H}/\partial \stackrel{\mathrm{\u0303}}{\mathit{p}}$. It follows from order1 homogeneity of the metric function ℱ_{*} that the surfacenormal slowness covector and this point velocity vector must always be conjugate, $\stackrel{\mathrm{\u0303}}{\mathit{p}}\cdot \mathit{v}=\mathrm{1}$.
Earlier, we asserted that surface points do not, in general, move in the surfacenormal direction, and now we have proof. Exploiting conjugacy, we can measure the angle ψ between the surfacenormal and the point velocity using their dot product $\mathrm{cos}\mathit{\psi}=\stackrel{\mathrm{\u0303}}{\mathit{p}}\cdot \mathit{v}/\left(pv\right)$. If the rate of erosion depends on surface tilt β, the corresponding metric function and Hamiltonian will both depend, in some nonlinear fashion, on the normal slowness components p_{x} and p_{z}, and so $\partial \mathcal{H}/\partial \stackrel{\mathrm{\u0303}}{\mathit{p}}$ and point velocity v will not in general be colinear with the surface normal. A gradientdependent erosion process is therefore anisotropic, and its degree of anisotropy is measured by the angle ψ.
The practical consequence of erosion driving anisotropic Hamiltonian flow lies in how it controls the propagation of information, in the sense of initial and boundary conditions, into a landscape. Each element of this Hamiltonian flow has both a point position a and a normal slowness $\stackrel{\mathrm{\u0303}}{\mathit{p}}$, i.e. each element contains information about the location and orientation of the surface and its reciprocal rate of erosion. Progression along the Hamiltonian flow occurs along successive point pairings; each pairing translates an element in space while carrying (and to some extent modifying) the surface information along with it. The angular disparity between the direction of information transfer (i.e. point velocity) and the intrinsic direction of surfacenormal motion is the anisotropy ψ.
1.2 Structure of the paper
The paper is organized into eight sections and a set of appendices. Section 2 summarizes how erosion in three dimensions (3D) can be tracked using implicit surfaces and level sets, makes a connection with the Hamilton–Jacobi equation, and demonstrates the natural link with Hamiltonian methods. Section 3 combines these concepts with those introduced in Sect. 1.1 and formulates a Hamiltonian theory of gradientdriven erosion (for a 2D slice of 3D space). It explores this theory using geometric mechanics and differential geometry and reveals how strong anisotropy lies at the heart of landscape surface evolution. Section 4 implements the geomorphic surface Hamiltonian using a particular model of gradientdriven erosion – an adaptation of the stream power incision model to handling erosion in the surfacenormal direction – and presents a nondimensionalization of the Hamiltonian and Hamilton's raytracing equations and a simple means of model solution. Section 5 shows how to use Hamiltonian ray tracing to obtain model surface solutions for a domain whose boundaries are subject to a constant vertical erosion rate. Section 6 discusses these numerical solutions and examines what they have to tell us about erosional anisotropy and the notion of two distinct directions of landscape erosion. It also relates model scales to realworld landscape time, space, and velocity scales. Section 7 looks at the broader implications of the Hamiltonian approach to erosion, and Sect. 8 draws some conclusions. Appendices A–F draw on disparate literature sources linked together here for the first time and use them to shed light on the theory presented in this paper.
2.1 Landscape as an implicit surface
In almost every model treatment of landscape erosion (Coulthard, 2001; Dietrich et al., 2003; Fowler, 2011; Pazzaglia, 2003; Tucker and Hancock, 2010; Tucker, 2015; van der Beek, 2013; Willgoose, 2005), the shape of the land surface in 3D space is written mathematically as a function of elevation h parameterized by the 2D horizontal coordinates {x,y} of points on the surface and by the time t at which the point elevations are assessed. In other words, the landscape is described by an explicit surface function $h(x,y;t)$.
An explicit surface description has advantages and disadvantages. On the plus side, theoretical development is relatively simple because it effectively involves only two spatial dimensions and because the numerical solution can be carried out on a 2D grid. On the negative side, the rate of erosion is only tracked in the vertical direction, through the partial derivative of elevation with time $\partial h/\partial t$: if there is any horizontal component of erosion, it is not tracked directly and has to be calculated indirectly using the lateral variation in elevation ∇h. Problems arise when the surface gradient becomes very steep, for example at knickpoints or channel banks, and any development of overhangs is obviously impossible.
If we instead describe the landscape using an implicit surface, many of these issues are eliminated. The price is greater complexity in the mathematics needed to formulate surface motion and to resolve it numerically. The extra cost is worth paying if it leads to greater insights into how landscapes form.
2.2 Landscape as the 2D zero contour of a 3D function
An implicit surface in 3D space is the set of points $\mathit{\left\{}x\right(t),y(t),z(t\left)\mathit{\right\}}$ that define the 2D “contour” or levelset surface of a function ϕ:
where ϕ is a nonlinear function defined at all points across the 3D domain of interest that varies with time and is nonlocal – in the sense that it can be a function of curvature or of values of itself at a distance. Put more simply, ϕ is a very flexible function that can be tailored to induce whatever surface motion is desired.
The term “implicit” is used because surface positions are not specified directly; instead, a surface is defined by “slicing” the function $\mathit{\varphi}(x,y,z)$ at some chosen value ϕ_{0} and finding positions $\mathit{\{}x,y,z\mathit{\}}$ for which $\mathit{\varphi}(x,y,z)={\mathit{\varphi}}_{\mathrm{0}}$. Think of how a visualization tool for a 3D scalar field, such as temperature, works: sequential slicing across a range of temperatures provides an animated view of its variation throughout a volume. This variation can be complex, revealing isolated blobs of high (or low) values that connect in topologically complicated ways as the slicing threshold temperature is changed. In this way, an implicit description of a surface can represent complex, multivalued geometry and topology without extra mathematical work.
Landscape evolution can be modelled with an implicit surface by writing an equation to drive evolution of the function ϕ and watching how its zero levelset $\mathit{\varphi}={\mathit{\varphi}}_{\mathrm{0}}=\mathrm{0}$ implicitly prescribes changes in surface positions $\mathit{\left\{}x\right(t),y(t),z(t\left)\mathit{\right\}}$ over time.
2.3 The levelset equation
Implicit surface motion in its most general form is described by the levelset equation (Gibou et al., 2018; Giga, 2006; Osher and Fedkiw, 2001, 2003; Sethian, 1999; Vladimirsky, 2001), in which $\mathit{\varphi}(x,y,z;t)$ is a 3D function constructed so as to evolve over time t with a velocity ξ, a vector function that in general varies with position and time, is possibly nonlocal, and only needs be defined where ϕ=0:
This is equivalent to holding the material derivative of the scalar field ϕ – driven to move by the vector field ξ – at zero along the zero contour of ϕ but otherwise allowing it to vary unconstrained.
Only the normal component ξ^{⟂} of the implicit surface velocity plays any role in driving motion: in the geomorphic case, this would be the surfacenormal erosion rate. So we can write the following equation:
The notation ξ^{⟂} is adapted from Osher and Merriman (1997).
This equation provides a very generic description of how a 2D surface evolves in 3D space, in the sense that it defers all description of processes into the formulation of the surfacechange rate function ξ^{⟂}. This function can readily treat topographic gradient and curvature and substrate erodibility; suitably provided with coupled process equations, it could also incorporate water flow depth and velocity, intermittent sediment cover, development of a vegetation layer, spatiotemporal precipitation, tectonic displacement, and so on. Such flexibility, however, is not our goal here. Instead, we seek geometric insights into the process of landscape erosion, which we can achieve if we limit the scope of this equation and make ξ^{⟂} a simplified function of local gradient and accumulated flow. A geomorphic levelset equation in this form makes it easier to tease out its fundamental behaviour.
2.4 Motion described by the Hamilton–Jacobi equation
If we restrict the surface velocity ξ to be a local function of position and time, Eq. (3) becomes the Hamilton–Jacobi equation (HJE):
where each vector r tracks a point as it moves from one zero level set of ϕ to another with velocity $\dot{\mathit{r}}=\mathrm{d}\mathit{r}/\mathrm{d}t$, while the front itself at that point moves in the direction ∇ϕ. These directions are not necessarily the same.
The HJE is a firstorder partial differential equation that plays a central role in classical mechanics (Arnold, 1989; Goldstein et al., 2000; Houchmandzadeh, 2020; Small and Lam, 2011; Whitham, 1999). Its driver is the Hamiltonian ℋ, which combines the surface velocity ξ with the gradient ∇ϕ in a way that lends it special properties.
The Hamiltonian in the HJE is required to be a local function – in the sense that it can depend on position r and instantaneous time t but cannot depend on the shape of the propagating surface at some distance away or on any historydependent quantities. Diffusive and quasidiffusive processes are not allowed either. However, viscosity solutions of the HJE (Crandall and Lions, 1981), which are the standard means of resolving profound mathematical challenges with this equation, ironically involve the addition of a weak, ultimately vanishing, secondorder term that can be considered a diffusive process at the subgrid scale.
2.5 Landscape as an erosion arrivaltime surface
If we wish to use the HJE to treat landscape evolution in terms of an implicit function, we need to consider how to write a Hamiltonian form of the erosion function driving that evolution. If this Hamiltonian is independent of (i.e. does not change with) time, it simplifies into a static HJE or eikonal equation ℋ(r,∇ϕ). The implicit surface function ϕ that solves this static HJE is a singlevalued, 3D function that defines the position and shape of arrival time surfaces. In other words, ϕ can be thought of as a first arrival time function $T(x,y,z)t$, where T defines the locus of surface points $\mathit{\{}x,y,z\mathit{\}}$ that satisfy at each time step t the equation
Another way to express this is to say that the contours of T are 2D isochrone surfaces embedded in 3D space that define the shape of the landscape as it changes.
In the eikonal equation, the Hamiltonian is a constant function of surface point position r and the gradient of the arrival time ∇T with the simple form:
Points on the surface move with velocity vector $\mathit{v}=\dot{\mathit{r}}$, while the surface itself moves with a slowness covector given by $\stackrel{\mathrm{\u0303}}{\mathit{p}}=\mathbf{\nabla}T$). It is important to emphasize that $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ is not a vector. Section 3.1 goes into more detail as to what is meant by the term “covector” and why the distinction is consequential.
Although both v and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ are both directional quantities describing surface motion, they only point in the same direction if the motion mechanism is isotropic. Measuring their angular disparity is the key to assessing the anisotropy. One of the aims of this study is use this measure to reveal the strong anisotropy of landscape erosion processes (see Sect. 3.18).
In this section, we formalize the ideas presented above into a Hamiltonian theory of erosion front motion. First, we provide a gentle introduction to the pivotal concept of a covector (Sect. 3.1) and show how useful it is for treating the direction and reciprocal speed of the propagating front. Then we show that the gradient of the surface arrival time is itself a covector (Sect. 3.2). Next, we make the case that the geomorphic processes driving erosional motion of a topographic surface can be represented by local functions (Sect. 3.3) parameterized by the surfacenormal covector, and how they constitute, broadly speaking, a form of geometric selfconstraint (Sect. 3.4). After imposing a gradientdependent form on the erosion function (Sect. 3.5), we show how the above ingredients lead, via the fundamental metric function, to a Hamiltonian description of erosion (Sects. 3.6–3.9). Next we delve into the connections between the fundamental function and erosional wavelets and use them to provide a graphic explanation of Huygens' principle as applied to erosion surface propagation (Sect. 3.10). We then express the equivalent Fermat's principle in terms of the variational path of least action (Sect. 3.11) to show that a point on the surface follows the path of least erosion time. This leads on to derivation of Hamilton's raytracing equations (Sects. 3.12–3.13) and a discussion of some of their properties (Sects. 3.14–3.15). A verification that the Lagrangian is constant (Sect. 3.16) follows. Then we discuss ray angles, their behaviour relative to surface tilt, and the existence of a critical tilt at which raypropagation behaviour changes (Sect. 3.17). This leads to an exploration of how the disparity between the two directions of erosion is a measure of erosional anisotropy (Sect. 3.18). Finally (Sect. 3.19), we look at the various ways the evolving surface tilt can be tracked in the model. Nondimensionalization is undertaken in Sect. 4.
Note that we use superscripts for contravariant tensor components (e.g. r^{x}) and subscripts for covariant tensor components (e.g. p_{z}); the Einstein summation convention (summing over similar tensor components) is adopted for brevity. Symbol usage is summarized in Table G.
3.1 Tracking erosion with covectors
Imagine a locally planar surface undergoing constant erosion (Fig. 3), where the surface tilt angle is β and the vector r takes values that lie along the erosion surface at a given time T(r). As time passes, erosion moves the surface progressively further into the substrate. Taking snapshots at regular intervals ΔT generates a uniformly spaced sequence of surfaces which we call erosional isochrones. These isochrones are level sets or contours of the arrival time function T. In Fig. 3, the time interval is chosen to be ΔT=1 year, and isochrones have been plotted for $T\left(\mathit{r}\right)=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\mathit{\}}$ years.
Let us fix the point of interest r at the location shown in Fig. 3. Here the surfacenormal rate or speed of erosion is ${\mathit{\xi}}^{\u27c2}=\mathrm{0.25}$ mm yr^{−1} and surface tilt is $\mathit{\beta}=\mathrm{60}{}^{\circ}$. Written as a vector, the erosion rate is as follows:
with a direction normal to the surface and an angle $\mathit{\beta}=\mathrm{60}{}^{\circ}$ to the vertical; its length or magnitude is the surfacenormal erosion rate:
Ideally, we should only have to compute the sine and cosine components to the erosion velocity vector ξ to get the horizontal and vertical rates of erosion. However, the vertical trigonometric component ξ^{z} does not equal the (negated) vertical rate of erosion ξ^{↓} (Eq. 10), nor does the horizontal trigonometric component ξ^{x} equal the horizontal rate of erosion ξ^{→}:
It seems almost too trivial to ask, but why does naive application of trigonometry let us down here? The answer lies in the fact that we have written the erosion rate as a vector: we should instead express it as a covector.
Consider $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ in Fig. 3, which can be written as a function with singlerow matrix form:
This scalar function takes as input a vector such as n and returns the number of isochrones crossed by that vector. Here n is the surfacenormal unit vector
Because we employ units of millimetres and years here, n has a length of $\left\mathit{n}\right=\mathrm{1}$ mm. Over this distance n crosses four 1year isochrones, and thus we obtain
Now consider the vertical component of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ (which is negative here) acting on n: counting downwards over a distance ${n}_{z}=\mathrm{1}/\mathrm{2}$ mm, we find one isochrone crossing:
The horizontal component of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ counts three isochrone crossings by the unit normal vector counting rightwards over a distance ${n}_{x}=\sqrt{\mathrm{3}}/\mathrm{2}$ mm:
These components can be added together because $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ is a linear function; this summation gives
which is the count of four we found by measuring along n directly. The count can of course take any real (fractional) value: for clarity, the example here has been constructed so as to yield round numbers.
The function $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ is called a “oneform” in the terminology of differential geometry, and instances of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ are called covectors. In general, a oneform operates on a vector and returns a scalar. Here, $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ takes in a unit vector and returns the slowness of erosion in the direction of that vector. In optics and seismology, $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ is known as the normal slowness; in classical mechanics it is called the generalized momentum. In a geomorphic context, this normal slowness can be interpreted as the maximum isochrone density, and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ the isochrone density covector, in that when applied to the unit normal vector n it calculates the maximum number of isochrones to be found in any direction from that point.
The slowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ is a more convenient measure of erosion rate because its sine and cosine components are the horizontal and vertical slownesses, which are (respectively) the reciprocal rates of erosion horizontally and vertically.
The magnitude of the covector here is the normal erosion slowness, i.e. the reciprocal erosion rate, and is given by
and surface slope is
In other words, by describing the rate of surface motion with an erosion slowness covector, instead of an erosion velocity vector, we can assess its variation with direction much more easily. Fundamentally, a covector is the correct way to represent motion of a surface at a given point, and a vector is the appropriate way to represent the position and motion of that point. See Appendix B for more details.
3.2 Gradient is a covector
The erosion slowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ has another facet: it is also the gradient of the arrivaltime function T. To see why, consider again Fig. 3 and its level sets of T at discrete intervals. These level sets are isochrones or contour surfaces of equal arrival time $T\left(\mathit{r}\right)=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots}\mathit{\}}$, which are represented schematically as simple straight lines in this figure. They successively increase in the direction of the normal vector n.
If we measure (in Fig. 3) the change in T in the x direction over a distance ${n}^{x}=\sqrt{\mathrm{3}}/\mathrm{2}$, we find that n^{x}dT_{x}=3. Similarly, if we measure the change in the −z direction over a distance ${n}^{z}=\mathrm{1}/\mathrm{2}$, we get n^{z}dT_{z}=1. In general terms,
and in this example we find
which is the normal slowness obtained in Eq. (13) written as a differential oneform. In other words, the rate of change dT(⋅) over a unit distance in the isochronenormal direction n is given by dT(n), and the isochrone or contour density dT(n) in the contournormal direction is the same as the covector magnitude p. We can now invoke the gradient operator ∇ and have
which says that the Euclidean gradient of the arrival time T of the erosion surface is the normal slowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$.
3.3 Modelling erosion in the surfacenormal direction
If we wish to frame a model of landscape evolution in terms of geometric mechanics, we need to employ the following three elements: (i) an implicit function to track the evolving landscape surface geometry; (ii) a surfacenormal erosion slowness covector, corresponding to the gradient of the implicit function, that encodes the reciprocal rate of motion of the surface; and (iii) an erosion model for the surfacenormal speed of erosion that can be parameterized using the slowness covector.
To supply the third element, we can write a generic model for the surfacenormal speed of geomorphic erosion that is a solely function of local fluxes and gradient:
Some erosion phenomena, such as quasidiffusive processes like rain splash, cannot be modelled under this local restriction, but this is a minor loss. Henceforth, the only flow we will consider is kinematic water flow resulting from spatially uniform rainfall runoff, and we will ignore complexities such as storm hydrograph cycles and the effects of sediment supply, transport, and cover.
A model in this form is not unambiguously local: its dependence on accumulated water flow presupposes a dependence on upstream catchment geometry; any change in catchment geometry, through motion of drainage divides, acts to change flow at distant points downstream. A fundamentally important assumption here is that divide motion is slow enough for the erosion equation to be considered effectively local. The validity of this assumption is discussed at the end of Sect. 7.1.
3.4 Erosion imposes a geometric selfconstraint
The process of landscape evolution represented by Eq. (23) is a kind of geometric straitjacket or geometric selfconstraint – in the sense that it essentially says the landscape obeys the following equation:
In other words, the shape of the landscape determines the patterns of surface flow and thereby the fluxes of material over the surface, and it mediates the effectiveness of these fluxes through its control of the gradients; these effects combine to set the rate at which the shape of the landscape changes: in short, change in landscape geometry is controlled by landscape geometry. This conclusion applies even if the erosion process is not spatially local.
The consequence of this geometric selfconstraint is that, at its heart, geomorphic erosion is driven by a particular kind of Hamiltonian. This Hamiltonian arises from how points on an erosion surface “see” (for want of a better term) their shortest path of erosion to the next set of surface points at little time later. The sections below explore this assertion in detail.
3.5 Separable, gradientdependent erosion rate model
The Hamiltonian approach developed here can in principle be applied to any erosion rate model, with the proviso that the bedrock surface can only undergo erosion, meaning that its motion must always be positive ${\mathit{\xi}}^{\u27c2}>\mathrm{0}$. If transient sediment deposition and bed cover are to be modelled, meaning that topographic elevation (in the bedrock reference frame) can rise as well as fall, alluvial geometry needs to be tracked as an additional model variable along with bedrock surface position. The resulting Hamiltonian would not be static, and the dimensionality of its phase space would be comparatively large. Such sophistication will eventually be needed, as models of this kind become the standard (e.g. Dietrich et al., 2003; Sklar and Dietrich, 2006; Zhang et al., 2015). However, in this introduction of geometric mechanics to the task of modelling erosion, we choose to avoid such complexity and instead settle on an erosion equation that (i) is a nonlinear (power) function of (space–time variable) rock surface gradient tan β(x,t); (ii) has a separable form, with spatial variables (constant in time) such as flow velocity and depth, sediment concentration, substrate erodibility, and the abrasion process itself aggregated into a separate multiplicative term φ(x); and (iii) describes the speed of erosion ${\mathit{\xi}}^{\u27c2}(\mathit{x},t)$ in the rocksurfacenormal direction:
Note that surface tilt relative to vertical is expressed as sin β rather than tan β because erosion rate is measured in the normal rather than the vertical direction. In a further simplification, we restrict the model to a 2D transect (Fig. 2).
3.6 The erosion equation in Hamiltonian coordinates
Covectors are an essential ingredient in the construction of a Hamiltonian framework for surface erosion. As we will show in the coming sections, the Hamiltonian endows each point on the surface at position r with an associated tangent covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ that represents the normal slowness of the surface at that point. The components of r and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ correspond to the axes of the phase space inhabited by the Hamiltonian.
Since our model is restricted here to a 2D transect of 3D Euclidean space, this Hamiltonian phase space is 4D; two of its four axes are spanned by the two components of the position vector, and the remaining two by the slowness covector components:
The Hamiltonian parameters $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ are coordinates in what, in mechanics, is usually called momentum phase space, and in differential geometry is called a cotangent bundle; we henceforth refer to this as the slowness phase space since momentum has no meaning in the current context. It has a dual, called the velocity space, or tangent bundle, where the Lagrangian corresponding to the geomorphic surface Hamiltonian is defined.
Reiterating Eq. (18) and reducing it to express the surface tilt angle β explicitly, we have
noting that p_{x}>0 and p_{z}<0 for the halfdomain shown in Fig. 2. Each point in phase space acts entirely independently.
The erosion equation (Eq. 25) is now easy to convert into a form parameterized by the components of r and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$:
This equation defines the surfacenormal reciprocal rate of erosion along a 2D profile, written in a form that neatly expresses the geometric selfconstraint inherent to the geomorphic erosion process. This selfconstraint is parameterized by vector position (r^{x},r^{z}) and covector normalslowness (p_{x},p_{z}), which respectively locate a particular point on the surface and encode the reciprocal speed of erosion orthogonal to the surface at that point.
3.7 The fundamental function
What we need to do now is reparameterize Eq. (28) to express the degree to which a coordinate $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ satisfies the geometric selfconstraint imposed by this equation. This is easily achieved using Okubo's technique (Antonelli et al., 1993; Bao et al., 2000; Shimada and Sabau, 2005; Yajima and Nagahama, 2009; Yajima et al., 2011), in which the covector parameter is scaled by a positive function ${\mathcal{F}}_{*}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$:
and substituted back in, rearranging to make ℱ_{*} the subject
The function ℱ_{*} is known as the fundamental (metric) function (see Appendix C; note that an asterisk in used in ℱ_{*} for reasons that will become clear in Sect. 3.9). It is also a Hamiltonian, and as such it is associated with a phase space defined by the four coordinate components $({r}^{x},{r}^{z},{p}_{x},{p}_{z})$. The subset of this 4D space whose locations satisfy the erosion equation given by Eq. (28) must meet the condition:
The power of a Hamiltonian comes from being able to trace a sequence of $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ across phase space for which this criterion holds continuously – a procedure otherwise known as solving Hamilton's equations – which yields the evolution over time of a single point on an erosion surface. However, for technical reasons (Sect. 3.8) it is best not to use ℱ_{*} directly as the geomorphic surface Hamiltonian; a little more work is needed.
To clarify the behaviour of ℱ_{*}, consider the combined meaning of Eqs. (30) and (31). The value of ℱ_{*} at a location in phase space with coordinates $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ is equal to the normal slowness $\sqrt{{p}_{x}^{\mathrm{2}}+{p}_{z}^{\mathrm{2}}}$ implied by that coordinate, i.e. its reciprocal erosion rate, multiplied by the erosion rate determined by the erosion process $\mathit{\phi}\left({r}^{x}\right)\phantom{\rule{0.125em}{0ex}}{{p}_{x}}^{\mathit{\eta}}/({p}_{x}^{\mathrm{2}}+{p}_{z}^{\mathrm{2}}{)}^{\mathit{\eta}/\mathrm{2}}$ acting at that coordinate. This product – of speed times slowness – is obviously equal to one for locations in phase space that represent geomorphically valid surface points in real space. All other locations of phase space are unphysical because at these values of $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ the erosion rate is not reciprocal to the erosion slowness, and this product is not equal to one.
3.8 The geomorphic surface Hamiltonian
The problem with using ℱ_{*} as a Hamiltonian is its order1 Euler homogeneity: functions of this type generate a metric tensor whose determinant is singular, meaning that the tensor cannot be inverted (e.g. Červený, 2002). This puts the Legendre transform and the Lagrangian out of reach. Fortunately there is a simple solution: just use the fundamental function in its squared form and define the geomorphic surface Hamiltonian as follows:
A prefactor of $\frac{\mathrm{1}}{\mathrm{2}}$ is included to make subsequent derivations tidier.
This quadraticform Hamiltonian has the advantage that it is order2 Euler homogeneous:
which makes its metric tensor nonsingular (if η≠1) and the Legendre transform feasible.
We know from Eq. (31) that ${\mathcal{F}}_{*}=\mathrm{1}$ for trajectories across slowness phase space that correspond to physically viable behaviour of surface points. So we can assert that the Hamiltonian is static and has the value
for solutions of the erosion equation. In more concrete terms, we can say that an arbitrary surface point located at r can only represent a point on an eroding surface if its associated orientation and slowness $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ satisfies this equation.
3.9 The geomorphic surface Lagrangian
The quadratic Hamiltonian $\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ has a dual quantity called the Lagrangian ℒ(r,v), which operates in a counterpart space spanned by coordinates giving the position r and velocity v of evolving points on the erosion surface. By symmetry, the Lagrangian is also the quadratic of a fundamental function, denoted ℱ. This function ℱ is the dual of ℱ_{*} and is similarly order1 homogeneous. Its quadratic ℒ is similarly order2 homogeneous:
To make the link between the spaces of ℋ and ℒ, we recognize that the normal slowness covector can be defined as the derivative of the Lagrangian with respect to the velocity coordinate
This is known as the “fibre derivative”.
Mapping from the Hamiltonian ℋ to the Lagrangian ℒ (and vice versa) exploits this property and is achieved with the Legendre transform:
A closed form for ℒ requires several more pieces of the puzzle before it can be derived, and the eventual equation is unwieldy. The contrasting simplicity of ℋ (Eq. 32) is why we prioritize the Hamiltonian over the Lagrangian in this paper.
In due course we will show that the dual fundamental function and the corresponding Lagrangian have constant values ℱ=1 and $\mathcal{L}=\frac{\mathrm{1}}{\mathrm{2}}$, in symmetry with ${\mathcal{F}}_{*}=\mathrm{1}$ and $\mathcal{H}=\frac{\mathrm{1}}{\mathrm{2}}$. Such constancy means that the Lagrangian does not vary with time and that the mutual variation of position r and erosion velocity v is tightly constrained.
3.10 Erosional wavelets and Huygens' principle
Geometric optics provides a way to visualize the Lagrangian and its relationship to the Hamiltonian (Figs. 4 and 5). Motion of an erosion front obeys Huygens' principle: we can imagine each point on the front generating a tiny erosional wavelet, and the coalescence of these wavelets forming the next erosion front. The shape of each erosional wavelet is defined by ℱ. Each shape is a velocity indicatrix giving the radial variation of ray velocity v at a point r or equivalently giving the distance that a point on the surface will erode in an infinitesimal interval.
Figure 4 visualizes a single erosional wavelet, its relationship both to the current erosion front at T=t and to the next at $T=t+\mathrm{\Delta}t$, the particular ray increment vector for which $\mathcal{H}=\mathcal{L}=\frac{\mathrm{1}}{\mathrm{2}}$, and the conjugate relationship of this vector to the frontnormal covector (see Sect. 3.15). Motion of the surface T(r)=t at point r over the interval Δt can be viewed in two mutually consistent ways: (i) the front moves a distance $\mathrm{\Delta}t/p$ in the surfacenormal direction given by $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ or (ii) the point moves a distance Δr=vΔt in the ray direction r. These directions are quite different because the erosion process is strongly anisotropic (Sect. 3.18).
Unconstrained, the point at r could be displaced onto any of the points along the erosional wavelet {Δr}. However, the only valid motion is onto the point r+Δr where the tangent to the wavelet curve is orthogonal to the front increment $\stackrel{\mathrm{\u0303}}{\mathit{p}}\mathrm{\Delta}t/{p}^{\mathrm{2}}$, i.e. the ray and front increments are conjugate to each other (Sect. 3.15).
When erosional wavelets at points along the surface are aggregated, moving T(r) onto T(r+Δr) as shown in Fig. 5, the result is anisotropic front motion that obeys Huygens' principle. The new front can also be found simply by propagating the old front a distance $\mathrm{\Delta}t/p$ in the direction $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ at each point r.
3.11 Fermat's principle as a least action integral
Huygens' principle emphasizes HJE solution in terms of propagation of a front. Fermat's principle, on the other hand, emphasizes solution in terms of tracing the trajectories of points along that front. These two principles are equivalent or dual (Holm, 2011; Houchmandzadeh, 2020; Small and Lam, 2011). Fermat's principle says that these trajectories are paths of stationary travel time: each trajectory obeys a variational principle which ensures its travel time is extremized; this extremal is almost always a minimum. The geomorphic equivalent is the principle that the path of erosion through a substrate from one point to another is the shortest route given the erodibility of the material and its anisotropy and inhomogeneity.
This principle is expressed mathematically by writing an action functional S_{γ}, in terms of the static Lagrangian ℒ, for the set of all possible paths {γ(t)} that a point on the erosion surface might take between two fixed points a=γ(t_{a}) and b=γ(t_{b}):
Note that the integrand ℒ is independent of time t and is a parametric function of positions along γ only. The path actually taken γ_{0} is the path for which the variation of the action is stationary:
For paths traced across the velocity space to which the geomorphic surface Lagrangian ℒ belongs, we can be sure that the action is minimized. Since ℒ is independent of t, we can deduce that γ_{0} is (locally) the path of the least erosion time. Such paths are known as geodesics.
In summary, by expressing a local erosion equation as a geomorphic surface Hamiltonian, converting it into its dual Lagrangian form, and writing the consequent variational principle as the minimization of an action functional for paths across velocity space, we can conclude that points on an erosion surface follow the shortest (in terms of erosion time) possible paths in real space. The next section derives Hamilton's raytracing equations from the Hamiltonian: integration of these rays across slowness phase space generates identical paths of shortest erosion time in real space.
3.12 Derivation of Hamilton's raytracing equations
The fundamental function ℱ_{*} generates a slowness phase space spanned by r and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ on which the geomorphic surface Hamiltonian $\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ operates, and we have a simple expression for ℋ given by Eq. (32). We inferred the existence of a dual fundamental function ℱ that generates a velocity space spanned by r and v on which a Lagrangian ℒ(r,v) operates, but we have yet to obtain expressions for ℱ and ℒ. We can nevertheless make use of the Lagrangian to derive equations of motion for the erosion surface that operate on the slowness phase space. These are called Hamilton's equations.
Our starting point is to examine the differentials of ℋ and ℒ and to compare them. The geomorphic surface Hamiltonian defined in Eq. (32) is static, meaning that it is constant over time, so its differential is
The differential of its counterpart Lagrangian ℒ(r,v) is
Substituting the “fibre derivative” form of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ in Eq. (36) into this equation, and adapting the terms in p_{i}, gives
Rearranging, we have an equation that contains the Legendre transform given in Eq. (37):
Consequently we have a second expression for the differential of ℋ:
Equating the terms in dℋ defined by this equation with those in Eq. (40), we obtain
The next step is subtle but important. Every coordinate (r,v) in velocity space is (potentially) an initial position and velocity for a point on some initial erosion surface. Similarly, every coordinate $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ in slowness phase space is (potentially) an initial position, surface orientation, and reciprocal surfacenormal erosion rate for a point on that initial erosion surface. However, most such phase space coordinates do not correspond to realworld points lying on physically viable paths {γ_{0}} that obey the principle of least erosion time established in Eq. (39). Conversely, for the locations in phase space that do lie on a paths of least action, we can write
Returning to the variation integral in Eq. (39), we can integrate by parts and simplify using the above result to get
The term in brackets [⋅] vanishes because a and b are fixed points associated with limit times t_{a} and t_{b} at which δr^{i}=0. The remaining integral gives the Euler–Lagrange equations for erosional surface motion:
Substituting Eqs. (36, and (46) into the two linking equations in Eq. (45), we obtain Hamilton's equations:
3.13 The meaning of Hamilton's equations
Hamilton's equations are coupled firstorder ordinary differential equations (ODEs) whose integration gives the motion of a single point on an erosion surface in terms of a trajectory across slowness phase space. Each point along the trajectory has phase space coordinates of position r^{i}=r (also the position in real space) and normal slowness covector ${p}_{i}=\stackrel{\mathrm{\u0303}}{\mathit{p}}$ (which encodes both the local tilt of the erosion surface and its reciprocal rate of erosion $p=\mathrm{1}/{\mathit{\xi}}^{\u27c2}$). If we aggregate the trajectories of a set of points from an initial surface we have the motion of the whole surface. This method of front tracking is called ray tracing.
The differential equations in Eq. (49) define the rates of change of the coordinates $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ in terms of the gradient components of the Hamiltonian. Since the Hamiltonian is a constant $\mathcal{H}=\frac{\mathrm{1}}{\mathrm{2}}$ along a ray or trajectory (Eq. 34), motion across the phase space must follow coordinates $(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ that trace a contour of ℋ. This is achieved by moving r in the direction $\partial \mathcal{H}/\partial {p}_{i}$ and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ in the direction $\partial \mathcal{H}/\partial {r}^{i}$, which is to say, orthogonal to the Hamiltonian gradient.
Hamilton's equations take concrete form if we substitute the expression for ℋ in Eq. (32) into Eq. (49). Since the model is limited here to 2D we have four coupled ODEs: two for the component rates of change of position,
and two for the component rates of change of normal slowness,
Raytracing solutions of Hamilton's equations are illustrated in Figs. 6, 12, 13, 14, and 17.
3.14 Constancy of the vertical erosion rate along a ray
The erosion model defined in Eq. (25) is independent of elevation. This makes the Hamiltonian ℋ independent of the vertical coordinate r^{z}, which leads to the zero element in $\dot{\stackrel{\mathrm{\u0303}}{\mathit{p}}}$ in Eq. (51), i.e. the vertical component of erosion slowness is constant:
This is a manifestation of Noether's theorem (Holm, 2011; Noether, 1971), which states that a continuous symmetry in the action implies a conservation law for the Euler–Lagrange equations. Here, we have symmetry with respect to r^{z} in ℋ, and therefore in ℒ, which implies a law of conservation of vertical slowness for the raytracing equations, i.e. that p_{z} must be conserved along a ray. Inasmuch as normal slowness can be crudely equated with the concept of momentum in classical mechanics, we have a “law of conservation of vertical momentum”. Similar conservation laws limited to particular coordinate directions arise in geometric optics (Holm, 2011).
This property simplifies the task of ray tracing by reducing the number of coupled ODEs in the numerical integration from four to three. Moreover, this constancy has the profound implication that the initial rate of vertical erosion ${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$ of a point is carried unchanged along its ray trajectory as the surface to which it is attached moves:
As such, each ray propagates information about the initial surface erosion rate upstream into the landscape until such time as it is destroyed at a cusp (which includes drainage divides, e.g. Fig. 6). Meanwhile the horizontal erosion rate can and does change along the ray because the horizontal component of the slowness covector p_{x} evolves as the surface erodes (Eq. 51).
3.15 Conjugacy of point velocity and front slowness
Hidden in the mathematics in previous sections is a simple relationship between the tangent velocity vector and cotangent normalslowness covector pair: they are conjugate to each other (Figs. 4 and 5), which is to say that their inner product is one. To prove this, consider the following property of an order2 homogeneous function like ℋ:
Combining Hamilton's equation for $\partial \mathcal{H}/\partial {p}_{i}$ (Eq. 49) with the definition of ray velocity ${v}^{i}=\mathrm{d}{r}^{i}/\mathrm{d}t$ and given the constant value of ${\mathcal{F}}_{*}=\mathrm{1}$ known from Eq. (30), this gives
which is the definition of conjugacy.
If the process of erosion were isotropic, conjugacy would obviously be true. Erosion velocity and normal slowness would be colinear, and since their magnitudes are mutually reciprocal, their product would be unity. However, the erosion process is manifestly not isotropic (see Sect. 3.18), which means that conjugacy also constrains the angular disparity between the ray and frontnormal directions.
3.16 Constancy of the Lagrangian
We can exploit conjugacy to reveal important behaviour of the fundamental function ℱ and the related Lagrangian ℒ. Since ℒ is (like ℋ) order2 homogeneous, it has the property
Using the fibre derivative form of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ in Eq. (36) and the definition of the Lagrangian in Eq. (35), we can deduce for physically valid ray trajectories that
In other words, the Lagrangian has the constant value of $\mathrm{1}/\mathrm{2}$ – just like the Hamiltonian (Eq. 34) – meaning that it is only those points with positions r and velocities v satisfying this equation that represent points on a moving erosion surface.
This shows that the geomorphic surface Lagrangian and Hamiltonian are both static given the model assumptions made here, such as constant external forcing and domain symmetry (Fig. 6): a more general theory that relaxes these restrictions would lead to nonconstancy of ℒ and ℋ.
3.17 Ray angle
An essential measure of ray direction is the angle α of the velocity vector v defined relative to horizontal as follows:
This definition, along with that for β given in Eq. (27), allows us to manipulate Hamilton's equations for the components of $\dot{\mathit{r}}$ (see Eq. 50) for which
into a relationship between the two angles (Fig. 7):
which inverts to give
where the choice of root depends on how far the point is along the ray trajectory (see below). By comparing α and β we can measure erosional anisotropy (see Sect. 3.18).
Examination of Eqs. (58) and (59) reveals an important property of the vertical motion of erosion rays and its dependence on η. Since p_{x}>0 and p_{z}<0 in the model halfspace and because v_{x}>0,
This switch in ray orientation as a function of slope scaling exponent η, which is illustrated in Fig. 12, echoes the observations in 1+1D of Weissel and Seidl (1998) and Royden and Perron (2013) of a change in upstream propagation with their gradient scaling exponent n. As their work has shown, this switch has important consequences for how and when knickpoints form (Stark and Stark, 2022).
The ray angle function (Eq. 60) has an extremum whose value is given by
This extremum represents a bound on permissible values of ray angle α. For η>1, the extremum is positive α_{ext}>0 and rays cannot point up more steeply than α<α_{ext}, while for η<1, the extremum is negative α_{ext}<0 and rays point down at negative angles limited by α>α_{ext}. The extremum is located at a critical value of β:
For $\mathit{\eta}=\frac{\mathrm{3}}{\mathrm{2}}$, the critical surface tilt is ${\mathit{\beta}}_{c}=\mathrm{50.77}{}^{\circ}$, while for $\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$ the critical tilt is ${\mathit{\beta}}_{c}=\mathrm{35.26}{}^{\circ}$ (see Fig. 7). At this critical angle the Lagrangian and the metric tensor are singular, which means that if the surface tilt reaches this angle, the link between ℋ and ℒ is broken, ℱ_{*} and ℱ are no longer metric functions, and the model space is no longer (pseudo)Finsler. What this means in practice is not yet clear; the critical angle may manifest as a transition in landscape geometric behaviour, but we can only speculate at this stage: further study is needed.
3.18 Erosional anisotropy
The difference between the erosion ray angle α and the erosion frontnormal angle β (rotated by 90^{∘} such that both angles are measured relative to horizontal) quantifies the anisotropy of the erosion process:
Defined in this way, $\mathit{\psi}={\mathrm{0}}^{\circ}$ for isotropic motion, and $\mathit{\psi}={\mathrm{90}}^{\circ}$ when anisotropy is so strong that rays and surface normal are orthogonal.
Figure 8 shows how ψ varies with surface tilt β when computed along a timeinvariant profile for $\mathit{\eta}=\frac{\mathrm{3}}{\mathrm{2}}$ and $\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$. As these plots demonstrate, the gradientdependent erosion process described by Eq. (25) is strongly anisotropic.
Figure 9 illustrates how anisotropy varies as a function of gradientscaling exponent η for a selection of ray angles α. As predicted in the previous section, the rays all point upwards (positive α) for η>1 and downwards (negative α) for η<1. Broadly speaking, anisotropy ψ reaches greater extremes for larger absolute values of $\mathit{\eta}\mathrm{1}$.
The physical relevance of anisotropy ψ is revealed by the following. The surfacenormal erosion rate can be computed from ray velocity by exploiting rayfront conjugacy (Eq. 55), which is equivalent to a dot product between ray vector and surfacenormal slowness
as well as by using the reciprocal relationship between erosion slowness and erosion speed $p=\mathrm{1}/{\mathit{\xi}}^{\u27c2}$ (Eq. 18), to get
While surface erosion takes place at a speed ξ^{⟂}, changes in external boundary conditions propagate much faster into the landscape along an erosion ray trajectory with a speed ξ^{⟂}sec ψ. The two are related by projecting the ray vector v onto the local unit surfacenormal vector, which lies at an relative angle ψ relative to the ray.
3.19 Measuring slope along the erosion front
Since the Hamiltonian tracks motion of the erosion front in a phase space spanned in part by the surfacenormal covector, solutions of front motion have the surface gradient encoded into them. Therefore, the gradient along the evolving topographic surface can be tracked in three distinct ways. One method is to take the ratio of the covector components:
A second method is to compute the topographic gradient:
In a numerical solution, this entails making a finitedifference approximation using values at nearest neighbour points. A third method is to construct a velocity triangle from the ray velocity components and the reciprocal covector slowness in the vertical direction, i.e. the vertical erosion rate:
Ideally, all three measurements of the topographic gradient should be equal. In practice, β_{ts} is computed nonlocally, whereas β_{p} and β_{vt} are strictly local but numerically different computations; we therefore expect the three estimates to be equal to within a precision set by choices such as ray density, time step, and interpolation method. A comparison of the methods is given in Fig. 10.
To keep development of a geomorphic Hamiltonian theory as simple as possible, the treatment so far (Sect. 3) has employed a somewhat abstract erosion model: it has assumed the erosion rate can be written as some combination of a power function of surface tilt and a spatially variable (but constant in time) function that encompasses flow rate, flow geometry, substrate erodibility, and so on. If we want to probe the behaviour of the geomorphic surface Hamiltonian and its implications for landscape erosion any further, we need to choose a particular form for the flow function component and to parameterize this spatial dependence. However, bear in mind that more general erosion models could also be transformed into Hamiltonian form and subjected to the analyses presented below.
4.1 A modified stream power incision model
Previous studies related to our work (Luke, 1972; Royden and Perron, 2013; Weissel and Seidl, 1998) have focused on the stream power incision model (SPIM) (e.g. Lague, 2014). In order maintain a clear conceptual link with these studies, and because SPIM can be adapted to satisfy the simplifying criteria adopted in Sect. 3.5, we use it here in a modified form. SPIM asserts that in channels,
where “slope” is the channel gradient tan β and upstream area, suitably scaled, is assumed to be a good composite proxy for the volumetric flow of water per contour width and its contributions to channel geometry, boundary flow, sediment transport, and rock surface abrasion. We modify this equation so that it instead tracks
where “slope” is now sin β. This model and classic SPIM coincide if η=1, since ${\mathit{\xi}}^{\downarrow}={\mathit{\xi}}^{\u27c2}/\mathrm{cos}\mathit{\beta}$ (Eq. 10), although they differ somewhat otherwise. Given this similarity, we can treat the slope and area exponents η⇔n and μ⇔m as roughly equivalent.
Our model domain is a 2D transect along a channel, which means we have to parameterize out catchment geometry and drainage accumulation into a function of distance downstream. If we consider upstream area to scale with an offset distance from the divide ${L}_{\mathrm{c}}(x\mathit{\epsilon})$, where ε is a very small regularization term, we can wrap this scaling into a power function form for the flow component of the erosion model:
In the numerical solutions presented in Sect. 6, the regularization term ε is given a nonzero value, but in the equations below it is ignored.
The surfacenormal channel erosion rate is then
In a similar manner to steady (constant erosion rate) solutions of SPIM (e.g. Lague, 2014), this model will generate channel profiles with the asymptotic slope–area scaling
assuming lowtomoderate slope angles where tan β≈sin β. To ensure that our numerical simulations all yield slope–area scaling consistent with that typically observed (e.g. Beeson and McCoy, 2020; Flint, 1974; Lague, 2014; Royden and Perron, 2013), we fix the exponent ratio (“concavity index”) at a constant $\mathit{\mu}/\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$.
4.2 Nondimensionalization
Before embarking on numerical solutions of the model, we nondimensionalize it. This is helpful in the following two ways: (i) it requires us to identify the characteristic length, time, erosion rate and slowness scales, which makes it easier to relate the model to realworld landscapes, and (ii) it makes generalization of model behaviour and solution geometries simpler.
An obvious length scale is the horizontal channel length L_{c}, i.e. the distance from the drainage divide x=L_{c} to the channel terminus x=0. The horizontal and vertical erosion rates at the terminus are
where ${\mathit{\xi}}^{{\to}_{\mathrm{0}}}/{\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}=\mathrm{tan}{\mathit{\beta}}_{\mathrm{0}}$ and the channel tilt angle at the terminus is
We choose ${\mathit{\xi}}^{{\to}_{\mathrm{0}}}$ as the characteristic velocity scale. The horizontal timescale is therefore
The vertical timescale is given by ${t}_{{\downarrow}_{\mathrm{0}}}={t}^{{\to}_{\mathrm{0}}}\mathrm{cot}{\mathit{\beta}}_{\mathrm{0}}$.
Now we can nondimensionalize the primary model variables as follows:
and the coordinate axes
Using them to rewrite the Hamiltonian we get
where we have defined the dimensionless number
We can think of 𝖢𝗂 as both an angle and a dimensionless erosion rate because when we nondimensionalize the vertical rate of erosion imposed at the boundary ${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$, we get the following equation:
Note that we can write
We can now rewrite Hamilton's equations in dimensionless form by rederiving them from Eq. (80). Alternatively, we can just substitute the nondimensionalized variables into Eqs. (50) and (51):
and so we get
and
Figure 11 provides a comparison of timeinvariant stream profiles for a selection of values of the dimensionless horizontal erosion rate $\mathsf{Ci}\in \mathit{\{}{\mathrm{0.1}}^{\circ},{\mathrm{1}}^{\circ},{\mathrm{4}}^{\circ}\mathit{\}}$. In all other figures illustrating numerical solutions, the value of this dimensionless number is set at $\mathsf{Ci}={\mathrm{4}}^{\circ}$.
4.3 Direct integration
For the simple scenario of a timeinvariant profile, the erosion equation (Eq. 73) can be directly integrated; more complex boundary and initial conditions do not allow it. The first step is to assume the vertical rate of erosion is constant everywhere ${\mathit{\xi}}^{\downarrow}={\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$, and thereby to manipulate Eq. (73) to expose its straightforward dependence on surface tilt β and position x (through φ(x)):
We can combine this equation with Eq. (68) to obtain a polynomial in surface gradient $\mathrm{tan}\mathit{\beta}=\mathrm{d}z/\mathrm{d}x$ and constrain it using the result (Eq. 53) that ${p}_{z}=\mathrm{1}/{\mathit{\xi}}^{\downarrow}=\mathrm{1}/{\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$ along the whole ray and thus everywhere along a timeinvariant profile. The resulting polynomial in surface gradient, in nondimensionalized form and for rational values of the gradient exponent such as $\mathit{\eta}=\frac{\mathrm{3}}{\mathrm{2}}$ or $\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$, is
We can use this function to compute the surface elevation as a 1+1D function $\widehat{z}(\widehat{x};\phantom{\rule{0.125em}{0ex}}\mathit{\eta},\mathit{\mu},\mathsf{Ci})$ as follows: (i) pick values of η, μ, and 𝖢𝗂; (ii) substitute these numbers into the above function to generate a polynomial in $\mathrm{d}\widehat{z}/\mathrm{d}\widehat{x}$ and $\widehat{x}$; (iii) define a set of sample positions $\mathrm{0}\le \mathit{\left\{}\widehat{x}\mathit{\right\}}<\mathrm{1}$ along the profile; (iv) at each $\widehat{x}$, find the positive, real root of this polynomial to infer the gradient $\mathrm{d}\widehat{z}/\mathrm{d}\widehat{x}$ at this position; and (v) use quadrature to integrate the gradient values along the profile to get $\widehat{z}\left(\widehat{x}\right)$.
Figure 11 shows a selection of nondimensionalized timeinvariant profiles obtained in this way. Notice how the profiles for the two different gradient exponents $\mathit{\eta}=\frac{\mathrm{3}}{\mathrm{2}}$ and $\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$ are essentially colinear for $\mathrm{0}\le \widehat{x}=x/{L}_{\mathrm{c}}<\mathrm{0.95}$. The practical upshot of this similarity is that it is unreasonable to expect to infer the scaling exponents η and μ from topography alone.
Direct integrations like this are also useful as a validation of the raytraced solutions: this is illustrated in Fig. 13, in which some examples of directly integrated timeinvariant profiles are shown to match those obtained by ray tracing.
The previous sections have shown how the geometric selfconstraint implicit in a broad class of erosion models can be transformed into a geomorphic surface Hamiltonian ℋ (Eqs. 32 and 80) and how this function can be used to derive Hamilton's equations of motion for points on an erosion surface (Eqs. 50, 51, 86, and 87). In this section we solve Hamilton's equations by numerical integration and use them to construct “steadystate”, timeinvariant surface profiles driven by constanterosionrate boundary conditions. In all solutions presented below, the dimensionless horizontal erosion rate is set at $\mathsf{Ci}={\mathrm{4}}^{\circ}$.
5.1 Model domain and boundary conditions
The domain is a vertical x–z transect (Fig. 2) along a stream profile that ranges from a drainage divide at x=L_{c} to a flow–exit boundary at x=0. Profile evolution is driven by a constant vertical erosion rate imposed at the exit, and evolution of the profile is tracked relative to the elevation of the exit. The drainage divide is pinned at a fixed horizontal position by mirroring (Fig. 6) the main profile with a symmetrical “image” profile spanning ${L}_{\mathrm{c}}\le x\le \mathrm{2}{L}_{\mathrm{c}}$; solution need only be performed over $\mathrm{0}\le x\le {L}_{\mathrm{c}}$. Although there is no need to invoke tectonic processes here, note that this model is geometrically equivalent to erosion of a (half) horst block whose uniform rock uplift is driven by constantrate vertical slip along a bounding normal fault, and whose topographic evolution is studied in the reference frame of the hanging wall.
5.2 Ray equations
In this model geometry, rays that initiate at x=0 (Fig. 12) and propagate in the positive x direction are annihilated at x=L_{c} when a paired ray, initiated at the same time at x=2L_{c}, arrives from the opposite direction. As such, the model induces a cusp to form at x=L_{c}, although its formation is not explicitly modelled here – instead, rays from x=0 are simply truncated at x=L_{c}.
Such ray tracing entails the numerical integration of Hamilton's equations in the form of four coupled, firstorder ODEs for ${\dot{r}}^{x}$ and ${\dot{r}}^{z}$ (Eq. 50) and ${\dot{p}}_{x}$ and ${\dot{p}}_{z}$ (Eq. 51). These are firstorder differential equations in time alone, so for each ray we need only supply four initial conditions, i.e. ${r}^{{x}_{\mathrm{0}}},{r}^{{z}_{\mathrm{0}}},{p}_{{x}_{\mathrm{0}}}$, and ${p}_{{z}_{\mathrm{0}}}$, one for each ray ODE. An oddity of ray tracing is that what would be boundary conditions in a partial differential equation (PDE) treatment become initial conditions for the rays, and what would be a separate Neumann velocity boundary condition for a PDE gets wrapped into those initial conditions.
Here we focus on obtaining the timeinvariant profile generated by a constant vertical velocity boundary condition ${\mathit{\xi}}^{\downarrow}={\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$ at x=0, for which we only need to perform ray tracing from x=0. We thus avoid having to generate rays along an initial topography and having to handle their transient interaction as the timeinvariant profile develops (a topic to be addressed in Stark and Stark, 2022).
The initial horizontal position for all rays is fixed at the stream terminus and location of the boundary condition ${r}^{{x}_{\mathrm{0}}}=x=\mathrm{0}$ (Fig. 2). The initial vertical position of a ray initiated at time t=t_{0} is given by simple integration of the vertical erosion rate: ${r}^{{z}_{\mathrm{0}}}={\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}{t}_{\mathrm{0}}$. The initial vertical component of the ray slowness covector must be consistent with this vertical velocity component, and thus we have ${p}_{{z}_{\mathrm{0}}}=\mathrm{1}/{\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$. Since ${\dot{p}}_{z}=\mathrm{0}$, this vertical covector component remains unchanged throughout ray propagation (see Sect. 3.14), and thus the number of coupled ODEs that need to be solved is effectively reduced from four to three.
The initial horizontal component of the slowness covector can be calculated if we realize that the topographic gradient at the boundary must be consistent with the orientation of the normal slowness, i.e. $\mathrm{tan}{\mathit{\beta}}_{\mathrm{0}}={p}_{{x}_{\mathrm{0}}}/{p}_{{z}_{\mathrm{0}}}$. As such, the initial value of the slowness covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ encodes the velocity boundary condition in both its direction and magnitude.
5.3 Numerical integration method
After some experimentation, the most accurate quadrature or numerical integration scheme for ray tracing with Eqs. (50) and (51) was found to be an implicit Runge–Kutta method designed for stiff ODEs: specifically, an implementation of the Radau IIA family of order 5 (see Hairer and Wanner, 2013, p. 72) provided by the Python package SciPy (Virtanen et al., 2020). Simpler and lowerorder Runge–Kutta quadrature methods also work well for most choices of model parameters, as does the highorder Runge–Kutta, dense output DOP853 method (see Hairer et al., 2008, p. 194).
All the numerical solutions presented here are reproducible using the following opensource software (split into two parts, both of which are needed for full operation): (i) the GME
package, which implements methods of geometric mechanics tailored to treating geomorphic erosion (v. 1.0, Stark, 2022a, c); and (2) a utilities library called GMPLib
(v. 1.0, Stark, 2022b, d).
5.4 Reference ray construction
Computation of the trajectory of a point on an erosion surface (and its normal slowness covector) is carried out by numerically integrating the coupled set of Hamilton's equations (dimensioned: Eqs. 50 and 51; nondimensionalized: Eqs. 86 and 87) with the boundary conditions described in Sect. 5.1. This constitutes the tracing of a single reference ray (Fig. 12), which suffices for construction of a timeinvariant topographic profile (see below). More rays need to be traced if we want to handle timevariable boundary conditions, evolution from an initial topography, or the transition between an initial surface and a slip boundary (Stark and Stark, 2022).
5.5 Synthesis of a timeinvariant profile
The following steps are required to construct a timeinvariant solution of the erosion equation akin to a faultdriven steadystate solution (Figs. 6 and 13):

choose values for the model parameters (notably the gradientscaling exponent η and upstream areascaling exponent μ);

specify the dimensionless vertical erosion rate at the boundary 𝖢𝗂;

generate a reference ray r_{ref}(t) by integrating Eqs. (50) and (51) (or their nondimensionalized equivalents Eqs. 86 and 87) from the boundary at $(\mathrm{0},{r}_{{z}_{\mathrm{0}}})$ and assign it an initiation time of t_{0}=0;

define the isochrone time T such that ${r}_{\text{ref}}^{x}\left(T\right)={L}_{\mathrm{c}}$;

generate a kth later ray ${\mathit{r}}_{k\mathrm{\Delta}{t}_{\mathrm{0}}}(t+k\mathrm{\Delta}{t}_{\mathrm{0}})$ with initiation time kΔt_{0} by making a copy of the reference ray, displacing it vertically by ${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}k\mathrm{\Delta}{t}_{\mathrm{0}}$, and pasting it at $(\mathrm{0},{r}_{{z}_{\mathrm{0}}}{\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}k\mathrm{\Delta}{t}_{\mathrm{0}})$;

truncate the copied ray at the point ${\mathit{r}}_{k\mathrm{\Delta}{t}_{\mathrm{0}}}(Tk\mathrm{\Delta}t)$;

repeat from step 4 until kΔt_{0}≥T;

collate the truncation points to generate a continuous curve T(r).
Some of these steps also entail interpolation and resampling.
This procedure generates the timeinvariant isochrone T(r) formed by the constant vertical velocity ${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$ boundary condition at x=0 (Figs. 6 and 13). Repetition of the procedure (or a simple copying of the solution), combined with a progressive offset of the initial ray location ${r}_{{z}_{\mathrm{0}}}$ at the boundary, simulates vertical normalfaultdriven erosion of a topographic profile at steady state in the reference frame of the (bedrock) substrate of the footwall bedrock (Fig. 14). Analysis of these composite results generates solutions for the alongprofile variations in the component erosion rates (Figs. 15c–e and 16c–e) and their anisotropy (Figs. 15a, 16a, and 17).
In all the solutions presented here, the areascaling exponent μ is chosen such that $\mathit{\mu}/\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$. The dimensionless rate of boundary erosion (Eq. 81) is fixed at $\mathsf{Ci}={\mathrm{4}}^{\circ}$ in all but Fig. 11.
In this section we present numerical solutions of timeinvariant topographic profiles in dimensionless form. These solutions help to validate the geomorphic surface Hamiltonian (Sects. 1.1–3); to test the inferences drawn from it (Sects. 3 and 4); to examine its nondimensionalization (Sect. 4.2) and the timescales, length scales, and velocity scales predicted by it (Sect. 6.1); to check how ray tracing by integrating Hamilton's equations performs as a means of modelling surface erosion and the propagation of boundarychange information (Sect. 5.4; Figs. 12–14), and to explore how erosional anisotropy ψ varies across a landscape.
Although the solutions here are limited to a 2D x–z transect, they provide a pilot test of elements needed to construct a fully 3D landscape evolution model around a geomorphic Hamiltonian: one in which (i) the denudation rate is defined as acting in the surfacenormal direction, rather than purely vertically, and (ii) topographic elevation is tracked as true geometric surface using an implicit “timeslice” function $T(x,y,z)$, instead of being modelled as a field using an explicit height function $h(x,y;t)$.
6.1 Scales
Tables 1 and 2 provide some example values of model parameters and their corresponding time, rate and vertical scales. For each example, the key choice is the dimensionless horizontal erosion rate 𝖢𝗂. This dimensionless number determines the dimensionless traversal time ${\widehat{t}}_{{L}_{\mathrm{c}}}^{\to}$, which is defined as the time it takes for a ray to travel from x=0 to x=0.95L_{c}, which is obtained by numerical ray tracing. Following this, by choosing the domain length L_{c} and the boundary rate of vertical erosion, dimensioned quantities can be computed. The parameters grad=tan β_{0}, ${\mathit{\xi}}^{{\to}_{\mathrm{0}}}$, and ${t}^{{\to}_{\mathrm{0}}}$ are derived exactly; the horizontal travel time ${t}_{{L}_{\mathrm{c}}}^{\to}$ and the profile height ${h}_{{L}_{\mathrm{c}}}$ close to the divide (at x=0.95L_{c}) are obtained by numerical solution. The values shown here are all rounded to one or two significant figures for clarity.
These tables demonstrate that boosting the imposed vertical erosion rate ${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$ linearly increases the consequent horizontal erosion rate ${\mathit{\xi}}^{{\to}_{\mathrm{0}}}$ and symmetrically decreases ${t}^{{\to}_{\mathrm{0}}}$ and ${t}_{{L}_{\mathrm{c}}}^{\to}$ but has no effect on the profile height ${h}_{{L}_{\mathrm{c}}}$. The most important result here is that by calculating the dimensionless traversal time ${\widehat{t}}_{{L}_{\mathrm{c}}}^{\to}$ we can estimate how long it takes for boundary change information to propagate into a landscape.
6.2 Timeinvariant solutions
Figures 13 and 14 illustrate raytraced timeinvariant solutions for two choices of the slope exponent $\mathit{\eta}\in \mathit{\{}\frac{\mathrm{3}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}}\mathit{\}}$ in the model equation (Eq. 73) for surfacenormal erosion rate ξ^{⟂}. Each raytraced isochrone T(r) is compared with an isochrone obtained by direct integration (Sect. 4.3), and in each case the match is excellent. Sequences of erosion surfaces resulting from similar timeinvariant solutions are shown in Fig. 14.
These solutions illustrate an important behaviour of the rays: for values of the slope exponent η>1 the ray velocities always have a positive vertical component ${\dot{r}}^{z}>\mathrm{0}$, whereas for η<1, the vertical component ${\dot{r}}^{z}$ always has a negative vertical component ${\dot{r}}^{z}<\mathrm{0}$.
6.3 Erosion rates
Figures 15 (for $\mathit{\eta}=\frac{\mathrm{3}}{\mathrm{2}}$) and 16 (for $\mathit{\eta}=\frac{\mathrm{1}}{\mathrm{2}}$) provide a sidebyside comparison of surface erosion rate components (ξ^{⟂}, ξ^{→}, ξ^{↓}) along raytraced timeinvariant profiles, together with some of the variables that contribute to their variation (anisotropy ψ and ray velocity components v^{x}, v^{z}). All plotted quantities are dimensionless.
As Figs. 15a–c and 16a–c show, the progressive upstream decrease in anisotropy ψ is reflected in upstream decreases in ray velocity (particularly the vertical component v^{z}) and the surfacenormal erosion rate ξ^{⟂}. The horizontal rate of erosion ξ^{→} decreases upstream in an apparently linear fashion, correlating with a similar behaviour in the horizontal component of the ray velocity v^{x}. The vertical rate of erosion ξ^{↓} is constant (to within the precision of the numerical solution), as expected for timeinvariant (“steadystate”) profiles.
Surfacenormal erosion rate is computed in two ways from the raytracing results (Figs. 15c and 16c). One way is to simply use the fact (Eq. 18) that normal speed is the reciprocal of normal slowness ${\mathit{\xi}}^{\u27c2}=\mathrm{1}/p$. The other is to project the ray velocity onto the surface normal (unit) vector using Eq. (66). The horizontal ξ^{→} and vertical ξ^{↓} erosion rate components are computed with Eqs. (9) and (10) using either of the estimates of ξ^{⟂}. Since ray tracing involves discrete sampling, values of ξ^{⟂} computed in these two ways are not numerically identical. The discrete sampling also entails having to generate interpolating functions so that erosion rate values can be calculated at arbitrary positions along the profile.
6.4 Anisotropy
Figure 17 provides a striking visualization of erosional anisotropy ψ(x) by plotting its variation with x along timeinvariant topographic profiles. The direction and magnitude of normal slowness covectors are represented with “fishbone” symbols (where the number of crosstick “bones” approximates slowness p), while arrows represent the ray velocity vectors. The colour attribute of each symbol visualizes the magnitude of the angular disparity ψ. The degree of anisotropy is evident in the strong angular disparity of r and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ for the same choices of $\mathit{\eta}\in \mathit{\{}\frac{\mathrm{3}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}}\mathit{\}}$. The strongest anisotropy is found downstream in the channels, where the channel tilt β is small, the normal covector points almost vertically downwards, and the ray velocity vector points almost horizontally upstream. Anisotropy decreases monotonically upstream as the normal covector rotates towards horizontal more rapidly than the ray vector angle. At the divide, the model erosion process is approximately isotropic; this limiting behaviour is moot, however, because the erosion model used here (Eq. 25) does not apply to steep channels.
7.1 Geometry controls (almost) everything
The main aspiration of this paper is to clarify what we mean when, in the context of landscape evolution, we speak of the direction of erosion. Our central mathematical tenet has been that while gradientdriven surface erosion takes place in the surfacenormal direction (Sect. 1.1; Eq. 3), points on successive erosion surfaces do not necessarily map in the same direction. Working from this premise and with the help of geometric mechanics, we have found unexpected complexity hidden in simple erosion models.
The concept of a covector is pivotal to our theory (Sect. 3.1). Once we realize that the surfacenormal erosion rate (imposed by the gradientdependent erosion model) can be written in terms of the normalslowness covector (which is the consequent motion of the surface), it takes only a few short steps to reach the geomorphic surface Hamiltonian. Hamilton's raytracing equations, the geomorphic surface Lagrangian, and the adherence to Huygens' and Fermat's principles all logically follow.
The essential ingredient of the theory is the realization that, at its core, the process of erosion is a geometric selfconstraint. If we disregard complexities such as sediment cover factors and external variations in forcing, a generic model of erosion is a statement about how a surface geometry (through its gradient and flow accumulation) determines the rate of change of that surface geometry. Reparameterizing this statement generates a fundamental function (and thus a Hamiltonian) that describes how to measure distance in the phase space of the erosion equation. The properties of this function reveal that landscape erosion is best described using Finsler geometry. This is important because it provides a fundamental explanation for why geomorphic erosion is anisotropic. As Fig. 17 demonstrates, this anisotropy is very strong.
Counterintuitively, the erosion rays point (obliquely) upwards if the scaling behaviour of slope in the erosion model has an exponent η>1: we might have expected points on an erosion surface to always move downwards since erosion is, after all, driving the surface downwards; this is indeed the case if η<1. Remember that the topographic profiles obtained by numerical solution here are timeinvariant solutions without an uplift term, as visualized as time slices in Fig. 14; upward motion of rays is therefore driven only by erosion and is not influenced by any tectonic motion.
The idea that surface erosion simultaneously drives two distinct motions – subvertically in the surfacenormal direction and subhorizontally in the ray direction – is an uncomfortable and apparently very abstract notion, but it has physical consequences. It means erosion drives information about boundary conditions upstream subhorizontally while also driving motion of the whole profile downwards subvertically. The timescale on which boundarycondition information propagates into the interior is the time it takes for a point to travel along a ray (the ray velocity is sometimes known as the signal velocity, which conveys the sense of information propagation well). This information may be the rate of erosion at the stream terminus, the equivalent slip rate on a boundary fault, or the erosion rate at a point on an initial profile.
For these timeinvariant solutions, all rays are identical and the boundary condition does not change. If the boundary erosion rate (or equivalent fault slip rate) were to change (e.g. Reinhardt et al., 2007), we would anticipate the ray paths to change, and we might expect them to intersect (depending on the value of η): this is one way that knickpoints form. Exploration of this topic is left for another paper (Stark and Stark, 2022). The only ray intersections presented here are those that implicitly take place at the drainage divide, as rays are imagined to approach symmetrically from a righthand halfdomain (Fig. 6). The crucial difference is that intersection at divides occurs when rays approach from opposite directions; knickpoints form where rays move in the same direction at different speeds, one overtaking the other.
Previous studies have considered knickpoint formation as the propagation and intersection of raylike characteristics (Luke, 1972; Royden and Perron, 2013; Weissel and Seidl, 1998) but always in terms of an explicit surface function and a onedimensional Hamilton–Jacobi equation describing elevation as a function of distance upstream and time. The parameter space traversed by these characteristics has no concept of the surfacenormal or of erosion slowness covectors, which prevents a direct comparison of the results of these studies with those presented here. However, they are broadly in agreement.
Perhaps the oddest outcome of the Hamiltonian theory, but one that is not surprising in retrospect, is that the vertical component of the erosion slowness covector is constant (Sect. 3.14; Eqs. 52–53). To be precise: as a surface point initiates at the boundary and moves along a ray into the interior, its vertical component of surface slowness is invariant, ${\dot{p}}_{z}=\mathrm{0}$, and thus the vertical component of the surface erosion rate is constant, ${\mathit{\xi}}^{\downarrow}\left(t\right)=\mathrm{1}/{p}_{z}\left(t\right)=\mathrm{1}/{p}_{{z}_{\mathrm{0}}}={\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$. For a timeinvariant (steadystate) profile, all rays are identical, meaning that there is only one ray solution; therefore, the surface at every point along the ray must be moving vertically at the same speed; all rays are independent; therefore, all rays must maintain constancy of the vertical component of the surface erosion rate $\mathrm{1}/{p}_{z}$. In this sense, a ray carries information of the boundary condition – the vertical slip rate – into the landscape until it is destroyed at a cusp.
The timescale of this information transfer is of crucial importance. If it is small relative to the timescale on which drainage divides move laterally and significantly changes accumulation areas and flows, then the assumption made in Sect. 3.5 is valid: namely, that the flow – at every point on the surface where flow influences erosion – can be parameterized by its surface geometry, i.e. its upstream area, in a manner constant with time (for the lifetime of a ray). This requirement can be weakened to allow for slow variation of the parameterization with time, in which case the Hamiltonian field would need to be recalculated periodically. This is not to say that the geomorphic surface Hamiltonian theory is invalidated if the timescale requirement is not met; rather, the theory would become nonlocal and more complicated. The degree to which such a step is necessary is a topic for future research.
On a side note, bear in mind that the following are all different ways of saying the same thing: (i) the directional pace (reciprocal rate) of erosiondriven surface motion, (ii) the surfacenormal slowness covector, (iii) the gradient of the erosionfront arrivaltime function, (iv) the directional density of erosion surface isochrones, and (5) the gradient of the geomorphic Hamilton action.
7.2 A geomorphic surface Hamiltonian in 3D
While the geomorphic surface Hamiltonian developed here is limited to erosiondriven motion of a linear front in 2D space, the goal is to construct a theory for surface evolution in 3D space. Several conceptual and computational hurdles will need to be overcome if this goal is to be met.
The main challenge for a 3D theory will be to find a way to treat channel formation that is consistent with the Hamiltonian methodology. It is tempting to want to resolve the channel shape itself, but this would entail having to add hydrodynamics, sediment transport, and abrasion processes to the mix; such a change would not only make the theory inordinately complex, it would run counter to our core premise that what matters is the geometric selfconstraint imposed by geomorphic processes, not the details of those processes. A parameterization of flow focusing in channels will be required: one that encapsulates channel crosssectional geometry without describing it explicitly.
Another challenge hinges on the assumption of locality. A firstcut 3D model can probably be framed with fixed catchment perimeters and static drainage divides; however, there will be a pressing need to generalize and allow for divide motion so that catchment shapes can selfform. The question of timescales raised in Sect. 7.1 will still apply in 3D: if the timescale of divide motion and catchment area change is large relative to the time it takes for erosion rays to traverse the catchment, we will probably be able to treat the flow component of the geomorphic surface Hamiltonian as approximately static, which will make it possible to derive Hamilton's equations for 3D ray tracing.
Numerical solution may require a change of approach, because ray tracing of a surface in 3D is much more cumbersome than for a line in 2D, particularly when dealing with ray intersections and cusp formation. An obvious alternative approach lies in the fact that the theory employs an implicit surface function to describe landscape geometry (Sect. 2): we can resolve erosion front motion on a regular grid and use a levelset method to solve the geomorphic HJE (Adalsteinsson and Sethian, 1995a; Mosaliganti et al., 2013; Sethian and Adalsteinsson, 1997). This will have to be done with some care, however, because of the Finsler nature of the geomorphic surface Hamiltonian and its inherent anisotropy. Onepass fast marching will not be possible, because even the most advanced algorithms for fast marching (Mirebeau, 2014b, 2019; Mirebeau and Portegies, 2019) are currently limited to metrics whose anisotropy is Riemannian (velocityindependent anisotropy) or Randers (velocitydependent anisotropy of a different type to that of the geomorphic Hamiltonian) (Mirebeau, 2014a). A further issue will be the nonconvexity of the geomorphic Hamiltonian for certain ranges of η and β (Appendix C): nonconvex Hamiltonians were addressed in the early literature on the levelset method (e.g. Adalsteinsson and Sethian, 1995b, c, 1997) and have been encountered in applications to nongeomorphic erosion (e.g. Radjenović et al., 2006a, b, 2010; Radjenović and RadmilovićRadjenović, 2009); recent methodological developments (Chow et al., 2018, 2019; Evans, 2014; Pinezich, 2019) may help. Methods developed in the field of seismology may also prove useful (e.g. Moser, 1991; Qian et al., 2003; Rawlinson et al., 2008; Wang et al., 2006).
7.3 The variational principle is not energy minimization
There is a substantial body of work founded on the idea that landscapes selforganize in order to minimize energy dissipation across their flow networks. Most of the literature developing this idea – broadly known as optimal channel network (OCN) theory – dates from the 1990s (IjjaszVasquez et al., 1993; Rigon et al., 1993; Rinaldo et al., 1992, 1998; RodriguezIturbe et al., 1992a, b) and is comprehensively reviewed in the book by RodriguezIturbe and Rinaldo (2001). In this section we compare and contrast OCN theory with our theory of the geometric mechanics of erosion.
OCN theory is framed in terms of the selfoptimization of a cost function. It identifies this cost function as the total rate of dissipation of mechanical potential energy released by water flowing down channels across the whole landscape. Initial development of the theory focused on the planform geometry and topology of channel networks; it was only later work that addressed the consequent formation of topography. Hillslopes were assumed to play no role. The strong geometric similarity between OCNs and natural stream networks, in particular their similar scaling behaviour, is often presented as a vindication of the theory.
Optimality is a commonly used concept in engineering; its cousin in physics is the notion that system behaviour arises through a variational principle that guarantees minimization of a key quantity. The most fundamental difference is that in engineering the optimization criterion is invoked as a design choice, whereas the variational principle arises as an expression of the underlying physics. It is from this difference that the following criticisms of OCN theory spring.
OCN theory arbitrarily requires minimization of the energy dissipated across the whole channel network; this stipulation is justified on the basis that many physical systems exhibit similar behaviour. Such a requirement implies the existence of a variational principle (Sinclair and Ball, 1996) guiding landscape evolution towards this optimal state, but OCN theory does not articulate this principle in words or mathematics. A corollary issue was the initial omission of a Hamiltonian, which was remedied to some extent in Rinaldo et al. (1998). The weakness of their Hamiltonian is that it cannot be used to derive equations for the timeevolution of the landscape: it constrains what shape the landscape must take, but it cannot explain how that shape comes about.
By framing an alternative theory in terms of geometric mechanics, these issues are avoided. The guiding variational principle is clearly articulated as follows (Sect. 3.11): topographic evolution obeys the principle of least erosion time. Adherence to this principle is not imposed; it arises geometrically from the way that geomorphic erosion propagates a topographic front and modifies the pattern of erosion rates. The correlative Hamiltonian (Sect. 3.8) generates equations that describe landscape evolution both in the form of Hamilton's equations and in their equivalent form of a Hamilton–Jacobi equation. Solution of these equations, for appropriate boundary conditions, evolves the topography to a timeinvariant shape, but this shape is the outcome of geometric interaction rather than a mechanism of energydissipation minimization.
Comparison of the two theories is a little premature because our theory needs further development if we want it to describe the evolution of a whole channel network. The current model also pins the drainage divide at a fixed position: as a result, the degrees of freedom present in a landscape evolving in 3D, notably those that permit different flow topologies and geometries, are absent from our model. It is these degrees of freedom that lead to the existence of many possible states of energy dissipation, i.e. many possible drainage network configurations; in our 2D theory, only one timeinvariant state (a simple linear profile), imposed by the model erosion (Eq. 25), is possible. Nevertheless, it will be interesting to see if a full 3D theory can throw light on what drives landscape selforganization and channel network formation: whether these phenomena arise primarily from the geometric selfconstraint imposed by geomorphic erosion and, if so, the extent to which the process of energy minimization is complementary.
When we say that the rate of erosion of a geomorphic surface is a function both of its tilt and of the fluxes passing over it, we are in essence saying that the rate of change of landscape geometry is a function of that geometry. Here we have shown how to express this geometric statement as a Hamiltonian and how to use this Hamiltonian to understand the meaning of the phrase “the direction of landscape erosion”.
Our foundational premise is that motion of an erosion surface intrinsically acts in the surfacenormal direction. On this basis, we can convert a gradientdependent erosion rate (i.e. speed) model into a model for the normal slowness (i.e. pace) of surface erosion, parameterized by surface tilt and upstream area and expressed as a covector. Using a simple mathematical trick (a scaling substitution), and by writing tilt in terms of slowness covector components, the model equation can be rearranged into what is called the fundamental function of its metric space; the square of this function is the geomorphic surface Hamiltonian.
This Hamiltonian is parameterized in terms of (i) the position of a single point on the surface, and (ii) the corresponding orientation and slowness of the surface at that point. The Hamiltonian thus occupies a sixdimensional phase space (which reduces to four if the model domain is restricted to a 2D slice). Although such extra dimensionality may seem to be just a mathematical abstraction, it provides real insight.
Study of the Hamiltonian and its phase space reveals that surface evolution simultaneously involves two distinct directions of motion: the surface (at a given point) moves in the surfacenormal direction, while the point itself moves in what may be an entirely different direction. The disparity between these two directions is a measure of the anisotropy of the process governing motion, and for the class of erosion models studied here such anisotropy is very strong.
This phenomenon is best explored using Hamilton's raytracing equations – derived from the Hamiltonian by simple differentiation – which express the motion of a surface point and its allied surfacenormal slowness in terms of ordinary differential equations (ODEs). They show that while changes in the surface erosion rate and direction are encoded in the normal slowness ODEs, information about boundary conditions and external changes is carried upstream by the ray ODEs.
There is an important dependence of ray tracing on surface tilt: if the model erosion rate scales faster than linearly with gradient (η>1), such rays always have a positive vertical component, i.e. they point upstream and obliquely upwards. However, if the model erosion rate scales sublinearly with gradient (η<1), erosion rays always have a negative vertical component, i.e. they point upstream but obliquely downwards along their trajectories.
We have shown how the phase space occupied by the geomorphic surface Hamiltonian is a metric space and how this leads us to deduce that the erosion rays traced by surface points are geodesics. In other words, they follow paths of the locally shortest erosion time: this is the variational principle that guides geomorphic surface erosion. It appears that energy dissipation need not be invoked and that instead all that matters is geometry.
A1 Geoscience applications of the HJE
The HJE has seen only sporadic use in the geosciences – except in the field of seismology, where its static or eikonal form has been found to be particularly useful. The eikonal equation is a good approximation for seismic wave propagation in the socalled “highfrequency limit” at which seismic wavelengths are very small compared to the scale of wave propagation (e.g. Červený, 1989, 2005, 2002; Dellinger, 1997; Mensch and Farra, 1999; Rawlinson et al., 2008; Slawinski, 2014; Virieux and Lambaré, 2007; Woodhouse and Deuss, 2007). From this approximation arises the convenient fiction of seismic rays, which are both the characteristics of the HJE and solutions of Hamilton's equations. Although there are disadvantages to its use in treating seismic wave propagation, e.g. dynamic interactions are not modelled and spectral information is lost, the Hamiltonian approach has proven insightful, particularly when dealing with anisotropic media (Antonelli et al., 2003a, b; Bóna and Slawinski, 2002, 2003; Bucataru and Slawinski, 2005; Červený, 2002; Klimeš, 2002; Yajima and Nagahama, 2009; Yajima et al., 2011).
An analogous form of the seismic Hamiltonian approach has been applied to studying the effects of anisotropy on fluid flow in porous media (Sieniutycz, 2000, 2007; Yajima and Nagahama, 2015).
A2 Applications of the HJE to geomorphology
In geomorphology, Luke (1972, 1974, 1976) pioneered application of the HJE to the modelling of fluvial knickpoints as shocks formed by kinematic waves (Lighthill and Whitham, 1955a, b; Whitham, 1999). Weissel and Seidl (1998) and Royden and Perron (2013) built on this approach to further understand the conditions under which knickpoints form and how they propagate. In all these studies, the HJE was deployed in an explicitsurface form, and its ability to model implicitsurface motion was not considered.
A3 Use of the eikonal equation in geomorphology
To our knowledge, only one previous study has attempted to model landscape evolution as an implicit surface moving according to an eikonal equation. Aronsson and Lindé (1982) did so in a treatment of weatheringlimited denudation of a rock cliff incised at its base by a river. By integrating the eikonal equation representing this erosion process, and by presenting levelset solutions as isochrones of the cliff transect, they demonstrated how variations in rock erodibility can lead to highly irregular surface geometry such as overhangs.
A4 Nongeomorphic erosion modelled with the HJE
There is a literature on erosion driven by nongeomorphic processes, and much of it is unfamiliar to the geomorphology community. The methods employed in some of these papers provide a partial foundation for our Hamiltonianbased approach. For example, both implicit surface motion and the HJE have been the basis for modelling erosion at microscopic scales in an engineering context.
Frank (1958) employed the concept of surfacemotion slowness as a means to model the anisotropic dissolution of crystal surfaces in 2D (although neither the HJE nor the concept of a covector were explicitly invoked). He later extended this approach to handling dissolution in 3D (Frank and Ives, 1960). His technique is widely cited in the crystallography literature (e.g. Frank and Ives, 1960; Ives, 1961; Osher and Merriman, 1997; Shemenski et al., 1965).
In materials science, the Frank (1958) method has been adapted to treat surface erosion at the micrometre scale driven by ion beam bombardment. Early work (Barber et al., 1973; Carter et al., 1971; Nobes et al., 1969) focused on amorphous substrates and isotropic erosion without mentioning the HJE. Subsequent advances introduced the HJE (Carter et al., 1984; Katardjiev, 1989; Katardjiev et al., 1989; Nobes et al., 1987; Smith et al., 1986; Witcomb, 1975) and the eikonal equation (Carter, 2001) and used them to address the issue of anisotropic erosion. Perhaps most relevant to our theoretical development is the review article by Smith et al. (1986), which is also notable for its invocation of an erosional Hamiltonian, and the papers by Carter et al. (1984), Katardjiev (1989) and Katardjiev et al. (1989), which connect the HJE and its Hamiltonian to Huygens' principle and the concept of erosional wavelets (see also Adalsteinsson and Sethian, 1995b, c, 1997; Sethian and Adalsteinsson, 1997).
A5 Front motion obeys Huygens' principle
Central to the ideas in the previous sections is Huygens' principle, one of the founding contributions to the field of optics. Using a graphical construction, the principle explains how a wavefront bends as it passes through media of varying resistance to motion (e.g. Arnold, 1989; Holm, 2011; Miller, 1991). At every instant, it pictures the front peppered with tiny wavelets. Each wavelet represents how far, if it were spreading in isolation, a point on the front would expand in the next instant to form its own microfront. Since the points are not isolated, they interfere to form a mutually tangential envelope, with each point moving to the location of its wavelet tangent. The set of successive of tangential envelopes constitutes the progressive motion of the front.
In wavepropagation terms, the wavelet represents the unit envelope of group velocity at the point of interest: its shape is called an indicatrix. There is a corresponding structure for phase velocity, known as the figuratrix, which is typically used in its reciprocal speed or slowness form. The velocity indicatrix and slowness figuratrix are linked through mutual conjugacy: as such, they contain the same information about front propagation but in different forms (Carathéodory, 1999; Perlick, 2000; Rider, 1926; Rund, 1959).
In other words, wavefront propagation can be tracked using either group information or phase information. For front propagation in general this equivalence translates into tracking using either (i) point velocities and their trajectories (ray paths) or (ii) pointwise frontnormal slownesses and their ensemble motions.
Huygens' principle is best known for explaining wave propagation in inhomogeneous but isotropic media, where the indicatrices and figuratrices are spherical but vary in size from place to place; isotropy ensures that the group and phasepropagation directions are the same. The principle is also often used to explain propagation in media whose anisotropy is symmetric but ellipsoidal (Arnold, 1989), where the group and phasepropagation directions are different. Recent efforts have further proved that the principle extends to asymmetric, nonellipsoidal indicatrices and figuratrices representing a generalized form of anisotropy (e.g. Dehkordi and Saa, 2019; Innami, 1995; Javaloyes et al., 2021; Markvorsen, 2016; Palmer, 2015) expressed in terms of something called Finsler geometry (see Appendices C and D).
A6 Wildfire spread and Finsler geometry
Several of the ideas discussed in previous sections have seen application in a totally different field – that of wildfire prediction – in the envelope model of fire spread. The earliest form of this 2D model was very simple (Van Wagner, 1969), postulating that winddriven fire growth can be approximated as a burn ellipse elongated and offset in the wind direction. Anderson et al. (1982) extended the model and deployed Huygens' principle to propagate a wildfire using elementary burn ellipses scattered along the fire front, each scaled and shaped according to the local fuel availability and wind direction.
These early efforts were purely graphical constructions (Sullivan, 2009). Subsequently, Richards (1990, 1995) formalized the firefrontpropagation process as a form of the HJE (without explicitly mentioning the equation by name). The model has subsequently evolved, and its most sophisticated version (Markvorsen, 2016) recognizes the elementary burn shapes as nonelliptical velocity indicatrices and frames the anisotropic motion in terms of Finsler structures. Finsler geometry is useful because it provides a convenient mathematical context in which to express the time it takes for a fire front to cover a given distance under the directional influence of wind and terrain. It is for similar reasons that Finsler geometry is important for understanding the anisotropy of geomorphic erosion, as Appendix D shows.
A7 Ray tracing the motion of a front
The rays of seismology and geometric optics are paths of least time, and they can be traced in two distinct ways: (i) by integrating Hamilton's equations, which are derived from the Hamiltonian contained in the HJE, or (ii) by transforming the Hamiltonian into (or writing directly) the corresponding Lagrangian, converting into the EulerLagrange equations, and integrating them (see Appendix E). In both cases, the essential step is to write a Hamiltonian version of the process governing motion. For simplicity, the derivation presented in this paper is limited to a 2D vertical slice of a landscape. A fully 3D treatment is the subject of ongoing research.
There is a connection between the Hamiltonian raytracing method developed here and the work of Luke (1972), Royden and Perron (2013), and Weissel and Seidl (1998). These previous approaches deployed the method of characteristics to solve a 1+1D form of HJE in which a 2D topographic profile is represented in an explicit fashion, and their results have some resemblance to those we obtain by full ray tracing (see Sect. 3.12). The main difference is the explicitly 1+1D form of the governing equation in these studies, which forces elevation to be a singlevalued function, and which coerces ray tracing into resolving horizontal motion only. If one were to write the Hamiltonian phase space covector coordinate (the direction and reciprocal speed of the surface at a point on the front) for these problems, it would take the reduced form of the slope patch variable of Royden and Perron (2013); this variable contains explicit information about horizontal motion of a surface patch (through its position), but vertical motion is implicit (see Royden and Perron, 2013, Eq. 15). As a result, the inherent anisotropy of the erosion process is hidden.
Slowness covectors and velocity vectors are different mathematical objects, and they live on different spaces, where “space” is meant in the abstract sense used in differential geometry. For each point r in the physical, Euclidean world we can create an allied tangent space that contains all the possible tangent velocity vectors (like ξ) at that point; we can also envisage a corresponding cotangent space to contain all the possible slowness covectors (like $\stackrel{\mathrm{\u0303}}{\mathit{p}}$) at that point. Bundled together, the tangent spaces for all points in real space constitute a tangent bundle or velocity space, while the union of cotangent spaces forms a cotangent bundle or slowness (classically called “momentum”) phase space. These two spaces are indispensable tools of geometric mechanics.
One way to see that vectors and covectors are different is to look at their tensor form: vectors are rank (1,0) contravariant tensors, whereas covectors are rank (0,1) covariant tensors (which is where the “co” prefix comes from). Tensors of different rank cannot be combined arithmetically; instead, operations such as contraction are needed to combine them. For example, in Eq. (13), the action of covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ on the unit vector n is a tensor contraction:
The expression p_{i}n^{i} here employs the Einstein summation convention: when an index (such as i) is shared by several terms, summation is automatically performed for those terms over all index elements (in this case, over $i\in \mathit{\{}x,z\mathit{\}}$). Upper indexes are used for contravariant tensor components; lower indexes are used for covariant tensor components.
The fundamental function ℱ_{*} has three key properties that are valid for a domain D of $\mathit{\{}{r}^{x},{r}^{z},{p}_{x},{p}_{x}\mathit{\}}$ phase space corresponding to physically reasonable values of surface tilt and erosion rate:

positive, order1 Euler homogeneity in the parameter $\stackrel{\mathrm{\u0303}}{\mathit{p}}$, i.e. if the covector $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ in ℱ_{*} is scaled by a positive scalar λ>0, the reparameterized function equals the original function scaled by λ,
$$\begin{array}{}\text{(C1)}& {\mathcal{F}}_{*}(\mathit{r},\mathit{\lambda}\stackrel{\mathrm{\u0303}}{\mathit{p}})=\mathit{\lambda}{\mathcal{F}}_{*}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})\phantom{\rule{2em}{0ex}}\text{for}\mathit{\lambda}\mathrm{0},\end{array}$$where the 1 in “order1” refers to the exponent in λ on the righthand side of this equation;

regularity, i.e. ℱ_{*} is smooth, in that it can be differentiated infinitely many times without encountering a discontinuity or undefined value;

the Hessian of ${\mathcal{F}}_{*}^{\mathrm{2}}$, i.e. the Hessian of the Hamiltonian ℋ, is
$$\begin{array}{}\text{(C2)}& {g}_{*}^{ij}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}:=\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial {\mathcal{F}}_{*}^{\mathrm{2}}}{\partial {p}_{x}\partial {p}_{z}}},\end{array}$$where for η>1 and ${p}_{x}/{p}_{z}=\mathrm{tan}{\mathit{\beta}}_{c}\ne \sqrt{\mathit{\eta}}$, both eigenvalues of ${g}_{*}^{ij}$ are real and positive, making ${g}_{*}^{ij}$ positive definite and ℱ_{*} strongly convex.
Given these properties, ℱ_{*} constitutes a type of Finsler metric (Bao et al., 2000; Shimada and Sabau, 2000). This means that ℱ_{*} provides a means of measuring distance and travel time between points in slowness phase space that is dependent on both position and direction of motion (Sect. A7 and Appendix D; see Bao, 2007). In other words, the shortest time path between two points in the corresponding real space may not be a straight line.
Strictly speaking, ℱ_{*} is a coFinsler metric on the cotangent space, and thus we are dealing with a coFinsler or Cartan geometric space (e.g. Miron et al., 2002; Yajima et al., 2011). The term “Finsler” is reserved for the counterpart tangent space and for the fundamental function ℱ, the dual of ℱ_{*}. Nevertheless, for brevity we use the term “Finsler” to apply to both spaces and metrics.
We also need to be cautious in generalizing about Finsler properties of ℱ_{*} and the geomorphic surface Hamiltonian. For η<1 and β<β_{c}, the Hessian of ℋ is not regular, g is indefinite with mixed signature (the eigenvalues are both positive and negative), ℱ_{*} is not convex (Beem, 1971; Červený, 2002; Giaquinta and Hildebrandt, 2004), and the Hamiltonian ℋ is nonconvex. For η>1 but β>β_{c}, the Hamiltonian is similarly nonconvex. Under these conditions, it is more appropriate to use the term pseudoFinsler for the metric and its phase space (see Asanov, 1985, pp. 21, 44, 266)
Having a Finsler, or at least pseudoFinsler, geometry is important for several reasons. The most immediate is the need to adopt a quadratic form of ℱ_{*} as a Hamiltonian (Sect. 3.8) because ℱ_{*} cannot be Legendre transformed directly (e.g. Červený, 2002; Giaquinta and Hildebrandt, 2004, p. 16). It also means that if we wish to solve erosion front motion as an HJE, we need to find an alternative to the fastmarching method because this algorithm is limited to Riemannian anisotropic metrics (Mirebeau, 2014b, 2019; Mirebeau and Portegies, 2019) and to a small subset of Finsler metrics whose velocitydependent, Randerstype anisotropy (Mirebeau, 2014a) differs from that of the geomorphic Hamiltonian.
In geomorphology, we are used to dealing with equations that operate in a flat Euclidean geometry where the space is spanned by Cartesian $\mathit{\{}x,y,z\mathit{\}}$ coordinates. In such a space, distances are measured directly using Pythagoras' theorem and the topology of curves across it is straightforward. Working in a flat space like this is fine for studies at the catchment scale and is a good approximation even at the orogen scale.
There are scales, however, where use of such a flat space is inadequate. For example, what if we are interested in processes on a global scale and need to account for the spherical geometry of Earth's surface? Switching to spheroidal coordinates is only half the battle because transport on a sphere is topologically different to that on a flat space: particles in locally straight motion follow looping paths, these paths are great circles, sets of great circles always converge and intersect, and so on. In this example, we need to understand the consequences of working in a curved space and its consequences if we want to understand physical phenomena acting at such scales.
The concept of curved spaces is relevant not just to processes on objects with topological curvature; in an abstract way, it can also apply to the space in which the governing equations operate. For some types of process, the governing equations can be mapped from Euclidean space into a nonflat phase space that both simplifies their solution and exposes their fundamental properties and behaviour.
The geomorphic surface Hamiltonian ℋ, which arises from the transformation of an erosion equation, operates in such a nonflat space. Distance and travel time are not Euclidean measures on this phase space because the fundamental metric function ℱ_{*} that defines ℋ has the properties described in Sect. 3.7.
The curved nature of nonEuclidean spaces lies in how distance is measured on them. The measurement of distance always requires a yardstick of some kind (on a metric space, this is a tensor derived from ℱ or ℱ_{*}, and an associated inner product), but the yardstick used in Finsler geometry is not the equivalent of a simple ruler. It is not an isotropic constant as it would be in a flat Euclidean space (where the metric tensor is a simple Kronecker delta), nor is it an anisotropic, possibly spatially variable, but otherwise static quantity as it would be in a curved Riemannian space (with a metric tensor whose variable elements are a function of position r alone). Instead, the yardstick is a function both of position and of the direction and magnitude of motion at that position, i.e. for ℱ_{*}, the metric tensor elements vary with both r and $\stackrel{\mathrm{\u0303}}{\mathit{p}}$, not just with r. Instead of measuring distance with a single inner product at each point, there is a family of inner products associated with each point (e.g. Shen, 2001).
This directional dependence of the “yardstick” or metric tensor is the defining characteristic of Finsler geometry (Bao et al., 2000; Chern, 1996; Holm, 2011; Shimada and Sabau, 2000). To be precise, ℱ_{*} specifies that the slowness phase space is a coFinsler or Cartan space, and its dual ℱ specifies that the velocity space is a Finsler space. The practical consequence is that the time taken to travel an infinitesimal distance across the space (Bao, 2007) at unit speed in a given direction is a function of that travel direction. Adding up such incremental times allows us to find the shortest path across the space, but the directional dependence of erosion time measurement makes this calculation nontrivial.
The fact that the geomorphic surface Hamiltonian operates in a Finsler geometry has profound consequences for the construction of erosiondriven equations of motion, the variational principle that underlies how landscape shape evolves, and the concept of erosional anisotropy. These consequences are explored in the next sections.
For further information on Finsler geometry, a good introduction is the nontechnical discussion in Gibbons and Warnick (2011), which also touches on several other topics important to this paper. A more mathematical but surprisingly approachable introduction can be found in Bao (2007), while more comprehensive treatments are provided by Antonelli et al. (1993); Antonelli (2000); Bao et al. (2000); Chern (1996); Giaquinta and Hildebrandt (2004); Mo (2006); Miron et al. (2002); Shen (2001, 2013); Shimada and Sabau (2000).
The Lagrangian method of ray tracing the motion of an erosion front is taken up in detail in Stark et al. (2022). This alternate approach is important because it demonstrates in practical terms how erosion rays are also geodesics, i.e. that they are solutions of the geodesic equation (Misner et al., 1973; Nolte, 2019) corresponding to the geomorphic surface Hamiltonian. Although geodesics have cropped up before in geomorphology as a means of delineating drainage on digital elevation models (Passalacqua et al., 2010a, b), their use in that context was a pragmatic means to an end, rather than a reflection of any underlying physics. In the our theory, however, the geodesic equation has physical meaning in that it is synonymous with the Euler–Lagrange equation of the geomorphic HJE; solutions to the geodesic equation follow the same paths of least time as solutions to Hamilton's raytracing equations derived from the geomorphic surface Hamiltonian. This assertion is proved in Stark et al. (2022).
Ray tracing through integration of Hamilton's equations is not the only way to solve for surface motion. In principle, we could instead use the geomorphic surface Hamiltonian $\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ in its HJE form and solve erosion front propagation using gridbased methods. In practice, numerical solution of this kind of eikonal equation is not straightforward (see Sect. 7.2 and Appendix C). The HJE is nevertheless instructive if we examine it in the context of some important concepts of classical mechanics. For example, Hamilton's principal function S(r,t), which is the Hamilton action S_{γ} (see Eq. 38) plus a constant, is
Use of the Legendre transform (Eq. 37) yields
The total derivative of S(r,t) with respect to time t has
Assuming the points {r} all lie on a path γ_{0} of least erosion time, we can write
Comparing this equation with Eq. (F2) leads to
such that the Hamiltonian $\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ can be written as
which is the standard form for the HJE. Now consider the arrival time function T(r), whose total time derivative is, given Eqs. (22), (55), and $\partial T/\partial t=\mathrm{0}$, as follows:
Integration here gives the abbreviated action; by choosing to integrate along a path of least action γ_{0} we obtain the shortest erosion time T(r):
Use of Eqs. (34) and (F2) connects S(r,t) with T(r), $\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$ and time t:
Differentiation gives
Substitution into the standard HJE in Eq. (F6) leads to
In this form, the HJE prescribes how the erosion front T(r) has a locus (a set of points {r}) that propagates such that the gradient ∇T (the directional density of T isochrones) satisfies the static Hamiltonian $\mathcal{H}=\frac{\mathrm{1}}{\mathrm{2}}$.
t  time 
L_{c}  distance from stream terminus to the drainage divide 
x  horizontal coordinate $\mathrm{0}\le x\le {L}_{\mathrm{c}}$ measured from stream terminus 
y  outofsection horizontal coordinate 
z  vertical coordinate: distance above the terminus 
{a},{b}  sets of points defining successive erosion front surfaces 
T_{a}, T_{b}  corresponding loci of erosion surfaces at successive times 
T(r)  isochrone of erosion surface at point r⇔ surface locus at T=t 
r  point vector on erosion front surface, i.e. point on erosion ray 
$\mathit{v}=\dot{\mathit{r}}$  tangent velocity vector of point moving along erosion ray 
$\stackrel{\mathrm{\u0303}}{\mathit{p}}$  covector of normal slowness of erosion front 
r^{x},r^{z}  horizontal, vertical components of ray point vector r 
${r}^{{x}_{\mathrm{0}}},{r}^{{z}_{\mathrm{0}}}$  boundary values of components of r 
v^{x},v^{z}  horizontal, vertical components of tangent ray velocity vector v 
p_{x},p_{z}  horizontal, vertical components of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ 
$p=\left\stackrel{\mathrm{\u0303}}{\mathit{p}}\right$  surface normal slowness, i.e. reciprocal erosion rate 
${p}_{{x}_{\mathrm{0}}},{p}_{{z}_{\mathrm{0}}}$  boundary values of components of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ 
α  ray angle, i.e. angle of v from horizontal 
α_{ext}  limit ray angle (maximum for η<1, minimum for η>1) 
β  angle of $\stackrel{\mathrm{\u0303}}{\mathit{p}}$ from vertical, also surface slope angle from horizontal 
β_{c}  critical surface slope angle 
β_{0}  boundary value of surface slope angle 
ℱ(r,v)  1homogeneous Finsler fundamental function 
${\mathcal{F}}_{*}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$  1homogeneous coFinsler (Cartan) fundamental function 
$\mathcal{H}(\mathit{r},\stackrel{\mathrm{\u0303}}{\mathit{p}})$  2homogeneous Hamiltonian 
ψ  erosional anisotropy $=\mathit{\alpha}\mathit{\beta}+\mathrm{90}{}^{\circ}$ 
h(x)  elevation as a 1+1D function of horizontal distance upstream 
$\mathit{\varphi}(x,y,z;t)$  levelset function 
ξ  erosion velocity vector; generic velocity function in levelset equation 
ξ^{x},ξ^{z}  horizontal, vertical components of erosion velocity vector ξ 
ξ^{⟂}  surfacenormal erosion rate (speed) 
ξ^{→}  horizontal (positive right) erosion rate 
ξ^{↓}  vertical (positive down) erosion rate 
${\mathit{\xi}}^{{\downarrow}_{\mathrm{0}}}$  boundary value of vertical (positive down) erosion rate 
n  surfacenormal unit vector 
n^{x},n^{z}  horizontal, vertical components of surfacenormal unit vector n 
η  gradientscaling exponent in surfacenormal erosion model 
μ  upstream areascaling exponent in surfacenormal erosion model 
λ  a real scalar 
ℒ(r,v)  2homogeneous Lagrangian 
$i,j\in \mathit{\{}x,z\mathit{\}}$  covariant (lower) or contravariant (upper) indices 
γ(t)  potential erosion ray path (parameterized by time t 
γ_{0}(t)  erosion ray path of least time 
S_{γ}  action functional for erosion ray paths 
${S}_{{\mathit{\gamma}}_{\mathrm{0}}}$  least action ⇔ (half) least erosion time 
n  gradientscaling exponent in vertical erosion model (SPIM) 
m  upstream areascaling exponent in vertical erosion model (SPIM) 
φ(x)  spatial component of rate of erosion at distance x upstream 
φ_{0}  base rate in flow component of erosion model 
ε  relative flow component rate at zero upstream area 
${t}^{{\to}_{\mathrm{0}}}$  horizontal erosion timescale for (computed from boundary rates) 
$\widehat{t},\widehat{x},\widehat{z}$  nondimensionalized coordinates 
$\widehat{\mathit{r}},{\widehat{r}}^{x},{\widehat{r}}^{z}$  nondimensionalized position variables 
$\widehat{\mathit{p}},{\widehat{p}}_{x},{\widehat{p}}_{z}$  nondimensionalized slowness variables 
𝖢𝗂  dimensionless boundary erosion rate 
k  ray index 
${t}_{{L}_{\mathrm{c}}}^{\to}$  timescale for erosion to traverse the domain 
${h}_{{L}_{\mathrm{c}}}$  height scale of timeinvariant topographic profile 
${g}_{*}^{ij}$  cometric tensor 
S  Hamilton's principal function 
The base utilities package is available at the GMPLib repository on GitHub (Stark, 2022b, https://geomorphysics.github.io/GMPLib) and archived at https://doi.org/10.5281/zenodo.6373127 (Stark, 2022d). The allied code and notebooks to solve the equations presented in this paper, both algebraically and numerically, are available at the GME repository on GitHub (Stark, 2022a, https://geomorphysics.github.io/GME) and archived at https://doi.org/10.5281/zenodo.6373103 (Stark, 2022c). GMPLib release version 1.0 and GME release version 1.0 were used to generate the results presented in this paper.
CPS led the research, wrote the code, and cowrote the paper. GJS collaborated on the research and in writing the paper.
The contact author has declared that neither they nor their coauthor has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Discussions with Petrus dos Anjos, Liam Reinhardt, and Mike Ellis are greatly appreciated.
This paper was edited by Eric Lajeunesse and reviewed by David J. Furbish and one anonymous referee.
Adalsteinsson, D. and Sethian, J. A.: A fast level set method for propagating interfaces, J. Comput. Phys., 118, 269–277, https://doi.org/10.1006/jcph.1995.1098, 1995a. a
Adalsteinsson, D. and Sethian, J. A.: A Level Set Approach to a Unified Model for Etching, Deposition, and Lithography I: Algorithms and TwoDimensional Simulations, J. Comput. Phys., 120, 128–144, https://doi.org/10.1006/jcph.1995.1153, 1995b. a, b
Adalsteinsson, D. and Sethian, J. A.: A Level Set Approach to a Unified Model for Etching, Deposition, and Lithography II: ThreeDimensional Simulations, J. Comput. Phys., 122, 348–366, https://doi.org/10.1006/jcph.1995.1221, 1995c. a, b
Adalsteinsson, D. and Sethian, J. A.: A Level Set Approach to a Unified Model for Etching, Deposition, and Lithography III: Redeposition, Reemission, Surface Diffusion, and Complex Simulations, J. Comput. Phys., 138, 193–223, https://doi.org/10.1006/jcph.1997.5817, 1997. a, b
Anderson, D. H., Catchpole, E. A., De Mestre, N. J., and Parkes, T.: Modelling the spread of grass fires, ANZIAM J., 23, 451–466, https://doi.org/10.1017/S0334270000000394, 1982. a
Antonelli, P. L.: Generalizations of Finsler geometry, Springer Science+Business Media, Dordrecht, https://doi.org/10.1007/9789401142359, 2000. a
Antonelli, P. L., Ingarden, R. S., and Matsumoto, M.: The theory of sprays and Finsler spaces with applications in physics and biology, SpringerScience+Business Media, B. Y., Dordrecht, https://doi.org/10.1007/9789401581943, 1993. a, b
Antonelli, P. L., Bóna, A., and Slawinski, M. A.: Seismic rays as Finsler geodesics, Nonlinear Anal.Real, 4, 711–722, https://doi.org/10.1016/S14681218(02)000731, 2003a. a
Antonelli, P. L., Rutz, S. F., and Slawinski, M. A.: A geometrical foundation for seismic ray theory based on modern Finsler geometry, in: Finsler and Lagrange Geometries, edited by: Anastasiei, M. and Antonelli, P. L., Springer Science+Business Media, 17–54, https://doi.org/10.1007/9789401704052, 2003b. a
Arnold, V. I.: Mathematical Methods of Classical Mechanics, 2nd edn., SpringerVerlag, https://doi.org/10.1007/9781475720631, 1989. a, b, c
Aronsson, G. and Lindé, K.: Grand Canyon – A quantitative approach to the erosion and weathering of a stratified bedrock, Earth Surf. Proc. Land., 7, 589–599, https://doi.org/10.1002/esp.3290070607, 1982. a
Asanov, G. S.: Finsler Geometry, Relativity and Gauge Theories, Springer Science+Business Media, https://doi.org/10.1007/9789400953291, 1985. a
Bao, D.: On two curvaturedriven problems in RiemannianFinsler geometry, in: Advanced Studies in Pure Mathematics, Finsler Geometry, Sapporo, 19–71, https://doi.org/10.2969/aspm/04810019, 2007. a, b, c
Bao, D., Chern, S. S., and Shen, Z.: An Introduction to RiemannFinsler Geometry, vol. 200 of Graduate Texts in Mathematics, Springer Science+Business Media, New York, NY, https://doi.org/10.1007/9781461212683, 2000. a, b, c, d
Barber, D. J., Frank, F. C., Moss, M., Steeds, J. W., and Tsong, I. S. T.: Prediction of ionbombarded surface topographies using Frank's kinematic theory of crystal dissolution, J. Mater. Sci., 8, 1030–1040, https://doi.org/10.1007/BF00756635, 1973. a
Beem, J. K.: Motions in two dimensional indefinite Finsler spaces, Indiana U. Math. J., 21, 551–555, https://www.jstor.org/stable/24890331 (last access: 21 April 2022), 1971. a
Beeson, H. W. and McCoy, S. W.: Geomorphic signatures of the transient fluvial response to tilting, Earth Surf. Dynam., 8, 123–159, https://doi.org/10.5194/esurf81232020, 2020. a
Bóna, A. and Slawinski, M. A.: Raypaths as parametric curves in anisotropic, nonuniform media: differentialgeometry approach, Nonlinear Anal.Theor., 51, 983–994, https://doi.org/10.1016/s0362546x(01)008732, 2002. a
Bóna, A. and Slawinski, M. A.: Fermat's principle for seismic rays in elastic media, J. Appl. Geophys., 54, 445–451, https://doi.org/10.1016/j.jappgeo.2003.08.019, 2003. a
Bucataru, I. and Slawinski, M. A.: Generalized orthogonality between rays and wavefronts in anisotropic inhomogeneous media, Nonlinear Anal.Real, 6, 111–121, https://doi.org/10.1016/j.nonrwa.2004.03.004, 2005. a
Carathéodory, C.: Calculus of variations and partial differential equations of the first order, AMS Chelsea Publishing, https://bookstore.ams.org/chel318/ (last access: 21 April 2022), 1999. a
Carter, G.: The physics and applications of ion beam erosion, J. Phys. D Appl. Phys., 34, R1–R22, https://doi.org/10.1088/00223727/34/3/201, 2001. a
Carter, G., Colligon, J. S., and Nobes, M. J.: The equilibrium topography of sputtered amorphous solids II, J. Mater. Sci., 6, 115–117, https://doi.org/10.1007/BF00550340, 1971. a
Carter, G., Nobes, M. J., and Cruz, S. A.: Surface morphology evolution of sputtered, moving substrates, J. Mater. Sci. Lett., 3, 523–527, https://doi.org/10.1007/BF00720990, 1984. a, b
Červený, V.: Ray tracing in factorized anisotropic inhomogeneous media, Geophys. J. Int., 99, 91–100, https://doi.org/10.1111/j.1365246x.1989.tb02017.x, 1989. a
Červený, V.: Fermat's variational principle for anisotropic inhomogeneous media, Stud. Geophys. Geod., 46, 567–588, https://doi.org/10.1023/A:1019599204028, 2002. a, b, c, d, e
Červený, V.: Seismic Ray Theory, Cambridge University Press, https://doi.org/10.1017/CBO9780511529399, 2005. a
Chern, S. S.: Finsler geometry is just Riemannian geometry without the quadratic equation, Not. Am. Math. Soc., 43, 959–963, https://doi.org/10.1090/conm/196/02429, 1996. a, b
Chow, Y. T., Darbon, J., Osher, S., and Yin, W.: Algorithm for overcoming the curse of dimensionality for certain nonconvex Hamilton–Jacobi equations, projections and differential games, Annals of Mathematical Sciences and Applications, 3, 369–403, https://doi.org/10.4310/AMSA.2018.v3.n2.a1, 2018. a
Chow, Y. T., Darbon, J., Osher, S., and Yin, W.: Algorithm for overcoming the curse of dimensionality for statedependent HamiltonJacobi equations, J. Comput. Phys., 387, 376–409, https://doi.org/10.1016/j.jcp.2019.01.051, 2019. a
Coulthard, T. J.: Landscape evolution models: a software review, Hydrol. Process., 15, 165–173, https://doi.org/10.1002/hyp.426, 2001. a
Crandall, M. G. and Lions, P. L.: Viscosity solutions of HamiltonJacobi equations., T. Am. Math. Soc., 277, 1–42, https://doi.org/10.2307/1999343, 1981. a
Dehkordi, H. R. and Saa, A.: Huygens' envelope principle in Finsler spaces and analogue gravity, Classical Quant. Grav., 36, 085008, https://doi.org/10.1088/13616382/ab0f03, 2019. a
Dellinger, J.: Anisotropic finite difference traveltimes using a HamiltonJacobi solver , in: SEG Technical Program Expanded Abstracts, Society of Exploration Geophysicists, 1786–1789, https://doi.org/10.1190/1.1885780, 1997. a
Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Stock, J. D., Heimsath, A. M., and Roering, J. J.: Geomorphic transport laws for predicting landscape form and dynamics, in: Prediction in Geomorphology, American Geophysical Union, Washington, D. C., 103–132, https://doi.org/10.1029/135GM09, 2003. a, b
Evans, L. C.: Envelopes and nonconvex Hamilton–Jacobi equations, Calc. Var. Partial Dif., 50, 257–282, https://doi.org/10.1007/s0052601306353, 2014. a
Flint, J. J.: Stream gradient as a function of order, magnitude, and discharge, Water Resour. Res., 10, 969–973, https://doi.org/10.1029/WR010i005p00969, 1974. a
Fowler, A.: Mathematical Geoscience, vol. 36 of Interdisciplinary Applied Mathematics, Springer London, London, https://doi.org/10.1007/9780857297211, 2011. a
Frank, F. C.: On the kinematic theory of crystal growth and dissolution processes, in: Growth and perfection of crystals: Proceedings of an International Conference on Crystal Growth held at Cooperstown, New York, 27–29 August 1958, edited by: Doremus, R. H., Roberts, B. W., and Turnbull, D., John Wiley, New York, 411–419, 1958. a, b
Frank, F. C. and Ives, M. B.: Orientationdependent dissolution of germanium, J. Appl. Phys., 31, 1996–1999, https://doi.org/10.1063/1.1735485, 1960. a, b
Giaquinta, M. and Hildebrandt, S.: Calculus of Variations II, vol. 311 of Grundlehren der mathematischen Wissenschaften, Springer Science+Business Media, Berlin, Heidelberg, https://doi.org/10.1007/9783662062012, 2004. a, b, c
Gibbons, G. W. and Warnick, C. M.: The geometry of sound rays in a wind, Contemp. Phys., 52, 197–209, https://doi.org/10.1080/00107514.2011.563515, 2011. a
Gibou, F., Fedkiw, R., and Osher, S.: A review of levelset methods and some recent applications, J. Comput. Phys., 353, 82–109, https://doi.org/10.1016/j.jcp.2017.10.006, 2018. a
Giga, Y.: Surface evolution equations: A level set approach, vol. 99 of Monographs in Mathematics, Birkhäuser Verlag, Basel, https://doi.org/10.1007/3764373911_2, 2006. a
Goldstein, H., Poole, C., and Safko, J.: Classical Mechanics, AddisonWesley, 3rd edn., https://doi.org/10.1119/1.1484149, 2000. a
Hairer, E. and Wanner, G.: Solving Ordinary Differential Equations II, vol. 14 of Stiff and Differential – Algebraic Problems, Springer Science+Business Media, Berlin, Heidelberg, https://doi.org/10.1007/9783662099476, 2013. a
Hairer, E., Nørsett, S. P., and Wanner, G.: Solving Ordinary Differential Equations I, Nonstiff Problems, Springer Science+Business Media, https://doi.org/10.1007/9783540788621, 2008. a
Harris, C. R., Millman, K. J., Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., Kerkwijk, M. H., Brett, M., Haldane, A., Río, J. F., Wiebe, M., Peterson, P., GérardMarchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array programming with NumPy, Nature, 585, 357–362, https://doi.org/10.1038/s4158602026492, 2020.
Holm, D. D.: Geometric Mechanics. Part I: Dynamics and Symmetry, 2nd edn., Imperial College Press, London, https://doi.org/10.1142/p801, 2011. a, b, c, d, e
Houchmandzadeh, B.: The Hamilton–Jacobi equation: An alternative approach, Am. J. Phys., 88, 1–8, https://doi.org/10.1119/10.0000781, 2020. a, b
IjjaszVasquez, E. J., Bras, R. L., RodríguezIturbe, I., Rigon, R., and Rinaldo, A.: Are river basins optimal channel networks?, Adv. Water Resour., 16, 69–79, https://doi.org/10.1016/03091708(93)90030J, 1993. a
Innami, N.: Generalized metrics for second order equations satisfying Huygens' principle, Nihonkai Mathematical Journal, 6, 5–23, https://projecteuclid.org/euclid.nihmj/1273780052, 1995. a
Ives, M. B.: Orientationdependent dissolution of lithium fluoride, J. Appl. Phys., 32, 1534–1535, https://doi.org/10.1063/1.1728391, 1961. a
Javaloyes, M. A., PendásRecondo, E., and Sánchez, M.: Applications of cone structures to the anisotropic rheonomic Huygens' principle, Nonlinear Anal.Theor., 209, 112337, https://doi.org/10.1016/j.na.2021.112337, 2021. a
Katardjiev, I. V.: A kinematic model of surface evolution during growth and erosion: Numerical analysis, J. Vac. Sci. Technol. A, 7, 3222–3232, https://doi.org/10.1116/1.576340, 1989. a, b
Katardjiev, I. V., Carter, G., and Nobes, M. J.: The application of the Huygens principle to surface evolution in inhomogeneous, anisotropic and timedependent systems, J. Phys. D Appl. Phys., 22, 1813–1824, https://doi.org/10.1088/00223727/22/12/003, 1989. a, b
Klimeš, L.: Relation of the wavepropagation metric tensor to the curvatures of the slowness and rayvelocity surfaces, Stud. Geophys. Geod., 46, 589–597, https://doi.org/10.1023/A:1019551320867, 2002. a
Kluyver, T., RaganKelley, B., Pérez, F., Granger, B., Bussonnier, M., Frederic, J., Kelley, K., Hamrick, J., Grout, J., Corlay, S., Ivanov, P., Avila, D., Abdalla, S., and Willing, C.: Jupyter Notebooks – a publishing format for reproducible computational workflows, in: Positioning and Power in Academic Publishing: Players, Agents and Agendas, edited by: Loizides, F. and Schmidt, B., IOS Press Ebooks, 87–90, https://doi.org/10.3233/978161499649187, 2016.
Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61, https://doi.org/10.1002/esp.3462, 2014. a, b, c
Lighthill, M. J. and Whitham, G. B.: On kinematic waves I. Flood movement in long rivers, P. Roy. Soc. Lond. A Mat., 229, 281–316, https://doi.org/10.1098/rspa.1955.0088, 1955a. a
Lighthill, M. J. and Whitham, G. B.: On kinematic waves II. A theory of traffic flow on long crowded roads, P. Roy. Soc. Lond. A Mat., 229, 317–345, https://doi.org/10.1098/rspa.1955.0089, 1955b. a
Luke, J. C.: Mathematical models for landscape evolution, J. Geophys. Res., 77, 2460–2464, https://doi.org/10.1029/JB077i014p02460, 1972. a, b, c, d
Luke, J. C.: Special solutions for nonlinear erosion problems, J. Geophys. Res., 79, 4035–4040, https://doi.org/10.1029/JB079i026p04035, 1974. a
Luke, J. C.: A note on the use of characteristics in slope evolution models, Z. Geomorphol. Supp., 25, 114–119, 1976. a
Markvorsen, S.: A Finsler geodesic spray paradigm for wildfire spread modelling, Nonlinear Anal.Real, 28, 208–228, https://doi.org/10.1016/j.nonrwa.2015.09.011, 2016. a, b
Mensch, T. and Farra, V.: Computation of qPwave rays, traveltimes and slowness vectors in orthorhombic media, Geophys. J. Int., 138, 244–256, https://doi.org/10.1046/j.1365246x.1999.00870.x, 1999. a
Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S., Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F., Pedregosa, F., Curry, M. J., Terrel, A. R., Roučka, Š., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A.: SymPy: symbolic computing in Python, PeerJ Computer Science, 3, e103–27, https://doi.org/10.7717/peerjcs.103, 2017.
Miller, D. A. B.: Huygens's wave propagation principle corrected, Opt. Lett., 16, 1370–1372, https://doi.org/10.1364/OL.16.001370, 1991. a
Mirebeau, J.M.: Anisotropic fastmarching on Cartesian grids using lattice basis reduction, SIAM J. Numer. Anal., 52, 1573–1599, https://doi.org/10.1137/120861667, 2014a. a, b
Mirebeau, J.M.: Efficient fast marching with Finsler metrics, Numer. Math., 126, 515–557, https://doi.org/10.1007/s0021101305713, 2014b. a, b
Mirebeau, J.M.: Riemannian fastmarching on Cartesian grids using Voronoi's first reduction of quadratic forms, SIAM J. Numer. Anal., 57, 2608–2655, https://doi.org/10.1137/17M1127466, 2019. a, b
Mirebeau, J.M. and Portegies, J.: Hamiltonian Fast Marching: A numerical solver for anisotropic and nonholonomic eikonal PDEs, Image Processing On Line, 9, 47–93, https://doi.org/10.5201/ipol.2019.227, 2019. a, b
Miron, R., Hrimiuc, D., Shimada, H., and Sabau, S. V.: The Geometry of Hamilton and Lagrange Spaces, Springer Science+Business Media, Dordrecht, https://doi.org/10.1007/0306471353, 2002. a, b
Misner, C. W., Thorne, K. S., and Wheeler, J. A.: Gravitation, W. H. Freeman, San Francisco, https://press.princeton.edu/books/hardcover/9780691177793/gravitation (last access: 21 April 2022), 1973. a
Mo, X.: An Introduction to Finsler Geometry, vol. 1 of Peking University Series in Mathematics, World Scientific, https://doi.org/10.1142/6095, 2006. a
Mosaliganti, K. R., Gelas, A., and Megason, S.: An efficient, scalable, and adaptable framework for solving generic systems of levelset PDEs, Front. Neuroinform., 7, 1–14, https://doi.org/10.3389/fninf.2013.00035, 2013. a
Moser, T. J.: Shortest path calculation of seismic rays, Geophysics, 56, 59–67, https://doi.org/10.1190/1.1442958, 1991. a
Nobes, M. J., Colligon, J. S., and Carter, G.: The equilibrium topography of sputtered amorphous solids, J. Mater. Sci., 4, 730–733, https://doi.org/10.1007/BF00742430, 1969. a
Nobes, M. J., Katardjiev, I. V., Carter, G., and Smith, R.: Analytic, geometric and computer techniques for the prediction of morphology evolution of solid surfaces from multiple processes, J. Phys. D Appl. Phys., 20, 870, https://doi.org/10.1088/00223727/20/7/008, 1987. a
Noether, E.: Invariant variation problems (translation from the German by M. A. Tavel; originally published as “Invariante Variationsprobleme”, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, 1918, p. 235–257), Transport Theor. Stat., 1, 186–207, https://doi.org/10.1080/00411457108231446, 1971. a
Nolte, D. D.: Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd edn., Oxford University Press, https://doi.org/10.1093/oso/9780198844624.001.0001, 2019. a
Osher, S. and Fedkiw, R. P.: Level set methods: An overview and some recent results, J. Comput. Phys., 169, 463–502, https://doi.org/10.1006/jcph.2000.6636, 2001. a
Osher, S. and Fedkiw, R. P.: Level set methods and dynamic implicit surfaces, vol. 153 of Applied mathematical sciences, Springer, NY, https://doi.org/10.1115/1.1760520, 2003. a
Osher, S. and Merriman, B.: The Wulff shape as the asymptotic limit of a growing crystalline interface, Asian J. Math., 1, 560–571, https://doi.org/10.4310/AJM.1997.v1.n3.a6, 1997. a, b
Palmer, B.: Anisotropic wavefronts and Laguerre geometry, J. Math. Phys., 56, 023503, https://doi.org/10.1063/1.4907215, 2015. a
Passalacqua, P., Do Trung, T., FoufoulaGeorgiou, E., Sapiro, G., and Dietrich, W. E.: A geometric framework for channel network extraction from lidar: Nonlinear diffusion and geodesic paths, J. Geophys. Res., 115, F01002–18, https://doi.org/10.1029/2009JF001254, 2010a. a
Passalacqua, P., Tarolli, P., and FoufoulaGeorgiou, E.: Testing spacescale methodologies for automatic geomorphic feature extraction from lidar in a complex mountainous landscape, Water Resour. Res., 46, W11535–17, https://doi.org/10.1029/2009WR008812, 2010b. a
Pazzaglia, F. J.: Landscape evolution models, in: The Quaternary Period in the United States, Elsevier, 247–274, https://doi.org/10.1016/S15710866(03)010121, 2003. a
Perlick, V.: Ray Optics, Fermat's Principle, and Applications to General Relativity, vol. 61 of Lecture Notes in Physics, SpringerVerlag, Berlin, Heidelberg, https://doi.org/10.1007/3540466622, 2000. a
Pinezich, J. D.: Propagation of singularities in nonconvex HamiltonJacobi problems: local structure in two dimensions, SIAM J. Math. Anal., 51, 3796–3818, https://doi.org/10.1137/18M1235570, 2019. a
Qian, J., Cheng, L.T., and Osher, S.: A level setbased Eulerian approach for anisotropic wave propagation, Wave Motion, 37, 365–379, https://doi.org/10.1016/S01652125(02)001014, 2003. a
Radjenović, B. and RadmilovićRadjenović, M.: 3D simulations of the profile evolution during anisotropic wet etching of silicon, Thin Solid Films, 517, 4233–4237, https://doi.org/10.1016/j.tsf.2009.02.007, 2009. a
Radjenović, B., Lee, J. K., and RadmilovićRadjenović, M.: Sparse field level set method for nonconvex Hamiltonians in 3D plasma etching profile simulations, Comput. Phys. Commun., 174, 127–132, https://doi.org/10.1016/j.cpc.2005.09.010, 2006a. a
Radjenović, B., RadmilovićRadjenović, M., and Mitrić, M.: Nonconvex Hamiltonians in three dimensional level set simulations of the wet etching of silicon, Appl. Phys. Lett., 89, 213102, https://doi.org/10.1063/1.2388860, 2006b. a
Radjenović, B., RadmilovićRadjenović, M., and Mitrić, M.: Level set approach to anisotropic wet etching of silicon, Sensors, 10, 4950–4967, https://doi.org/10.3390/s100504950, 2010. a
Rawlinson, N., Hauser, J., and Sambridge, M.: Seismic ray tracing and wavefront tracking in laterally heterogeneous media, in: Advances in Geophysics, vol. 49, Elsevier, 203–273, https://doi.org/10.1016/S00652687(07)490033, 2008. a, b
Reinhardt, L. J., Bishop, P., Hoey, T. B., Dempster, T. J., and Sanderson, D. C. W.: Quantification of the transient response to baselevel fall in a small mountain catchment: Sierra Nevada, southern Spain, J. Geophys. Res., 112, F03S05, https://doi.org/10.1029/2006JF000524, 2007. a
Richards, G. D.: An elliptical growth model of forest fire fronts and its numerical solution, Geophys. Res., 30, 1163–1179, https://doi.org/10.1002/nme.1620300606, 1990. a
Richards, G. D.: A general mathematical framework for modeling twodimensional wildland fire spread, Int. J. Wildland Fire, 5, 63–72, https://doi.org/10.1071/WF9950063, 1995. a
Rider, P. R.: The figuratrix in the calculus of variations, Ann. Math., 28, 640–563, https://doi.org/10.2307/1989068, 1926. a
Rigon, R., Rinaldo, A., RodriguezIturbe, I., IjjaszVasquez, E., and Bras, R. L.: Optimal channel networks: a framework for the study of river basin morphology, Water Resour. Res., 29, 1635–1346, https://doi.org/10.1029/92WR02985, 1993. a
Rinaldo, A., RodriguezIturbe, I., Rigon, R., Bras, R. L., IjjaszVasquez, E., and Marani, A.: Minimum energy and fractal structures of drainage networks, Water Resour. Res., 28, 2183–2195, https://doi.org/10.1029/92WR00801, 1992. a
Rinaldo, A., RodríguezIturbe, I., and Rigon, R.: Channel networks, Annu. Rev. Earth Pl. Sc., 26, 289–327, https://doi.org/10.1146/annurev.earth.26.1.289, 1998. a, b
RodriguezIturbe, I. and Rinaldo, A.: Fractal River Basins: Chance and SelfOrganization, Cambridge University Press, ISBN 9780521004053, 2001. a
RodriguezIturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., IjjaszVasquez, E., and Marani, A.: Fractal structures as least energy patterns: The case of river networks, Geophys. Res. Lett., 19, 889–892, https://doi.org/10.1029/92GL00938, 1992a. a
RodriguezIturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., Marani, A., and IjjaszVasquez, E.: Energy dissipation, runoff production, and the threedimensional structure of river basins, Water Resour. Res., 28, 1095–1103, https://doi.org/10.1029/91WR03034, 1992b. a
Royden, L. and Perron, J. T.: Solutions of the stream power equation and application to the evolution of river longitudinal profiles, J. Geophys. Res.Earth, 118, 497–518, https://doi.org/10.1002/jgrf.20031, 2013. a, b, c, d, e, f, g, h
Rund, H.: The Differential Geometry of Finsler Spaces, Springer Science+Business Media, Berlin, Heidelberg, https://doi.org/10.1007/9783642516108, 1959. a
Sethian, J. A.: Level Set Methods and Fast Marching Methods, Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, 2nd edn., Cambridge University Press, https://www.cambridge.org/jp/academic/subjects/mathematics/ (last access: 21 April 2022), 1999. a
Sethian, J. A. and Adalsteinsson, D.: An overview of level set methods for etching, deposition, and lithography development, IEEE T. Semiconduct. M., 10, 167–184, https://doi.org/10.1109/66.554505, 1997. a, b
Shemenski, R. M., Beck, F. H., and Fontana, M. G.: Orientationdependent dissolution of iron whiskers, J. Appl. Phys., 36, 3909–3916, https://doi.org/10.1063/1.1713969, 1965. a
Shen, Z.: Lectures on Finsler Geometry, World Scientific, https://doi.org/10.1142/4619, 2001. a, b
Shen, Z.: Differential Geometry of Spray and Finsler Spaces, Springer Science+Business Media, https://doi.org/10.1007/9789401597272, 2013. a
Shimada, H. and Sabau, S. V.: Finsler geometry, in: Generalizations of Finsler geometry, edited by: Antonelli, P. L., Kluwer, Dordrecht, 15–24, https://doi.org/10.1007/9789401142359, 2000. a, b, c
Shimada, H. and Sabau, S. V.: Introduction to Matsumoto metric, Nonlinear Anal.Theor., 63, e165–e168, https://doi.org/10.1016/j.na.2005.02.062, 2005. a
Sieniutycz, S.: Dynamic programming approach to a Fermat type principle for heat flow, Int. J. Heat Mass Tran., 43, 3453–3468, https://doi.org/10.1016/S00179310(00)000326, 2000. a
Sieniutycz, S.: A variational theory for frictional flow of fluids in inhomogeneous porous systems, Int. J. Heat Mass Tran., 50, 1278–1287, https://doi.org/10.1016/j.ijheatmasstransfer.2006.09.014, 2007. a
Sinclair, K. and Ball, R. C.: Mechanism for global optimization of river networks from local erosion rules, Phys. Rev. E, 76, 3360–3363, https://doi.org/10.1103/physrevlett.76.3360, 1996. a
Sklar, L. S. and Dietrich, W. E.: The role of sediment in controlling steadystate bedrock channel slope: Implications of the saltationabrasion incision model, Geomorphology, 82, 58–83, https://doi.org/10.1016/j.geomorph.2005.08.019, 2006. a
Slawinski, M. A.: Waves and rays in elastic continua, 3rd edn., World Scientific Publishing Company, https://doi.org/10.1142/7486, 2014. a
Small, A. and Lam, K. S.: Simple derivations of the Hamilton–Jacobi equation and the eikonal equation without the use of canonical transformations, Am. J. Phys., 79, 678–681, https://doi.org/10.1119/1.3553462, 2011. a, b
Smith, R., Carter, G., and Nobes, M. J.: The theory of surface erosion by ion bombardment, P. Roy. Soc. Lond. A Mat., 407, 405–433, https://doi.org/10.1098/rspa.1986.0103, 1986. a, b
Stark, C. P.: Geometric Mechanics of Erosion software package (GME
), https://geomorphysics.github.io/GME (last access: 21 April 2022), 2022a. a, b
Stark, C. P.: Geomorphysics Python library (GMPLib
), https://geomorphysics.github.io/GMPLib (last access: 21 April 2022), 2022b. a, b
Stark, C. P.: GME: Geometric Mechanics of Erosion (v1.0), Zenodo [code], https://doi.org/10.5281/zenodo.6373103, 2022c. a, b
Stark, C. P.: GMPLib: Geomorphysics Python Library (v1.0), Zenodo [code], https://doi.org/10.5281/zenodo.6373127, 2022d. a, b
Stark, C. P. and Stark, G. J.: Knickpoints, cusps, and the geomorphic surface Hamiltonian, in preparation, 2022. a, b, c, d
Stark, C. P., Stark, G. J., and dos Anjos, P. H. R.: Erosion and the geodesic equation, in preparation, 2022. a, b
Sullivan, A. L.: Wildland surface fire spread modelling, 1990–2007. 3: Simulation and mathematical analogue models, Int. J. Wildland Fire, 18, 387–403, https://doi.org/10.1071/WF06144, 2009. a
Tucker, G. E.: Landscape evolution, in: Crustal and Lithosphere Dynamics: Treatise on Geophysics, edited by: Watts, A. B., Elsevier, 593–630, https://doi.org/10.1016/B9780444538024.00124X, 2015. a
Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50, https://doi.org/10.1002/esp.1952, 2010. a
van der Beek, P.: Modelling landscape evolution, in: Environmental Modelling, John Wiley & Sons, Ltd, Chichester, UK, 309–331, https://doi.org/10.1002/9781118351475.ch19, 2013. a
Van Wagner, C. E.: A simple firegrowth model, Forestry Chron., 103–104, https://doi.org/10.5558/tfc451032, 1969. a
Virieux, J. and Lambaré, G.: Theory and observations – Body waves: Ray methods and finite frequency effects, in: Seismology and Structure of the Earth: Treatise on Geophysics, edited by: Romanowicz, B. A. and Dziewonski, A. M., Elsevier, 127–155, https://doi.org/10.1016/B9780444527486.000043, 2007. a
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Jarrod Millman, K., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I. L., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nat. Methods, 17, 261–272, https://doi.org/10.1038/s4159201906862, 2020. a
Vladimirsky, A. B.: Fast methods for static HamiltonJacobi Partial Differential Equations, PhD thesis, University of California, Berkeley, https://escholarship.org/uc/item/8k28k9t9 (last access: 21 April 2022), 2001. a
Wang, Y., Nemeth, T., and Langan, R. T.: An expandingwavefront method for solving the eikonal equations in general anisotropic media, Geophysics, 71, T129–T135, https://doi.org/10.1190/1.2235563, 2006. a
Weissel, J. K. and Seidl, M. A.: Inland propagation of erosional escarpments and river profile evolution across the southeast Australian passive continental margin, in: Prediction in Geomorphology, American Geophysical Union, Washington, DC, 189–206, https://doi.org/10.1029/GM107p0189, 1998. a, b, c, d, e
Whitham, G. B.: Linear and Nonlinear Waves, Whitham/Linear, John Wiley & Sons, Inc., Hoboken, NJ, USA, https://doi.org/10.1002/9781118032954, 1999. a, b
Willgoose, G.: Mathematical modeling of whole landscape evolution, Annu. Rev. Earth Pl. Sc., 33, 443–459, https://doi.org/10.1146/annurev.earth.33.092203.122610, 2005. a
Witcomb, M. J.: Frank's kinematic theory of crystal dissolution applied to the prediction of apex angles of conical ionbombardment structures, J. Mater. Sci., 10, 669–682, https://doi.org/10.1007/BF00566576, 1975. a
Woodhouse, J. H. and Deuss, A.: Theory and observations – Earth's free oscillations, in: Seismology and Structure of the Earth: Treatise on Geophysics, edited by: Romanowicz, B. A. and Dziewonski, A. M., Elsevier, 31–65, https://doi.org/10.1016/B9780444527486.00002X, 2007. a
Yajima, T. and Nagahama, H.: Finsler geometry of seismic ray path in anisotropic media, P. Roy. Soc. AMath. Phy., 465, 1763–1777, https://doi.org/10.1098/rspa.2008.0453, 2009. a, b
Yajima, T. and Nagahama, H.: Finsler geometry for nonlinear path of fluids flow through inhomogeneous media, Nonlinear Anal.Real, 25, 1–8, https://doi.org/10.1016/j.nonrwa.2015.02.009, 2015. a
Yajima, T., Yamasaki, K., and Nagahama, H.: Finsler metric and elastic constants for weak anisotropic media, Nonlinear Anal.Real, 12, 3177–3184, https://doi.org/10.1016/j.nonrwa.2011.05.018, 2011. a, b, c
Zhang, L., Parker, G., Stark, C. P., Inoue, T., Viparelli, E., Fu, X., and Izumi, N.: Macroroughness model of bedrock–alluvial river morphodynamics, Earth Surf. Dynam., 3, 113–138, https://doi.org/10.5194/esurf31132015, 2015. a
 Abstract
 Introduction
 Core principles
 Theory
 Implementation
 Raytracing solutions
 Results
 Discussion
 Conclusions
 Appendix A: Related studies
 Appendix B: Phase spaces and tensors
 Appendix C: ℱ_{*} is a metric function
 Appendix D: Finsler geometry and curved space
 Appendix E: Lagrangian and geodesics
 Appendix F: HJE and Hamilton action
 Appendix G: Notation
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References
 Abstract
 Introduction
 Core principles
 Theory
 Implementation
 Raytracing solutions
 Results
 Discussion
 Conclusions
 Appendix A: Related studies
 Appendix B: Phase spaces and tensors
 Appendix C: ℱ_{*} is a metric function
 Appendix D: Finsler geometry and curved space
 Appendix E: Lagrangian and geodesics
 Appendix F: HJE and Hamilton action
 Appendix G: Notation
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References