Simulating the evolution of the topography–climate coupled system

Landscape evolution models simulate the long-term variation of topography under given rainfall scenarios. In reality, local rainfall is largely affected by topography, implying that surface topography and local climate evolve together. Herein, we develop a numerical simulation model for the evolution of the topography–climate coupled system. We investigate how simulated topography and rain field vary between “no-feedback” and “co-evolution” simulations. Co-evolution simulations produced results significantly different from those of no-feedback simulations, as illustrated by transects and time evolution in rainfall excess among others. We show that the evolving system keeps climatic and geomorphic footprints in asymmetric transects and local relief. We investigate the roles of the wind speed and the time lags between hydrometeor formation and rainfall (called the delay time) in the co-evolution. While their combined effects were thought to be represented by the non-dimensional delay time, we demonstrate that the evolution of the coupled system can be more complicated than previously thought. The channel concavity on the windward side becomes lower as the imposed wind speed or the delay time grows. This tendency is explained with the effect of generated spatial rainfall distribution on the area–runoff relationship.

1 Introduction

Mutual influences between topography and climate have been widely recognized (e.g., Molnar and England, 1990; Masek et al., 1994; Willett, 1999; Roe et al., 2008). The spatial distribution of precipitation leads to spatial variability in surface runoff and erosion rates (e.g., Reiners et al., 2003; Moon et al., 2011; Bookhagen and Strecker, 2012; Ferrier et al., 2013), affecting landscape formation in the long term (e.g., Anders et al., 2008; Colberg and Anders, 2014). On the other hand, topography greatly affects precipitation distribution. Encountering upslope terrain on the path of moist air movement, wind raises an air parcel to a higher altitude, leading to cooling, condensation, and formation of hydrometeors (Smith, 1979). Such an orographic effect gives rise to the pattern in precipitation fields correlated with topography, which has been reported around the world (Puvaneswaran and Smithson, 1991; Park and Singh, 1996; Bookhagen and Burbank, 2006; Falvey and Garreaud, 2007, among others). The interaction between surface topography and precipitation field has a far-reaching implication; i.e., landscape and local climate evolve together as a single entity. Roe et al. (2008) suggested an analytical framework to describe the feedbacks among orographic precipitation, fluvial erosion, and critical wedge orogeny. Field-based evidences of the co-evolution have also been reported (e.g., Norton and Schlunegger, 2011; Champagnac et al., 2012). Nevertheless, the co-evolutionary dynamics and the resultant patterns in topography as well as local climate remain largely unexplored.

Mathematical modeling can be a sensible approach for the deeper investigation of feedbacks and controls in the evolutionary dynamics of the topography–climate coupled system. This requires a sophisticated numerical simulation model which links various surface processes and local rainfall generation processes. Since the 1990s, whole landscape evolution models have been popularly utilized for quantitative analysis and theoretical understanding of terrestrial processes (e.g., Willgoose et al., 1991a; Howard, 1994; Densmore et al., 1998; Zhang et al., 2016). Readers may refer to Temme et al. (2017) for a comprehensive review of mathematical models. Most early simulation studies assumed spatially uniform rainfall. Some studies have imposed spatially varying overland flow and evaluated its contribution to landform evolution (e.g., Pelletier, 2009; Wobus et al., 2010). Han et al. (2014) considered simple spatial rainfall distributions in simulating landform evolution of the island of Hawaíi, using the landscape evolution model CHILD (Tucker et al., 2001). Ward and Galewsky (2014) studied the impact of spatial rainfall distribution on stream incision using a 1D model. While these studies adopted spatial variability in rainfall or runoff, their spatial distribution was treated as invariant over time, ignoring topographic feedbacks on the precipitation field.

Co-evolution modeling with the consideration of topographic feedbacks was initiated by Masek et al. (1994), who coupled a simple orographic rainfall model with the cellular-automata surface process model of Chase (1992). Braun and Sambridge (1997) simulated the co-evolution of topography and local rain field with their landscape evolution model CASCADE. They adopted a simple orographic rainfall function where the rainfall amount becomes proportional to the surface elevation. Anders et al. (2008) incorporated the physically based orographic rainfall model of Smith and Barstad (2004) (referred to as the SB model in this paper) into the CASCADE model. Adopting a similar approach, Han et al. (2015) coupled the SB model with the CHILD model and investigated the imprint of orographic rainfall on the steady-state morphology, focusing on longitudinal profiles and river network organization. It was suggested that rainfall gradients perpendicular to a mountain range produce narrower catchments. Colberg and Anders (2014) adopted a simple Gaussian function for the spatial rainfall distribution wherein the local surface elevation is used as a variable. They ran co-evolution simulations by imposing the rainfall function into the CASCADE model and claimed that the location of precipitation maxima can determine whether the passive margin escarpment retreats and becomes gentler or steepens over time.

Despite aforementioned progresses, co-evolution modeling is a relatively new area, and many open questions remain. In this study, we aimed to advance our understanding of two fundamental subjects: (1) topographic changes that co-evolutionary dynamics bring in and (2) the respective role of meteorological variables in co-evolution. First, we explored how the relief and transect of topography evolve under the co-evolutionary setting. In this context, we investigated the feedbacks among local climate, surface processes, and the given tectonic uplift condition. Particularly, we questioned how the channel profile evolves under the co-evolutionary setting. The effect of orographic rainfall on the channel profile was investigated by Roe et al. (2002) for a 1D single corridor. For the river network within a landscape, Han et al. (2015) discussed the effect of orographic rainfall on channel concavity, as “high rainfall rates at the ridge top lead to mainstem channels that have relatively low concavity”. We sought a framework, which generalizes this description, from theoretical investigation in the present study.

Second, we investigated how meteorological variables control co-evolution. The most up-to-date strategy for simulating co-evolution can be found in the studies of Anders et al. (2008) and Han et al. (2015). In those studies, a rain field is continuously updated with evolving terrain by coupling the SB and landscape evolution models, and therefore, feedbacks between topography and climate were simulated. The SB model generates the rainfall field with many meteorological variables such as the wind vector v, the conversion time t_c from cloud water to hydrometeors, and the fallout time of hydrometeors t_f. The sum of the latter two quantities is referred to as the total delay time t_d (Smith and Evans, 2007), i.e., $t_{\mathrm{d}}=t_{\mathrm{c}}+t_{\mathrm{f}}$. A greater $\Vert \mathit{v}\Vert$ or t_d locates rainfall further downward. In this sense, their combined effects on an orographic rain field can be given as the non-dimensional delay time (or non-dimensional cloud drift time) (Barstad and Smith, 2005), i.e.,

$$\begin{array}{}\text{(1)}& t_{*}=\Vert \mathit{v}\Vert t_{\mathrm{d}}/L,\end{array}$$

where the length scale L is given as the mountain half-width. Considering the long-term effect of an orographic rain field on the topography formation, Anders et al. (2008) claimed that t_* “has a profound impact on topography (and precipitation patterns)”. However, we question whether t_* can represent the effect of local climate on long-term topography evolution. While both $\Vert \mathit{v}\Vert$ and t_d contribute horizontal displacement of the rain field, they may play oppositely in producing overall rainfall amount.

We aimed to address these questions by improving our understanding of the co-evolutionary system. We approached the proposed agenda with mathematical simulations. We investigated a hypothetical landscape and local climate converging to a quasi-steady state. The rest of this paper is organized as follows: we present the proposed co-evolution model structure in Sect. 2; numerical simulation results along with scientific insights in Sect. 3; in-depth discussions about implications of simulation results and current model limitations in Sect. 4; and conclusions in Sect. 5.

2 Coupling orographic rainfall generation and topography evolution

Modeling orographically induced precipitation is complex and can be attempted in various ways depending on the level of detailed physics incorporated (Barros and Lettenmaier, 1994). Among those, vertically integrated analytical models provide a reasonable balance between simulation accuracy and model complexity. In these models, the vertically averaged wind is expressed as the advection wind vector v:

$$\begin{array}{}\text{(2)}& \mathit{v}=U\mathit{i}+V\mathit{j},\end{array}$$

where U and V are wind speed components, and i and j are unit vectors along the x and y directions (we adopt the Cartesian coordinate system), respectively. Note that U or V is a mean speed, while the real wind exhibits a velocity profile often expressed as the logarithmic (Rossby, 1932) or power (Archibald, 1884) function of altitude. Therefore, the U or V value imposed in the SB model is considered greater than the ground wind speed component.

One of the well-known models of this type is the upslope model (e.g., Collier, 1975), which expresses cloud water S at a coordinate (x, y) in a raster domain as the inner product of v and the gradient ∇z of surface elevation field z (Smith and Barstad, 2004):

$$\begin{array}{}\text{(3)}& S(x,y)=C_{\mathrm{w}}\mathit{v}\cdot \mathrm{\nabla }z+S_{\mathrm{o}}.\end{array}$$

Note that both v and ∇z are vectors evaluated at the coordinate (x, y). Here, S_o is the background condensation rate, driven by the large-scale (synoptic) vertical wind component. If the orographic effect is nil (either v or ∇z is zero), S=S_o. In Eq. (3), C_w is the thermodynamic uplift sensitivity factor, given as

$$\begin{array}{}\text{(4)}& C_{\mathrm{w}}={\mathit{\rho }}_{\mathit{\nu }}{\mathrm{\Gamma }}_{\mathrm{m}}/\mathit{\gamma },\end{array}$$

where ρ_ν, Γ_m, and γ are the saturation vapor density, moist adiabatic lapse rate, and environmental lapse rate, respectively.

Assuming that condensed water immediately falls on the ground, the precipitation rate P=S (Eq. 3). Smith (2003) relaxed the assumption of immediate fallout by introducing the hydrometeor conversion time t_c and fallout time t_f. Accordingly, the Fourier transform of P(x,y) was presented as

$$\begin{array}{}\text{(5)}& \widehat{P}(k,l)={\displaystyle \frac{\widehat{S}(k,l)}{\left(\mathrm{1}+i\mathit{\omega }{t}_{\mathrm{c}}\right)\left(\mathrm{1}+i\mathit{\omega }{t}_{\mathrm{f}}\right)}},\end{array}$$

where $\widehat{S}(k,l)$, i, and ω are the Fourier transform of S(x,y), the imaginary unit ($\sqrt{-\mathrm{1}}$), and the intrinsic frequency, respectively. Here, k and l are the components of the horizontal wave number vector and $\mathit{\omega }=Uk+Vl$. Then, precipitation distribution can be obtained via its inverse transformation as

$$\begin{array}{}\text{(6)}& P(x,y)=\int \int \widehat{P}(k,l)e^{i(kx+ly)}dkdl+P_{\mathrm{o}}.\end{array}$$

Here, P_o is an optional term of the background precipitation rate, i.e., the rate without the presence of any orographic effect. In the absence of terrain, P_o=S_o.

Smith and Barstad (2004) improved the calculation of cloud water by considering airflow dynamics as

$$\begin{array}{}\text{(7)}& \widehat{S}(k,l)={\displaystyle \frac{C_{\mathrm{w}}i\mathit{\omega }\widehat{z}(k,l)}{\mathrm{1}-i{H}_{\mathrm{w}}m}},\end{array}$$

where $\widehat{z}(k,l)$ is the Fourier transform of the elevation field z, H_w is the moist layer height, and m is the vertical wave number. H_w is given as

$$\begin{array}{}\text{(8)}& H_{\mathrm{w}}=-{\displaystyle \frac{R_{\mathrm{v}}{T}^{\mathrm{2}}}{\mathrm{\Lambda }\mathit{\gamma }}},\end{array}$$

where R_v, T, and Λ are the gas constant for vapor (461 J kg⁻¹ K⁻¹), the near-surface air temperature, and the latent heat of vaporization (2501 kJ kg⁻¹ at 273.15 K), respectively. Instead of keeping Λ a constant, its variation with T is considered in the proposed model, using

$$\begin{array}{}\text{(9)}& \mathrm{\Lambda }=\mathrm{1918.46}{\left({\displaystyle \frac{T}{T-\mathrm{33.91}}}\right)}^{\mathrm{2}}\left(\mathrm{in}\mathrm{kJ}{\mathrm{kg}}^{-\mathrm{1}}\right),\end{array}$$

which was suggested by Henderson-Sellers (1984). In the hydrostatic limit ($N_{\mathrm{m}}^{\mathrm{2}}\gg {\mathit{\omega }}^{\mathrm{2}}$), m is given as

$$\begin{array}{}\text{(10)}& m(k,l)=N_{\mathrm{m}}\sqrt{k^{\mathrm{2}}+l^{\mathrm{2}}}\mathrm{sgn}\left(\mathit{\omega }\right)/\mathit{\omega },\end{array}$$

where the effective moist static stability (Fraser et al., 1973) is approximated as

$$\begin{array}{}\text{(11)}& N_{\mathrm{m}}=\sqrt{g\left(\mathit{\gamma }-{\mathrm{\Gamma }}_{\mathrm{m}}\right)/T}.\end{array}$$

Here, g is the gravitational acceleration. Considering the time lags (Eq. 5) and airflow dynamics (Eq. 7), a better estimation of the precipitation field is given (Smith and Barstad, 2004) as

$$\begin{array}{}\text{(12)}& \widehat{P}(k,l)={\displaystyle \frac{C_{\mathrm{w}}i\mathit{\omega }\widehat{z}(k,l)}{\left(\mathrm{1}-i{H}_{\mathrm{w}}m\right)\left(\mathrm{1}+i\mathit{\omega }{t}_{\mathrm{c}}\right)\left(\mathrm{1}+i\mathit{\omega }{t}_{\mathrm{f}}\right)}},\end{array}$$

which is the model mentioned earlier in this paper as the SB model. The orographic precipitation can be obtained from its inverse transformation (Eq. 6). Precipitation can be in various phases such as rain and snow. Our scope of precipitation phase is the rainfall at the time when the drop touches the ground.

The SB model is coupled with the surface process model LEGS (Landscape Evolution model using Global Search) (Paik, 2012). LEGS is a physically based mathematical model implemented on the same raster domain as the SB model. For this study, we revised and upgraded the source code of LEGS, adopting many improved features. One of them is the replacement of the GD8 (global eight-direction) method (Paik, 2008) with the improved GD8 method (Shin and Paik, 2017). LEGS adopted the GD8 method for surface flow path extraction. While this helps obtain reasonable results in terms of evolutionary speed and characteristics of simulated topography, some technical issues have been reported with the GD8 method, which are resolved by the improved GD8 method. The improved GD8 method requires the assignment of a start cell in a digital elevation model. This requires sorting cells according to their elevation values. The merge sort algorithm, invented by John von Neumann (see, e.g., Knuth, 1987), is adopted for this task.

In a real landscape, transport-limited and detachment-limited conditions are mixed, where the latter often appears upstream. In most theoretical modeling studies, however, the entire landscape has been assumed to be one of these conditions. Such a simple assumption eases the interpretation of simulation results by excluding higher-order complexities. The choice between the two conditions depends on the scope and purpose of each study. For example, Willgoose et al. (1991a) adopted the transport-limited condition. Earlier co-evolution studies of Anders et al. (2008) and Han et al. (2015) used the detachment-limited condition. The dominance between the two conditions can be related to timescales. In this study, the simulations ran for a long term (5 Ma) initiated from a flat topography. Temme et al. (2017) noted that “the assumption that there is always a large supply of transportable material is more likely true when considering very long (Ma) timescales – but on smaller timescales, temporary lack of transportable material can occur”. Following this argument, a transport-limited alluvial landscape is considered in this study.

P calculated from the SB model is used to generate the overland flow. There is no evapotranspiration or infiltration in our simulations. In this sense, we name P the rainfall excess rate. Then, the flow discharge per unit channel width q is calculated for every cell considering its upslope area A and P given in the area. The sediment transport rate per unit channel width q_s (mass/time/length) is calculated for each cell along a flow path as a function of q and the local slope $\Vert \mathrm{\nabla }z\Vert$, i.e.,

$$\begin{array}{}\text{(13)}& q_{\mathrm{s}}=C\Vert \mathrm{\nabla }z{\Vert }^{\mathit{\alpha }}{q}^{\mathit{\beta }},\end{array}$$

where C is an empirical coefficient. Equation (13) has been widely used in many landscape evolution studies (e.g., Smith and Bretherton, 1972; Willgoose et al., 1991a; Tucker and Slingerland, 1994; Braun and Sambridge, 1997), and its origin is described in Appendix A. Shallow landslide is considered in the simple way, previously applied by Tucker and Slingerland (1994); any slopes greater than the angle of repose Θ shall fail, translating upslope mass at the rate q_f to downhill until the local slope reaches Θ. This contribution is combined with the fluvial sediment transport to form the total sediment flux $\mathit{\xi }(=q_{\mathrm{s}}+q_{\mathrm{f}})$ at every cell. With the cell-to-cell mass transfer calculated above, the elevation field is updated over time t on the basis of the mass conservation along the flow path p, given as

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial z}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-\mathit{\lambda }){\mathit{\rho }}_{\mathrm{s}}}}{\displaystyle \frac{\partial \mathit{\xi }}{\partial p}}+u,\end{array}$$

where λ, ρ_s, and u are the porosity, density of a sediment particle, and tectonic uplift rate, respectively.

In the coupled model proposed in this study, the rain field generated by the SB model is fed into the LEGS model. Then, the updated surface topography obtained from the LEGS model becomes input to the SB model to generate the rain field for the next time step. The SB model involves the Fourier transform of the elevation field (Eq. 12) and the inverse transform of the rain field (Eq. 6). We implemented the full coupling between the LEGS and SB models; i.e., the topography and local rain field feed into each other at every iteration. To lighten any computation load in these operations, Kim (2014) adopted the fast Fourier transform method (Cooley and Tukey, 1965) and developed an efficient code, which we utilized in this study.

3 Applications

We applied the developed model to theoretical cases of mountain range formation, uplifting from an initially (almost) flat surface at the sea level. Simulations were implemented on a square domain composed of 200×200 cells, wherein each cell has an area of 0.5×0.5 km². In numerical simulations, spatial and temporal resolutions should be related to each other. For the selected calculation time interval (described below) and given simulation conditions, optimal Δx and Δy were found to be 500 m, which maximizes the numerical stability. In the given domain, 10 rows each at northern and southern boundaries are regarded as ocean. Hence, the effective terrestrial domain is composed of 200×180 cells. Eastern and western boundaries are regarded as extended mountain ranges wherein the periodic boundary condition is applied. This setting induces the development of a mountain range parallel to the coastline. In co-evolution simulations, wind is consistently blown from the north. We found that the open water surface (10 rows near boundaries) in front of the terrain on the course of airflow helps improve simulation results of the SB model and generate a smooth rain field over the terrestrial domain. Sediments eroded over the landscape, entering the surrounding ocean, are considered permanently lost; i.e., sea floor is assumed to be of an indefinite depth.

Whole landscape evolution modeling is typically implemented over geologic timescales. Accordingly, the calculation time interval is usually greater than a year. However, simulations should be implemented at the timescale of an individual storm event to consider orographic storms. As a storm event lasts at best a few days, to accommodate the effect of each storm event on erosion and deposition, the simulation might be implemented at a daily interval. Nevertheless, it is computationally exhausting to run the model at a daily time step over Ma timescales. To deal with these mismatching timescales, we adopted two concepts of (1) the dominant discharge and (2) dual time steps.

We first noted that not all storm events are geomorphologically effective. The streamflow corresponding to an event largely responsible for alluvial channel formation is referred to as the dominant discharge (Leopold et al., 1964). The dominant discharge is often considered similar to the effective discharge as well as the bankfull discharge and corresponds to the return period of between 1 and 2 years (Knighton, 1998). In accordance with this notion, we only impose a single representative storm event (referred to as the dominant storm) within the return period Δt₁ while taking Δt₁ as the primary calculation time interval (a similar approach with the bankfull discharge was shown in Byun and Kim, 2011). The flow generated by the given storm is considered the dominant discharge. We adopted Δt₁=2 years, which is close to the typical return period of a dominant storm reported in the literature.

The dominant discharge is of a significant amount, and so is ξ. If such a large amount is imposed in a single time step, we may encounter over-erosion and the formation of depression zones, the area which has elevation z lower than any neighboring cell. This occasion is considered an artifact due to the limited resolution of a numerical approach. Let us take an example of two adjacent cells where the upstream cell is to be eroded and the eroded material is supposed to be deposited on the downstream cell. As the erosion and deposition continue, the local slope between them gradually reduces. Because the sediment flux is a function of the local slope (Eq. 13), the sediment flux accordingly decreases and becomes nil if the gradient reaches zero. However, in the modeling, the local slope is invariant during the calculation time interval. Therefore, the sediment flux from the upstream cell is likely to be over-estimated. If the over-estimation is excessive, the elevation of the cell can become even below the downstream cell. To prevent or minimize over-erosion, the calculation time interval needs to be reduced. This can be done by dividing the storm duration into multiple time steps. In this respect, we adopted another time step of Δt₂, smaller than the duration of the dominant storm event. With this, Eq. (14) is discretized as

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }z=u\mathrm{\Delta }{t}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-\mathit{\lambda }){\mathit{\rho }}_{\mathrm{s}}}}\sum _{n=\mathrm{1}}^D\left({\mathit{\xi }}_{<\mathrm{in}>n}-{\mathit{\xi }}_{<\mathrm{out}>n}\right){\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{2}}}{\mathrm{\Delta }p}},\end{array}$$

where ${\mathit{\xi }}_{<\mathrm{in}>n}$ and ${\mathit{\xi }}_{<\mathrm{out}>n}$ are incoming and outgoing total sediment fluxes, respectively, at the time step n. Note that ${\mathit{\xi }}_{<\mathrm{in}>n}$ is the sum of fluxes originating from all neighboring cells contributing flow towards the given cell. Meanwhile, ${\mathit{\xi }}_{<\mathrm{out}>n}$ is calculated as the mass flux from the given cell towards the downstream cell. Δp is the distance between the given cell and its downstream cell. If the local flow direction is cardinal, $\mathrm{\Delta }p=\mathrm{\Delta }x=\mathrm{\Delta }y$, otherwise (diagonal flow direction), $\mathrm{\Delta }p=\sqrt{\mathrm{2}}\mathrm{\Delta }x$. In this study, we adopted Δt₂=1 d. The dominant storm event may last for D days, and Eq. (15) is solved for each day (D iterations per Δt₁). By combining the concepts of the dominant discharge and dual time intervals, we can implement virtually daily simulations while saving computational resources substantially.

For the purpose of this study, an ideal simulation target can be a high and active mountain range (such as the Andes and Himalaya) where orographic effects appear clearly. However, we considered the physical limit in the applicable peak altitude of the single-layer SB model. The peak altitude of the Olympic Mountains, USA, for which the SB model was evaluated in a previous study (Smith and Barstad, 2004), is about 2.4 km. Taking this as a reference, the tectonic uplift rate u=1 mm yr⁻¹ is used for all simulations, which yields a peak elevation of about 1.7 km (in most simulations) with a maximum of 2.7 km. In this way, we sought a reasonable compromise between the need to generate a high mountain range and the model limitation. This uplift rate is imposed uniformly over the domain and constantly throughout the entire simulation duration of 5 Ma. The tectonic response to unloading patterns induced by spatial variations in precipitation is ignored. The other model parameter values are summarized in Table 1.

Table 1 Model parameters commonly used in all simulations.

Most landscape evolution modeling studies adopt random perturbations to reproduce dendritic stream network formation. Randomness can be imposed in various media, such as the initial topography (e.g., Han et al., 2015) and surface material (e.g., Paik and Kumar, 2008). We adopted the former approach, and the random noise, following a uniform distribution (not spatially auto-correlated), is given in the initial topography. The initial topography is nearly flat with a very mild slope ($\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$) imposed to make sure surface water flows towards the ocean. The given random perturbation ($\pm \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m) is tiny enough relative to the given initial elevation gap between two nearby cells (10⁻² m), and so no depression zone forms in the initial topography. We implemented simulations initiated with the common setting described above but with varying model complexities. Two classes of simulations, namely “no feedback” and “co-evolution”, are presented below.

3.1 No-feedback simulations

No-feedback simulation results are obtained by running the LEGS model with spatially uniform rainfall, invariant over time. In no-feedback simulations, we focus on the effect of different rainfall excess rates P on long-term topography evolution. In the SB model, this can be seen as varying P_o in Eq. (6) with zero wind velocity, which will negate the orographic effect (i.e., P=P_o). We applied three scenarios of P=50, 100, and 150 mm d⁻¹. We assumed these storms continue for 4 d (D=4); hence, the total rainfall depth for each event is 200, 400, and 600 mm, respectively. Each event is considered to correspond to a dominant storm occurring every 2 years (Δt₁=2 years). With this Δt₁, 5 Ma simulation time is discretized into 2.5 million time steps. Actual numerical simulations are implemented at daily resolution, i.e., four iterations per every Δt₁. Hence, the computational load is 10 million iterations for 5 Ma simulation time in each scenario.

As the constant uplift is continuously imposed, the topographic relief increases over time in all cases, starting from zero (the sea level), and approaches a quasi-steady state (Fig. 1a). In this paper, the topographic relief or briefly relief refers to the gap between the peak elevation in the domain and the base level (which is zero as it refers to the sea level). As the greater P is imposed, the topography converges to the lower relief. This agrees with results from the earlier experimental study of Bonnet and Crave (2003).

Figure 1 Results of no-feedback simulations. Time variations of (a) the topographic relief and (b) the sediment yield from the entire domain are shown.

The sediment yield rapidly increases along with the topographic relief during the early stage of evolution (Fig. 1b). Rapid stream network formation at the early stage (Paik and Kumar, 2008) accelerates the sediment yield. As simulations proceed, the sediment yield converges to the rate of tectonic mass increase by the uplift (shown as a red horizon in Fig. 1b), regardless of P. This implies that the steady-state sediment yield is determined by the tectonic uplift rate, not rainfall. In contrast to the clear convergence of sediment yield (Fig. 1b), the convergence of relief is unclear as it exhibits continuous minor variations, after rapid early growth (Fig. 1a). This implies that our simulations produce quasi-steady, instead of purely steady, landscapes. The strict definition of the steady-state landscape requires no cell in the domain to show any elevation change over time. Referring to Eq. (14), this can be achieved if “the erosion rate is areally uniform and balances the rate of tectonic uplift” (Howard, 1994).

The balancing condition between the sediment yield and the uplift rate at the quasi-steady state helps explain distinct levels of quasi-steady reliefs with different P (Fig. 1a). To generate the same sediment yield at steady states, a high rainfall excess scenario is expected to pair with a low relief, as the sediment flux is given as the product of the runoff and gradient (Eq. 13). This argument is also in accordance with an earlier study about the detachment-limited condition (Whipple and Tucker, 1999). Note that the inverse relationship between the imposed P and the resultant relief is found to be nonlinear (as P linearly grows from 50 to 100 and 150 mm d⁻¹, relief reduces nonlinearly, approximately from 2.7 to 1.5 and 1 km). This is due to the nonlinearity in Eq. (13).

One may notice that, in the scenario of the greatest rainfall (P=150 mm d⁻¹), the sediment yield trend exhibits sharp fluctuations over time (Fig. 1b). The major mechanism causing these fluctuations is the formation and filling of depression zones. There is no depression zone initially, but they may appear during simulations as a result of over-erosion. Once formed, it only receives sediment inputs and generates no sediment output, and so the sediment yield from the domain drops. As a depression zone is filled with sediments, it disappears over time (Temme et al., 2006). To reduce over-erosion and depression zone formation, it would be needed to limit P and so ξ or to improve the time resolution, adopting a smaller Δt₂. Enhancing resolution will further increase the already large number of iterations or computational load. Here, we found that Δt₂=1 d is the resolution sufficiently fine to preclude over-erosion in all scenarios except the highest P(=150 mm d⁻¹) scenario. Orographic rainfall generated in our co-evolution simulations (to be shown below) is less than P=150 mm d⁻¹ as well, except some locally concentrated areas. Hence, we adopt Δt₂=1 d throughout this study to save the computational load. Note that episodic landslides along coastlines can lead to sudden increase in sediment yield, although these events are undetected in our simulations.

Figure 2 Schematic showing orographic rainfall over the 3D bird view of topography. Results after 5 Ma from co-evolution simulations with V=16 m s⁻¹ and t_d=1200 s are shown (transect shown in Fig. 3b). Horizontal displacement of raindrop from generation to falling location is given as t_d×V.

3.2 Co-evolution simulations

We evaluated the long-term consequences of considering the orographic rainfall effect by coupling the LEGS and SB models. In co-evolution simulations, the rainfall field spontaneously evolves over time and space. We followed basic simulation settings of no-feedback simulations: dominant storm events lasting for 4 d (D=4, Δt₂=1 d), considered per every 2 years (Δt₁=2 years), and uniform tectonic uplift rate u=1 mm yr⁻¹ imposed continuously for the entire simulation time of 5 Ma.

A number of simulation sets were designed with various atmospheric conditions. We evaluated the combined effects of the wind v, the delay times (t_c and t_f), and the temperature T on the rain field generation and long-term topography (Fig. 2). The wind direction is considered consistently north, i.e., U=0. As such, $\Vert \mathit{v}\Vert =V$ and a wide range of V were applied. Delay times are dependent on T. Anders et al. (2008) reported that $t_{\mathrm{d}}(=t_{\mathrm{c}}+t_{\mathrm{f}})$ (in seconds) is inversely related to T (in Kelvin), which can be approximated as

$$\begin{array}{}\text{(16)}& t_{\mathrm{d}}=\mathrm{58}\mathrm{000}-\mathrm{200}T.\end{array}$$

According to Eq. (16), three pairs of t_d and T are chosen for simulations: (1) t_d=600 s, T=287 K; (2) t_d=1200 s, T=284 K; and (3) t_d=1800 s, T=281 K. It has been known that the relative contribution of t_c or t_f in t_d is insignificant for the generation of orographic rainfall (Smith and Evans, 2007). Therefore, we simply set $t_{\mathrm{c}}=t_{\mathrm{f}}(=t_{\mathrm{d}}/\mathrm{2})$. For each pair, three levels of wind speed (V=8, 16, 24 m s⁻¹) were tested, yielding nine scenarios in total. We imposed the base level rainfall P_o=20 mm d⁻¹ in Eq. (6). Each of these nine scenarios is considered a dominant storm.

Finally, the environmental lapse rate γ was carefully chosen to consider the physical limitation of the SB model. With T, γ determines the moist layer height (Eq. 8). Because no vertical variability is allowed in the SB model, applying the SB model to a topography higher than the moist layer is questionable. As described earlier, the maximum topographic height within which the SB model has been applied is 2.4 km. Considering this, we set γ=−6.3 K km⁻¹, so that the range of the moist layer depth among nine scenarios is between 2.32 and 2.44 km. The moist adiabatic lapse rate Γ_m in Eq. (4) is given as −6.5 K km⁻¹ in our theoretical simulations, referring to examples of Smith and Barstad (2004).

Topography evolves in a dramatically different manner between no-feedback and co-evolution simulations, implying significant topography–climate feedbacks. For example, two simulations of “no feedback” with P=100 mm d⁻¹ and “co-evolution” with V=16 m s⁻¹ and t_d=1200 s are compared in Fig. 3. These are chosen in that rainfall depths accumulated over the total simulation period are similar. Transects from no-feedback simulations are symmetric, whereas co-evolution simulations result in asymmetric transects. Co-evolution simulations resulted in the peak being pushed toward the lee side. Such peak migration matches the postulation of Smith (2006), the physical experiments of Bonnet (2009), and earlier simulation studies with the detachment-limited condition (Anders et al., 2008; Han et al., 2015). This will be further discussed in the next section. The formation of asymmetric transects is directly related to the asymmetric rainfall distribution (Willett, 1999). The SB model generates the rainfall as a function of wind vector, local topographic gradient, delay times, temperature, and other meteorological variables. Accordingly, the rain fields generated here differ from those from simpler models which assume rainfall as a monotonic function of elevation; the generated rainfall exhibits sharp increase from the northern coast, where the moist air first confronts the slope, and then gradual reduction as the moist air penetrates deeper inland (Fig. 3b).

Figure 3 North–south transects of landscapes and rainfall fields after 5 Ma simulation times. Results from (a) no-feedback (P=100 mm d⁻¹) and (b) co-evolution (V=16 m s⁻¹, t_d=1200 s, $t_{*}=\mathrm{0.427}$) simulations are shown.

4 Discussion

4.1 Asymmetry and relief in evolving topography

The formation of asymmetric topography is related to the time evolution of the rain field. On the initial flat surface, rainfall is given uniformly as P_o without any orographic effect. Channel incision initiates at coasts (outlets) and then migrates upstream. This process should occur symmetrically across the initial transect. As the topography uplifts, the uniformity in the generated rain field is broken and the rain field evolves in an asymmetric fashion (Fig. 4a). The windward side receives more rainfall and so more runoff as well as sediment yield than the lee side. As a result, the channels on the windward side propagate upstream at a faster rate compared to those on the lee side (Fig. 4b). This leads to the downward peak migration and a steeper slope on the lee side (Fig. 4c).

Figure 4 Transects of evolving landscapes and rainfall field from a co-evolution simulation after (a) 0.2, (b) 0.6, and (c) 1 Ma of simulation times. Simulation is implemented with V=16 m s⁻¹ and t_d=1200 s. The results at 5 Ma are shown in Figs. 2, 3b, and 6.

As described above, the landscape experiences asymmetric erosion between the two sides (Fig. 5). On the other hand, the tectonic uplift is kept spatially uniform. Topography at each time step is determined through the trade-off between asymmetric erosion and uniform uplift. An interesting outcome of this complex adjustment is the development of the distinct peaks, which can be shown as a high ridge-valley relief at the lee side (Fig. 6a). These sharp peaks have been created by the interference between different evolutionary processes on the windward and lee sides. Whereas the windward side experiences high erosion due to enhanced rainfall during the evolution, many peaks remain intact on the lee side due to much less rainfall while being elevated (Fig. 4b). This results in a sharp contrast in elevation along the ridge line (Fig. 6a). Therefore, such features are a remnant of paleo-topography, similar to the monadnock.

Figure 5 Time variation of sediment yield from co-evolution simulation. Simulation is implemented with V=16 m s⁻¹ and t_d=1200 s. Red horizontal line corresponds to the mass balance with the tectonic uplift rate.

Figure 6 Contour map of simulated (a) surface topography and (b) rainfall field from the co-evolution simulation with V=16 m s⁻¹ and t_d=1200 s after 5 Ma (its transect shown in Fig. 3b). Numbers in the legend indicate the surface elevation (m) and rainfall excess rate (mm d⁻¹), respectively.

Catchments on windward and lee sides develop alternately (Fig. 6a). For a large catchment developed on the windward side, its channel head deeply penetrates beyond the domain center and toward the lee side (e.g., cross section A–A in Fig. 6a). The space between upper areas of two adjacent catchments from the windward side is taken as a niche where a major lee-side catchment develops. Along the valley line of a major catchment on the lee side (e.g., cross section B–B in Fig. 6a), the head of the lee-side channel reaches close to the domain center. Through the topographic feedback on the local climate, generated rainfall also exhibits alternate spatial pattern (Fig. 6b).

In their study with the detachment-limited condition, Roe et al. (2002) claimed that the orographic rainfall reduces relief. In their model, rainfall increases monotonically with elevation and reaches a maximum at the peak elevation, which leads to a large peak incision and hence lower relief. However, the real rainfall pattern induced by the orographic effect is often far from such monotonic increase with elevation (e.g., Bookhagen and Burbank, 2006). The SB model generates the rainfall distribution much more complexly and diversely, in comparison with Roe et al. (2002). As a result, our simulations generate relief comparable to, or even greater (higher peak elevation) than, the relief obtained from no-feedback simulations for similar accumulated rainfall depths (e.g., Fig. 3). Depending on local atmospheric and topographic conditions, rainfall concentrates on different altitudes, influencing relief in different ways.

4.2 Roles of wind speed and total delay time in co-evolution

In the co-evolution modeling, the domain-averaged rainfall amount increases over time from the initial rate of P_o (Fig. 7a and b). Such increase is in accordance with the relief growth and consequent feedback on rainfall. As shown in Eq. (3), the wind speed together with the topographic gradient is the primary factor that controls the rainfall. Essentially all rainfall depth generated beyond P_o is due to the orographic effect (Eq. 6), which becomes activated only if a moist air parcel is pushed by wind and climbs uphill. Clearly, orographic rainfall shall be nil when V=0, and the generated rainfall amount tends to increase as a greater wind speed is imposed (Fig. 7a).

Figure 7 Variations of (a, b) domain-averaged rainfall excess, (c, d) topographic relief, and (e, f) mean elevation during co-evolution simulations with different (a, c, e) V (t_d fixed as 1200 s) or (b, d, f) t_d (V fixed as 16 m s⁻¹).

For a greater wind speed V, rainfall increases at a faster rate and converges to a greater amount. Such a relationship between V and P, though, was found to be nonlinear. A physical upper limit appears to exist as V greater than a certain value contributes no further increase in P. With t_d=1200 s and T=284 K, V=16 m s⁻¹ seems to be the limit (Fig. 7a), although a greater rainfall is possible with a smaller t_d and warmer weather. The generation of P depends on both climatic (v) and topographic (∇z) factors (Eq. 3). Hence, to understand this nonlinear relationship and the limit in P, we need to understand the long-term effect of V on topography through rainfall and the topographic feedback on P, which is discussed below.

The resulting relief tends to be lower with a greater V (Fig. 7c). A similar tendency can be claimed for the mean elevation (Fig. 7e). These can be explained with the same reasoning described for no-feedback simulations (Fig. 1); i.e., greater wind speed is associated with more rainfall, and so it must be paired with a lower gradient to yield a sediment flux matching the mass generated at the imposed uplift rate. Interestingly, those two runs with V=16 and 24 m s⁻¹, while producing similar domain-average rainfall (Fig. 7a), show a relief difference (Fig. 7c). This can be explained by another mechanism: a greater wind speed carries rainfall further uphill, toward the lee side (compare Figs. 3b and 8a). This results in a greater erosion over the peak, which also contributes to lower relief.

Figure 8 Transects of landscapes and rainfall fields from co-evolution simulations after 5 Ma simulation times. Simulation settings are (a) V=8 m s⁻¹, t_d=1200 s and (b) V=16 m s⁻¹, t_d=600 s.

Despite the overall positive correlation between V and P (Fig. 7a), growing wind speed may act negatively on P. As discussed, with increasing V the relief tends to decrease. This gives a negative feedback on P. Hence, as V grows, there is a trade-off between the direct positive feedback on P and indirect (via a smaller relief) negative feedback on P. In two cases of V=16 and 24 m s⁻¹ shown in Fig. 7a, these effects have coincidently become identical. In detail, the run with V=24 m s⁻¹ has a greater positive feedback on P compared to the run with V=16 m s⁻¹. However, this advantage is fully compensated by its greater negative feedback on P due to the lower relief, resulting in the same P.

In contrast to the V versus P relationship, P decreases as t_d grows (Fig. 7b). This is because a part of hydrometeors generated within the domain falls out of the domain in proportion to t_d. For example, the scenario with V=24 m s⁻¹ and t_d=1800 s embeds 43 km horizontal displacement between locations of hydrometeor formation and rainfall (Fig. 2). Because a shorter delay time is related to a greater rainfall, one may also expect a robust relationship between t_d and the steady-state relief. Indeed, the robust relationship was found between t_d and the mean elevation (Fig. 7f). However, a similar relationship was not consistently obtained between t_d and relief (Fig. 7d). Although we find the inverse relationship between t_d and the rainfall amount (domain averaged) (Fig. 7b), a high rainfall induced by a short t_d concentrates downhill (compare Figs. 3b and 8b), which gives little effect on the peak elevation. Hence, no general relationship between t_d and relief can be established.

As discussed, both wind and t_d play primary roles in rain field generation. In fact, the combined effect of $\Vert \mathit{v}\Vert$ and t_d, rather than their individual effects, controls the co-evolutionary system. In this regard, the non-dimensional delay time t_* (Eq. 1) was claimed as an important measure in shaping the resultant topography (Anders et al., 2008). However, our co-evolution modeling shows that steady-state topography can substantially differ even for an identical t_*. For example, Fig. 8 illustrates two simulations with the same t_* of 0.213 but very different resultant topography. The dramatic difference is due to multiple effects of V and t_d. While they give a similar effect on the horizontal displacement of rainfall, they act oppositely in producing overall rainfall amount: increased V mostly induces greater rainfall (domain average), whereas increased t_d reduces rainfall (Fig. 7).

As described previously, the pure steady state can be claimed if every cell in the domain exhibits no elevation change over time. This would be referred to as the dynamic equilibrium. In a macro perspective, the dynamic equilibrium can be inferred when the relief becomes steady (Hack, 1960; Willgoose, 1994). In co-evolution simulations, generated rainfall amount converges to a constant, which was claimed as another indicator of a steady state (Han et al., 2015). In our simulations, the rainfall (Fig. 7a and b), relief (Fig. 7c and d), mean elevation (Fig. 7e and f), and sediment yield (Fig. 5) all exhibit converging tendency, implying that any of these can be a potential indicator of the dynamic equilibrium. However, the convergence times all differ among these criteria. Further, the relief shows no strict convergence but continuous minor variations. These imply that the pure steady states have not been achieved.

4.3 Profile concavity

Along a natural stream, the local slope $\Vert \mathrm{\nabla }z\Vert$ and drainage area A are associated as (Flint, 1974)

$$\begin{array}{}\text{(17)}& \Vert \mathrm{\nabla }z\Vert \propto A^{-\mathit{\theta }},\end{array}$$

where θ is the concavity index. For alluvial channels, θ is mostly found to be positive (Howard, 1980; Tucker and Whipple, 2002; Byun and Paik, 2017), which indicates a concave profile. No-feedback simulations have well-reproduced concave profiles (e.g., Willgoose et al., 1991b; Willgoose, 1994; Paik, 2012). Here, a basic question is raised about the role of orographic rainfall in the profile formation. Roe et al. (2002) reported the impact of orographic rainfall on channel concavity in the 1D context. Below, we discuss the formation of channel profiles in a more general 3D landscape context.

Channel profiles at quasi-steady states, simulated from our co-evolutionary model, are mostly concave. Interestingly, we found that the profile concavity on the windward side becomes lower as V or t_d increases (Fig. 9). We can explain such variation with the scaling relationship below (Willgoose et al., 1991b):

$$\begin{array}{}\text{(18)}& \mathit{\theta }={\mathit{\alpha }}^{-\mathrm{1}}[\mathit{\varphi }\mathit{\beta }-\mathrm{1}].\end{array}$$

In Eq. (18), ϕ is the exponent of the power-law relationship between the streamflow Q and drainage area A, i.e., Q∝A^ϕ. A more comprehensive equation, considering the hydraulic geometry (Leopold and Maddock, 1953), Hack's law (Hack, 1957), and downstream fining (Brush, 1961) relationships, was derived by Byun and Paik (2017). Here, we use a briefer Eq. (18) for a simpler explanation. With fixed sediment transport parameters α and β (used in Eq. 13), Eq. (18) implies that the concavity θ is proportional to ϕ. Considering that the flow within a catchment propagates through a stream network downstream, if rainfall concentrates downstream, high ϕ is expected (significant downstream increase in streamflow). As V or t_d increases, the rainfall distribution moves upstream, which will result in a lower ϕ. From Eq. (18), hence, a greater V or t_d can lead to a decreased concavity. If ϕ is excessively low (and hence ϕβ<1), even a convex profile may form.

Figure 9 Channel profiles after 5 Ma from simulations with (a) different V (t_d=1200 s) and (b) different t_d (V=16 m s⁻¹). The profile of the main channel of the largest catchment in each simulation is shown. For comparison, profiles are normalized with maximum elevation and lengths. Blue straight line is drawn to help visually identify the profile concavity.

4.4 Limitations

It is important to note the limitations of the SB model in reproducing complex rainfall patterns in mountainous terrain. A critical weakness, relevant to our modeling approach, is that the SB model considers only single values of T and V (i.e., spatially uniform) for the entire domain. Ideally, a 2D gridded or even 3D weather model can be coupled with the surface process model. However, running such a model for a geologic timescale currently faces a computational challenge (we ran 10 million iterations for each scenario).

We have demonstrated that generated rainfall decreases as the delay time grows (Fig. 7b). This accounts for missing hydrometeors generated inside but falling out of the domain. However, there are also hydrometeors formed out of the domain and blown into the domain. Such an amount is yet to be counted for more realistic simulations. To count such an effect, we need to deeply investigate to what spatial extent such an effect shall be considered.

The scope of the present study is the bare soil landscape. In reality, vegetation has played a significant role in shaping the Earth's terrestrial surface. Incorporating vegetation dynamics could make co-evolution results much more complicated, as their feedback mechanisms have a complex link with water, solar radiation, nutrient, carbon, sediment (soil), topography, etc. Modeling vegetation dynamics in the whole landscape context has gradually progressed in the last decade (e.g., McGrath et al., 2012; Yetemen et al., 2015). These efforts promise that vegetation dynamics can be eventually incorporated into the proposed co-evolution framework.

Finally, we face a basic question on whether simulation results are testable. Documented co-evolution records are absolutely lacking. While this is a fundamental issue in the surface process modeling community, we need cooperative efforts to address it. We indeed witness new lights such as thermochronology in topography reconstruction and paleo-climate simulations using general circulation models.

5 Conclusions

As the close linkage between topography and climate has been acknowledged, numerical modeling of their co-evolution has become an important subject. In this study, a sophisticated numerical model is developed for the given purpose, and results from “no-feedback” and “co-evolution” simulations are compared. Detailed simulation results provide insights into the feedbacks and controls in co-evolution. Major scientific contributions from our theoretical study can be summarized below.

• Evolving topography carries climatic and geomorphic footprints as shown in asymmetric transects and high valley relief on the lee side. Orographic rainfall affects resulting topographic relief in complex ways as spatial rainfall distribution greatly varies with meteorological conditions.

• Evolved topography under the same non-dimensional delay time t_* can significantly differ because of different roles of the wind speed and total delay time.

• Concavity in the channel profile on the windward side becomes lower as the wind speed or delay time increases. This is explained with the known scaling relationship between the streamflow and upslope area, which is associated with the spatial rainfall distribution, driven by wind and topographic gradients.

In most landscape evolution modeling studies, rainfall has been considered as an input. On the basis of results shown in the present study, we suggest that the modeling community reconsider this practice: we should incorporate local climate dynamics into landscape evolution. Probably, this is the time to revise the terminology of the whole landscape evolution modeling as modeling of the evolution of the topography–climate coupled system. It is desired to further explore the roles of model (particularly meteorological) parameters and their combinations. In this sense, simulation results presented here are examples calling for future research needs in co-evolution modeling.

Appendix A

The sediment transport equation (Eq. 13) adopted in this study follows the classical formula, widely applied over a century. The first literature on this functional form traces back to du Boys (1879), whose sediment transport equation follows

$$\begin{array}{}\text{(A1)}& q_{\mathrm{s}}\propto {\mathit{\tau }}_{\mathrm{b}}\left({\mathit{\tau }}_{\mathrm{b}}-{\mathit{\tau }}_{\mathrm{c}}\right),\end{array}$$

where τ_b is the bed shear stress and τ_c is to encounter the initiation of motion (similar to but not exactly the same as the critical shear stress).

The bed shear stress can be estimated as ${\mathit{\tau }}_{\mathrm{b}}=\mathit{\rho }gH\Vert \mathrm{\nabla }z\Vert$, where ρ is the density of the fluid and H is the flow depth. In analogy, we can suggest ${\mathit{\tau }}_{\mathrm{c}}=\mathit{\rho }g{H}_{\mathrm{c}}\Vert \mathrm{\nabla }z\Vert$, where H_c is the threshold depth corresponding to the initiation of the motion for the given slope $\Vert \mathrm{\nabla }z\Vert$. Substituting these into Eq. (A1) gives

$$\begin{array}{}\text{(A2)}& q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{2}}\left(H^{\mathrm{2}}-H_{\mathrm{c}}H\right).\end{array}$$

If we take the Chezy formula, $H\propto \Vert \mathrm{\nabla }z{\Vert }^{-\mathrm{1}/\mathrm{3}}{q}^{\mathrm{2}/\mathrm{3}}$. Similarly, $H_{\mathrm{c}}\propto \Vert \mathrm{\nabla }z{\Vert }^{-\mathrm{1}/\mathrm{3}}{q}_{\mathrm{c}}^{\mathrm{2}/\mathrm{3}}$, where q_c is the threshold value of q corresponding to H_c and τ_c. Then, Eq. (A2) can be rewritten as

$$\begin{array}{}\text{(A3)}& q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{1.33}}\left(q^{\mathrm{1.33}}-q^{\mathrm{0.67}}{q}_{\mathrm{c}}^{\mathrm{0.67}}\right).\end{array}$$

Instead of Chezy, the Manning formula can also be used, i.e., $H\propto \Vert \mathrm{\nabla }z{\Vert }^{-\mathrm{0.3}}{q}^{\mathrm{0.6}}$ and $H_{\mathrm{c}}\propto \Vert \mathrm{\nabla }z{\Vert }^{-\mathrm{0.3}}{q}_{\mathrm{c}}^{\mathrm{0.6}}$. Substituting these into Eq. (A2) gives

$$\begin{array}{}\text{(A4)}& q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{1.4}}\left(q^{\mathrm{1.2}}-q^{\mathrm{0.6}}{q}_{\mathrm{c}}^{\mathrm{0.6}}\right).\end{array}$$

Equation (A2) states q_s as a function of H. By transforming Eq. (A2) into Eq. (A3) or Eq. (A4), q_s is expressed as a function of q. This is practically advantageous in that the evaluation of flow discharge is easier than the flow depth, particularly in the context of landscape evolution modeling.

Austrian engineer Schoklitsch implemented laboratory experiments to evaluate this type of equations with the grain size between 1 and 2 mm. The equation of Schoklitsch (1934) published in 1934 (introduced in English by Shulits, 1935) is

$$\begin{array}{}\text{(A5)}& q_{\mathrm{s}}=C_{\mathrm{o}}{d}^{-\mathrm{0.5}}\Vert \mathrm{\nabla }z{\Vert }^{\mathrm{1.5}}\left(q-q_{\mathrm{o}}\right),\end{array}$$

where q_o is the threshold value of q for the initiation of motion and d is the grain size. The coefficient C_o has a complex dimension. In English units (q_s in lb s⁻¹ ft⁻¹, d in inches, q in ft³ s⁻¹ ft⁻¹), C_o=86.7 lb in^0.5 ft⁻³ was suggested from flume experiments (see Shulits, 1935), which is about 7000 kg mm^0.5 m⁻³ in metric units (q_s in kg s⁻¹ m⁻¹, d in millimeters, q in m³ s⁻¹ m⁻¹).

The above Eqs. (A3)–(A5), if the threshold term is neglected, commonly express the sediment transport rate for non-cohesive bed streams as the function of local slope and flow discharge, suggesting the form of Eq. (13). Quantifying the initiation of motion (the threshold term) is difficult, especially for small grains such as sands (Byun and Paik, 2017). Movement of sand particles has been observed even at very weak flow in flume experiments (e.g., Shvidchenko and Pender, 2000). Considering such a large uncertainty associated with it, the threshold term is often neglected when applying these equations to the scale of landscape evolution (Willgoose et al., 1991a; Tucker and Slingerland, 1994; Dietrich et al., 2003).

Other sediment transport equations such as Shields, Einstein–Brown, and Kalinske can also be expressed in the form of Eq. (13) (Raudkivi, 1967). Taking the Chezy or Manning formula, the Shields equation can be expressed as

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{1.67}}{q}^{\mathrm{1.67}}\left(\mathrm{using}\mathrm{Chezy}\right)\mathrm{or}\\ & q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{1.7}}{q}^{\mathrm{1.6}}\left(\mathrm{using}\mathrm{Manning}\right).\end{array}\end{array}$$

Derivation from the Einstein–Brown equation is given in Henderson (1966) and later by Willgoose et al. (1991a) in more detail as

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{2}}{q}^{\mathrm{2}}\left(\mathrm{using}\mathrm{Chezy}\right)\mathrm{or}\\ & q_{\mathrm{s}}\propto \Vert \mathrm{\nabla }z{\Vert }^{\mathrm{2.1}}{q}^{\mathrm{1.8}}\left(\mathrm{using}\mathrm{Manning}\right).\end{array}\end{array}$$

As shown in this Appendix, Eq. (13) can be derived from many existing sediment transport formulae. It has also been directly supported by flume experiments. While all these approaches result in the form of Eq. (13), they differ in the values of exponents α and β. For example, on the basis of Eq. (A3), $\mathit{\alpha }=\mathit{\beta }=\mathrm{1.33}$, while Eq. (A4) implies α=1.4 and β=1.2. Even greater values are shown in Eq. (A6) or Eq. (A7). Another experimental study by MacDougall (1933) is also noteworthy, where the same functional form of Eq. (13) was suggested with a range of α between 1.25 and 2, while β is kept as unity. From the literature reviewed above, we can see that the ranges of α and β are $\mathrm{1}<\mathit{\alpha }<\mathrm{2.1}$ and $\mathrm{1}<\mathit{\beta }<\mathrm{2}$. Details of these ranges have been studied by Prosser and Rustomji (2000). Comparing Eqs. (13) and (A5), the coefficient C in Eq. (13) is interpreted as $C=C_{\mathrm{o}}{d}^{-\mathrm{0.5}}$ given α=1.5 and β=1. As α and β vary, the value and dimension of C in Eq. (13) are subject to variation.