A coupled soilscape-landform evolution model: Model formulation

22 This paper describes the coupling of the State Space Soil Production and Assessment Model 23 (SSSPAM) soilscape evolution model with a landform evolution model to integrate soil profile 24 dynamics and landform evolution. SSSPAM is a computationally efficient soil evolution model which 25 was formulated by generalising the mARM3D modelling framework to further explore the soil profile 26 self-organization in space and time, and its dynamic evolution. The landform evolution was integrated 27 into SSSPAM by incorporating the processes of deposition and elevation changes resulting from 28 erosion and deposition. The complexities of the physically based process equations were simplified by 29 introducing state-space matrix methodology that allows efficient simulation of mechanistically linked 30 landscape and pedogenesis processes for catena spatial scales. The modelling approach and the 31 physics underpinning the modelled processes are described in detail. SSSPAM explicitly describes the 32 particle size grading of the entire soil profile at different soil depths, tracks the sediment grading of 33 the flow, and calculates the elevation difference caused by erosion and deposition at every point in the 34 soilscape at each time step. The landform evolution model allows the landform to change in response 35 to (1) erosion and deposition, and (2) spatial organisation of the co-evolving soils. This allows 36 comprehensive analysis of soil landform interactions and soil self-organization. SSSPAM simulates 37 fluvial erosion, armouring, physical weathering, and sediment deposition. The modular nature of the 38 SSSPAM framework allows integration of other pedogenesis processes in follow-on research projects. 39 This paper presents the initial results of soil profile evolution on a dynamic landform. These 40 simulations were carried out on a simple linear hillslope to understand the relationships between soil 41 characteristics and the geomorphic attributes (e.g. slope, area). Process interactions which lead to such 42 relationships were also identified. The influence of the depth dependent weathering function on 43 soilscape and landform evolution was also explored. These simulations show that the balance between 44 erosion rate and sediment load in the flow accounts for the variability in spatial soil characteristics 45 while the depth dependent weathering function has a major influence on soil formation and landform 46 evolution. 47


Abstract 22
This paper describes the coupling of the State Space Soil Production and Assessment Model 23 (SSSPAM) soilscape evolution model with a landform evolution model to integrate soil profile 24 dynamics and landform evolution. SSSPAM is a computationally efficient soil evolution model which 25 was formulated by generalising the mARM3D modelling framework to further explore the soil profile 26 self-organization in space and time, and its dynamic evolution. The landform evolution was integrated 27 into SSSPAM by incorporating the processes of deposition and elevation changes resulting from erosion 28 and deposition. The complexities of the physically based process equations were simplified by 29 introducing state-space matrix methodology that allows efficient simulation of mechanistically linked 30 landscape and pedogenesis processes for catena spatial scales. SSSPAM explicitly describes the particle 31 size grading of the entire soil profile at different soil depths, tracks the sediment grading of the flow, 32 and calculates the elevation difference caused by erosion and deposition at every point in the soilscape 33 at each time step. The landform evolution model allows the landform to change in response to (1) 34 erosion and deposition, and (2) spatial organisation of the co-evolving soils. This allows comprehensive 35 analysis of soil landform interactions and soil self-organization. SSSPAM simulates fluvial erosion, 36 armouring, physical weathering, and sediment deposition. The modular nature of the SSSPAM 37 framework allows integration of other pedogenesis processes to be easily incorporated. This paper 38 presents the initial results of soil profile evolution on a dynamic landform. These simulations were 39 carried out on a simple linear hillslope to understand the relationships between soil characteristics and 40 the geomorphic attributes (e.g. slope, area). Process interactions which lead to such relationships were 41 also identified. The influence of the depth dependent weathering function on soilscape and landform 42 evolution was also explored. These simulations show that the balance between erosion rate and 43 sediment load in the flow accounts for the variability in spatial soil characteristics while the depth 44 dependent weathering function has a major influence on soil formation and landform evolution. The 45 results demonstrate the ability of SSSPAM to explore hillslope and catchment scale soil and landscape 46 evolution in a coupled framework. Soil is one of the most important substances found on planet Earth. As the uppermost 55 layer of the earth surface, soil supports all the terrestrial organisms ranging from microbes to 56 plants to humans and provides the substrate for terrestrial life [Lin, 2011]. Soil provides a 57 transport and a storage medium for water and gases (e.g. carbon dioxide which influences the 58 global climate) [Strahler and Strahler, 2006]. The nature of the soil heavily influences both 59 geomorphological and hydrological processes [Bryan, 2000]. In addition to the importance of 60 soil from an environmental standpoint, it provides a basis for human civilization and played an 61 important role in its advancement through the means of agricultural development [Jenny, 62 1941]. Understanding the formation and the global distribution of soil (and its functional 63 properties) is imperative in the quest for sustainable use of this resource. 64 Characterization of soil properties at a global scale by sampling and analysis is time 65 consuming and prohibitively expensive due to the dynamic nature of the soil system and its 66 complexity [Hillel, 1982]. However over the years researchers have found strong links between 67 different soil properties and the geomorphology of the landform on which they reside [Gessler 68 et al. 2000,1995]. Working on this relationship several statistical methods have been 69 developed to determine and map various soil properties depending on other soil properties and 70 geomorphology such as pedotransfer functions, geostatistical approaches, and state-factor (e.g. 71 clorpt) approaches [Behrens and Scholten, 2006]. Pedotransfer functions (PTFs) use easily 72 measurable soil attributes such as particle size distribution, amount of organic matter, and clay 73 content to predict hard to measure soil properties such as soil water content. Although very 74 useful, PTFs need a large database of spatially distributed soil property data and require site 75 specific calibration [Benites et al., 2007]. Geostatistical methods use a finite number of field 76 samples to interpolate the soil property distribution over a large area. Developing soil property 77 maps using geostatistical methods is possible for smaller spatial scales, however soil sampling 78 and mapping soil attributes can be prohibitively expensive and time consuming for larger 79 spatial domains [Scull et al., 2003]. State-factor methods, such as scorpan (developed by 80 introducing existing soil types and geographical position to clorpt framework) use digitized 81 existing soil maps and easily measurable soil attribute data to generate spatially distributed soil 82 property data using mathematical concepts such as fuzzy set theory, artificial neural network 83 or decision tree methods [McBratney et al., 2003]. However these techniques also suffer from 84 scalability issues and the typical need for site specific calibration. 85 While spatial mapping of soil properties is important, understanding the evolution of 86 these soil properties and processes responsible for observed spatial variability of soil properties 87 is also important. In order to quantify these processes and predict the soil characteristics 88 evolution through time, dynamic process based models are required [Hoosbeek and Bryant,89 1992]. These mechanistic process models predict soil properties using both geomorphological 90 attributes and various physical processes such as weathering, erosion, and bioturbation 91 [Minasny and McBratney, 1999] Since only three layers were used in MILESD, the representation of the particle size 139 distribution down the soil profile was limited. Although LORICA incorporated additional soil 140 layers into the MILESD modelling framework, detailed exploration of soil profile evolution or 141 interactions between landform evolution and soil profile evolution has not yet been done with 142 this model. Importantly, particle size distribution of the soil can be used as a proxy for various 143 soil attributes such as the soil moisture content [Arya and Paris, 1981;Schaap et al., 2001].

144
The main objective of this paper is to present a new soilscape evolution model capable of 145 predicting the particle size distribution of the entire soil profile by integrating a previously 146 developed soil grading evolution model in to a landscape evolution model.

147
Here we present the methodology for incorporating sediment transport, deposition and

156
The introduction of a landform into the SSSPAM framework is done using a digital 157 elevation model. The structure of the landform evolution model follows that for transport-158 limited erosion [Willgoose et al., 1991] but modified so as to facilitate its coupling with the 159 soilscape soil grading evolution model SSSPAM described in [Welivitiya et al., 2016]. Here a 160 regular square grid digital elevation model was used and converted into a two dimensional 161 array which can be easily processed and analysed in the Python/Cython programing language.

162
Using the "steepest-slope" criterion [Tarboton, 1997]   The soil profile evolution of each pixel is determined using the interactions between 167 the soil profile and the flowing water at the surface. Figure 1 shows these layers and their 168 potential interactions. This is similar to the schematic for the standalone soil grading evolution 169 model but is different in that the erosion/deposition at the surface is a result of the imbalance 170 between upslope and downslope sediment transport. The water layer acts as the medium in 171 which soil particle entrainment or deposition occurs depending on the transport capacity of the 172 water at that pixel. The water provides the lateral coupling across the landform, by the sediment 173 transport process. The soil profile is modelled as several layers to reflect the fact that the soil 174 grading changes with soil depth depending on the weathering characteristics of soil. Erosion where r is runoff excess generation (m 3 /s/m 2 ) and is contributing area (m 2 ) of that pixel.

Sediment deposition 267
If the total mass of incoming sediment is higher than the transport capacity of the 268 sediment transport capacity at the pixel (Table 1 shows the need to have this adjustment vector and the method we used to calculate it.

285
Consider the example values given in Table 2. The total mass of the incoming sediments value (which is 24) and rescaling these values with D (total deposition mass), we can calculate 291 the masses of sediment which need to be deposited from each grading class. In some cases 292 (when the total deposition D is higher than the ∑ , value) the mass of material which 293 needs to be deposited can be larger than the available sediments in that particular size class. In  Depending on the position of the sediment particle entering into the pixel with respect 332 to critical immersion depth, whether or not that particle will deposit in that pixel can be 333 determined. Particles entering to the pixel below the critical immersion depth will settle within 334 the current pixel, while particles entering above the critical immersion depth will stay in 335 suspension and exit the current pixel. The critical immersion depth is greater for larger (or 336 denser) particles and less for smaller (or less dense) particles. For sediment particles in larger 337 size classes, the critical immersion depth can be larger than the flow depth (m) (thickness 338 of the water column). That means all the particles in that particle size class will settle in the 339 pixel. Using the critical immersion depth and the flow depth we can define the diagonal 340 elements , of the deposition transition matrix in following manner.
Note the deposition transition matrix is a diagonal matrix which contains only  The following discussion briefly describes the methodology used to calculate the above 347 variables. The average settling velocity of all the particle sizes classes can be calculated for 348 typical sediment sizes using Stoke's Law [Lerman, 1979].  Figure 4 shows three different cases that can occur during the deposition process. In

395
Case 1 (Figure 4(b)) the deposition height Δℎ is less than the surface armour thickness . The methodology used for simulating weathering within the soil profile is detailed by 412 Welivitiya et al. [2016]. It uses a physical fragmentation mechanism where a parent particle 413 disintegrates into n number of daughter particles with a single daughter particle retaining 414 fraction of the parent particle by volume and the remaining n-1daughter particles retaining 415 fraction 1 − of the parent particle volume. By changing n and we can simulate a wide

464
However as the number of model layers increases, the time required for the each simulation 465 also increases. During our initial testing, we found that the soil depth rarely increased beyond 466 10 m and decided to set 10.05 m as the maximum soil depth for this scenario.

467
Two soil grading data sets (  Figure 6 shows six outputs at different times during hillslope and soil profile evolution.

486
The upper section in each of the panels in Figure 6 is the cross-section median diameter

491
The soil depth is the depth below the surface at which d50 reaches the maximum possible 492 particle size (i.e. the bedrock grading). Figure 6(a) shows the initial condition for the soilscape:   in to upper footslope and lower toeslope. Generally the soil grading distribution is fine at the 563 summit, coarsens from the summit to the shoulder and backslope followed by fining from 564 backslope to the base [Birkeland, 1984]. Furthermore, the soil depth is typically high in summit

568
In order to better understand the dynamics of soilscape evolution we also plotted the 569 elevation, slope, rate of erosion (and/or deposition), surface d50, soil depth and profile d50 for 570 four sites (Figure 6(a)). The first two sites (sites 1 and 2) are either side of the plateau-hillslope  When the erosion front crosses site 1, the gradient increases as does the erosion rate (at 591 around 20,000 years). During this phase of increasing erosion the surface d50 also increased.

592
However, the surface d50 stabilizes around 14 mm before the erosion rate reaches its maximum 593 value. This is because once total armouring occurs, the erosion is reduced to a very low value.

594
Although the erosion is low, the slope of the site 1 continues to increase until it reaches a weathering modelling work done byKirkby [Kirkby, 1977;1985;. However we will be 682 incorporating a physically based chemical weathering model described by Willgoose [2018] 683 into SSSAPAM in the future. All available evidence suggests that in order to effectively model The deposition model of SSSPAM is designed in such a way that the difference between 689 the transport capacity and the sediment load of the flow is always deposited regardless of the 690 settling velocities. This is done to prevent the flow from being over the transport capacity. weathering mechanism which leads to the heterogeneity of soilscape properties in a region. A 768 future paper will discuss how this work can be extended to include the impact of chemical 769 weathering into soilscape evolution.