The linear feedback precipitation model (LFPM 1.0) – a simple and efficient model for orographic precipitation in the context of landform evolution modeling
- 1Institut für Geo- und Umweltnaturwissenschaften, Albert-Ludwigs-Universität Freiburg, Freiburg, Germany
- 2Fachbereich Geographie und Geologie, Paris Lodron Universität Salzburg, Salzburg, Austria
Correspondence: Stefan Hergarten (firstname.lastname@example.org)
The influence of climate on landform evolution has attracted great interest over the past decades. While many studies aim at determining erosion rates or parameters of erosion models, feedbacks between tectonics, climate, and landform evolution have been discussed but addressed quantitatively only in a few modeling studies. One of the problems in this field is that coupling a large-scale landform evolution model with a regional climate model would dramatically increase the theoretical and numerical complexity. Only a few simple models have been made available so far that allow efficient numerical coupling between topography-controlled precipitation and erosion. This paper fills this gap by introducing a quite simple approach involving two vertically integrated moisture components (vapor and cloud water). The interaction between the two components is linear and depends on altitude. This model structure is in principle the simplest approach that is able to predict both orographic precipitation at small scales and a large-scale decrease in precipitation over continental areas without introducing additional assumptions. Even in combination with transversal dispersion and elevation-dependent evapotranspiration, the model is of linear time complexity and increases the computing effort of efficient large-scale landform evolution models only moderately. Simple numerical experiments applying such a coupled landform evolution model show the strong impact of spatial precipitation gradients on mountain range geometry including steepness and peak elevation, position of the principal drainage divide, and drainage network properties.
The redistribution of moisture from the oceans towards continental domains governs the global erosion engine. Spatial variability in precipitation and hence in the availability of water or ice as principal agents of erosion controls the shape of landforms (e.g., Ellis et al., 1999; Willett, 1999; Anders et al., 2008; Bonnet, 2009; Menking et al., 2013; Colberg and Anders, 2014; Goren et al., 2014; Chen et al., 2019; Han et al., 2015; Paik and Kim, 2021). However, feedbacks between topography, precipitation, and erosion may even make it difficult to distinguish between cause and effect (Molnar and England, 1990).
Long-term fluvial erosion is a field in which simple numerical models have been applied with great success for some decades. The simplest model in this context is often referred to as the stream-power incision model (SPIM) and is the key component of several models of long-term fluvial landform evolution (for an overview, see, e.g., Willgoose, 2005; Wobus et al., 2006). The SPIM considers rivers to be linear elements (so without explicitly accounting for the width and the cross-sectional shape) and predicts the erosion rate E as a function of the upstream catchment size A and the channel slope S in the form
While the exponents m and n are kept constant, all site-specific influences on erosion are subsumed in a single lumped parameter K, called erodibility. The SPIM implements the concept of detachment-limited erosion (Howard, 1994) in the sense that all particles entrained by the river are immediately swept out of the system. The applicability of this concept even to bedrock rivers in high mountain regions has been questioned (e.g., Turowski, 2012). However, several extensions of the SPIM by sediment transport were proposed (e.g., Whipple and Tucker, 2002; Davy and Lague, 2009; Hergarten, 2020), with efficient numerical schemes having recently become available (Yuan et al., 2019; Hergarten, 2020). Even extensions towards glacial erosion were recently proposed (Deal and Prasicek, 2021; Hergarten, 2021a).
The SPIM and its derivates are well-suited for problems of tectonic geomorphology, e.g., variations in uplift rate or contrasts in lithology. In turn, the occurrence of the erodibility as a single, lumped parameter is a serious limitation concerning the influence of precipitation.
A framework for extending the SPIM to spatial variation in precipitation was already applied in several studies (e.g., Yanites and Ehlers, 2012; Goren et al., 2014; Garcia-Castellanos and Jiménez-Munt, 2015; Salles, 2016; Yuan et al., 2019). The idea is that erosion rates should depend on discharge rather than on catchment size, although the SPIM (Eq. 1) is typically written in terms of catchment size for historical reasons. Let P be the effective precipitation for the moment, i.e., the part of the precipitation that contributes to discharge. If we assume that the actual erodibility K refers to a given uniform reference precipitation P0 and thus to a reference discharge q0=P0A, the catchment size A can be replaced by in Eq. (1). It is then assumed that Eq. (1) holds for any discharge q if A is replaced by . Since the equations become somewhat cumbersome if q is replaced by the integral of P over the upstream catchment, Hergarten (2021a) defined
as the catchment-size equivalent of the discharge. It defines the catchment size needed to generate the actual discharge q at the reference precipitation P0. The advantage of using this terminology is that all relations in the context of erosion keep their simplicity, just with Aeq instead of A.
One might think of using scenarios of a regional climate model for computing precipitation, e.g., the Weather Research and Forecasting (WRF) model (Skamarock et al., 2021). The recent version of this model can in principle be run on PCs at spatial resolutions of a few kilometers, which allows for a consideration of orographic effects at the catchment scale. However, there would still be a huge imbalance between the complexity of the precipitation model and the simplicity of the erosion model. This imbalance would not only concern the computing effort, but also the level of detail of the prediction. While we might think of long-term mean precipitation rates in the extension of the SPIM by nonuniform precipitation described above, individual large rainstorms and related floods contribute much to landform evolution in reality. So even the question of whether an ensemble average or rather some kind of maximum over scenarios of a regional climate model yields a better input for landform evolution modeling is nontrivial.
Preserving the simplicity of the SPIM and its derivates requires simple models focusing on orographic effects on the relative precipitation , which allows for computing Aeq (Eq. 2). The main challenge is finding a level of complexity much below that of regional climate models that still provides new insights into landform evolution. On a qualitative level, reproducing an increased precipitation rate at the windward side of orogens and a rain shadow behind mountains would be some minimum requirement. Taking into account scales larger than the width of individual orogens, it may also be desirable to reproduce the overall decrease in precipitation with increasing distance from the reservoir of moisture (typically an ocean).
On a fundamental level, even extremely simple approaches have been used. Goren et al. (2014) distinguished between the windward side and the leeward side of a mountain belt just by the main drainage divide and assigned increased relative precipitation to the windward region. This extremely simple model turned out to be sufficient for explaining a shift and an asymmetry in the drainage divide.
In turn, the models proposed by Roe et al. (2003), Smith and Barstad (2004), and Garcia-Castellanos (2007) use the concept of vertically integrated water contents and the respective fluxes per unit width. Assuming steady-state conditions, precipitation is derived from the negative divergence of the flux per unit width. All these models bring the topography into play by a thermodynamic equilibrium that depends on altitude via temperature.
The earliest among these models (Roe et al., 2003), however, does not model any fluxes explicitly but directly proposes an equation for the divergence and thus for precipitation. The model predicts the rate of precipitation explicitly as a function of local surface elevation and slope in the wind direction. As the only nonlocal component of the model, a Gaussian smoothing in the upwind direction was used in order to reduce effects of surface roughness. Due to these properties, the model is able to reproduce increased precipitation at the windward side compared to the leeward side of a mountain belt, but it fails to describe the large-scale shadow in a plane behind the mountain range or the decrease in precipitation with increasing distance to the ocean.
The two other models (Smith and Barstad, 2004; Garcia-Castellanos, 2007) consider spatially variable water contents and the respective fluxes, with transport at a given wind velocity assumed. The model proposed by Smith and Barstad (2004) defines two components, interpreted as cloud water and hydrometeors. This model focuses on condensation and fallout at small scales, while it cannot predict transport over long distances (see Sect. 6). It therefore requires a refilling from an additional reservoir and is, similarly to the model of Roe et al. (2003), not able to predict large-scale precipitation patterns. In turn, the model of Garcia-Castellanos (2007) describes the vertically integrated water content by a single variable. Using a quite ingenious approach for describing deviations from equilibrium, it is able to capture the increase in precipitation with elevation as well the slow decrease in precipitation with increasing distance from the ocean. In turn, it requires an artificial smoothing at small scales, similarly to the model of Roe et al. (2003).
It seems that the model of Smith and Barstad (2004) (SB model in the following) received the most attention from the landform evolution modeling community among these models. It was adopted by some other authors in the context of co-evolution of topography and climate (e.g., Anders et al., 2008; Han et al., 2015; Paik and Kim, 2021), although the model of Garcia-Castellanos (2007) has some advantages (see Sect. 6). The main advantage of the SB model seems to be that it can be implemented numerically on a regular grid using a forward and backward Fourier transform without the need to carry additional variables and to think about numerical stability and efficiency.
The goal of this study is developing a model that captures both the direct response of precipitation to changes in topography and large-scale precipitation patterns without the need for ad hoc assumptions such as an additional reservoir or smoothing. Beyond this, the numerical complexity should be not much higher than in the existing models. In particular, the linear time complexity (i.e., that the computing effort increases only linearly with the grid size) achieved by contemporary fluvial landform evolution models (Hergarten and Neugebauer, 2001; Braun and Willett, 2013; Yuan et al., 2019; Hergarten, 2020) should be preserved.
The model developed in the following is inspired by the concepts of Smith and Barstad (2004) and Garcia-Castellanos (2007). Similarly to these models, we describe the distribution of water in the atmosphere in terms of vertically integrated water contents measured in meters, which can be interpreted as water column heights. Following the ideas of Smith and Barstad (2004), we use two components, while the model of Garcia-Castellanos (2007) uses a single component and thus seems to be simpler at first sight. However, we will see in Sect. 4.1 that the effort of using two components pays off.
2.1 The governing equations
Let Qv be the content of vapor and Qc be the content of cloud water, both vertically integrated and measured as the height of a water column. Following the concepts of Smith and Barstad (2004) and Garcia-Castellanos (2007), we assume that advection with a given velocity is the predominant transport mechanism. If is the respective velocity of advection, the advective flux per unit width is
(measured in square meters per second), where the subscript “v/c” means that the relation holds for either vapor (v) or cloud water (c).
In a general formulation, and would be vectors in the direction of advection. Let us, for simplicity, assume that the coordinate system is aligned in such a way that advection acts in the x direction. In addition, we assume dispersion in the y direction, i.e., in the direction normal to the advection. Then the vertically integrated moisture balance for each of the components reads
where is a source term (measured in meters per second) describing the interaction between the components and the loss by precipitation.
The dispersion term used in Eq. (4) is a specific form of a diffusion term with a diffusivity , where Ld is the dispersion length. Dispersion terms in advection equations typically arise from a spatial variability in velocity that is not resolved by the large-scale description of the flow field. Assuming a constant dispersion length Ld reflects the idea that the fluctuations in velocity and thus the diffusivity are directly proportional to the large-scale velocity. However, assuming a constant dispersion length is not essential for the model developed here. Similarly to assuming the same dispersion length for both components, it is just a convenient choice that keeps the equations simple.
Dispersion in the longitudinal direction is not taken into account, although there is no reason why it should be smaller than the transversal direction. The reason for including only transversal dispersion is that it has a larger effect on the properties of the model. Transversal dispersion is the only process that links points with different y values. Without transversal dispersion, the precipitation pattern would fall into a set of individual lines parallel to the x axis. So transversal dispersion is an essential component of the approach in combination with two-dimensional landform evolution models. In turn, we will see in Sect. 4.1 that longitudinal dispersion is not essential for the properties of the model, while it would make the numerical treatment more complicated (Sect. 3).
Since the timescales of processes in the atmosphere are much shorter than the timescales involved in landform evolution, steady-state conditions can be assumed in Eq. (4). If we furthermore assume that the velocities are constant, Eq. (4) can be written conveniently in terms of the fluxes per unit width (Eq. 3):
Following the concepts of Smith and Barstad (2004), we assume that condensation (from Qv to Qc) and precipitation (from Qc) are linear processes with given time constants τc and τf, respectively. In contrast to this model and also to the model of Garcia-Castellanos (2007), we do not introduce an equilibrium water content explicitly. Instead, we start from a more fundamental level by considering condensation of vapor and re-evaporation of cloud water (e.g., Roe, 2005) as competing processes in the form
The nondimensional coefficient α defines the dynamic equilibrium between the two processes. An equilibrium between Qc and Qv is achieved if . Rewriting Eq. (6) in terms of fluxes per unit width yields
with the length scale of condensation Lc=vvτc and the modified coefficient (still nondimensional).
Since the extension of the conversion of vapor into cloud water by a (negative) linear feedback term (evaporation) is the key idea behind our approach, the model is called the linear feedback precipitation model (LFPM) in the following.
The rate of precipitation (measured in meters per second) can also be expressed in terms of the flux per unit width according to
with the length scale Lf=vcτf. Then the source term of Qc is
and the full system of differential equations for the two fluxes reads
2.2 The effect of topography
Orographic precipitation is related to a dependence of the equilibrium on altitude (e.g., Roe, 2005). Since the re-evaporation of cloud water requires energy, altitude has an immediate effect here. The rate of re-evaporation should decrease with decreasing temperature and thus with increasing altitude. As the simplest approach, we consider only this effect and assume that the length scales of condensation (Lc) and fallout (Lf) are constant. The Arrhenius relation,
with a constant a, provides the simplest model for the dependence of β on the temperature T, where both a and T are measured in Kelvin. Using a linear decrease in temperature with altitude H,
where T0 is the temperature at sea level and Γ the lapse rate (measured in Kelvin per meter), Eq. (12) can be written in the form
where β refers to the altitude H and β0 to sea level. Defining , Eq. (14) can be written in the form
For simplicity, we assume that Eq. (15) also holds for the vertically integrated cloud water content with H as the surface elevation and neglect the term . The latter is a first-order approximation concerning H, which requires that the decrease in temperature ΓH is small compared to the absolute temperature T0 at sea level. Using these approximations, Eq. (15) reduces to
This relation describes the decrease in β by a single lumped parameter H0, which defines a vertical length scale and describes the elevation at which β has decreased by a factor e compared to sea level. While the description of the height dependence by a single, lumped parameter is convenient, it is not an essential part of the LFPM. Any other relation, e.g., the more elaborate version used by Garcia-Castellanos (2007), which does not rely on the two approximations introduced above, could be used as well.
2.3 Boundary conditions
Since the system of differential equations defined by Eqs. (10) and (11) is of first order in x and of second order in y, it is a parabolic system. Finding a unique solution in a rectangular domain (, ) requires boundary conditions at x=0, y=0, and y=ymax (but not at x=xmax).
Since moisture is coming in at x=0, it is straightforward to define Fv and Fc there. Then the integral of the total influx over this boundary defines the total amount of water available for precipitation in the domain. However, the question of how to distribute a given total influx F to Fv and Fc is not trivial and requires more knowledge about the properties of the model. It will be addressed in Sect. 4.1.
All types of boundary conditions could be used at y=0 and y=ymax. Neumann boundary conditions or periodic boundary conditions are more useful than Dirichlet boundary conditions here since the Fv and Fc are fluxes along these boundaries, and it is not trivial to define reasonable prescribed values for Fv and Fc. Homogeneous Neumann boundary conditions define , which means that there is no transversal dispersion across these boundaries. The implementation in the landform evolution OpenLEM presented in the following section uses periodic boundary conditions in the y direction by default, which are convenient in many applications.
Taking into account advection only along one of the coordinate axes and neglecting longitudinal dispersion considerably facilitates the numerical implementation of the model. Let us first rewrite Eqs. (10) and (11) in matrix form:
Let us further assume unit grid spacing in both directions, , for simplicity. This means that the length scales Ld, Lc, and Lf must be measured in terms of the grid spacing in the following. If we use a left-hand (so upwind) difference quotient for the advection term, the discretized form of Eq. (17) can be written in the form
where the indices x and y correspond to the positions. So the values can be computed from the values by solving a one-dimensional problem (in the y direction). The respective linear equation has a tridiagonal structure of 2×2 blocks and can be written in the form
with the 2×2 matrix
and the 2×2 identity matrix 1. This equation system can be solved, e.g., by the direct Gaussian scheme based on 2×2 blocks.
The examples shown in the following section are computed using the open-source landform evolution model OpenLEM. This model already contains up-to-date implementations of fluvial erosion such as the shared stream-power model (Hergarten, 2020), which will be used in Sect. 8. All components of OpenLEM are of linear time complexity at arbitrary time step lengths, which means that the numerical effort increases only linearly with the size of the lattice. The computation of the precipitation proposed here preserves this property. Independently of the size of the lattice, we found an increase in computing time by a factor of about 2.4 compared to the simplest form of the SPIM and a factor of about 2.2 compared to the shared stream-power model, which includes sediment transport. This increase is owing to taking into account transversal dispersion, for which Neumann boundary conditions would be cheaper than the periodic boundary conditions used in OpenLEM.
4.1 Characteristic length scales
Let us for the moment consider the model only in the longitudinal direction, i.e., without the dispersion term, and let us assume a constant elevation for the moment. Then the set of parameters consists of two horizontal length scales Lc and Lf as well as a nondimensional parameter β. In this section, it is shown that the relevant length scales that characterize the properties of the model differ from Lc and Lf.
In this situation, Eq. (17) reduces to a linear system of two ordinary differential equations,
The behavior of the solutions is determined by the eigenvalues of the matrix A defined in Eq. (18). These are found by solving the characteristic equation of A,
Since β and ϕ are nondimensional, the eigenvalues λ± are also nondimensional. The eigenvalues describe exponentially decaying solutions of the form . These solutions can also be written in the form and , respectively, where
describe the respective length scales of the decay. Since and for all values of β and ϕ,
for β>0. In turn, , and thus
As a consequence,
for β>0. So the approach based on the dynamic equilibrium creates two characteristic length scales outside the range between Lc and Lf. These scales differ even if we assume Lc=Lf as suggested by Smith and Barstad (2004).
The longer length scale Ll describes the ability to transport moisture over large distances. Smith and Barstad (2004) suggested timescales of 200 to 2000 s for the conversion of cloud water and for fallout, corresponding to length scales Lc and Lf of 10 to 100 km at wind speeds of 50 m s−1. If we, e.g., assume β=10, Ll is in the range between 119 and 1191 km. So the transport range may be considerably larger than the length scales Lc and Lf of the involved processes. This property is essential for simulating long-range transport over large continental areas with a closed water balance.
Interestingly, Makarieva et al. (2009) indeed found an exponential decay of mean precipitation rates with increasing distance from the ocean at nonforested areas in a worldwide analysis. They found a decay length of about 600 km, which is 6 to 60 times larger than the reasonable range of Lc and Lf that even refers to quite high wind speeds. This result supports the idea behind the LFPM and provides an idea about the order of magnitude of Ll.
The mode of long-range transport is obtained by inserting λ− into the first row of the eigenvalue equation AF=λF,
According to Eq. (7), this ratio is in equilibrium. So the vapor content is slightly above its equilibrium value in the long-range transport mode, which results in a low net rate of condensation.
Figure 1 illustrates the long-range transport for a boxcar-shaped topography with a height H=H0. All properties are considered nondimensional values. The parameter values are and β0=10. The incoming fluxes are Fv=10 and Fc=0 at x=0. According to Eqs. (24) and (26), the length scale of long-range transport is Ll≈11.9 at sea level (H=0). Both fluxes and the rate of precipitation decrease exponentially with this length scale for x<5 and x>10, except for the beginning of the ranges.
The plateau (H=H0) is characterized by a lower value of according to Eq. (7), resulting in a lower length scale Ll≈5.5. So the precipitation is higher at the plateau but in turn decreases more rapidly with x. This difference is reflected in a lower ratio and therefore in a lower ability to keep moisture in form of vapor.
In turn, the short length scale Ls describes the adjustment if the ratio of Fc and Fv deviates from Eq. (30). Since these deviations predominantly arise from changes in topography (via the elevation-dependence of β), the length scale Ls can be considered the length scale of orographic precipitation.
Three transition zones characterized by Ls occur in Fig. 1, starting at x=0, x=5, and x=10. The length scale of the first and the third transition (H=0) is Ls≈0.08, while it is Ls≈0.18 for the second transition (H=H0). In general, Ls increases with elevation, while the length scale of long-range transport Ll decreases with elevation. Their product is constant according to Eq. (27).
While the second and the third transitions arise from changes in topography, the first transition occurs because it is assumed that the influx only contains vapor (Fv=10) but no cloud water (Fc=0), which is far off from equilibrium and from the long-range transport mode. It is therefore useful to adjust the boundary condition in such a way that the incoming fluxes are in the long-range transport mode described by Eq. (30). Then the incoming flux of cloud water must be
where F is the total incoming flux, and . This modified boundary condition is used in all subsequent examples throughout this study.
The non-instantaneous reaction to abrupt changes in topography is a central property of the model. Without this, the small-scale roughness of topography would directly affect the precipitation pattern so that an additional smoothing procedure would be required. Longitudinal dispersion could also used for smoothing, but as the LFPM generates a scale of smoothing on its own, taking into account longitudinal dispersion is not urgently required. Taking this result into account, our approach based on two components of water content with a two-way conversion is some kind of minimum model that is able to capture both continentality (a slow decrease in precipitation at large scales) and a delayed reaction to small-scale changes in topography.
While the model in its original form involves two longitudinal length scales Lc and Lf, a transversal length scale Ld, a vertical length scale H0, and a nondimensional parameter β0 (referring to sea level), it is also possible to replace β0 by the length scale of long-range transport Ll (alternatively also by Ls, but that would be less useful). The value of β can be computed conveniently from Eq. (23):
While this relation is valid for the respective values of Ll and β at any elevation, it is particularly useful for computing β0 from Lc, Lf, and Ll at sea level.
If we use Ll instead of β0, Eq. (27) reveals that the short length scale Ls only depends on Ll and on the product LcLf, while the individual values of Lc and Lf are not relevant for Ls at sea level. This is, however, not the case for the elevation dependence. Figure 2 illustrates the relevance of the individual values of Lc and Lf in combination with topography. The default scenario (solid blue line) is the same as in Fig. 1, except that the fluxes at the boundary were adjusted according to Eq. (31). In all scenarios, Ll (at sea level) was kept constant (so not β0). As expected, the behavior at sea level remains the same for Lc≠Lf (red lines) as long as the product LcLf is constant. This is, however, not true at H>0, where the increase in precipitation with elevation becomes distinctly weaker for Lc≠Lf (red vs. blue lines), regardless which of the values is greater.
The decrease in Ll with elevation and thus the respective increase in precipitation can be computed from Eq. (16) according to
The remaining derivative can be computed from Eq. (23),
Inserting this result into Eq. (34) yields
after some basic transformations, and finally after inserting Eq. (33),
It is easily recognized that the second factor is always lower than Ll. So the decrease in Ll and thus the increase in precipitation with elevation are always smaller than the decrease in β, which is characterized by the vertical scale H0. The relation is symmetric concerning Lc and Lf, and the elevation dependence is strongest for .
At least for topographies with moderate relief (in relation to H0), a difference between Lc and Lf can be replaced by an increased value of H0. If we define
the behavior of the model essentially remains the same for moderate elevations. This result is illustrated by the green lines in Fig. 2. The precipitation obtained for with an increased reference elevation H0=1.2 (Eq. 38) is close to those for Lc=3 and with H0=1 for the topography with . For the higher topography with H=1, however, the remaining deviation is larger.
Keeping in mind that the definition of H0 in Sect. 2.2 already required some approximations, it is not a problem that H0 has to be increased artificially if we replace different values of Lc and Lf by the same value . If we accept that there is a residual overestimation of the effect of topography that increases with elevation, we can assume Lc=Lf without losing much of the model's fundamental capabilities.
4.2 The influence of transversal dispersion
As mentioned above, transversal dispersion is the only component that prevents the model from falling into a set of independent one-dimensional models. In contrast to the other length scales of the model, however, the dispersion length cannot be interpreted directly as a spatial scale. It rather links longitudinal and transversal length scales of the moisture pattern.
Let us for the moment assume that condensation and fallout are switched off () and that the topography is flat (H=0). Let us further assume that the incoming flux (at x=0) has some transversal variation in water content (either in Fv or Fc or in both) according to
So small-scale transversal patterns decay much faster in the direction of advection than large-scale patterns. The length scale of the decay, Lx, decreases quadratically with the scale Ly of the transversal pattern. The length scale of dispersion, Ld, describes the strength of dispersion so that an increase in Ld reduces Lx.
Figure 3 illustrates the effect of dispersion for an obstacle of a width Ly=1 and a height H=1, where the parameter values are the same as in the previous examples. Without transversal dispersion (Ld=0), the higher precipitation falling on the obstacle causes an infinite precipitation shadow. For Ld=0.01, the longitudinal scale of decay is Lx≈10 (Eq. 41). The precipitation shadow has become considerably weaker at this distance behind the obstacle but is still visible. Finally, the scale of decay decreases to Lx≈1 for Ld=0.1 so that the shadow vanishes rapidly behind the obstacle.
Evaporation including the transpiration by plants, called evapotranspiration, plays a major part in the water balance. While the potential rate of evapotranspiration mainly depends on the climatic conditions and on vegetation, the actual rate is often much lower due to limited availability of water at the surface and in the shallow subsurface.
However, the concept for including variations in precipitation in large-scale landform evolution models is not able to predict the availability of water. The water balance is stated in terms of fluxes, which are not directly related to amounts of stored water. Estimating the amount of stored water would require a model for the flow velocity and would introduce additional complexity. Garcia-Castellanos (2007) presented a first step in this direction by distinguishing lake areas and assigning a rate of evaporation to these areas.
Here we propose a simpler idea by assuming that the rate of evapotranspiration is proportional to the rate of precipitation instead of the amount of stored water. This means that a given fraction ϵ of the precipitation evaporates immediately. This leads to one additional term in Eq. (10) so that the system of differential equations turns into
While the total precipitation is still , the effective precipitation that contributes to runoff is then
In order to understand the effect of evapotranspiration, Eqs. (42) and (43) can be brought to the same form as Eqs. (10) and (11) by introducing an increased coefficient for re-evaporation in the atmosphere,
and an increased fallout length,
Then ϕ (Eq. 18) changes to
It is easily recognized that . So only the last term ϕ in Eq. (24) is affected by evapotranspiration. The smaller eigenvalue λ− comes closer to zero then, while the greater eigenvalue λ+ does not change much. So Ll increases considerably, while Ls remains almost constant, and thus
according to Eq. (27).
While these results suggest that the effect of evapotranspiration could be mimicked by modifying the parameters β and Lf (Eqs. 45 and 46), it must be kept in mind that this only holds at constant elevation. Otherwise, β depends on H and thus also the effect of ϵ on Ll. Since β describes the re-evaporation of water in the atmosphere, it makes sense to assume that the rate of evapotranspiration has the same dependence on H as β,
although ϵ refers to the surface and β to vertically integrated properties.
Then the dependence of Eq. (45) remains valid for all elevations. However, as ϵ decreases with elevation, the increase in the effective Lf (Eq. 46) and thus also in the effective Ll becomes weaker at large elevations. So mimicking evapotranspiration by adjusting β and Lf (Eqs. 45 and 46) would overestimate the effect of evaporation at large elevations or, in turn, underestimate the effect of elevation. Figure 4 illustrates this result for the example from Fig. 1, where an evaporation ratio ϵ0=0.5 is assumed at sea level. The results are compared to the model without evaporation but with modified parameter values (Eq. 45) and (Eq. 46). The effective precipitation in the mountain region and the effective length scale of transport differ by almost a factor of 2 for H=H0. So mimicking evapotranspiration by adjusted parameters is only useful for small topography.
As discussed in Sect. 1, the models of Smith and Barstad (2004) (SB model) and Garcia-Castellanos (2007) also use vertically integrated water contents and advective transport at a given wind velocity. In spirit, the LFPM is somewhat similar to these models.
It may seem at first that taking into account transversal dispersion was the major progress of our approach. However, the numerical scheme proposed in Sect. 3 could in principle also be used for including transversal dispersion in the two other models. The fundamental differences are hidden in details of the model structure that have a bigger effect than it seems at first.
The SB model assumes only a one-way coupling between the two moisture components. In our terminology, this would be β=0 in Eqs. (10) and (11). Following the considerations of Sect. 4.1, the length scale of long-range transport of moisture is then, which is between about 10 and 100 km at rather high wind speeds of 50 m s−1. So assuming a one-way conversion practically removes the ability to transport moisture over large distances of several hundred kilometers. Therefore, the SB model requires a permanent refilling of the water storages at some given background rate. When considering a large plain, the precipitation rate will always approach this prescribed background rate, regardless of the topography in front of the plain. So this model focuses on the behavior at intermediate scales but cannot capture large-scale precipitation patterns. This is presumably the reason why the two moisture components are interpreted as cloud water and hydrometeors in the SB model instead of vapor and cloud water in the LFPM.
Beyond this, the feedback parameter β carries the information about the elevation dependence in the LFPM. So the effect of topography must be included in another way if this feedback is not taken into account. Smith and Barstad (2004) assumed that the rate of conversion is not proportional to the absolute value of the cloud water content (Qv in the LFPM) but to the difference of this content towards an elevation-dependent equilibrium content. This results in an additional source or sink term in the equations, which carries all information about the topography.
In the earliest version of the SB model (Smith, 2003), it was assumed that the source term is directly proportional to the topographic slope in flow direction, . So upslope flow introduces a positive source term, which is converted into precipitation with some lag and smoothing. Smith and Barstad (2004) proposed a more elaborate source term taking into account airflow dynamics in more detail.
In the following, we compare the two versions of the SB model to the LFPM in a one-dimensional example that describes the rise to a plateau. In contrast to the similar example considered by Smith and Barstad (2004, Sect. 3c), we do not use a smooth arctangent function but a ramp with a constant slope between two horizontal planes for clarity. Four different topographies are considered in Fig. 5, where nondimensional coordinates are used. Similarly to Fig. 1, Lc=Lf is used as the horizontal length scale, while H0 defines the vertical length scale. The four scenarios refer to plateau elevations of H=H0 and H=4 H0 and to ramp lengths of 2.5 Lc and 10 Lc. All precipitation values are arbitrarily scaled but by constant factors for each model throughout all scenarios. So the absolute precipitation values cannot be compared among the SB model and the LFPM but among different scenarios for the same model. The background rate is zero in all results of the SB model so that the values shown in Fig. 5 should not be interpreted as absolute values but as differences towards a prescribed background precipitation rate for the SB model.
The upslope version of the SB model shows the simplest behavior. The precipitation rate is zero (or equal to the prescribed background rate) in front of the ramp as well as at the plateau far behind the ramp. Since the ramp introduces a constant source term, a constant precipitation rate is also approached at the ramp if the ramp is sufficiently long (Fig. 5b and d). Otherwise, the increase in precipitation ceases at the end of the ramp before a constant rate is approached, which results in a distinct maximum in precipitation at the transition to the plateau.
For the version of the SB model with the more elaborate airflow dynamics, the hydrostatic limit was considered. This model version involves one additional nondimensional parameter (Smith and Barstad, 2004, Sect. 3b); we use suggested as a typical value there. Since the simple upslope version corresponds to , the effect of values could be estimated qualitatively from the curves. The most important effect of the more elaborate airflow dynamics is that the source term no longer depends only on the local slope at the considered location, but also on the slopes in a larger part of the domain. As a consequence, the precipitation at a given point not only depends on the topography in the upwind range, but to some extent also on the topography of the downwind region. This becomes obvious at the plain in front of the ramp, where the precipitation rate already increases at a considerable distance to the ramp. Next, the peak in precipitation is shifted upwind and therefore towards lower elevations at the ramp. Finally, the values become negative at the plateau, which means that the precipitation rate will be lower than the background rate.
Overall, the more elaborate airflow dynamics increase the spatial range over which the topography affects precipitation compared to the simple upslope version of the SB model. However, the precipitation approaches the same value (the background rate) in both directions far away from the ramp for all versions of the SB model, independent of the elevation of the plateau. So none of the versions of the SB model can predict precipitation rates at a high plateau that are much lower than those in the plain in front of the plateau, as occurs, e.g., in the Tibetan Plateau.
The results obtained from the LFPM are readily understood from the considerations made in the previous sections, with the depletion of moisture by precipitation as the most important phenomenon. If the plateau is low (Fig. 5a and b), the increase in precipitation is moderate. Since the depletion is also moderate then, the plateau is still exposed to considerable precipitation. In turn, precipitation decreases rapidly along the high plateau. The amount of moisture that enters the plateau depends on the length of the ramp, while the further depletion along the plateau is related to the reduced value of β at large elevations. For the long and high ramp (Fig. 5d), the depletion at the ramp is so strong that the entire plateau is quite dry.
Restricted to the region close to the ramp, however, the LFPM and the SB model can be adjusted to yield similar results, although for different reasons. In the LFPM, it is the combination of the elevation-dependent conversion process and the limited amount of moisture, while much of the behavior depends on the model used for the source term in the SB model. The linearity of the SB model also deserves attention in this context. While both models are linear with regard to the water contents, the SB model is also linear concerning the topography. As a consequence, precipitation patterns for the low and high ramps in Fig. 5 (a, c and b, d) are the same, with only the absolute values scaling like the elevation of the plateau. Owing to the linearity, precipitation patterns depend only on the lateral structure of the topography but not on the absolute elevation in all versions of the SB model. While this can be fixed to some extent by adjusting the model parameters for individual scenarios, it becomes a serious limitation in the context of co-evolution of topography and climate, e.g., if the same uplift pattern is considered at different absolute rates. As will be shown in Sect. 8, the LFPM is more powerful here.
The linearity of the SB model also affects the precipitation at the leeward side of mountains. The precipitation of a declining ramp would just be opposite to that of an increasing ramp. So if we added a declining ramp behind a large plateau in Fig. 5, the precipitation rates would be negative there. Taking into account the fact that a background rate has to be added to the precipitation rates shown in Fig. 5, this might not be a crucial problem in some of the scenarios, but in particular for short and high ramps (Fig. 5c), the negative rates would be too high to be compensated for by the background rate. In order to overcome this problem, Smith and Barstad (2004) suggested truncating the precipitation term explicitly at negative values. However, the SB model still predicts extremely dry leeward sides, and improving this behavior was obviously one of the motivations for extending the model by multiple layers proposed by Barstad and Schüller (2011). In contrast, the LFPM does not tend towards extremely dry leeward regions, and it is guaranteed that neither Fv nor Fc can become negative. So the LFPM requires no artificial measures at the leeward side, in particular no truncation in order to avoid negative precipitation rates.
In this sense, including the feedback by re-evaporation in the LFPM may look more complicated first, but it is the key to capturing large-scale precipitation patterns and avoids the need for taking additional measures at the leeward side of mountains.
The model of Garcia-Castellanos (2007) uses a single moisture component and thus seems to be simpler than our approach. The fundamental structure of this model can be explained by considering the limit Lc→0. This means that the contents of vapor and cloud water immediately achieve an equilibrium , defined by the elevation-dependent value β. Adding Eqs. (10) and (11) then yields
where is the total flux per unit width. This is the fundamental structure of the model of Garcia-Castellanos (2007) except for the dispersion term. The precipitation rate is the ratio of F and an elevation-dependent value (1+β)Lf. The length scale of transport is then, which is much larger than Lf if β is sufficiently large. So replacing the coupling between the two flux components by an equilibrium preserves the ability to capture large-scale precipitation dynamics.
However, the short length scale Ls is lost if we assume instantaneous equilibrium. Precipitation reacts immediately to changes in topography then so that the precipitation pattern becomes sensitive to the small-scale roughness. In order to overcome this problem, Garcia-Castellanos (2007) introduced an additional smoothing by applying a convolution with half of a Gaussian curve in the upwind direction. As illustrated in Fig. 6, the behavior of the one-component version with additional smoothing is indeed similar to our two-component approach in 1D. Here, a smoothing length Δx=2 Ls appears to be a reasonable choice at H=0. However, it should be emphasized that the concept of smoothing including the choice of a smoothing length is an ad hoc assumption, while smoothing automatically emerges in the LFPM. Beyond this, applying the convolution with half of a Gaussian function is numerically more expensive than the treatment of two moisture components.
In this sense, considering two moisture components is indeed the simplest choice if we want to combine long-range transport and a smooth response to sharp topographic gradients in a linear model. This result is directly related to the occurrence of two eigenvalues, i.e., horizontal length scales discussed in Sect. 4.1. A single-component model seems to be simpler at first sight, but we have to pay for this simplicity later as soon as we need the second length scale.
Moreover, the model of Garcia-Castellanos (2007) contains a nonlinear component, which is not included in the LFPM. While the precipitation term in Eq. (50) was obtained by simplifying our two-component model, Garcia-Castellanos (2007) introduced an expression with this structure directly in the form
where P0 is some reference precipitation, and Qmax(H) is an elevation-dependent maximum water content. While this expression is still equivalent to the precipitation term in Eq. (50), it was extended in such a way that Eq. (51) is only applied for Q<Qmax(H). Otherwise, it was assumed that the excess water content Q−Qmax(H) is immediately converted into precipitation. Transferred to the formalism of Eq. (50), this occurs if the ratio exceeds a given threshold.
The red curve in Fig. 6 illustrates this effect; it was assumed that the air at the left-hand boundary is just at the limit Q=Qmax(0). While Q<Qmax(0) in the foreland, Qmax(H) (via 1+β in our formalism) decreases suddenly at x=5 due to the sudden increase in H so that Q>Qmax(H). This would even result in a sharp peak in the precipitation curve without smoothing. Smoothing in the upstream direction reduces the height of the peak and widens it into the plateau.
As a main effect of this extension, the model of Garcia-Castellanos (2007) becomes able to predict some kind of overshooting in precipitation at the windward side of mountains. For a high plateau, the original version predicts a gentle increase in precipitation to the level at the plateau. A decrease in precipitation only occurs due to the decrease in water content or, if we consider a mountain range, due to decreasing elevation at the leeward side.
Such an effect might be useful, although the physical basis is not trivial since a dynamic equilibrium between vapor and cloud water is already included in the linear model. The LFPM could be extended by nonlinearity in several ways. The precipitation process may be a candidate here since coagulation plays a part in the growth of hydrometeors, and the rate of coagulation increases rather quadratically than linearly with concentration. However, following the concept of parsimony, we do not follow ideas of nonlinearity further in this study.
This section presents an application of the LFPM to real topography. It should, however, be seen as an illustration rather than as a validation. While the attempt to validate the SB model by Barstad and Smith (2005) suffered from the availability of data at a sufficient spatial resolution, we must keep in mind that all models discussed in this paper were not tailored for reproducing precipitation patterns exactly. As discussed in Sect. 6, predicting the effects of changes in topography on precipitation is more important and more challenging in the context of landform evolution than, e.g., adjusting a model to predict the precipitation rate at the front of a mountain range close to the coast. Thus, a serious validation would have to be based on several locations and scenarios in which the LFPM and other models would have to be tested against real-world data and regional climate models.
Figure 7 compares the precipitation pattern of the India–Asia collision zone modeled with the LFPM to the annual precipitation pattern of the TRMM2b31 dataset (Bookhagen and Burbank, 2010) and the WorldClim2.1 precipitation data (Fick and Hijmans, 2017). The TRMM2b31 dataset is an outcome of the Tropical Rainfall Measurement Mission, for which the average precipitation rates are based on a 12-year (1998 to 2009) time series. The dataset has a spatial resolution of roughly 5 km and covers the region between 36∘ S and 36∘ N. The WorldClim2.1 dataset is a global dataset of climate variables with spatial resolutions between 10 arcmin (≈18 km) and 30 arcsec (≈0.9 km). Monthly precipitation data are averaged for the years 1970 to 2000 and based on a large number of weather stations. Integrating the average precipitation from January to December results in the annual precipitation dataset shown in Fig. 7d.
Simulations of the LFPM were performed on a regular grid with 1 km mesh width for two different directions of atmospheric flow (from south to north and from southwest to northeast). Ocean areas were considered to be boundaries, where a uniform influx was assumed. The domain was extended in such a way that each flow line (either from south to north or from southwest to northeast) starts from a point in the ocean.
Several simulations with different parameter values were conducted. However, it is immediately recognized in Fig. 7 that the flow direction has a strong influence on the precipitation pattern. Compared to this influence, the effect of the parameter values is much weaker. In particular, similar results in terms of the root mean square (rms) deviation were found for different combinations of the parameter values. So a formal fit would not be very useful. Instead, km and Ll=500 km were chosen as reasonable values according to the considerations in the previous sections. A rather high evaporation ratio ϵ0=0.75 was assumed, while lower values would yield similar results with greater values of Ll according to the findings of Sect. 5. Finally, H0=2 km and Ld=25 km were assumed. The influx at the boundaries was adjusted automatically in such a way that the average precipitation over the domain matches the respective average TRMM2b31 precipitation (0.99 m yr−1), which is very close to the average WorldClim2.1 precipitation (1.00 m yr−1).
While some large-scale properties of the precipitation pattern such as wet regions at the western coast of India (Ghats escarpment) and at the orographic front of the Himalaya and the dry Tibetan Plateau are at least qualitatively reproduced by the LFPM, deviations in precipitation rate between the different datasets are distinct. The root mean square (rms) deviation towards the TRMM2b31 dataset for the given parameters with a wind direction from south and from southwest amount to 0.81 m yr−1 and 0.84 m yr−1, respectively. The rms deviations towards the Worldclim2.1 dataset are slightly lower (0.77 and 0.80 m yr−1).
The rms deviations can be reduced to about 0.6 m yr−1 by tuning the parameters. In particular, increasing H0 goes along with higher precipitation rates in the Tibetan Plateau (closer to the measured data) and lower precipitation rates at the windward side of topographic barriers. While this tuning reduces the rms deviation, it does not necessarily result in a better reproduction of the most striking precipitation features such as high orographic precipitation at the windward side and rain shadows at the leeward side. Furthermore, large values of H0 would lead to very small orographic effects in other much lower mountain ranges.
In general, we should be careful not to compensate for principal limitations of the model with potentially unrealistic parameter values. In particular, this applies to the pre-defined uniform wind direction in the LFPM, which cannot describe the atmospheric circulation pattern of the entire region sufficiently well. Furthermore, obvious differences also occur between the TRMM2b31 and the Worldclim2.1 dataset. Compared with the WorldClim2.1 data, the TRMM2b31 data indicate much higher precipitation rates at the western Himalaya but distinctly drier conditions at the eastern Tibetan Plateau. These differences suggest that the precipitation rates still involve a considerable uncertainty. The rms deviation between these two datasets is 0.42 m yr−1, which is not far below the deviation of our “best-fit” parameter set.
Similarly to the approaches of Roe et al. (2003), Smith and Barstad (2004), and Garcia-Castellanos (2007), the scope of the model developed in this study is not a precise prediction of precipitation rates but its combination with long-term landform evolution. This section provides some examples exploring the effect of orographic precipitation and continentality on fluvial landform evolution and the feedback of the resulting topography on the precipitation pattern.
As described in Sect. 1, the SPIM and its derivates can easily be extended by variable effective precipitation and thus be coupled with the LFPM. In the following, we use a model that is not restricted to pure bedrock incision but also takes into account sediment transport. While the idea behind this model dates back to Whipple and Tucker (2002) and Davy and Lague (2009) or partly even to older studies (Howard, 1994; Kooi and Beaumont, 1994), it is used here in the most recent formulation, the so-called shared stream-power model (Hergarten, 2020). The constitutive equation of the shared stream-power model reads
where Q is the sediment flux (volume per time; not to be confused with the atmospheric water contents Qv and Qc). This model contains two erodibilities, where Kd describes the erodibility in the absence of transported sediment, while Kt characterizes the ability to transport sediment at zero erosion. For a deeper insight into the properties of the shared stream-power model and the meaning of its parameters, the reader is referred to Hergarten (2021b).
Following Robl et al. (2017), we use n=1 (i.e., the linear version of the model), m=0.5, and K=2.5 Myr−1. Studies of natural and experimental river profiles at the transition to a foreland by Guerit et al. (2019) suggest a ratio (G in their notation) for n=1, which leads to Kd≈6.5 Myr−1 and Kt≈4.1 Myr−1.
The concept of expressing discharges as their catchment-size equivalent (referring to a hypothetic uniform reference precipitation rate) is particularly useful in the context of the shared stream-power model. If both occurrences of A in Eq. (52) are interpreted as catchment-size equivalents of the actual discharge, Eq. (52) remains formally the same, including the values and the units of Kd and Kt.
We use a regular mesh of 2000×2000 nodes for a domain of 500 km linear size, corresponding to a spatial resolution of 250 m. Although rather coarse, hillslope processes are still relevant at this scale. Using a purely fluvial model would lead to artificially steep slopes and thus increased elevations at the drainage divides, which may affect the precipitation pattern despite the robustness of the model against the small-scale roughness of the topography. In order to avoid this, we use an approach that was brought into play in the context of debris flows in steep valleys by Stock and Dietrich (2003) and developed further by Hergarten et al. (2016). This approach replaces the term Aθ in Eq. (1) (or here in Eq. 52) by , where and Ac is a given constant. This modification acts like an increased catchment size or like an increased discharge and thus avoids the occurrence of extremely steep slopes at small catchment sizes. We use Ac=1 pixel ≈ 0.06 km2 here, which is roughly consistent with the estimate Ac=0.05 km2 obtained by Hergarten et al. (2016) for the topography of Taiwan.
In the following section, we show how the decrease in precipitation rate with distance to the source of moisture (i.e., continentality) controls the shape and the height of large mountain ranges. Then we illustrate topographic patterns resulting from feedbacks between rock uplift, orographic precipitation, and fluvial erosion.
8.1 Impact of continentality on landform evolution
The large-scale precipitation pattern over continental areas is controlled by the length scale of long-range moisture transport Ll. If the extension of a mountain belt in the wind direction reaches the order of magnitude of Ll, the precipitation pattern may have a strong influence on its height and shape even without any immediate effect of elevation on precipitation. In terms of the model parameters, this situation is described by a reference elevation H0 much larger than the surface elevation. Then the precipitation pattern reflects increasing continentality with an exponential drop in precipitation rate from the moisture source towards the continental inland. The precipitation rate is solely controlled by the influx Fin and the length scale Ll, which is further stretched by considering evaporation.
Figure 8 shows the effect of continentality on topography. The considered mountain range is 300 km wide and uplifted at a rate of 0.25 km Myr−1 in all three examples. The two 100 km wide foreland regions are tectonically inactive. Moisture enters the model domain at the southern boundary and is advected towards north. While the geometry, the length scales km, the elevation-independent evaporation ratio ϵ=0.5, and the amount of moisture entering the southern boundary are the same in all three scenarios, the length scale Ll of the long-range transport varies from 50 to 600 km. Only steady-state topographies are considered, with the term steady state used in a sloppy way here since the drainage pattern in an inactive foreland permanently reorganizes (Yuan et al., 2019). Since this reorganization has a minor effect on the mountain topography, the mountain range comes close to a steady state with some fluctuations.
Two effects of continentality must be distinguished. First, Ll has an influence on the amount of precipitation that the mountain range receives in total. This amount should affect the total height of the steady-state mountain range. It is small if Ll is small since most of the moisture is already lost in the foreland, but it is also small if Ll is large since most of the moisture passes the domain without much precipitation then. So there must be a length Ll at which the amount of precipitation on the mountain range becomes maximal.
For Ll=600 km (Fig. 8a and b), evaporation with ϵ=0.5 effectively stretches Ll to about 1200 km following the considerations of Sect. 5. Since this is 4 times the width of the mountain range, the precipitation rate varies by a factor of e¼≈1.28 over the mountain range. Differences in discharge are, however, smaller. If the drainage divide is in the middle, the difference in total precipitation differs only by a factor of e⅛≈1.13 between the windward and the leeward half of the mountain range. So the discharges of the big rivers only differ by this factor at the edges of the mountain range.
This difference is visible in the swath profile (Fig. 8b). The distribution of maximum elevations across the mountain range is already quite asymmetric since small catchment sizes at ridge lines and hillslopes cause the erosion rate to be directly related to the local precipitation rate. Hence, the highest domains become increasingly steeper towards the leeward side. In contrast, the minimum elevation of the swath profiles describes the large rivers, which do not differ so much in their discharge. So the profile of the minimum elevation is still quite symmetric here. This also implies that the relative incision of the rivers in relation to the hillslopes is deeper at the leeward side than at the windward side.
According to the findings of Sect. 5, the length scale Ll=100 km considered in Fig. 8c and d is stretched by evaporation to 206 km. This value is close to the length scale with which the mountain range receives the maximum amount of precipitation in total, km. Consequently, the overall height of the mountain range is quite low here. Since precipitation varies by a factor of here, the topography becomes strongly asymmetric. Again, this asymmetry mainly concerns the maximum elevation in the swath profile (Fig. 8d), while the asymmetry of the minimum elevation referring to the largest rivers is smaller. Nevertheless, the asymmetry in the minimum elevation is clearly visible here, and it goes along with a shift of the principal drainage divide towards the leeward side. While the total area drained by the leeward part of the mountain range was 48 % of the total area of the mountain range for Ll=600 km, this fraction has decreased to 42 % now. Despite this moderate shift, the highest peaks are already separated from the main drainage divide. For Ll=50 km (Fig. 8e and f, effectively 114 km with evaporation), the total amount of precipitation on the orogen decreases. This results in an increasing overall surface elevation. More importantly, the topography becomes extremely asymmetric since the precipitation rate varies by a factor of over the mountain range. This variation is strong enough to make the river profiles (minimum elevation in Fig. 8f) strongly asymmetric. The overall asymmetry is so strong that it also dominates the mean elevation. Apart from very high peaks close to the leeward border of the mountain range, the highest mean elevation is also achieved there, while the version with Ll=100 km featured an almost constant mean elevation over the leeward part of the mountain range. In turn, the shift of the main drainage divide is only moderate compared to the previous scenario. The leeward fraction of the total drained area has decreased from 42 % to 39 %. As a consequence, the large massifs that have formed in the northern part drain almost entirely towards the leeward side. So rivers starting from high regions partly drain towards the south first but then change their flow direction towards the large north-trending valleys.
8.2 Orographic precipitation controlling mountain range geometry
We finally consider the effect of topography on the precipitation pattern and the resulting feedback on landform evolution. The overall geometry and the parameter values are the same as before, except for a fixed length scale Ll=500 km, which was chosen in such a way that the effect of continentality over the mountain range is rather weak. In contrast to the previous examples, transversal dispersion is relevant here, where Ld=5 km was chosen.
While the vertical length scale could be defined arbitrarily based on the erodibility and the uplift rate in the previous examples, it is defined here by the reference elevation H0 that describes the decrease in β and ϵ with elevation. We use a fixed value H0=1 km and consider scenarios with different uplift rates.
The results shown in Fig. 9 reveal a distinct difference not only in mountain height, but also in mountain range asymmetry and spatial gradients in precipitation rate. As expected, mean and peak elevations increase with uplift rate. In strong contrast to scenarios of uniform precipitation, the highest mean and peak elevations are shifted towards the leeward side of the mountain range. The observed asymmetry with a gentle increase in elevation on the windward side and a strong decrease on the leeward side increases with uplift rate. In all scenarios, the highest rates of effective precipitation occur at the windward side of the mountain range, but these spatial gradients in precipitation rate increase distinctly with uplift rate. At an uplift rate of 1.5 km Myr−1, the average effective precipitation decreases by a factor of about 5 from the windward to the leeward mountain front (Fig. 9i). The principal drainage divide is shifted towards the leeward side with increasing precipitation gradient and uplift rates but to a much lesser extent than the distribution of mountain heights would suggest.
As illustrated in Fig. 10, the relationship between uplift rate and mountain height becomes nonlinear in contrast to simple scenarios assuming a uniform precipitation rate. At small uplift rates, the behavior is dominated by the increase in precipitation with increasing topography, which results in higher erosion rates. As a consequence, the topography still increases with uplift rate, but the increase is weaker than linear. For the scenario considered here, this holds for uplift rates of up to about 0.5 km Myr−1.
The concavity of the relation between uplift rate and topography is lost at higher uplift rates. The limited amount of moisture supplied from the boundary plays a central part here. As a consequence, the mean effective precipitation rate approaches a constant value for large uplift rates. However, even a constant mean precipitation rate does not imply a linear relation between uplift rate and topography since the spatial distribution of the precipitation becomes increasingly inhomogeneous. This effect is recognized in the maximum precipitation rate in Fig. 10, which continuously increases with increasing uplift.
Figure 11 illustrates the relation between topography and precipitation for the three considered uplift rates. A bimodal distribution is found at low elevations. As shown in Fig. 9c, f, and i, low elevations occur along big valleys close to the boundaries of the mountain range. Since precipitation decreases systematically from the windward side to the leeward side, the low-elevation range splits up into a rather wet windward domain and a rather dry leeward domain. This distinction is, however, lost with increasing surface elevation since intermediate elevations are distributed over the entire mountain range. This goes along with rapidly increasing variability in precipitation at given elevation.
While some decline of the increase in precipitation with elevation is already visible for U=0.5 km Myr−1, it even turns into an absolute decrease at large elevation. The highest precipitation rates are found at H≈2 km for U=1 km Myr−1 as well as for U=1.5 km Myr−1 and decrease above this elevation. This decrease is not an immediate effect of the elevation since the model itself predicts a continuous increase in precipitation with elevation at a given moisture content. As discussed above, it arises from the limited amount of moisture supplied from the windward boundary. So the occurrence of the highest precipitation rates at H≈2 km is not only related to the parameters of the precipitation model, but also to the properties of the erosion model. As a further consequence, an extreme variation in precipitation occurs at H≈2 km with high rates at the windward side and very dry regions in the shadow or high mountains where most of the available moisture has already been consumed.
The decrease in precipitation at large elevations results in a strong interaction with landform evolution. Since the erosion rate at a given channel slope decreases then, an equilibrium between uplift and erosion can only be achieved by increasing the channel slopes. Since this also requires increasing elevations, a positive feedback occurs. This feedback is visible as a convex relation between uplift rate and elevation in Fig. 10, which means that elevation increases stronger than linearly with uplift rate. This effect is particularly strong in the maximum elevation. It corresponds to the formation of very high peaks close to the leeward boundary of the mountain range (Fig. 9e and i), while the major valleys between these peaks are not particularly high.
While the precipitation pattern explains several properties of the resulting topography, we should keep in mind that the erosion rate depends on the discharge and not on the local precipitation rate. A consequence of this difference is recognized in the leeward foreland in Fig. 9. While the precipitation rate is overall low here, the topography becomes highly variable for U=1.5 km Myr−1.
Huge alluvial fans form behind the highest regions (at x≈200 km and at x≈400 km). Their occurrence is related to the low precipitation in the northern foreland. Figure 12a reveals that the respective catchments are quite small and completely in the precipitation shadow of the large massifs. This results in a very low discharge and thus in a limited ability to transport the sediment coming from the high mountain region. In equilibrium, the low discharge must be compensated for by a high channel slope, which leads to the formation of the huge alluvial fans.
The regions between these fans feature rivers with large catchments, which reach deep into the mountain range up to the principal drainage divide. Since some of the moisture coming from the windward boundary passes the principal drainage divide, parts of these catchments are exposed to high precipitation. As a consequence, the discharges are rather high, although the precipitation rate in the leeward foreland is overall low. Therefore, the respective rivers are able to carry the sediment coming from the mountain region without forming large alluvial fans.
Figure 12c shows that the windward foreland also features rivers with strongly different ratios of discharge and catchment size. While the reason for this variation is basically the same as in the leeward foreland, the variation is less pronounced here since the windward foreland region is exposed to higher precipitation than the leeward foreland. Concerning the resulting topography, however, the main difference between the two foreland regions consists of the presence of low-discharge rivers bringing sediments from the mountain range to the leeward foreland. In the windward foreland, rivers with a low ratio of discharge to catchment size typically have their source in the foreland itself and thus carry only a small amount of sediment. Therefore, the windward foreland features no big alluvial fans.
As shown in Fig. 12d, the difference in the discharge characteristics between the leeward side and the windward side is not restricted to the foreland regions. The discharge at any given catchment size A≥10 km2 varies by less than a factor of 2 at the windward side, while the variation at the leeward side is higher.
Following individual rivers downstream, the change in the discharge characteristics is opposite in both domains. It is recognized in Fig. 12c that rivers originating close to the principal drainage divide start with similar ratios of discharge to catchment size at both sides. However, tributaries at the windward side are exposed to higher precipitation so that the ratio of discharge to catchment size typically increases downstream at the windward side until the rivers reach the foreland. In turn, tributaries at the leeward side are typically quite dry and thus cause a downstream decrease in the ratio of discharge to catchment size. As a consequence, equilibrium profiles of tributaries at the leeward side are rather steep compared to the main rivers, which corresponds to a deeper incision of the large river valleys than at the windward side.
After performing some simple comparisons of the LFPM with existing models and exploring its potential in the context of landform evolution modeling, we now recapitulate what the LFPM is and what it is not.
As discussed in Sect. 6, the LFPM has some advantages over the models proposed by Smith and Barstad (2004) and by Garcia-Castellanos (2007). However, all these models share the restriction to pre-defined atmospheric flow patterns. This restriction reduces not only the numerical complexity in comparison to regional climate models, but also the predictive power.
At the actual level, models of this type are particularly useful for theoretical considerations of the co-evolution of topography and climate. There are numerous open questions not only about the overall asymmetry of mountain ranges, but also about asymmetries at smaller scales (particularly individual drainage divides, e.g., Trost et al., 2020) or the longevity of large plateaus. In a first step, such theoretical studies, in which artificial topographies are typically used under well-defined conditions, aim at a principal understanding. The second step, however, is finding out what topographic signatures can tell us about the climatic conditions in the past.
When leaving the realm of theoretical studies and sensitivity experiments and looking at real orogens in the geologic past, additional challenges arise. Large changes in topography not only affect the precipitation pattern within an orogen, but may even change the global climate. As an example, Takahashi and Battisti (2007) found that the Andes play a central part for the climatic north–south asymmetry around the Equator. So at least the input of moisture at the windward boundary and probably also the pre-defined direction of the atmospheric fluxes would have to be taken into account, e.g., by coupling the model to a general circulation model (GCM).
While general circulation models have played a central part in paleoclimate reconstructions, the recent study of Mutz and Ehlers (2019) focuses on properties that are particularly relevant for Earth surface processes. This direction of research can be seen as a first attempt to approach the co-evolution of topography and climate from large scales. So it is somehow complementary to theoretical studies and sensitivity experiments with landform evolution models. Coupling GCMs with landform evolution models and simple models of orographic precipitation such as the LFPM might become a point at which the two complementary approaches meet in the future.
This study presents a new model for orographic precipitation for use in large-scale landform evolution models such as the SPIM or the shared stream-power model. The goal was to arrive at a model that clearly goes beyond the simplest concepts, such as predicting the precipitation rate directly from surface elevation or local slope, but to stay at a level of complexity consistent with simple landform evolution models. In particular, the numerical complexity should not be much higher than that of the respective landform evolution models.
The linear feedback precipitation model (LFPM) developed in this study describes two moisture components, which are interpreted as vapor and cloud water. In contrast to previous models used in this context, a two-way conversion between the two components was assumed without considering a thermodynamic equilibrium explicitly. While this concept, in which an equilibrium evolves dynamically, seems to be more complicated first, it helps to navigate around some problems and requires few further assumptions.
As a key property, the LFPM captures a decrease in precipitation with increasing distance from the ocean (or any other source of moisture). This decrease is very slow over large continental areas with little topography but becomes faster if orographic precipitation at large mountain ranges consumes a considerable part of the available moisture. While precipitation overall increases with elevation in the model, it may decrease again at high elevations due to the limited amount of moisture. As a second important property, precipitation responds to changes in topography not instantaneously but with a finite length scale. Therefore, the model is not sensitive to the small-scale roughness of the topography and can be operated without any additional smoothing.
The length scale of the decrease in precipitation due to the limited amount of water and the length scale of the response to changes in topography can be computed from the velocity of transport in the atmosphere and the timescales of the conversion between vapor and cloud water and of the fallout of precipitation. The model structure proposed here is in principle the minimum model that is able to reproduce long-range transport and a response to changes in topography with a finite length scale.
The model also includes dispersion of the moisture fluxes in a direction perpendicular to the main transport direction. This component of the model is particularly useful in combination with two-dimensional landform evolution models since precipitation shadows of infinite length would occur behind individual peaks otherwise. Numerically, the dispersion is the most expensive part of the model. However, it can be implemented as a series of one-dimensional diffusion problems as long as the main direction of transport follows one of the principal coordinate axes. The numerical complexity is still linear then, which means that the computing effort increases only linearly with the total number of nodes of the grid. Since contemporary, fully implicit numerical schemes for the respective erosion models are also of linear complexity, this property is essential for preserving the high numerical efficiency of these models.
The model can easily be extended by a simple model of evapotranspiration wherein an elevation-dependent fraction of the precipitation is returned to the atmosphere. While this extension increases the length scale of long-range transport further, it does not change the properties of the model fundamentally.
Even the simple examples presented in this study show the remarkable impact of continentality and orographic precipitation on mountain range geometry and on the co-evolution of topography and precipitation pattern. Future studies can use this numerically efficient approach to address a wide range of research questions in the field of landscape evolution for which the assumption of uniform precipitation is too simple to explain landscape metrics and topographic patterns.
Codes and results of the simulations are available at the FreiDok data repository 219131 (https://doi.org/10.6094/UNIFR/219131, Hergarten and Robl, 2021). The open-source landform evolution model OpenLEM is freely available at http://hergarten.at/openlem (Hergarten, 2022), which also contains an implementation of the LFPM. Interested users are advised to download the most recent version of OpenLEM or the respective stand-alone version of the LFPM. The authors are happy to assist interested readers in reproducing the results and performing subsequent research.
The Video Supplement includes the respective time-dependent simulations of Sect. 8.2 and a reduced version of Fig. 12 optimized for color-blind readers. The supplement related to this article is available online at: https://doi.org/10.5194/gmd-15-2063-2022-supplement.
SH developed the theory and the main parts of the codes. JR performed the simulations and analyzed the data. Both authors wrote the paper.
The contact author has declared that neither they nor their co-author has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This open-access publication was funded by the University of Freiburg.
This paper was edited by Travis O'Brien and reviewed by Kyungrock Paik and Sebastian G. Mutz.
Barstad, I. and Schüller, F.: An extension of Smith's linear theory of orographic precipitation: introduction of vertical layers, J. Atmos. Sci., 68, 2695–2709, https://doi.org/10.1175/JAS-D-10-05016.1, 2011. a
Bookhagen, B. and Burbank, D. W.: Toward a complete Himalayan hydrological budget: Spatiotemporal distribution of snowmelt and rainfall and their impact on river discharge, J. Geophys. Res.-Earth, 115, F03019, https://doi.org/10.1029/2009JF001426, 2010. a, b
Braun, J. and Willett, S. D.: A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution, Geomorphology, 180–181, 170–179, https://doi.org/10.1016/j.geomorph.2012.10.008, 2013. a
Colberg, J. S. and Anders, A. M.: Numerical modeling of spatially-variable precipitation and passive margin escarpment evolution, Geomorphology, 207, 203–212, https://doi.org/10.1016/j.geomorph.2013.11.006, 2014. a
Deal, E. and Prasicek, G.: The sliding ice incision model: A new approach to understanding glacial landscape evolution, Geophys. Res. Lett., 48, e2020GL089263, https://doi.org/10.1029/2020GL089263, 2021. a
Garcia-Castellanos, D.: The role of climate during high plateau formation. Insights from numerical experiments, Earth Planet. Sc. Lett., 257, 372–390, https://doi.org/10.1016/j.epsl.2007.02.039, 2007. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t
Garcia-Castellanos, D. and Jiménez-Munt, I.: Topographic evolution and climate aridification during continental collision: Insights from computer simulations, PLoS ONE, 10, e0132252, https://doi.org/10.1371/journal.pone.0132252, 2015. a
Goren, L., Willett, S. D., Herman, F., and Braun, J.: Coupled numerical–analytical approach to landscape evolution modeling, Earth Surf. Proc. Land., 39, 522–545, https://doi.org/10.1002/esp.3514, 2014. a, b, c
Guerit, L., Yuan, X. P., Carretier, S., Bonnet, S., Rohais, S., Braun, J., and Rouby, D.: Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes, Geology, 47, 853–856, https://doi.org/10.1130/G46356.1, 2019. a
Han, J., Gasparini, N. M., and Johnson, J. P. L.: Measuring the imprint of orographic rainfall gradients on the morphology of steady-state numerical fluvial landscapes, Earth Surf. Proc. Land., 40, 1334–1350, https://doi.org/10.1002/esp.3723, 2015. a, b
Hergarten, S., Robl, J., and Stüwe, K.: Tectonic geomorphology at small catchment sizes – extensions of the stream-power approach and the χ method, Earth Surf. Dynam., 4, 1–9, https://doi.org/10.5194/esurf-4-1-2016, 2016. a, b
Kooi, H. and Beaumont, C.: Escarpment evolution on high-elevation rifted margins: insights derived from a surface process model that combines diffusion, advection and reaction, J. Geophys. Res., 99, 12191–12209, 1994. a
Makarieva, A. M., Gorshkov, V. G., and Li, B.-L.: Precipitation on land versus distance from the ocean: evidence for a forest pump of atmospheric moisture, Ecol. Complex., 6, 302–307, https://doi.org/10.1016/j.ecocom.2008.11.004, 2009. a
Menking, J. A., Han, J., M.Gasparini, N., and Johnson, J. P.: The effects of precipitation gradients on river profile evolution on the Big Island of Hawai'i, GSA Bull., 125, 594–608, https://doi.org/10.1130/B30625.1, 2013. a
Molnar, P. and England, P.: Late Cenozoic uplift of mountain ranges and global climate change: chicken or egg?, Nature, 346, 29–34, 1990. a
Mutz, S. G. and Ehlers, T. A.: Detection and explanation of spatiotemporal patterns in Late Cenozoic palaeoclimate change relevant to Earth surface processes, Earth Surf. Dynam., 7, 663–679, https://doi.org/10.5194/esurf-7-663-2019, 2019. a
Roe, G. H.: Orographic precipitation and the relief of mountain ranges, Annu. Rev. Earth Planet. Sci., 33, 645–671, https://doi.org/10.1146/annurev.earth.33.092203.122541, 2005. a, b
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D., and Huang, X.-Y.: A Description of the Advanced Research WRF Version 4.3, Tech. Rep. NCAR/TN-556+STR, National Center for Atmospheric Research, Boulder, Colorado, https://doi.org/10.5065/1dfh-6p97, 2021. a
Smith, R. B. and Barstad, I.: A linear theory of orographic precipitation, J. Atmos. Sci., 61, 1377–1391, https://doi.org/10.1175/1520-0469(2004)061<1377:ALTOOP>2.0.CO;2, 2004. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r
Takahashi, K. and Battisti, D. S.: Processes controlling the mean tropical pacific precipitation pattern. Part I: The Andes and the Eastern Pacific ITCZA, J. Climate, 20, 3434–3451, https://doi.org/10.1175/JCLI4198.1, 2007. a
Trost, G., Robl, J., Hergarten, S., and Neubauer, F.: The destiny of orogen-parallel streams in the Eastern Alps: the Salzach–Enns drainage system, Earth Surf. Dynam., 8, 69–85, https://doi.org/10.5194/esurf-8-69-2020, 2020. a
Turowski, J. M.: Semi-alluvial channels and sediment-flux-driven bedrock erosion, in: Gravel-bed Rivers, chap. 29, edited by: Church, M., Biron, P., and Roy, A., John Wiley & Sons, Ltd, 399–418, https://doi.org/10.1002/9781119952497.ch29, 2012. a
Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., and Tian, D.: The Generic Mapping Tools version 6, Geochem. Geophy. Geosy., 20, 5556–5564, https://doi.org/10.1029/2019GC008515, 2019. a
Willgoose, G.: Mathematical modeling of whole landscape evolution, Annu. Rev. Earth Planet. Sci., 33, 443–459, https://doi.org/10.1146/annurev.earth.33.092203.122610, 2005. a
Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, in: Tectonics, Climate, and Landscape Evolution, vol. 398 of GSA Special Papers, edited by: Willett, S. D., Hovius, N., Brandon, M. T., and Fisher, D. M., Geological Society of America, Boulder, Washington, DC, 55–74, https://doi.org/10.1130/2006.2398(04), 2006. a
Yanites, B. J. and Ehlers, T. A.: Global climate and tectonic controls on the denudation of glaciated mountains, Earth Planet. Sc. Lett., 325–326, 63–75, https://doi.org/10.1016/j.epsl.2012.01.030, 2012. a
Yuan, X. P., Braun, J., Guerit, L., Rouby, D., and Cordonnier, G.: A new efficient method to solve the stream power law model taking into account sediment deposition, J. Geophys. Res.-Earth, 124, 1346–1365, https://doi.org/10.1029/2018JF004867, 2019. a, b, c, d