By Milad Hooshyar, Shashank Kumar Anand, and Amilcare Porporato
Earth’s surface evolves in response to hydro-climatic and tectonic forcings. Modeling the formation and evolution of topographic features, such as channels and hillslopes, is of great importance in hydrology and geomorphology due to their direct effect on sediment and water transport, as well as ecosystem regulation. Quantifying the response of the land surface elevation to climate change and human activities is key to preventing or reversing disastrous events like land degradation in South Carolina’s Calhoun Critical Zone, where intensive cotton cultivation has led to high erosion and soil loss [3].
Precipitation over the land surface either infiltrates through the soil, evaporates, or flows on the surface. The latter erodes the top layer of soil by mobilizing soil particles. Movement of soil particles is still possible in the absence of water flow through what is called “soil creep.” At the same time, tectonic activities and weathering processes at the top of the bedrock produce soil, which—in the steady-state condition—should be balanced by erosion and creep.
The processes involved in landscape evolution—from precipitation to surface-flow generation, erosion, and tectonic activities—are very complicated and in some cases, not fully understood. The study of such a complex system is possible through investigation of natural basins, experimental simulations, and mathematical modeling. The data from natural basins has revealed several universal laws in river networks concerning properties such as drainage density, junction angles, and channel network structure [8]. In some cases, interpreting the relationships between external forcings (e.g., climate) and properties of the natural basin is difficult due to time variability of the forcings and the complex feedback between subsystems like climate, vegetation, and erosion. Alternatively, we can use a minimalist system of partial differential equations to model the dynamics of the land surface elevation \(h\),
\[ \frac {\partial h}{\partial t}=D \nabla^2h-Ka^m|\nabla h|^n+U, \tag1\]
where \(D\) is the diffusion coefficient, \(K\) is an erosion coefficient, \(m\) and \(n\) are model parameters, \(U\) is the uplift rate, and \(a\) is the specific drainage area.
Figure 2. Schematic definition of the specific drainage area \(a\). Here, \(A\) is the area confined between two streamlines and a contourline with length \(w\). \(a\) is the ratio of \(A\) to \(w\) as \(w\) approaches zero.
Equation \((1)\) describes the temporal change in surface elevation as a result of the processes depicted on the right-hand side (RHS) of the equation. The landscape reshapes as a combined effect of the diffusive soil-creep (first term on the RHS), erosion by the surface runoff (second term on the RHS), and tectonic uplift (third term on the RHS).
The schematic definition of \(a\) is presented in Figure 2, which is the horizontal projection of the upstream length that passes through a point. Mathematically, \(a\) is linked to \(h\) [2, 5],
\[ \nabla \cdot \left( a \frac{\nabla h}{|\nabla h|} \right) +1 = 0. \tag 2 \]
In reference [4], equations \((1)\) and \((2)\) are non-dimensionalized to obtain the "Channelization Index" \((\mathcal{C_I})\), thus capturing the intrinsic property of the systems with only one dominant length scale (e.g., squared or semi-infinite domain), as
\[ \mathcal{C_I} = \frac {Kl^{m+n}}{D^nU^{1-n}}, \tag3 \]
where \(I\) is the domain's length scale. In numerical landscape evolution models (LEMs), instead of explicitly solving equation \((2)\) over the entire domain, a flow-direction (e.g., \(\textrm{D}8\) and \(\textrm{D}\infty\)) method approximates the specific drainage area based on the given elevation field [7, 9]. Equation \((1)\) is then solved to update the elevation using a numerical scheme until the steady-state is reached.
Figure 3. Simulation results for the square domain (side length \(=I=100\textrm{m}\)) for \(\mathcal{C_I} = 20\), with smooth landscape profile (3a) and \(\mathcal{C_I} = 200\) (3b) for the dissected landscape.
The nonlinear term \((\)the erosion term in equation \((1))\) increases the computational cost of any numerical solver. Especially for an implicit solver—where the current state of the system is used to update the elevation field—solving the coupled nonlinear equations is computationally expensive. To this aim, we use an efficient implicit solver that provides accurate steady-state solutions for wide values of \(\mathcal{C_I}\) [1].
Numerical simulations of the coupled system of equations \((1)\) and \((2)\) show that as parameter \(\mathcal{C_I}\) increases (e.g., higher erosion), the steady-state solution transitions from a smooth and unchannelized surface towards a dissected one with branching channels (see Figure 5). The linear stability analysis for a semi-infinite domain with \(m=n=1\) predicts the transition to occur at \(\mathcal{C_I} \approx 37\) [4].
Figure 4. Numerical solution for the rectangular domain (width = 100m, length = 500m) with model parameters \(m=n= 1.0\) and \(\mathcal{C_I}=500\), starting from the flat surface until the topographic equilibrium is reached.
In the dissected surfaces, the branches follow a hierarchical pattern where small channels connect to form bigger ones. Overall, this behavior resembles the control of the Reynolds number \((Re)\) on the flow condition in turbulence, where increasing \(Re\) first initiates the instability. Furthermore, a higher \(Re\) creates a hierarchy of vortices in which bigger vortices transfer energy to smaller ones until flow viscosity dissipates the energy [4].
Figure 5. Increasing \(\mathcal{C_I}\) results in progressively dissected surfaces with branching channels. The colors in 5a-5c show the value of specific catchment area \(a\), where higher values correspond to channels. The lowering of the mean elevation profile, computed by averaging elevation \(h\) along the \(x\)-axis, and the formation of the logarithmic region with increasing \(\mathcal{C_I}\) are depicted in 5d and 5e. The mean elevation and distance from the boundary are appropriately non-dimensionalized in 5e with average elevation at \(y = 50\textrm{m}\), denoted by \(\bar{h}_{max}\) and the parameters in equation \((1)\). The logarithmic region corresponds to the linear segments in Figure 5e’s semi-log space.
Increasing the value of \(\mathcal{C_I}\) also results in a lower mean-elevation profile, where the mean-elevation profile is simply computed by averaging the elevation field along the \(x\)-axis in Figure 5. In the asymptotic condition when the value of \(\mathcal{C_I}\) is relatively high, the emerging mean-elevation profile of the surface contains a logarithmic scaling at the intermediate region of the boundary, very similar to the universal logarithmic scaling in wall-bounded turbulence. Such a logarithmic region hints to a common behavior of bounded complex systems, where the propagation of the patterns is restricted by the system boundary [6].
Amilcare Porporato presented this work during a minisymposium presentation at the 2019 SIAM Conference on Applications of Dynamical Systems, which took place in May in Snowbird, Utah.
References
[1] Anand, S.K., Hooshyar, M., & Porporato, A. (2019). Linear layout of multiple flow-direction networks for landscape-evolution simulations. Preprint, arXiv:1909.03176.
[2] Bonetti, S., Bragg, A.D., & Porporato, A. (2018). On the theory of drainage area for regular and non-regular points. Proceed. Royal Soc. A: Math., Phys. Eng. Sci., 474(2211), 20170693.
[3] Bonetti, S., Richter, D.D., & Porporato, A. (2019). The effect of accelerated soil erosion on hillslope morphology. Earth Surf. Process. Land.
[4] Bonetti, S., Hooshyar, M., & Camporeale, C., & Porporato, A. (2019). Channelization cascade. Preprint, arXiv:1812.03696.
[5] Gallant, J.C., & Hutchinson, M.F. (2011). A differential equation for specific catchment area. Wat. Res. Research, 47(5).
[6] Hooshyar, M., Bonetti, S., Singh, A., Foufoula-Georgiou, E., & Porporato, A. (2019). From turbulence to landscapes: universality of logarithmic mean profiles in bounded complex systems. Preprint, arXiv:1909.13367.
[7] O’Callaghan, J.F., & Mark, D.M. (1984). The extraction of drainage networks from digital elevation data. Comp. Vis. Graph. Im. Process., 28(3), 323-344.
[8] Rodriguez-Iturbe, I., & Rinaldo, A. (2001). Fractal River Basins: Chance and Self-organization. New York, NY: Cambridge University Press.
[9] Tarboton, D.G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Wat. Res. Research, 33(2), 309-319.
 |
Milad Hooshyar is a postdoctoral fellow in the Department of Civil and Environmental Engineering, the Princeton Institute for International and Regional Studies, and the Princeton Environmental Institute at Princeton University. |
 |
Shashank Kumar Anand is a graduate student in the Department of Civil and Environmental Engineering at Princeton University |
 |
Amilcare Porporato is Thomas J. Wu '94 Professor in the Department of Civil and Environmental Engineering and the Princeton Environmental Institute at Princeton University. |