The Direction of Erosion

Abstract. The rate of erosion of a geomorphic surface depends on its local gradient and on the material fluxes over it. 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 self-constraint to convert an 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 surface-normal slowness. By rewriting surface tilt in terms of normal slowness components, and by 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 of solving for the evolution of an erosion surface: here we use it to derive Hamilton's ray tracing equations, which describe both the velocity of a surface point and the rate of change of the surface-normal slowness at that point. In this context, erosion involves two distinct directions: (i) the surface-normal 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. Analysis of the Hamiltonian shows that these rays carry boundary-condition information upstream, and that they are geodesics, meaning that erosion takes the path of least erosion time. This constitutes a definition of the variational principle governing landscape evolution. In contrast with previous studies of network self-organization, neither energy nor energy dissipation is invoked in this variational principle, only geometry.



Introduction
When geomorphologists describe the evolution of a landform, a direction of erosion is often invoked: for example, we speak of a bank cutting laterally, or a cliff retreating, or a knickpoint eroding upstream, or a river channel incising down into bedrock. 25 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 sub-vertically downwards, hewing closely to gravity, except at knickpoints where it occurs sub-horizontally upstream, and along the channel walls where it acts sub-horizontally 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 30 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 apparently many directions involved in the driving evolution of an erosion surface, and it is not immediately clear how to unify them.
In other physical contexts, the logical step would be to deploy tracers to establish motion. But an erosion surface, by defini-35 tion, destroys itself as erosion proceeds. If we were to tag a set of points on some initial surface, these points would immediately vanish as the surface evolves. The idea of persistent points moving with an erosion surface is an abstract one, and its utility is not obvious.
So, what does it really mean when we speak of the direction of erosion? Does it make sense to talk of erosion in terms of motion of surfaces? And can we learn anything fundamentally useful and insightful by doing so? The goal of this paper is to 40 answer these questions with an emphatic "yes", and to justify this conclusion in a rigorous fashion using the mathematics of differential geometry and classical mechanics.

Structure of the paper
The paper is organized into seven sections (including this introduction). Section 2 provides some background on geometric mechanics and its potential applications to geomorphology: it discusses how surface motion in three dimensions (3D) can 45 be described mathematically, makes some connections with geometric optics, and cites related applications in fields such as seismology, crystallography, materials science, and wildfire prediction. Section 3 develops these ideas in a more formal fashion, formulates a Hamiltonian theory of erosion front propagation (albeit limited to a 2D slice of 3D space), and explores aspects of this theory from the points of view of geometric mechanics and differential geometry. Section 4 explains how to generate numerical solutions of erosion front motion using ray tracing. Section 5 presents solutions for one class of problem: 50 a time-invariant erosion profile driven by slip on a bounding fault. These results are analyzed to reveal how erosion defined in the surface-normal direction connects with vertical and horizontal erosion rates and with erosion ray velocities; further analysis reveals how strong anisotropy lies at the heart of the erosion process. Section 6 looks at the broader implications of this Hamiltonian approach to erosion, and Sect. 7 draws some conclusions. Further discussion of more advanced mathematical aspects of the theory are provided in the appendices, Implicit surface motion in its most general form is described by the level-set equation (Gibou et al., 2018;Giga, 2006;Fedkiw, 2001, 2003;Sethian, 1999;Vladimirsky, 2001), in which φ(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 non-local, and only need 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 surface-normal erosion rate. So we can write ∂φ ∂t + ξ ⊥ |∇φ| = 0 (3) The notation ξ ⊥ is adapted from Osher and Merriman (1997).
With this equation we can describe how a 2D surface evolves in 3D space as a result of both erosion and deposition; the effects of tectonic displacement are easy to incorporate, as are spatiotemporal changes in weathering rates, vegetation and precipitation, 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 by limiting the scope of this equation. The simpler, more constrained form of the 110 geomorphic level-set equation makes it easier to tease out its fundamental behaviour.

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, or HJE: H (r, ∇φ; t) = ∂φ ∂t (4) 115 where each vector r tracks a point as it moves from one zero level-set of φ to another with velocityṙ = dr dt, while the front itself at that point moves in the direction ∇φ. These directions are not necessarily the same.

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 so-called "high frequency limit" at which seismic wavelengths are very small compared to the scale of wave propagation (e.g., Cervený, 1989Cervený, , 2005Cervený, , 2002Dellinger, 1997;Mensch and Farra, 1999;Rawlinson et al., 2008;Slawinski, 2014;Virieux and 130 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 modeled and spectral information is lost, the Hamiltonian approach has proven insightful, particularly when dealing with anisotropic media Antonelli et al. (2003); Slawinski (2002, 2003); Bucataru and Slawinski (2005);Červený (2002); Klimeš (2002); Yajima and Nagahama (2009); Yajima et al. 135 (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, 2007;Yajima and Nagahama, 2015).

Applications of the HJE to geomorphology
In geomorphology, Luke (1972Luke ( , 1974Luke ( , 1976) pioneered application of the HJE to the modelling of fluvial knickpoints as 140 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 explicit-surface form, and its ability to model implicit-surface motion was not considered.

Landscape as an erosion arrival-time surface 145
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 H(r, ∇φ). Imposing such time constancy in the Hamiltonian is not as onerous as it sounds. For example, if we use the eikonal equation to model landscape erosion driven by slip along bounding faults, its ray-tracing solution can still incorporate changes in fault slip rate over time (this topic is covered in Stark 150 and Stark (2021a)).
The implicit surface function φ that solves a static, eikonal form of Eq. (4) is a single-valued, 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 {x, y, z} that satisfy at each time step t the equation T (x, y, z) − t = 0 (5) 155 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: H (r, ∇T ) = const (6) 160 Points on the surface move with velocity vector v =ṙ, while the surface itself moves with a slowness covector ∇T . It is important to emphasize that ∇T is not a vector, which will likely come as a surprise to those not schooled in differential geometry. 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 ∇T 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 165 this study is use this measure to reveal the strong anisotropy of landscape erosion processes (see Sect. 5.3).

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 weathering-limited denudation of a rock cliff incised at its base by a river. By integrating the eikonal equation representing this erosion process, and by presenting level-set solutions 170 as isochrones of the cliff transect, they demonstrated how variations in rock erodibility can lead to highly irregular surface geometry such as overhangs.

Non-geomorphic erosion modeled with the HJE
There is a literature on erosion driven by non-geomorphic 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 Hamiltonian-based approach.

175
For example, both implicit surface motion and the HJE have been the basis for modeling erosion at microscopic scales in an engineering context. Frank (1958) employed the concept of surface-motion 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 method has been adapted to treat surface erosion at the micron 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 ( 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 185 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 paper by Katardjiev et al. (1989), which connects the HJE and its Hamiltonian to Huygens' principle and the concept of erosional wavelets.

Front motion obeys Hugyens' principle
Central to the ideas in the previous sections is Huygens' principle, one of the founding contributions to the field of optics.

190
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 195 of the front.
In wave propagation 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 (Perlick, 2000;Rider, 1926;Rund, 200 1959).
In other words, wavefront propagation can be tracked using either phase information or group information. For front propagation in general this equivalence translates into tracking using either (i) point velocities and their trajectories (ray paths), or (ii) point-wise front-normal slownesses and their ensemble motions.
Huygens' principle is best known for explaining wave propagation in inhomogeneous but isotropic media, where the indica-205 trices and figuratrices are spherical but vary in size from place to place; isotropy ensures that the group and phase propagation directions are the same. The principle is also often used to explain propagation in media whose anisotropy is ellipsoidal (Arnold, 1989), where the group and phase propagation directions are different. Recent efforts have further proved that the principle extends to non-ellipsoidal indicatrices and figuratrices representing a generalized form of anisotropy (e.g., Dehkordi and Saa, 2019;Innami, 1995;Markvorsen, 2016) expressed in terms of something called Finsler geometry (see Appendices B and C).

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 wind-driven 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 215 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 (1990Richards ( , 1995 formalized the fire front propagation 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 non-elliptical velocity indicatrices and frames the anisotropic motion in terms of Finsler structures. Finsler geometry is useful 220 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 C shows.

Fermat's principle of least travel time
Huygens' principle emphasizes HJE solution in terms of propagation of a front; Fermat's principle, on the other hand, empha-225 sizes 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 least travel time.
To be more precise, it states that each path obeys a variational principle which ensures its travel time is extremized; this extremal is almost always a minimum.

Ray tracing the motion of a front 230
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 Euler-Lagrange equations, and integrating them (see Appendix D). 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 235 is the subject of ongoing research.
There is a connection between the Hamiltonian ray tracing method developed here and the work of Luke (1972), Royden andPerron (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 240 governing equation in these studies, which forces elevation to be a single-valued 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 245 of the erosion process is hidden.

Variational principle governing erosion patterns
The key lesson to be learned here is that each point on an erosion surface follows the path of least erosion travel time. We therefore propose that the variational principle driving landscape formation is the principle of least erosion time. Since the Hamiltonian derived here (Sect. 3.8) is limited to a 2D slice ( Fig. 1), for now we can only be fully confident in making this 250 claim for 2D topographic profiles. Nevertheless, this variational principle will likely hold for 3D landscapes as well. If true, this assertion casts serious doubt on the idea that a principle of energy (dissipation) minimization governs landscape selforganization Rigon et al., 1993;Rinaldo et al., 1992Rinaldo et al., , 1998Rodriguez-Iturbe et al., 1992a, b; Rodriguez-Iturbe and Rinaldo, 2001). Sections 3.11, 3.18 and 6.3 develop this idea further.

255
In this section, the ideas presented above are formalized 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 (connections with the concept of phase spaces and with tensor calculus are given in Appendix A). 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 260 (Sect. 3.3) parameterized by the surface-normal covector, and how they constitute, broadly speaking, a form of geometric selfconstraint (Sect. 3.4). After specifying a separable form for such an erosion function (Sect. 3.5), we then show how the above ingredients lead, via the fundamental metric function (Appendix B), to a Hamiltonian description of erosion (Sects. 3.6-3.9; a supplementary discussion of how the geomorphic surface Hamiltonian inhabits a Finsler phase space, and a brief introduction to some pertinent concepts of differential geometry, are provided in Appendix C). Next we delve into the connections between the 265 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 ray tracing equations (Sects. 3.12-3.13) and discussion of some of their properties (Sects. 3.14-3.15). Next, the ray velocity angle is used to establish a measure of erosional anisotropy along with a closed-form for the Lagrangian (Sect. in preparation for undertaking numerical simulations. Note: we use superscripts for contravariant tensor components (e.g., r x ), and subscripts for covariant tensor components 275 (e.g., p z ); the Einstein summation convention (summing over similar tensor components) is adopted for brevity. Symbol usage is summarized in Table A1.

Tracking erosion with covectors
Imagine a locally planar surface undergoing constant erosion (Fig. 2), 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 280 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. 2, the time interval is chosen to be ∆T = 1 y and isochrones have been plotted for T (r) = {0, 1, 2, 3, 4, 5} y.
Let's fix the point of interest r at the location shown in Fig. 2. Here the surface-normal rate or speed of erosion is ξ ⊥ = 0.25 mm y −1 and surface tilt is β = 60 • . Written as a vector, the erosion rate is: with a direction normal to the surface and an angle β = 60 • to the vertical; its length or magnitude is the surface-normal erosion rate: Ideally, we should only have to compute the sine and cosine components to the erosion velocity vector ξ to get the horizontal 290 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 295 that we have written the erosion rate as a vector: we should instead express it as a covector.
Consider p in Fig. 2, which can be written as a function with single-row 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 surface-normal unit vector Because we employ here units of millimetres and years, n has a length of |n| = 1 mm. Over this distance n crosses four 1-year isochrones, so we obtain Now consider the vertical component of p (which is negative here) acting on n: counting downwards over a distance n z = 305 1/2 mm, we find one isochrone crossing, so: The horizontal component of p counts three isochrone crossings by the unit normal vector counting rightwards over a distance These components can be added together because p is a linear function; this summation gives p x (n) + p z (n) = 3 + 1 = 4 = p(n) 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. Normal slowness here is p = p(n) = 4 y mm −1 (Eqs. 13, 18) for a surface tilted at β = 60 • corresponding to a surface-normal erosion rate of ξ ⊥ = 1/4 mm y −1 (Eq. 8). Simple trigonometry applied to p gives the vertical and horizontal slownesses (Eq. 17), and their reciprocals are the vertical ξ ↓ and horizontal ξ → erosion rates (Eqs. 9, 10, and 17). The front covector is also the gradient of the arrival times, or isochrone density, given by p = ∇T , which counts the number of isochrones crossed in unit time in the front-normal direction (Eq. 21).
The function p is called a one-form in the terminology of differential geometry, and instances of p are called covectors. In 315 general, a one-form operates on a vector and returns a scalar. Here, p takes in a unit vector and returns the slowness of erosion in the direction of that vector. In optics and seismology, 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 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.

320
The slowness covector 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 aka the reciprocal erosion rate, and is given by 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 A for more details.

Gradient is a covector
The erosion slowness covector p has another facet: it is also the gradient of the arrival-time function T . To see why, consider again Fig If we measure (in Fig. 2) the change in T in the x direction over a distance n x = √ 3/2, we find that n x dT x = 3. Similarly, if we measure the change in the −z direction over a distance n z = 1/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). 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 contour-normal direction is the same 345 as the covector magnitude p. We can now invoke the gradient operator ∇ and state that which says that the gradient of the arrival time function T of the erosion surface is the same as its normal slowness covector p.
Note that equations 19 and 21 are arguably an abuse of notation in the way they mix operators and differentials; this problem could be resolved using Dirac or bra-ket notation (e.g., Cohen-Tannoudji et al., 2020, chapter 2).

A surface-normal erosion rate model
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 surface-normal 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 surface-normal speed of erosion that can be parameterized using the slowness covector.

355
To supply the third element, we can write a generic model for the surface-normal speed of geomorphic erosion that is a solely function of local fluxes and gradient: surface-normal erosion rate ∼ func(flow, gradient) 13 https://doi.org/10.5194/esurf-2021-59 Preprint. Discussion started: 28 July 2021 c Author(s) 2021. CC BY 4.0 License.
Some erosion phenomena, such as quasi-diffusive processes like rain splash, cannot be modeled 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 360 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 365 equation to be considered effectively local. The validity of this assumption is discussed at the end of Sect. 6.1.

Erosion imposes a geometric self-constraint
The process of landscape evolution represented by Eq. 22 is a kind of geometric straitjacket, or geometric self-constraint -in the sense that it essentially says the landscape obeys: 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 self-constraint is that, at its heart, geomorphic erosion is driven by a particular kind of 375 Hamiltonian. This Hamiltonian arises from how points on an erosion surface (for want of a better term) "see" their shortest path of erosion to the next set of surface points at little time later. The sections below explore this assertion in detail.

A specific form for the erosion equation
To see what form the geomorphic Hamiltonian may take, we need to make some concrete choices about the erosion model and its domain. As Fig. 1 illustrates, we restrict the model domain to a slice of 3D space -which is geometrically distinct 380 from slicing 2+1D space (e.g., Luke, 1972;Royden and Perron, 2013;Weissel and Seidl, 1998). This 2D transect is defined as a half-space with coordinates (x, z), aligned along the mainstream of a catchment with simple off-axis geometry, a drainage divide at x = x 1 , and predetermined variation of vertical position z with time at x = 0. The model catchment is shaped such that its upstream area is proportional to the square of the distance downstream from the divide (x 1 − x) 2 .
For the purposes of simplicity and consistency with past practice, we choose to work with a generalization of the stream 385 power model (e.g., van der Beek, 2013) tailored to describe motion of an implicit erosion surface (Sect. 2.1). The model defines the surface-normal speed of erosion ξ ⊥ as a separable function of a flow component ϕ(x) and a surface tilt component | sin β| η : a positive exponent. The flow component ϕ(x) encapsulates catchment geometry and the downstream accumulation of water and tools, while the compound effects of gradient on flow velocity, impact wear rate, etc. are wrapped into the |sin β| η term.
The exponent η here is approximately equivalent to the slope exponent n in the stream power model, which instead expresses the vertical erosion rate ξ ↓ in terms of |tan β| n -noting that ξ ↓ = ξ ⊥ / cos β as given by Eq. (10), and that for small surface angles sin β ≈ tan β. As such, we expect the slope exponent to lie in the range 1/4 ≤ η ≤ 2 (e.g., Royden and Perron, 2013).

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 p that represents the normal slowness of the surface at that point. The components of r and p correspond to the axes of the phase space inhabited by the Hamiltonian.

400
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 (r, 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 405 has no meaning in the current context. It has a dual, called the velocity phase 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: 1 noting that p x > 0 and p z < 0 for the half-domain shown in Fig. 1. Each point in phase space acts entirely independently.

410
The erosion equation given in Eq. (24) is now easy to convert into a form parameterized by the components of r and p: This equation defines the surface-normal reciprocal rate of erosion along a 2D profile, written in a form that neatly expresses the geometric self-constraint inherent to the geomorphic erosion process. This self-constraint is parameterized by vector position (r x , r z ) and covector normal-slowness (p x , p z ), which respectively locate a particular point on the surface and encode the 415 reciprocal speed of erosion orthogonal to the surface at that point.

The fundamental function
What we need to do now is reparameterize Eq. (27) to express the degree to which a coordinate (r, p) satisfies the geometric self-constraint imposed by the equation. This is easily achieved using Okubo's technique (Antonelli et al., 1993;Bao et al., by a positive function F * (r, p): and substituted back in, rearranging to make F * the subject: The function F * is a Hamiltonian, although for problems of this type it is better known as a fundamental function (see Ap-425 pendix B; note that an asterisk in used in F * for reasons that will become clear in Section 3.9). As a Hamiltonian 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. (27) must meet the condition: The power of a Hamiltonian comes from being able to trace a sequence of (r, p) across phase space for which this criterion 430 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 consider F * as the geomorphic surface Hamiltonian; a little more work is needed.
To clarify the behaviour of F * , consider the combined meaning of Eqs. (29) and (30). The value of F * at a location in phase space with coordinates (r, p) is equal to the normal slowness p 2 x + p 2 z implied by that coordinate, aka its reciprocal erosion 435 rate, multiplied by the erosion rate determined by the erosion process ϕ(r x ) p x η (p 2 x + p 2 z ) η/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 (r, p) the erosion rate is not reciprocal to the erosion slowness, and this product is not equal to one. Note that in physics (quantum field theory) jargon the physically valid coordinates are called "on shell" whereas the invalid coordinates are called "off shell": although 440 this terminology has no physical meaning here, we deploy it occasionally to clarify which parts of phase space we are talking about.

The geomorphic surface Hamiltonian
The problem with using F * as a Hamiltonian is its order-1 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 trans-445 form, 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: A prefactor of 1 2 is included to make subsequent derivations tidier.
which makes its metric tensor non-singular (if η = 1), and which puts the Legendre transform within reach.
We know from Eq. (30) that F * = 1 for trajectories across slowness phase space that correspond to physically viable behavior of surface points (aka "on shell"). So we can assert that for solutions of the erosion equation.

The geomorphic surface Lagrangian
The quadratic Hamiltonian H(r, p) has a dual quantity called the Lagrangian L(r, v), which operates in a counterpart phase 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 F . This function F is the dual of F * , and is similarly 460 order-1 homogeneous. Its quadratic L is similarly order-2 homogeneous: To make the link between the phase spaces of H and L, 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 H to the Lagrangian L (and vice versa) exploits this property and is achieved with the Legendre transform: A closed form for L requires several more pieces of the puzzle before it can be derived, and the eventual equation is very long 470 and unwieldy. The contrasting simplicity of H (Eq. 31) 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 F = 1 and L = 1 2 , in symmetry with F * = 1 and H = 1 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.

475
The preceding discussion is rather abstract, so we return to Huygens' principle to help visualize some of the ideas. Compare the following: in Sect. 2.10, we argued that stepwise motion of an erosion front can be imagined as the propagation of an ensemble of tiny erosional wavelets; in Sect. 2.11, we connected this idea to an extension of Hugyens' principle in which its elliptical wavelets (for the 2D model domain here) can take on a more general shape; in Sect.s 3.7 and 3.9, and Appendix B, we introduced the concept of a fundamental (metric) function and defined such a function F (r, v) for the velocity phase space 480 (which for the geomorphic surface Lagrangian is a Finsler tangent space). These three concepts are all the same thing: the shape of each erosional wavelet is also the shape implicitly defined by F , and these shapes are the elementary wavelets used by Huygens' principle to visualize how front motion proceeds. 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.

Figures 3 and 4 provide a graphic explanation.
485 Figure 3 visualizes a single erosional wavelet, its relationship both to the current erosion front at T = t and to the next at T = t + ∆t, the particular ray increment vector for which H = L = 1 2 , and the conjugate relationship of this vector to the front normal 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 p∆t/p 2 in the surface-normal direction given by p = ∇T ; (ii) the point moves a distance ∆r = v∆t in the ray direction r (Sect. 5.3). These directions are quite different, because the erosion process is 490 strongly anisotropic (Sect. 2.7).
Unconstrained, the point at r could be displaced onto any of the points along the "erosional wavelet" {∆r} (the locus of 2D indentation produced by erosion at that point alone). However, the only valid motion is onto the point r + ∆r where the tangent to the wavelet curve is orthogonal to the front increment p∆t/p 2 , i.e., the ray and front increments are conjugate to each other (Sect. 3.15).

495
When erosional wavelets at points along the surface are aggregated, moving T (r) onto T (r + ∆r) as shown in Fig. 4, the result is anisotropic front motion that obeys Huygens' principle. This new front can also be found simply by propagating the old front a distance ∆t/p in the direction p at each point r.

Fermat's principle as a least action integral
As first mentioned in Sect. 2.12, Fermat's principle (Holm, 2011) says that the path of a light ray between two points is such 500 that its optical length is stationary against any small deflection onto a nearby path (where optical distance is defined as the product of geometric distance and of the refractive index of the medium). This is equivalent to saying that the travel time is stationary, and, for all practical purposes, 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.

505
This principle is expressed mathematically by writing an action functional S γ , in terms of the static Lagrangian L, 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 L 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 phase space to which the geomorphic surface Lagrangian L belongs, we can be sure that the action is minimized (on other types of phase space the stationary action may be a maximum). Since L is independent of t, we can deduce that γ 0 is the path of (locally) 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 La-  The fixed divide is enforced by mirroring the profile at x = x1 and by imagining a symmetrical profile to x = 2x1, such that symmetrically generated rays annihilate each other at a cusp formed at x = x1. The boundary condition imposed at x = 0 (and simulated at x = 2x1) is a constant vertical erosion rate ξ ↓ 0 , mimicking the behavior of a vertical normal fault slipping at a constant rate ξ ↓ 0 at the boundary. The initial value of the front slowness covector p at x = 0 is chosen such that the surface tilt β0 and vertical slowness pz 0 are consistent with this rate. The model therefore simulates a horst block undergoing constant uplift and consequent erosion. Model topography is solved for by constructing surface isochrones {T (r)} from the rays. Since rays are traced only from the boundary, and none from an initial surface, the isochrones are time-invariant (more complex initial and boundary conditions are explored in Stark and Stark (2021a, b)). The standard term for such topography is "steady state", but the Hamiltonian dynamical system has no stable point and the term is somewhat misleading here.

520
The fundamental function F * generates a slowness phase space spanned by r and p on which the geomorphic surface Hamiltonian H(r, p) operates, and we have a simple expression for H given by Eq. (31). We inferred the existence of a dual fundamental function F that generates a velocity phase space spanned by r and v on which a Lagrangian L(r, v) operates, but we have yet to obtain expressions for F and L. 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.

525
Our starting point is to examine the differentials of H and L and to compare them. The geomorphic surface Hamiltonian defined in Eq. (31) is static, meaning that it is constant over time, so its differential is: The differential of its counterpart Lagrangian L(r, v) is: Rearranging, we have an equation that contains the Legendre transform given in Eq. (36), Consequently we have a second expression for the differential of H: Equating the terms in dH defined by this equation with those in Eq. (39), we obtain: The next step is subtle but important. Every coordinate (r, v) in velocity phase space is (potentially) an initial position and velocity for a point on some initial erosion surface. Similarly, every coordinate (r, p) in slowness phase space is (potentially) 540 an initial position, surface orientation, and reciprocal surface-normal erosion rate for a point on that initial erosion surface.
However, most such phase space coordinates do not correspond to real-world points lying on physically viable paths {γ 0 } that obey the principle of least erosion time established in Eq. (38). Conversely, for the locations in phase space that do lie on a paths of least action, we can write:

545
Returning to the variation integral in Eq. (38), 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:

The meaning of Hamilton's equations
Hamilton's equations are coupled first-order ordinary differential equations (ODEs) whose integration gives the motion of a 555 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 = p (which encodes both the local tilt of the erosion surface and its reciprocal rate of erosion p = 1/ξ ⊥ ). 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. (48) define the rates of change of the coordinates (r, p) in terms of the gradient components 560 of the Hamiltonian. Since the Hamiltonian is a constant H = 1 2 along a ray or trajectory (Eq. 33), motion across the phase space must follow coordinates (r, p) that trace a contour of H. This is achieved by moving r in the direction ∂H ∂p i and p in the direction −∂H ∂r i , which is to say, orthogonal to the Hamiltonian gradient.
Hamilton's equations take concrete form if we substitute the expression for H in Eq. (31) into Eq. (48). 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, 3.14 Constancy of the vertical erosion rate along a ray The erosion model defined in Eq. (24) is independent of elevation. This makes the Hamiltonian H independent of the vertical coordinate r z , which leads to the zero element in˙ p in Eq. (50), 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 H, and therefore in L, which implies a law of conservation of vertical slowness for the ray tracing 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 ξ ↓0 of a point is carried unchanged along its ray trajectory as the surface to which it is attached moves:

585
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. 5). 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. 50).

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 590 normal-slowness covector pair: they are conjugate to each other (Figs. 3 and 4), which is to say, their inner product is one. To prove this, consider the following property of an order-2 homogeneous function like H: 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. 5.3), which means that conjugacy also constrains the angular disparity between the ray and 600 front-normal directions.
where the root is chosen to conform with the requirement that β > 0 for small α. By comparing α and β we can measure This switch in ray orientation as a function of slope scaling exponent η, which is illustrated in Fig. 6, 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, 2021a).

Constancy of the Lagrangian
We can exploit conjugacy to reveal important behaviour of the fundamental function F and the related Lagrangian L. Since L is (like H) order-2 homogeneous, it has the property Using the fibre derivative form of p in Eq. (35) and the definition of the Lagrangian in Eq. (34), we can deduce that i.e., the Lagrangian has the constant value of 1 2 just like the Hamiltonian (Eq. 33) for physically realistic aka "on shell" states. 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. 5): a more general theory that relaxes these restrictions would lead to non-constancy of L and H.

HJE and Hamilton action
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 H(r, p) in its HJE form and solve erosion front propagation using grid-based methods. In practice, numerical solution of this kind of eikonal equation is not straightforward (see Sect. 6.2 and Appendix B).
The HJE is nevertheless instructive if we examine it in the context of some important concepts of classical mechanics. For 635 example, Hamilton's principal function S(r, t), which is the Hamilton action S γ (see Eq. 37) plus a constant, is Use of the Legendre transform (Eq. 36) yields Assuming the points {r} all lie on a path γ 0 of least erosion time ("on shell"), we can write which is the standard form for the HJE. Now consider the arrival time function T (r), whose total time derivative is, given Eqs. (21), (54), and ∂T ∂t = 0: 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 Eq. (33) and Eq. (61) connects S(r, t) with T (r), H(r, p) and time t: 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 660 gradient ∇T (the directional density of T isochrones) satisfies the static Hamiltonian H = 1 2 .

Measuring slope along the erosion front
Since the Hamiltonian tracks motion of the erosion front in a phase space spanned in part by the surface-normal 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 finite-difference 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 670 direction, aka the vertical erosion rate: Ideally, all three measurements of the topographic gradient should be equal. In practice, β ts is computed non-locally while β 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 illustrated 675 in Fig. 7.

Formulation of erosion model
The analysis presented above is quite general: the only restriction on the erosion model has been (Eq. 24) that it take a separable form in flow and slope as ξ ⊥ := ϕ(x) | sin β| η . We can further generalize and note that the slope component can Previous studies related to our work (Luke, 1972;Royden and Perron, 2013;Weissel and Seidl, 1998) have focused on the stream power incision model or SPIM (e.g., Lague, 2014), and because of its relative simplicity we make a broadly similar choice here. The SPIM approach asserts that, in channels, vertical erosion rate ∝ (area) m (gradient) n (74) 685 under the assumption that upstream area, suitably scaled, is a good proxy for the rate contribution of the volumetric flow of water per contour width and where gradient is defined as tan β. In a modest modification, we instead assert that, along a channel transect, surface-normal erosion rate ∼ (area) µ (slope) η (75) Figure 7. Estimation of the surface-normal angle from vertical β, aka the angle of the surface from horizontal for (a) η = 3 2 and (b) η = 1 2 , and with µ/η = 1 2 . This angle can be computed in three ways; their mutual consistency shown here provides a partial validation of the ray tracing method.
where "slope" is sin β. Note that, since ξ ↓ = ξ ⊥ / cos β (Eq. 10), the two models are the same when η = 1. More broadly, we 690 can make the approximate asymptotic associations η ⇔ n and µ ⇔ m. In concrete terms, we set where x 1 is the transect length, such that (1 − x/x 1 ) is dimensionless distance downstream, ϕ 0 is a rate constant, and ϕ 0 ε is the value of ϕ(x) at the divide. This model generates a channel profile with the asymptotic scaling at constant erosion rate (see, e.g., Lague, 2014) and for low-to-moderate slope angles for which tan β ≈ sin β. To keep complexity to a minimum, in the remainder of the paper we fix the exponent ratio (aka "concavity index") at a constant µ/η = 1 2 such that the resulting slope-area scaling is close to that typically observed (e.g., Beeson and McCoy, 2020;Flint, 1974;Lague, 2014;Royden and Perron, 2013).

Model domain
The model domain is limited to a vertical transect (Fig. 1) in the theoretical development by Section 3. The particular setup implemented here is a long-stream profile ranging from a fault-slip, flow-exit boundary at x = 0 to a fixed drainage divide at x = x 1 . The divide is pinned by mirroring (Fig. 5) the main profile with a symmetrical "image" profile spanning x ∈ [x 1 , 2x 1 ] and bounded by a paired fault at x = 2x 1 .

710
In this model geometry, rays that initiate at x = 0 and propagate in the positive x direction are annihilated at x = x 1 when a paired ray, initiated at the same time at x = 2x 1 , arrives from the opposite direction. As such, the model induces a cusp to form at x = x 1 , although its formation is not modeled here (instead, rays from x = 0 are simply truncated at x = x 1 ). Cusp self-formation and propagation, which can be modeled explicitly by locating and tracking ray intersections (Stark and Stark, 2021a).

Boundary and initial conditions
Ray tracing across the 2D transect here (Fig. 6) involves the numerical integration of Hamilton's equations in the form of four coupled, first-order ODEs forṙ x andṙ z (Eq. 49),ṗ x andṗ z (Eq. 50). These are first-order differential equations in time alone, so for each ray we need only supply four initial conditions -i.e., r x0 , r z0 , p x0 , p z0 -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 720 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 solving the formation of a time-invariant profile generated by constant slip at x = 0, and so we need only perform ray tracing from the boundary at x = 0; this avoids the complexity of the supplying boundary conditions along an initial topography. For the model 2D transect (Fig. 1), the initial horizontal position for all rays is fixed at the fault r x0 = x = 0.

725
The vertical fault slip rate is mimicked by forcing the boundary rate of vertical erosion ξ ↓0 to match it. Consequently, the initial vertical position of a ray initiated at time t = t 0 is given by simple integration of this slip (aka erosion) rate: Correspondingly, the initial vertical component of the ray slowness covector must be consistent with the vertical slip (aka erosion) rate, requiring p z0 = −1/ξ ↓0 . Sinceṗ z = 0, this covector component remains unchanged throughout the propagation of the ray, as Sect. 3.14 has discussed, and so the number of coupled ODEs that need to be solved is effectively reduced to 730 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., tan β 0 = −p x0 /p z0 . As such, the initial value of the slowness covector p encodes the velocity boundary condition in both its direction and magnitude.

Direct integration 735
For the simple scenario of a time-invariant profile, the erosion equation (Eq. 24) 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 ξ ↓ = ξ ↓0 , and thereby to manipulate Eq. (24) to expose its straightforward dependence on surface tilt β and position x (through ϕ(x)): which, for sin β > 0, can be converted into We can combine this equation with Eq. (72) to obtain a polynomial in surface gradient tan β = dz/dx, and constrain it using 745 the result (Eq. 52) that −p z = 1/ξ ↓ = 1/ξ ↓0 along the whole ray and thus everywhere along a time-invariant profile. The resulting polynomial in dz/dx takes the form, for η = 3 2 , and for η = 1 2 , By finding at each x the roots of the respective polynomial, and by integrating over x the appropriate root dz/dx, we obtain the surface elevation z(x) as a function that combines horizontal distance upstream x, the flow component function ϕ(x), and the slip rate at the boundary ξ ↓0 : In order to obtain values of z(x), numerical root finding is required for all but the simplest choices of η and ϕ(x). These direct 755 integrations are useful, because they serve to validate the ray-traced solutions, as illustrated in Fig. 8.

Numerical integration method
After some experimentation, the most accurate quadrature or numerical integration scheme for ray tracing with Eqs. (49) and (50) 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.,760 2020). Simpler, and lower-order Runge-Kutta quadrature methods also work well for most choices of model parameters, as does the high-order Runge-Kutta, dense output, DOP853 method (see Hairer et al., 2008, p. 194).

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 (Eqs. 49, 50) with the above boundary conditions (Sect. 4.2). This constitutes the tracing of a single reference ray (Fig. 6), which suffices for construction of a time-invariant topographic profile (see below). More rays need to be traced if we want to handle time-variable boundary conditions, evolution from an initial 770 topography, or the transition between an initial surface and a slip boundary. 3. generate a k th later ray r k∆t0 (t + k∆t 0 ) with initiation time k∆t 0 by making a copy of the reference ray, displacing it vertically by −ξ ↓0 k∆t 0 , and pasting it at (0, r z0 − ξ ↓0 k∆t 0 ); 4. truncate the copied ray at the point r k∆t0 (T − k∆t); 780 5. repeat step (ii) until k∆t 0 ≥ T ; 6. collate the truncation points to generate a continuous curve T (r).
Some of these steps also entail interpolation and resampling.
This procedure generates the time-invariant isochrone T (r) formed by the constant fault-slip ξ ↓0 boundary condition at x = 0 (Figs. 5, 8). Repetition of the procedure (or a simple copying of the solution), combined with a progressive offset of the 785 initial ray location r z0 at the boundary, simulates vertical normal-fault-driven erosion of a topographic profile at steady state in the reference frame of the (bedrock) substrate of the footwall (Fig. 9). Analysis of the shape of the time-invariant profile gives the component erosion rates (Figs. 10, 11) and their anisotropy (Figs. 13, 14).

Results
Numerical solutions of time-invariant topographic profiles (Sect. 4) provide a way to validate the geomorphic surface Hamil-790 tonian, test the inferences drawn from it (Sect. 3), and check that ray tracing by integrating Hamilton's equations is a viable approach to modeling surface erosion. More broadly, they demonstrate the practicality of modeling landscape evolution in true 3D (albeit with a 2D transect) with an equation that describes the erosion rate in the surface-normal direction, rather than in 2+1D (which would be reduced to 1+1D here) with erosion constrained to act in the vertical direction only. Figure 8 illustrates ray-traced time-invariant solutions for two choices of the slope exponent η ∈ 3 2 , 1 2 in the model equation (Eq. 24) for surface-normal erosion rate ξ ⊥ . Each ray-traced isochrone T (r) is compared with an isochrone obtained by directly integration (Sect. 4.3), and in each case the match is excellent. Sequences of erosion surfaces resulting from similar time-invariant solutions are shown in Fig. 9.
These solutions illustrate an important behaviour of the rays: for values of the slope exponent η > 1 the ray velocities 800 always have a positive vertical componentṙ z > 0, whereas for η < 1, the vertical componentṙ z always has a negative vertical componentṙ z < 0.

Erosion rates
The component rates of surface erosion along a time-invariant profile behave in some obvious and some non-obvious ways. (Figs. 10, 11).

805
There is no fault-driven uplift in the model -there is only fault slip at the left-boundary simulated by a vertical velocity condition at x = 0. Therefore, the time-invariant profile as a whole must move downwards, without changing shape and with a constant rate of vertical motion ξ ↓ . If we were to track the profile in a reference frame moving with the slipping fault, which would be equivalent to adding a vertical uplift rate, the profile would remain static, and all surface isochrones would be colinear.
Such a solution would constitute what is better known as a steady-state topographic profile.

810
The surface-normal erosion rate can be computed in two ways from the ray-tracing results (Figs. 10a and 11a). We can simply use the fact (Eq. 18) that normal slowness is the reciprocal of the normal speed ξ ⊥ = 1/p. Or we can invoke ray-front conjugacy and deduce that: As Figs. 10a,b and 11a,b show, the surface-normal erosion rate decreases progressively towards the divide, and this trend 815 clearly originates in the angular disparity (α − β + π 2 ) in the above equation (see Sect. 5.3). 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. It turns out that ξ ⊥ = 1/p is the more numerically stable estimator (Figs. 10a and 11a),

Anisotropy
The relative directions of motion of the erosion front and of a point on the front provide an insightful measure of erosional anisotropy. The front moves with a normal slowness covector tilted at an angle β to the vertical, while a point on the front moves with a ray vector tilted at an angle α to the horizontal (Fig. 12). The angular disparity (α − β), with 90 • added to put both angles in the same frame of reference, is the quantity to focus on.    ber of cross-tick "bones" approximates slowness p), while arrows represent the ray velocity vectors. The colour attribute of each symbol indicates its angular disparity (α − β + 90 • ). The degree of anisotropy is evident in the strong angular disparity of r and p for the same choices of η ∈ 3 2 , 1 2 . 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 840 the ray vector angle. At the divide, the erosion process (for this simplistic erosion model) is approximately isotropic.

Geometry controls (almost) everything
A key aspiration of this paper is to clarify what we mean when we speak of the direction of erosion. Our central tenet has been that, from a mathematical point of view, it is only meaningful to consider erosion of a static object in its surface-normal 845 direction (Eq. 3). Working from this premise, with the help of geometric mechanics, we have found unexpected complexity in this seemingly benign issue.
The concept of a covector is pivotal to our theory (Sect. 3.1). Once we realize that the rate of surface-normal erosion (which is imposed by the model erosion process) can be written in terms of the normal-slowness covector (which is the consequent motion of the surface), it takes only a few short steps to reach the geomorphic surface Hamiltonian. Hamilton's ray tracing 850 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 self-constraint.
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 855 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. 14 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, 860 driving the surface downwards; this is indeed the case if η < 1. Remember that the topographic profiles obtained by numerical solution here are time-invariant solutions without an uplift term, as visualized as time slices in Fig. 9; 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 surface-normal direction and subhorizontally in the ray direction -is an uncomfortable and apparently very abstract notion, but it has physical consequences.

865
It means erosion drives information about boundary conditions upstream subhorizontally (generally with an upward tilt steeper than the channel) while also driving motion of the whole profile downwards subvertically. The time scale on which boundary condition information propagates into the interior is the time it takes for a point to travel along a ray. This information may be Figure 14. Anisotropy of erosion for time-invariant profiles with η = 3 2 and η = 1 2 , and with µ/η = 1 2 . Ray velocity vectors v are represented by arrows (where arrow length approximates speed v); normal-slowness covectors p are represented by fishbone symbols (where the number of cross-tick "bones" approximates slowness p). The degree of anisotropy is evident both in the divergence of the r and p directions, and in the colour attribute used to visualize their angular disparity (α − β). Anisotropy progressively decreases upstream. the instantaneous slip rate on the boundary fault (Sect. 4.2), or the erosion rate at a point on an initial profile (Stark and Stark, 2021a).

870
For these time-invariant solutions, all rays are identical and the fault slip condition does not change. If the 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, 2021a). 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 right-hand half-domain (Fig. 5). The crucial difference is that intersection 875 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 ray-like 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 880 traversed by these characteristics has no concept of the surface-normal 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;. 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,ṗ z = 0, and thus the 885 vertical component of the surface erosion rate is constant, ξ ↓ (t) = −1 p z (t) = −1 p z0 = ξ ↓0 . For a time-invariant (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 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 (Royden and Perron (2013) came to a similar conclusion).

890
The time scale of this information transfer is of crucial importance. If it is small relative to the time scale on which drainage divides move laterally and significantly change accumulation areas and flows, then the assumption made in Sect. 3.3 is valid: namely, that the flow -at every point on the surface where flow influences erosion -can be parameterized by its surface geometry, aka 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 895 periodically. This is not to say that the geomorphic surface Hamiltonian theory is invalidated if the time scale 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.

A geomorphic surface Hamiltonian in 3D
While the geomorphic surface Hamiltonian developed here is limited to erosion-driven motion of a linear front in 2D space, 900 the ultimate goal is to construct a theory for surface evolution in 3D space. Several conceptual as well as 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 905 would run counter to our core premise that what matters is the geometric self-constraint imposed by geomorphic processes, not the details of those processes. A parameterization of flow focusing in channels will be required: one that encapsulates channel cross-sectional geometry without describing it explicitly.
Another challenge hinges on the assumption of locality. A first-cut 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 910 that catchment shapes can self-form. The question of time scales raised in Sect. 6.1 will still apply in 3D: if the time scale 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 915 than for a line in 2D, particularly when dealing with ray intersections and cusp formation. An obvious alternative lies in the fact that the theory employs an implicit surface function to describe landscape geometry: we can resolve erosion front motion on a regular grid and use a level-set method such as fast marching to solve the geomorphic HJE. This will have to be done with some care, however, because of the Finsler nature of the geomorphic surface Hamiltonian and its inherent anisotropy. Direct application of one-pass fast marching will not be possible, because even the most advanced algorithms for fast marching are 920 limited to metrics whose anisotropy is velocity-independent or Riemannian (Mirebeau, 2014a(Mirebeau, , 2019Mirebeau and Portegies, 2019), and to a small subset of Finsler metrics (Mirebeau, 2014b). It may nevertheless be possible to use Riemannian fast marching in an incremental, level-set propagation fashion to approximate integration of a Finsler HJE; other methods may also be viable (Moser, 1991;Qian et al., 2003;Rawlinson et al., 2008;Wang et al., 2006).

The variational principle is not energy minimization 925
There is a substantial body of work founded on the idea that landscapes self-organize 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 Rigon et al., 1993;Rinaldo et al., 1992;Rodriguez-Iturbe et al., 1992a, b;Rinaldo et al., 1998) and is comprehensively reviewed in the book by Rodriguez-Iturbe and Rinaldo (2001). In this section we compare and contrast OCN theory with our theory of the geometric mechanics of erosion.

930
OCN theory is framed in terms of the self-optimization 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 935 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.

940
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 time-evolution 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 (sects. 2.12, 2.14, and 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 950 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 time-invariant shape, but this shape is the outcome of geometric interaction rather than a mechanism of energy-dissipation minimization.
Comparison of the two theories is a little premature, because our theory needs further development if we want it to describe 955 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 aka many possible drainage network configurations; in our 2D theory, only one time-invariant state (a simple linear profile), imposed by the model erosion (Eq. 24), is possible. Nevertheless, it will be interesting to see if a full 3D theory can throw light 960 on what drives landscape self-organization and channel network formation: whether these phenomena arise primarily from the geometric self-constraint imposed by geomorphic erosion, and if so, the extent to which the process of energy minimization is complimentary.

Conclusions
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 965 are in essence saying that the rate of change of landscape geometry is a function of that landscape 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 erosion".
A key assumption in the development of this theory is the understanding that the rate of surface erosion is a velocity normal to the surface: differential geometry tells us that no other direction has any meaning. With this in mind, we can write a generic 970 model of erosion -albeit one that does not explicitly treat the effects of sediment cover, weathering, or any hydrodynamic effects -as the normal slowness of the surface as a function of tilt angle and upstream area. Using a simple mathematical trick (a scaling substitution), and by writing tilt in terms of slowness 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 975 orientation and slowness of the surface at that point. The Hamiltonian therefore occupies a six-dimensional phase space.
Although such extra dimensionality may seem to be a distracting mathematical abstraction, it provides some real insights.
Study of the Hamiltonian and its phase space makes it clear that surface evolution simultaneously involves two distinct directions of motion: the surface (at a given point) moves in the surface-orthogonal 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 980 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 ray tracing equations -derived from the Hamiltonian by simple differentiation -which express the motion of a surface point and its allied surface-normal 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 are carried upstream along the rays.

985
There is an important dependence of ray tracing on surface tilt: if the model erosion rate 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 990 us to deduce that the erosion rays traced by surface points are geodesics, In other words, they follow paths of locally shortest erosion time: this is the variational principle that governs geomorphic surface erosion. It appears that energy dissipation need not be invoked, and that instead all that matters is geometry.
Code availability. Software to solve and visualize the model has been developed as a pair of platform-independent Python3 packages (built around the SymPy, NumPy and SciPy libraries (Harris et al., 2020;Meurer et al., 2017;Virtanen et al., 2020)) and a set of allied IPython/Jupyter notebooks (Kluyver et al., 2016). The base utilities package is available at the "Geomorphysics Python library" (GMPLib) repository on GitHub at Stark (2021a) and archived on Zenodo at Stark (2021d) . The allied code and notebooks to solve the equations presented in this paper, both algebraically and numerically, is available at the "Geometric Mechanics of Erosion" (GME) repository on GitHub (Stark, 2021b) and archived on Zenodo at Stark (2021c). GMPLib release version 1.0 and GME release version 1.0 were used to generate the results presented in this paper.

1000
Appendix A: Phase spaces and tensors 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 p) at that point. Bundled together, the tangent 1005 spaces for all points in real space constitute a tangent bundle or velocity phase space, while the union of cotangent spaces forms a cotangent bundle or slowness (classically called "momentum") phase space. These two phase 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 1010 rank cannot be combined arithmetically; instead, operations such as contraction are needed to combine them. For example, in Eq. (13), the action of covector p on the unit vector n is a tensor contraction: p(n) = p i n i = i∈{x,z} p i n i = p x n x + p z n z (A1) 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 ∈ {x, z}). Upper indexes 1015 are used for contravariant tensor components; lower indexes are used for covariant tensor components. The utility of tensor calculus, specifically the use of metric tensors to express the direction of erosion, is explored in more detail in Stark and Stark (2021b).

Appendix B: F * is a metric function
The fundamental function F * has three key properties that are valid for a domain D of {r x , r z , p x , p x } phase space correspond-1020 ing to physically reasonable values of surface tilt and erosion rate: 1. Positive, order-1 Euler homogeneity in the parameter p: if the covector p in F * is scaled by a positive scalar λ > 0, the reparameterized function equals the original function scaled by λ, F * (r, λ p) = λF * (r, p) for λ > 0 (B1) where the 1 in "order-1" refers to the exponent in λ on the right hand side of this equation 1025 2. Regularity: F * is smooth, in that it can be differentiated infinitely many times without encountering a discontinuity or undefined value.
3. For η > 1, the Hessian of (F * ) 2 , i.e., the Hessian of the Hamiltonian H, is regular, which means that all the values of the matrix g ij are positive, where g ij := 1 2 ∂(F * ) 2 ∂p x ∂p z (B2) 1030 We deduce that g is positive-definite (all eigenvalues are real and positive), and that F * is strongly convex.
Given these properties, F * constitutes a type of Finsler metric (Bao et al., 2000;Shimada and Sabau, 2000). This means that F * 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. 2.13 and Appendix C; see Bao, 2007). In other words, the shortest time path between two points in the corresponding real space may not be a straight line.
In geomorphology, we are used to dealing with equations that operate in a flat Euclidean geometry where the space is spanned by Cartesian {x, y, z} coordinates. In such a space, distances can be 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 1055 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.

1060
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 non-flat phase space that both simplifies their solution and exposes their fundamental properties and behaviour.
The geomorphic surface Hamiltonian H, which arises from the transformation of an erosion equation, operates in such 1065 a non-flat space. Distance and travel time are not Euclidean measures on this phase space, because the fundamental metric function F * that defines H has the properties described in Sect. 3.7.
The curved nature of non-Euclidean 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 F or F * , 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 1070 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 F * , the metric tensor elements vary with both r and 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).

1075
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, F * specifies that the slowness phase space is a co-Finsler or Cartan space, and its dual F specifies that the velocity phase 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 1080 directional dependence of erosion time measurement makes this calculation non-trivial.
The fact that the geomorphic surface Hamiltonian operates in a Finsler geometry has profound consequences for the construction of erosion-driven equations of motion, for the variational principle that underlies how landscape shape evolves, and for the concept of erosional anisotropy. These consequences are explored in the next sections.

1090
The Lagrangian method of ray tracing the motion of an erosion front is taken up in detail in Stark and Stark (2021b). 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 DEMs (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 1095 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 ray tracing equations derived from the geomorphic surface Hamiltonian. This assertion is proved in Stark and Stark (2021b).