1.2.1 Numerical Experimentation Mode

In this mode models are used in the same way as field and laboratory experiments. Experiments are done where components of the physics are changed to explore changes in behaviour. These changes may be changes in rate parameters or parameters that change the interactions and feedbacks, or may involve complete substitution of one set of physics for another. In this way we can quantitatively discern the impact of changes in process properties and how they change system response. The main problem with this is that, just as with field and laboratory experiments, it can be difficult to determine whether the observed changes in system response occur only for the parameters adopted for those experiments (or a particular site for field experiments). The main advantage that numerical experiments have is that it is possible to completely control the system. The user has the ability to hold some parameters steady while changing others, something that can be difficult to do for field catchments and even sometimes in laboratory experiments. By doing this it may be possible to ascertain how to design an experiment that can distinguish and test competing hypotheses for field behaviour. For instance, if changing a process (parameters, mathematical formulation and so on) does not change some characteristic we can observe and measure in the field, then, at least for this process, there is no value in performing such an experiment because the field experiment will also probably not show any significant differences. In this sense using models for numerical experimental laboratories is no different from the way science was done before computers, it is just that we can now do preliminary tests of what to expect in the field and hopefully design more robust experiments. This does not negate the value of purely curiosity-driven experiments to ‘see what happens’. However, it does impose rigour into experiments where quantitative conclusions might be possible about cause and effect. Without the model it is not possible to know if the proposed processes, when quantified, actually lead to the quantitative effect that is observed. Quantification is key.

The type of model that is well suited to numerical experimentation is one where it is easy to change the construction and/or parameters of the model. Large complex models, while having the kitchen sink of all processes built into them, tend to be less well suited to this problem because it can be hard to be confident that you fully understand all the assumptions built into the model. Also, large complex models more commonly have problems with equifinality (e.g. Hancock et al., Reference Hancock, Coulthard and Lowry2016). Equifinality (Beven, Reference Beven, Rhoads and Thorn1996) is where (1) two processes or (2) multiple combinations of parameters lead to similar responses so it is difficult to discern which process/parameters actually cause the observed behaviour.

1.2.2 Prediction Mode

In this mode models are used for making predictions of behaviour (1) at other sites or (2) for future situations. These predictions may then be compared with some societal, environmental or regulatory requirements to see if the modified system’s behaviour is satisfactory. In principle our numerical experimentation models can be used for this problem, but there are a number of subtle differences. The most important difference is that the model must be capable of predicting system behaviour for the criteria that are used for assessment and ideally provide a measure of how reliable that prediction is. In some cases this acceptance criteria may be relative to unmodified behaviour (e.g. how much better will the system be after modification, expressed as, say, a percentage or absolute change relative to no modification).

It is common to observe that the exact deterministic results from a landform evolution model are sensitive to many model inputs. While small changes may change the average behaviour of the model only a small amount (e.g. average erosion rate across a catchment), these same changes may dramatically change the location of specific features (e.g. the two simulations may have a completely different location for a gully). Thus if the acceptance criterion is a limit on the maximum erosion at a specific location, it may be difficult to demonstrate that the modified behaviour at a specific location is better, even though the average for the entire landform is, in fact, better. Moreover, if small changes in model inputs have a dramatic effect on the location of a gully, it may be impossible to convince regulators that the model predictions are useful indicators of the performance of proposed field modifications.

Finally, if project proponents (e.g. a low-level nuclear waste repository that needs to be shown to be safe for 10,000 years – current requirements in the United States and Australia for radioactive mining waste and low-level nuclear waste) are to spend large amounts of money, they need to know that the quantitative predictions of the model are accurate, or at least justifiable to regulators and defensible in court. Regulators and courts, for instance, will often want to see short-term predictions of behaviour that can be tested quantitatively against observed behaviour to assess whether the model is ‘correct’ or not.

1.2.3 Model Calibration, Validation and Uncertainty Analysis

Finally, we need to be able to test our models. This is challenging because we are dealing with systems that evolve over thousands to millions of years, and because it is unusual to have situations where we can do independent replicates of landform evolution. Thus we are left with the difficult situation that our models will not match observed data perfectly (since they are approximations of the field) and we don’t know how much difference between the models and the field is acceptable. Moreover, we may be able to calibrate our models to a specific field case (i.e. change the physics and its parameters until we get a satisfactory fit), but this demonstrates only that the model is feasible (i.e. with appropriate parameters the model can fit the data), not that the model is correct. Indeed, the hypothesis-testing literature (and a computer model is simply a hypothesis written in a computer language) indicates that we can never prove a model correct, we can only prove it incorrect. It may fit all the current data, but the data to disprove it may be just around the corner. This ability to reject a model is called falsifiability and is a key component of model testing (Popper, Reference Popper1959). A related aspect of this is that a good model makes predictions that can potentially be falsified. If the model cannot make predictions that can be falsified, then it is not possible to test it (Pfister and Kirchner, Reference Pfister and Kirchner2017).

Nor can we blindly adopt similar physics from other types of models that have been tested for different applications. For instance, there are numerous well-validated and tested agricultural erosion models (e.g. RUSLE, CREAMS, WEPP). They are based on one of the broad range of shear stress–driven detachment and sediment transport models (see Chapter 4). They have all been empirically fitted to field data so that they all satisfactorily predict erosion over relatively short time scales of a few years for paddock-scale problems. However, when applied over long time scales of millennia or more, the specific parameters used in the different transport models yield qualitatively different landforms, even though over agricultural time scales of a few years the differences between the models are small (e.g. Willgoose and Gyasi-Agyei, Reference Willgoose and Gyasi-Agyei1995). The reason for this is the long-term feedbacks between changes in topography and erosion, something that is not modelled in the agricultural erosion models, so differences become apparent only over the longer term.

One unique aspect of landform evolution models (and it is not yet clear whether this is also true of soilscape evolution models; Phillips, Reference Phillips1993) is they show sensitivity to initial conditions and external inputs so that the deterministic output (e.g. the exact values of elevation at every node in the model) can vary dramatically with only small changes in inputs or parameters (e.g. Willgoose and Gyasi-Agyei, Reference Willgoose and Gyasi-Agyei1995; Willgoose et al., Reference Willgoose, Hancock, Kuczera, Wilcock and Iverson2003). Channel and valley locations can move significant distances laterally for only small changes in model inputs. Thus it is not possible to make a direct comparison of the exact values at each location of the model and the field data. Rather it is necessary to identify statistics of the landforms that will be the same if the physics is modelled correctly and then compare these statistics. One might then ask, What is the set of essential statistics that is required to completely describe a landscape? This set of essential statistics can then form the focus of any model testing.

Key outputs needed for model testing are confidence limits on model predictions or confidence limits on field measurements. For instance, the t-test is commonly used for testing the differences between models and experiments. The t-statistic is the difference between the two experiments divided by the variability of the experiments, and this variability can be determined from the confidence limits. As we have mentioned, field replication of landscapes is difficult, so it is rarely possible to calculate confidence limits on field measurements. An alternative is to calculate confidence limits on model output. In principle this is simple if we have accuracy estimates on model parameters and exogenic inputs. We can use Monte Carlo simulation with the model. This process is to repeatedly run the model with different parameters, where the different parameters are determined by sampling the parameter values from the probability distribution of the model parameters. In this way we can generate a number of output simulations, and we can derive the probability distribution of the outputs for the model and thus determine confidence limits on the predictions (e.g. Willgoose et al., Reference Willgoose, Hancock, Kuczera, Wilcock and Iverson2003). Likewise we can explore the effect of unknown climate by running simulations with randomly generated climate data to explore the effect of climate variability. The problem with this process is that it can be extremely computer intensive. Depending on how many parameters need to be varied, it is common for thousands of simulations to be necessary.

This problem of excessive computer time is even worse for satisfying many regulatory requirements in engineering applications because it is common that acceptance criteria are expressed as a small probability of some criterion or criteria being exceeded (e.g. one in one million chance of failure in any one year is common for high-consequence industrial applications). Thus not only is it necessary to do Monte Carlo simulation to obtain the probability distribution of the output, but the small probability required for acceptance means that the results we are trying to determine are on the tail of the probability distribution. The further out we are on the tail of the probability distribution, the greater is the number of simulations that are required to adequately sample enough rare events to define the tail of the probability distribution.

The key conclusion from this is that to be able to provide confidence limits on model simulations it is crucial to have fast models. Thus it is not simply enough to be able to simulate some processes, but it is also necessary to be able to simulate these processes efficiently so that confidence limits can be generated in a reasonable compute time. The book will return to this issue of algorithmic efficiency repeatedly.

1.3.1 Continental Scale

Figure 1.1 overviews the important processes at the continental scale. At this length scale (10–100 km) the weight of the land on the crust has significant interactions with the tectonics of mountain building. The loss of mass by erosion lightens the load on the crust and the crust rises in response by isostacy (Section 12.2). If mountains are being created by convergence of two plates, then this unloading by erosion also has significant feedbacks on the rate and nature of the interaction between the two plates underneath the emerging mountain range (Section 12.3). An important interaction is that the topography of the mountain range changes the spatial distribution of rainfall (Section 3.1.1) with higher rainfall on the windward side of mountains than the leeward side. However, not only does this increased rainfall increase runoff (and everything else being equal increases erosion), but it also increases the vegetation density (increased density decreases erosion), and it is the competing effects of increased runoff and increased vegetation that determines whether the erosion rate is higher on the rainy side or the drier side (Section 16.1.2).

Figure 1.1: A schematic of the processes that occur at the continental scale including mountain building, isostatic rebound, river erosion, spatially variable rainfall and deep marine deposition. The heavy grey arrows are process fluxes and their predominant direction of action. The two circular arrows represent the feedbacks between (left) erosional reduction of elevation and rising in the crust, and (right) deep-sea deposition and depression of the crust.

The rivers (and the hillslopes linked to them) are the transport mechanism to remove sediment from the mountain ranges to the ocean (Chapter 4); these river network and hillslopes will be discussed below.

While this book will not cover marine processes, the fate of the eroded sediment is shown in Figure 1.1. Most of the particulate and dissolved load is deposited either on the continental shelf (either in the delta at the river mouth or farther offshore) or in basins in the deep ocean beyond the shelf. The extra load from this deposited material results in the crust being depressed downward; that is, not only does erosion on the land interact with the crust and mantle so the crust rises, but the deposited sediment in the ocean also interacts with the crust and mantle, so the crust falls.

The continental shelf is an important transition region between terrestrial processes and marine processes. The depth below current sea level of the ocean edge of continental shelves worldwide is of the order of 140 m, which coincides with the level of the ocean during the ice ages of the last million years. Ice core data indicate that about 80% of the last million years has been ice age, while the current, interglacial, conditions (since about 10,000 years ago) have prevailed for less than 20% of this period. Thus the topography of the continental shelf regions around many continents will largely reflect terrestrial processes, with about 10,000 years of recent marine processes overprinted on them.

Some final points to note about Figure 1.1 are the following:

1. 1. The arrows for uplift of the crust under the mountain range are not vertical, because mountain building involves convergence of the crust as well as vertical uplift. Moreover there is the possibility of relative lateral motion horizontally between the converging plates (e.g. the San Andreas fault in California).

2. 2. The rainfall is spatially and temporally variable to reflect the effect of the changing topography on rainfall. Note that this spatial distribution of rainfall may vary seasonally, and year to year, as a result of changes in the prevailing wind directions. Some signature research sites have broadly consistent prevailing wind directions (e.g. Taiwan, New Zealand, Chile), while for other sites wind directions even today vary in direction seasonally as a result of seasonal monsoons (e.g. the Himalayas).

1.3.2 Catchment and River Network Scale

Figure 1.2 overviews the important processes at the catchment and river network scale. At the continental scale the main transport process for sediment is the channel network. Most of the sediment that is transported by the river network comes from the hillslopes, either by (1) mass movement (e.g. landslides) in parts of the river network that are steep or (2) by fluvial erosion and soil creep in flatter areas (Figure 1.2a).

Figure 1.2: A schematic of the processes and interactions in channel networks at the catchment scale: (a) a plan of the river network showing the ‘string and beads’ arrangement of channels and floodplains, (b) a long section along a river showing the periodic change in channel type from bedrock to alluvial and back again, (c) a typical alluvial channel/floodplain cross section at low flow showing the lateral movement of the channel (left to right in this figure) within the floodplain as a result of meandering and (d) a typical alluvial channel/floodplain cross section at high flow showing the excavation of the channel cross section down to and exposing the underlying bedrock, and the deposition of particulate matter on the floodplain.

Broadly speaking the river network consists of two types of channels (Sections 4.4 and 4.5): (1) bedrock channels where the river flows across bare exposed bedrock and (2) alluvial channels where the river flows through deposited sediment. These two channel types are end members in a transition from bedrock channels to alluvial channels as the sediment load in the river increases, but this distinction suffices for the moment. Alluvial channels occur in regions of deposited sediment called floodplains and typically meander from side to side within this floodplain. Bedrock channels are typically constrained laterally by the erosion-resistant bedrock channel sides and so tend to be straight rather than meandering. The floodplain is a result of net deposition of sediment on the floodplain by the alluvial channel (typically at high flows), and these sediments are remobilised as the river meanders from side to side (typically at lower flows). The residence time of the sediment within the floodplain is measured in many thousands of years. Recent work has also indicated that floodplains are a major storage site for organic carbon (Section 10.7).

Thus the process of erosion of the mountains is one of delivery of sediments from the hillslopes to the channel. This sediment (both mineral and organic matter) then passes through a mixture of bedrock and alluvial channels on its way to the ocean. In the alluvial channels there is an exchange of sediment between the channel and the floodplain deposits with some of the channel sediment being deposited in the floodplain, while some of the floodplain materials are remobilised from the floodplain deposits (Sections 4.5 and 10.7). Typically the finest, lightest and most mobile sediments will pass through the floodplain without being exchanged with the floodplain sediments (Section 4.2.1.3). This leads to what is commonly referred to as a ‘strings and beads’ plan of the river network where the strings are channel reaches without adjacent floodplains and the beads are river reaches with their adjacent floodplains.

Figure 1.2b shows a schematic of a river long profile showing the transition from a bedrock channel reach to an alluvial channel/floodplain reach and back again. For low sediment loads the bottom of the bedrock channel will be clear of sediment so bedrock incision will erode the channel. For high sediment loads the bottom of the channel in the bedrock reach will be periodically covered by sediment waves and pulses that move down the reach. Thus at high sediment loads the bottom of the bedrock channel will be partially protected by sediment. Incision can occur only when the bedrock channel base is not covered by sediment. Thus the bedrock channel will incise into the bedrock periodically, with the percentage of time in the bedrock incision mode being a function of the percentage of time that the bedrock bed is not covered by sediment. Likewise for the alluvial channel there will localised stretches where the bedrock underneath the floodplain is exposed. Moreover, during high flows the channel is scoured and can be excavated down to bedrock (Figure 1.2d). Thus periodically an alluvial channel can be incising bedrock (and thus the bedrock under the floodplain is lowered), though for most of the time and for most locations the channel is largely in balance with the meandering channel excavating the floodplain sediment on the outside of the meander bend and depositing inside the bend, so that there is an exchange of sediment and organic matter between the channel and the floodplain (Figure 1.2c). At high flow, sediment and organic matter are deposited on the floodplain (Figure 1.2d). The residence time in the floodplain deposits is then a function of the rate of floodplain excavation by meandering relative to the width of the floodplain. The residence time is a key parameter for organic matter sequestration in the floodplain because it determines how much of the deposited organic matter decomposes in situ versus being remobilised by bank erosion at some time after its original deposition.

1.3.3 Hillslope Scale

Figure 1.3 overviews the important processes at the hillslope scale. The four panels overview different components of the hillslope processes: (a) mineral transport, (b) soil organic matter and carbon transport, (c) water balance and (d) vegetation and fire. The hillslopes deliver the material that is then transported to the coast by the rivers. Most of the material transported by the rivers has been sourced in some fashion from the hillslopes, so the characteristics of the material eroded from the hillslope determine the material in the river and the rate of river transport.

Figure 1.3: A schematic showing two hillslopes draining down to a central channel. The four panels each show the processes operating for the four main components of the hillslope system: (a) soil mineral matter, (b) soil organic matter, (c) the water balance and (d) vegetation and fire. The heavy grey arrows show process fluxes and their predominant direction of action. The aspect label indicates those processes and properties that are different between polar- and equatorial-facing hillslopes for mid-latitude catchments.

Starting with the mineral matter (Figure 1.3a) the main transport mechanisms are fluvial erosion (Chapter 4), soil creep (Chapter 9) and mass movement (Chapter 13). The fluvial erosion is largely a function of runoff generation from rainfall (Chapter 3). The main modifier of the fluvial erosion rate is the groundcover provided by vegetation and its leaf litter (Chapter 14) though the canopy also provides some protection, but it is less important than groundcover. The dynamics of soil production and weathering influence the characteristics of the sediment on the surface of the soil. Soil production is the process of conversion of the underlying saprolite into soil, and its rate is a function of the thickness of the soil (Chapter 6). Within the soil profile a range of weathering processes are active that are (1) breaking the soil fragments into smaller and smaller fragments (Chapter 7) and (2) chemically transforming the soil rock fragments into secondary minerals (mostly clays) and leachate (that is transported out of the hillslope by water as dissolved load) (Chapter 8). The soil is vertically mixed by biological processes called bioturbation (Sections 7.3.5, 8.3.4). The roots of vegetation provide much of the strength of the hillslope to resist landsliding and other forms of mass movement (Chapters 13 and 14). Finally at mid-latitudes the aspect of the hillslopes can make a significant difference to vegetation cover, so many of the transport processes above vary depending on whether the hillslope is polar-facing or equatorial-facing.

Figure 1.3b shows the organic matter cycle on the hillslope. The main processes are vegetation providing organic matter, in the form of leaf litter, which is incorporated into the soil profile (Chapters 10 and 14). Much organic matter is generated within the profile from plants roots (the rhizosphere), mycorrhizae hyphae and microbiology. Mycorrhizae hyphae (i.e. fungi roots) are important in developing soil structure and increasing soil infiltration rates, and by mass are the most significant component of soil organic matter (SOM). In the same way that bioturbation mixes the mineral matter from top to bottom in the soil profile, so bioturbation mixes the SOM. Within the profile the SOM decomposes over time, and the decomposition rate is initially fast for fresh SOM as the most reactive components decompose. The decomposition rate slows over time since the remaining SOM contains a greater percentage of components more resistant to decomposition (Chapter 10). SOM is preferentially eroded (i.e. sediment is enriched with SOM) at the surface relative to the mineral matter because it is less dense than the mineral matter. Finally the decomposition of the SOM (as well as respiration from vegetation roots; see below) generates carbon dioxide within the soil profile. Carbon dioxide levels within the soil profile can be very high (near 100%), and when this dissolves in water it generates acid that drives chemical weathering of the soil mineral matter (Chapter 8). The main impacts of aspect on the organic matter cycle are in changing vegetation density (and thus the generation rate for SOM source material) and soil temperature (the SOM decomposition rate increases strongly with soil temperature).

The water cycle is illustrated in Figure 1.3c. For most of the processes in this book the main characteristic of importance is the soil moisture because it influences (1) the growth rate for plants, (2) the amount of soil organic matter stored in the soil (wetter soils biodegrade faster), (3) the rate of chemical weathering in the soils and (4) the initial wetness of the soil when rainfall occurs so that wetter soil generates more runoff. The main sources that increase the soil moisture in the water cycle are rainfall and infiltration. The main sinks that decrease soil moisture are (1) bare soil evaporation, (2) plant transpiration (a function of plant growth rate, Chapter 14), (3) surface water runoff (a driver of fluvial erosion, Chapter 4) and (4) groundwater flow and leachate transport (Chapter 8). Higher soil moisture is associated with higher plant growth rate, surface runoff, creep and landsliding rates and lower fire occurrence. The main impacts of aspect are (1) the infiltration rate is changed by changes in soil structure due to SOM and soil grading and (2) higher transpiration from vegetation on equatorial-facing slopes.

Finally, the vegetation and fire cycles are summarised in Figure 1.3d. The main feedbacks from vegetation are in (1) plant roots, which increase the strength of the soil making it more resistant to landsliding, (2) groundcover and leaf litter, which provide soil surface protection against fluvial erosion, and (3) root respiration, which generates carbon dioxide within the soil profile, a source of acid active in chemical weathering. The main impacts of fire on landscape evolution are by modifying the impacts of vegetation. Fire removes groundcover vegetation and leaf litter layer and, if intense enough, the canopy vegetation as well. The loss of ground protection by the removal of groundcover and leaf litter results in fluvial erosion increasing dramatically until the vegetation recovers in a year or two. A further fire impact is that the SOM in the top of the soil is volatised so that the soil structure created by SOM (mostly mycorrhizae hyphae and polysaccharides) is destroyed. This results in reduced infiltration, increased runoff and increased erosion. For those trees killed by fire their root strength is reduced so that the resistance to mass movement and creep provided by plant roots is reduced, and mass movement rates increase for a decade or so after fire. The charcoal that is created by fire is incorporated into the soil, and since it is decomposition resistant, it is effectively sequestered forever within the soil. The main impact of aspect on vegetation is that the differences in vegetation and soil moisture lead to different fire recurrence rates and intensity. The denser the vegetation the higher the fire intensity, while the lower the soil moisture the dryer the fuel (e.g. leaf litter), leading to higher fire intensity.

What should be clear from the complexity of Figure 1.3 is that many factors on the hillslope interact and may result in either positive or negative feedbacks. Thus modelling the evolution of only the mineral matter processes (as was common for the early landform evolution models) independently of the soil organic matter, vegetation, fire and water is modelling only one part of a highly connected story. The challenge is that some parts of the processes in Figures 1.1–1.3 are quite difficult to quantify, so modelling them is difficult. This difficulty is no reason to ignore them, however, and the aim of this book is to provide an overarching framework within which we can quantify these processes and test their impacts on landscape evolution.