glacial erosion and its implementation in large-scale landform-evolution models

Modeling glacial landform evolution is more challenging than modeling fluvial landform evolution. While several numerical models of large-scale fluvial erosion are available, there are only a few models of glacial erosion, and their application over long time spans requires a high numerical effort. In this paper, a simple formulation of glacial erosion which is similar to the fluvial stream-power model is presented. The model reproduces the occurrence of overdeepenings, hanging valleys, and steps at confluences at least qualitatively. Beyond this, it allows for a seamless coupling to fluvial erosion and 5 sediment transport. The recently published direct numerical scheme for fluvial erosion and sediment transport can be applied to the entire domain, where the numerical effort is only moderately higher than for a purely fluvial system. Simulations over several million years on lattices of several million nodes can be performed on standard PCs. An open-source implementation is freely available as a part of the landform evolution model OpenLEM.


75
The product of ice thickness and total velocity, h(v s + v d ), defines the ice flux per unit width. Inserting this property into the mass balance of the glacier yields a differential equation for the thickness h. If the topography of the ice surface (H + h, where H is the topography of the bedrock) is considered as the variable, it is a diffusion equation with a high and strongly variable diffusivity. This property makes its numerical treatment challenging and allows only for small time increments.

95
According to Eqs. (5 : 4) and (4 : 5), the ratio of v d and v s is proportional to h 2 . Thus, the relative contribution of deformation to the total flux converges to zero in the limit of a thin layer, provided that the ice is not frozen to the bedrock. So let us assume for the moment that the entire ice flux is dominated by sliding. If we consider a rectangular cross section of a width w, the total ice flux (volume per time, not per unit width) is 100 The thickness h can be eliminated by combining Eqs. (5) and (7), which yields v s ∼ q i w ψ−1 ψ S.
If we assume that a given erodibility K refers to a uniform reference precipitation p 0 , the long-term mean discharge from a 155 catchment of size A is q = p 0 A. If we further take into account that erosion rather depends on discharge than on catchment size, :::: Then : Eq.
where p is the total precipitation, p i is the part of p that is converted into ice, and s is the pixel area of the considered grid cell.
The sum extends over all neighbors that deliver their discharge to the considered site, called donors in the following.
This concept is widely used in modeling fluvial erosion. However, the spatial scales are different here. At least small rivers are usually narrower than the mesh width of the lattice. As this may cause problems for detachment-limited erosion in combination 210 with hillslope processes, approaches taking into account that rivers cover only a part of the pixels of the grid were developed (Howard, 1994;Perron et al., 2008;Pelletier, 2010). Glaciers are, however, often wider than typical mesh widths of some tens to hundreds of meters. Thus, glacial erosion acts over an area around to the cardinal flow path, called swath in the following. The scaling arguments developed in the previous section suggest that the width of the swath is a function of the ice flux according to Eq. (12 :: 13). This concept is in principle the same as following the cardinal flow path with a pen of a variable width w.
Let c 0 be a point on the cardinal flow path, and c i for i ≥ 1 the upstream points on the cardinal flow path, defined by the condition that c i+1 is the biggest donor of c i . Let u be any donor of c 0 , which is not on the cardinal flow path (u = c 1 ). The point u is added to the swath and assigned to c 0 if the distance between u and c i is not greater than wi 2 for any i ≥ 0, where w i is the width according to Eq. (12 :: 13) applied to the ice flux of the point c i . After selecting the donors of c 0 that satisfy 250 this condition, the same procedure is applied to all donors of these points. The procedure continues until no more donors that satisfy the condition are found.
However, there is no straightforward definition of cardinal flow lines on an absolute scale since dendritic flow patterns may cover a wide range of scales. Therefore, the procedure described above is applied to all points c 0 with a nonzero ice flux.
Finally, the largest value of A i (resulting from different starting points c 0 ) is taken for each point.

255
An example of the obtained swaths and the respective fluxes is shown in Fig. 1 : 2. The topography used here is a fluvial equilibrium topography with K = 1, m = 0.5, and n = 1 under uniform uplift U = 1 on a grid of 5000 × 5000 nodes. This topography is shown in Fig. 2 : 3. It was already used by Hergarten (2020, Fig. 1a), where the only difference is a shift of the periodic eastern and western boundaries in such a way that the three biggest rivers do not cross the boundaries.
As expected, a strong downstream increase in valley width is the main difference towards the fluvial equilibrium topography.
Due to the U-shape of the big valleys, the parts of the tributaries that are captured by the swath are almost flat, which is 295 equivalent to lowering their base level compared to the fluvial topography. As a consequence, the surface heights of high regions decrease, although the same erodibility as for fluvial erosion was assumed.
Both the occurrence of wide, U-shaped valleys and the reduced elevation are also visible in the topographic profiles shown A simple, linear model for the part of the total precipitation that is converted into ice, 310 is used in this study, where H e is the ELA and H f is the altitude where the entire precipitation is converted into ice. The surface height H cannot be included easily in the fully implicit scheme for erosion, so it is evaluated at the beginning of the time step. This will, however, not strongly affect the stability of the implicit scheme and is a minor restriction of the maximum time increment compared to treating the flow directions in an explicit manner (Hergarten, 2020). A is already more then three times greater than A i for the considered flow paths. This means that the discharge of water is considerably higher than the ice flux for a major part of the flow path.
Assuming that all sites with A i > 0 are eroded glacially, while fluvial erosion only affects ice-free sites, is the simplest 320 concept of a coupled model. While this concept was used by Prasicek et al. (2020), the following example illustrates that it may cause problems when applied over long time spans.
The χ transform introduced by Perron and Royden (2013) provides a simple way to analyze longitudinal river profiles quantitatively. It transforms the upstream coordinate x to a new coordinate The χ transform eliminates the inherent concavity of river profiles arising from the upstream decrease in catchment size.
Equilibrium profiles under spatially uniform conditions turn into straight lines with a slope of U K 1 n . In our example (K = 1, The upper solid lines in Fig. 7 depict the χ-transformed profiles. The plot reveals that the profiles are even steeper than the 345 fluvial part on average almost up to the ELA. So the lower ice flux compared to the total discharge even shadows the higher glacial erodibility here. This effect could be compensated by increasing the glacial erodibility further, but the rather unrealistic steep fronts in the bedrock topography would persist.
Topography obtained under the same conditions as in Fig. 6, but assuming that erosion by meltwater follows the same relation as fluvial erosion with the same erodibility.

405
While the technical implementation would be simple if a model for the factor ω is given, it is not considered as an essential part of the model proposed in this study and is therefore not considered further.
Since the flow of meltwater should be confined to narrow channels, the best approach would be to consider the two sediment fluxes Q g and Q f as separate variables without any mixing. The numerical scheme proposed by Hergarten (2020) could in 435 principle be extended accordingly. However, this would require additional theoretical and numerical effort. Beyond this, the two sediment fluxes merge anyway, e.g., if material is deposited by the glacier due to decreasing transport capacity when approaching the glacier terminus and eroded again by the meltwater stream. Therefore, a simpler approach using a single sediment flux is suggested. According to Eqs. (30) and (31), the total erosion rate E = E g + E f follows the relation

440
where the definition of the effective erodibility for incision, is the same as in Eq. (27). If the total sediment flux Q = Q g + Q f is given, the fraction of the stream power spent for sediment transport depends on how Q is distributed. Let us assume an optimized distribution in the sense that this fraction is minimized.
It is easily recognized that this is the case if either the glacier or the meltwater carries the entire load, and thus This equation can be written in the same form as the original shared stream-power model (Eq. 3), with

450
Figures 9 and 10 :: 10 ::: and ::: 11 : provide a numerical example. Since sediment transport is more interesting in transient states than in a steady state, a declining ELA was chosen here. The ELA starts from H e = 16 (slightly higher than the maximum surface height of the fluvial equilibrium topography) and decreases at a rate of 4 (four times faster than the uplift), so that an ELA of H e = 8 (and H f = 10) is reached at t = 2.
The fluvial erodibilities are set to K d,f = 2.6 and K t,f = 1.625 ::: for ::::: rivers ::: and ::::::::: meltwater. Their ratio is Kd,f Kt,f = 1.6. This ratio is 455 the sediment deposition coefficient in the notation of Davy and Lague (2009), where a recent analysis of steady-state topographies suggested 1.6 as a realistic value for n = 1 (Guerit et al., 2019, data supplement). Furthermore, these values satisfy the with K f = 1, which ensures that the fluvial equilibrium topography computed for the detachment-limited model is also in equilibrium here (Hergarten, 2020). For glacial erosion, it is assumed that bedrock incision is not more efficient as fluvial bedrock incision (K d,g = K d,f = 2.6), while the glacier transports sediment 10 times more efficiently than water (K t,g = 10K t,f = 16.25).
While the topography ( Fig. 9 :: 10) is similar to the topography shown in Fig. 8 : 9, the longitudinal profiles ( Fig. 10 :: 11) show distinct differences to the detachment-limited version. First, the fluvial parts of the profiles have become steeper compared to 465 the initial fluvial topography. This increase in steepness is due to the high sediment flux from the glaciated part. As the erosion rate in the glacial part is higher than in equilibrium, the glacier brings more sediment than the upstream part of the river brought under fluvial conditions. This increased sediment flux reduces the ability to erode the bedrock according to Eq. (3), and the fluvial erosion cannot follow the uplift any more.
The glaciated parts of the profiles are less steep than the initial fluvial profiles, but still steeper than it would be expected in 470 equilibrium. Their slope in the χ plot is about 0.8, while Eq. (37) predicts an equilibrium slope of 1 Kd,g + 1 Kt,g = 0.45 with the glacial erodibilities assumed here. So the glaciated part is even steeper than the fluvial part if both are considered in relation to their equilibrium slope. The reason for the increased steepness is basically the same as for the fluvial part. Converting a Vshaped valley into a U-shaped valley goes along with a high erosion rate and thus yields a high amount of sediment. Although the efficiency of the glacier in transporting sediment was assumed to be ten times higher than for the rivers, the large amount 475 of sediments limits the ability of the glacier to erode the bedrock during the conversion of the valley shape.

Finite ice thickness
The approach developed in Sect. 2 considers the limit of zero ice thickness. The finite thickness, however, has a strong influence on glacial landform evolution. Overdeepened valleys would be neither possible in the stream-power model proposed here nor in the original shallow-ice approximation in the limit of zero thickness. Beyond this, the ice surface defines a base level for the 480 tributaries and may thus play an important part in the formation of hanging valleys.
A first estimate of the thickness h can be obtained by combining Eqs. (7) and (8) in the form and thus

485
In combination with Eqs :: Eq. (??)and (15) ::: 13), this yields Inserting the channel slope at the bedrock surface would lead to an unreasonable thickness h → ∞ at locations where S → 0.
Three effects of the finite thickness can be included in the model. First, the ice surface H +h can be taken into account in the 515 glacial mass balance (Eqs. 21 and 23), as it is usually done in glacial models. Since already H has to be taken at the beginning of the time step in Eq. (23), h is treated the same way without affecting the stability of the scheme strongly.
The most important effect of a finite ice thickness is that the channel slope S should refer to the ice surface H + h instead of the bedrock surface H . : as :::::::: assumed :: in ::: the ::::::: previous :::::::: examples. ::::: With ::: the :::::::::::::: parameterization :::::: h ∼ w, ::: the :::::: implicit ::::::: scheme ::: for ::: the :::::: erosion ::: part :::: can :: be :::::::::: maintained :::: with ::: the :::: slope :: of ::: the ::: ice ::::::: surface. Technically, this can be easily implemented by adding h to H 520 before performing the erosion step and subtracting it afterwards. Here the scheme for drawing the swath around the cardinal flow path described in Sect. 3 is particularly useful. It ensures that all points of the swath draining to a given point on the cardinal flow path, except for those located upstream on the cardinal path, have the same prolonged ::::::: extended : ice flux and thus also the same ice thickness. So the finite ice thickness has only a strong effect within the swath close to the glacier terminus While the slopes along the cardinal flow path are similar to the scenario with zero ice thickness, the glaciers advance further into the valleys here. This difference is related to taking into account the ice surface instead of the bedrock surface in Eq. (23), which defines the part of the precipitation that is converted into ice. This part has increased, which results in a higher ice flux and thus in an advancing of the glaciers.

535
The upper parts of the glaciated valleys feature smooth segments, which are interrupted by distinct steps. These steps occur at major confluences. Due to the abrupt increase in ice flux at confluences, the ice thickness also increases abruptly. Over long times, however, erosion will smoothen the ice surface, so that the step is transferred from the ice surface to the bedrock topography.
The overdeepening of the valley floor in the lower part of the glacier can be recognized well in the χ-transformed profiles.

540
The longitudinal shape of the bedrock surface looks quite irregular, with several short segments of strong overdeepening. This The third point where the model could be improved by the finite ice thickness is the separation of the total flow velocity into a deformation velocity v d and a sliding velocity v s . While the ice flux (Eq. 7) depends on the sum of both velocities, the erosion rate (Eq. 6) depends only on v s . As discussed in the beginning of Sect. 2, the ratio of the two velocities is proportional to h 2 .
In the 1-D model of Prasicek et al. (2018), this relation was used for eliminating v d from the equations. Equation (7) then turns with a parameter β depending on the rheology of the ice (for details, see Prasicek et al., 2018). The thickness h cannot be eliminated easily then. If we, however, assume that the last term in Eq. (41) is only a small correction, we can insert the estimate of h developed in this section there. Then all subsequent relations remain the same, but where :: the ::::::::::::: catchment-size 555 :::::::: equivalent : A i :: of ::: the ::: ice ::: flux : is replaced by with another parameter ξ. So the influence of the ice flux on erosion and sediment transport decreases if deformation becomes relevant. While Eq. (42) could easily be implemented, it introduces an additional parameter and will presumably not yield fundamentally different results. This extension is therefore not considered further in this study.

Numerical performance
The concepts developed in the previous sections were implemented in the open-source landform evolution model OpenLEM. A fully implicit scheme for the fluvial shared stream-power model was already available in OpenLEM before (Hergarten, 2020).
The behavior concerning the time increment δt is basically the same as for the fluvial version. Practically, the accuracy is not limited by the accuracy of the implicit scheme itself, but by changes in the flow pattern, which are treated in an explicit The time-dependent simulations with the declining ELA and the finite ice thickness (Sect. 6) were performed with δt = 10 −2 , δt = 10 −3 , and δt = 10 −4 . On a visual level (Figs. 11 and 12 :: 12 ::: and ::: 13), the results were almost indistinguishable for δt = 10 −3 570 and δt = 10 −4 , while a small difference was observed for δt = 10 −2 . However, the difference is rather a small shift on the time axis than a principal difference in the shape of the glaciers and the resulting landforms. So δt = 10 −3 should be a safe choice, while larger values should also be possible if required, e.g., in long-term simulations. If we, e.g., assume a fluvial erodibility of K = 2.5 Myr −1 (Robl et al., 2017), a unit of nondimensional time corresponds to 400,000 yr. So time increments of some hundred years appear to be safe, while some thousand years will also yield reasonable results.

575
However, the fluvio-glacial model requires a higher numerical effort per time step than the purely fluvial version. On the 5000 × 5000 grid, an increase in CPU time by a factor of about 1.7 was found for the fluvio-glacial shared stream-power model compared to the respective fluvial version even if the ice flux is zero everywhere. This factor increases to about 2.7 if the topography is completely glaciated. This factor, however, depends on the width of the glaciers. An increase of the factor of proportionality in Eq. (12) will result in an increasing numerical effort. In principle, the factors may also increase slightly for 580 larger grids than the 5000 × 5000 nodes considered here since the width of the largest glaciers may also increase then. In turn, the implementation in OpenLEM used here is still in a preliminary state and leaves room for further optimization.
So the increase in numerical effort compared to the purely fluvial model is moderate, and simulations even over several million years are possible on standard PCs.
Regarding the implementation in a large-scale landform evolution model, the main difference towards fluvial erosion is that glaciers are usually wider than the grid spacing in contrast to rivers. Including the finite width of glaciers is the main challenge in the implementation.
While the first model formulation assumes an infinitely thin ice layer, a finite thickness can also be included, where further 595 approximations are necessary. With this extension, overdeepenings, hanging valleys, and steps at confluences can be simulated in an at least qualitatively reasonable way.
The implementation in the open-source landform evolution model OpenLEM uses the fully implicit scheme for erosion and sediment transport proposed by Hergarten (2020) in the context of fluvial landform evolution. This scheme allows for arbitrary time increments in principle, where changes in the flow pattern practically define an upper limit. The numerical effort 600 is moderately higher than for the purely fluvial version and should be some orders of magnitude lower than for models based on the shallow-ice equations such as ICE-CASCADE and iSOSIA. Simulations even over several million years can be performed on standard PCs.
As a main limitation, the model presented here requires empirical relations for the width and the thickness of glaciers as a function of their ice flux. In contrast, models that implement the shallow-ice approximation directly are able to adjust the 605 geometry of the cross section according to the initial geometry of the valley, its slope, and the parameters of the ice flux. A detailed benchmarking against the iSOSIA model as a reference with regard to the efficiency and to the question how well the stream-power-based model captures glacial erosion will be subject of a subsequent study.
Code and data availability. The open-source landform evolution model OpenLEM including the extensions presented here is freely available at http://hergarten.at/openlem. An additional package that contains all codes and simulated data is available at http://hergarten.at/openlem/esurf-610 2021-1.zip (preliminary location during the review phase). The author is happy to assist interested readers in reproducing the results and performing subsequent research.

Author contributions. N/A
Competing interests. The author declares that there is no conflict of interest.