Interactive comment on “ Exploring the sensitivity on a soil area-slope-grading relationship to changes in process parameters using a pedogenesis model ” by

The authors present a model (SSSPAM5D) to investigate the relation between soil grading, surface erosion and weathering. Minor changes to an existing model (mARM3D) are done and a number a 2D simulations have been performed at the soilscape level in the framework of a sensitivity analysis. The main conclusion is that an earlier published area-slope-grading relationship holds when model parameters are altered. Interestingly, the sensitivity analysis shows a dependency of soil grading on the ratio between area and slope exponents, similar to what has been found for other earth surface processes.


Introduction
Soil is a product of various physical processes acting on earth's crust.Weathering is a major contributor to soil production, along with transport processes that transport new material away and bring new material into a point.Weathering is a general term used to describe all the processes which cause rocks or rock fragments to disintegrate or alter through physical (Ollier, 1984;Wells et al., 2006Wells et al., , 2008;;Yokoyama and Matsukura, 2006), chemical (Green et al., 2006;Ollier, 1984) or biological means (Strahler and Strahler, 2006).Disintegration of rock material through physical weathering can occur by (1) unloading, (2) expansion and contraction of rock through heating and cooling cycles, (3) stress developing in rock fractures due to freezing water, (4) salt crystal growth or tree root intrusions, and (5) abrasion of rock by harder materials transported by flowing water or glaciers (Thornbury, 1969).Physical weathering where larger soil particles are broken down into smaller particles is dominant in the surface layer of material where it is more exposed.Weathering also occurs underneath the surface and the weathering rate at these subsurface layers can be modelled with depthdependent weathering functions.
Spatial redistribution of soil can occur due to different processes such as soil creep and erosion.Soil creep is the process of downslope movement of soil over a low grade slope with a substantial soil mantle under the force of gravity and friction (Ollier and Pain, 1996).Although soil creep can have a significant influence on some soil properties on some landforms (Braun et al., 2001;Roering et al., 2007;West et al., 2014) on landforms with interlocking rock fragments, its influence is not significant.On the other hand erosion can occur in all landforms in one form or another.Erosion is a term used for removal of material from an existing soil profile.Erosion can occur due to a number of processes such as (1) surface water flow (fluvial erosion), (2) wind (aeolian erosion), (3) flow of glaciers (glacial erosion) and (4) animal or plant activity (biological erosion) and others.Fluvial and Aeolian erosion tend to create an "Armour" on the soil surface.Depending on the energy of the erosion medium (water or air), transportable fine particles are preferentially entrained and transported from the surface soil layer.This process coarsens the remaining surface soil layer enriching it with coarser, less mobile, material.With time, if the energy of the transport medium remains constant, an armoured layer is formed with all the transportable material removed.At this time the sediment transport reaches zero.This armour, where all the materials are larger than the largest grains which the transport medium can entrain, prevents erosion of material from the subsurface.If the energy of the transport medium increases, the existing armour can be disrupted, and a newer stable armour with coarser material can be formed (Sharmeen and Willgoose, 2006).Armouring in river beds has been widely understood and studied extensively for mostly streams and rivers (Gessler, 1970;Gomez, 1983;Lisle and Madej, 1992;Little and Mayer, 1976;Parker and Klingeman, 1982).
The importance of soil as an agricultural and commercial resource, and as an influencing factor on environmental processes such as climate regulation, is well established (Jenny, 1941;Bryan, 2000;Strahler and Strahler, 2006;Lin, 2011).However, spatially distributed quantification of soil properties is difficult because of the complexity and dynamic nature of the soil system itself (Hillel, 1982).The necessity for quantified and spatially distributed soil functional properties is clear (Behrens and Scholten, 2006;McBratney et al., 2003).Moreover, explicit soil representation in models of environmental processes and systems (e.g.landform evolution, and hydrology models) has increased rapidly in the last few decades.For accurate prediction these physically based and spatially explicit models demand high-quality spatially distributed soil attributes such as hydraulic conductivity (McBratney et al., 2003).
The need for improved soil data arises in two main areas: (1) better mapping of the description of the soil (e.g.particle size distribution, soils classification), and (2) improved representation of soil functional properties (e.g.hydraulic conductivity, water holding capacity).For most environmental models the soil functional properties are of greatest interest since they determine the pathways and rates of environmental process.Accordingly this paper is focussed on a soil representation that can underpin the derivation of functional properties.Pedotransfer functions exist (albeit with large uncertainty bounds) to then relate these soil descriptions to functional properties.The existence of these pedotransfer functions intellectually underpins the rationale of the work in this paper.While these techniques are not the focus of this paper, some discussion of them is pertinent so that the importance of the scaling relationship discussed in this paper can be fully appreciated.
Traditional soil mapping typically uses field sampling and classifies soils into different categories based on a mixture of quantitative (e.g.pH) and qualitative features (e.g.colour).It does not directly provide the functional soil properties required by environmental models.Several techniques have been introduced to tackle this lack of functional description such as pedotransfer functions, geostatistical approaches, and state-factor (Clorpt) approaches (Behrens and Scholten, 2006).Pedotransfer functions (PTFs) have been developed to predict functional soil properties using easily measurable soil properties such as particle size grading, organic content, and clay content.Although PTFs are very useful, they are limited because they need spatially distributed soil descriptions and, in many cases, site-specific calibration (Benites et al., 2007).Geostatistical approaches interpolate field data to create soilattribute maps.Clorpt or Scorpan approaches (McBratney et al., 2003) use regression or fuzzy-set theory to create soilattribute or soil-class maps (Behrens and Scholten, 2006).
Geostatistical digital soil mapping using field sampling of soil is possible for a specific site where the area is small (Scull et al., 2003).However, it can be prohibitively expensive and time consuming for larger sites.Soil mapping techniques, such as Clorpt or Scorpan, use digitisation of existing soil maps.They generate soil classes through decision tree methods and artificial neural networks using easily measurable soil attributes (similar to PTFs) to generate the digital soil maps (McBratney et al., 2003).Although much work has been carried out these methods also suffer the need for site-specific calibration.
Remote-sensing technologies such as gamma ray spectroscopy have introduced novel methods of characterizing soil properties and developing digital soil maps (Triantafilis et al., 2013;Wilford, 2012).The digital soil maps produced by gamma ray spectroscopy are relatively coarse and their spatial coverage is limited while their links with functional properties remain uncertain (McBratney et al., 2003).Developments in geographic information systems (GIS) have enabled fast and efficient characterization and analysis of large amounts of spatial and non-spatial data (Scull et al., 2003).The ease of use of GIS has revolutionised modelling by making distributed modelling easier to do and interpret (Singh and Woolhiser, 2002).This is the rationale for the Global-SoilMap initiative, which aims to provide a global 90 m map of soil properties for the world (Sanchez et al., 2009).
Many researches have reported strong relationships between terrain attributes and soil.For example, using field measurements Moore et al. (1993) found significant correlations between terrain attributes and soil properties such as soil wetness index and soil organic carbon content.Poesen et al. (1998) reported a strong relationship between the slope gradient and the rock fragment size on the soil surface.Statistical (Gessler et al., 1995(Gessler et al., , 2000) ) and process-based models (Govers et al., 2006) have been proposed to predict these re-lationships.They have been implemented to predict the soil attributes data using terrain attributes as a proxy.ARMOUR (Sharmeen and Willgoose, 2006) is a physically based model that simulates (1) rainfall-runoff event overland flow, (2) erosion and the selective entrainment of fine sediments that creates armouring of the soil surface, and (3) weathering of the particles on the surface that breaks down the armour.AR-MOUR simulates the evolution of the soil surface on a hillslope.However, the very high computing resources and long run times (1000s of years at minute time resolution) of the physically based modelling prevented the coupling of AR-MOUR with a hillslope evolution model.Cohen et al. (2009) developed a state-space matrix soils model, mARM, and calibrated it to output from ARMOUR.mARM was significantly more computationally efficient than ARMOUR, and was able to simulate more complex hillslope geometries.It was sufficiently fast that it could be used to simulate the spatial distribution of the soil profile as well the surface properties.By incorporating the weathering characteristics of soil profile into mARM, Cohen et al. (2010) developed mARM3D, which was able to explore the evolution of soil profiles for small catchments.Cohen et al. (2009) was the first to identify using pedogenic processes the relationship between the hillslope soil grading, and the hillslope gradient that this paper further investigates.However, it was only tested for a small number of cases, and for one set of climate and pedogenic data.Cohen et al. (2010) showed the robustness of the relationship with changes in in-profile weathering relationship but did not investigate the full range of parameter values.This paper generalises the mARM3D formulation and extends its numerics to allow us to test the relationship for more general conditions.We present the results and insights obtained by the new modelling framework, State Space Soilscape Production and Assessment Model (SSSPAM).The state-space based model we developed using the SSSPAM framework simulates soil evolution in two horizontal dimensions (i.e.x and y), depth down the soil profile, time, and the soil particle size distribution with depth.

Modelling approaches
The combined effect of armouring and weathering on the soil evolution on hillslopes was first explored by Sharmeen and Willgoose (2006).They investigated interactions between particle weathering and surface armouring and its effect on erosion using a physically based one-dimensional hillslope soil erosion model called ARMOUR.To carry out their simulations they used surface soil grading data from two mine sites: (1) Ranger Uranium Mine (Northern Territory, Australia), and (2) Northparkes Gold Mine (New South Wales, Australia).They demonstrated that the influence of weathering was significant in the armouring process, sediment flux, and erosion rate.ARMOUR could also modify the armour properties, and even prevent the development of armour, by rapidly disintegrating the coarse material.If the amount of sediment generated by weathering is large it can be stored on the surface during times when the transport capacity is not large enough to entrain all the material.They called this a "transport-limited" regime.On the other hand for low weathering rates the armour will build up and prevent the subsurface material from eroding which was called "weatheringlimited" regime.In between these two extremes they identified an equilibrium region where the erosion and weathering balance each other and where only the fine fraction generated by weathering is removed form the surface leading to a stable armour layer.The grading of the armour layer was found to be different to the underlying soil grading (Sharmeen and Willgoose, 2006).Using ARMOUR they demonstrated the feasibility of using a physically based model to represent soil evolution to study geomorphological evolution and as a simple model for pedogenesis.The main drawback of the numerical approximation used in ARMOUR model was its high computational complexity and very long run times which prevented it from being used for more complex geometries such as 2-D catchments (Cohen et al., 2009), or its coupling with a landform evolution model.Cohen et al. (2009) simplified ARMOUR by reformulating it as a state-space matrix model, mARM, where the complex nonlinear physical processes of particle entrainment in ARMOUR were modelled using transition matrices.By doing so Cohen et al. (2009) was able to reduce the numerical complexity of ARMOUR and significantly reduce runtimes.The computational efficiency of mARM allowed Cohen et al. (2009) to explore (1) time-and space-varying relationships between erosion and physical weathering rates at the hillslope scale, (2) more complex planar drainage geometries, and (3) interactions between the soil profile and the soil surface properties.They found that for erosion-dominated slopes the surface coarsens over time, while for weatheringdominated slopes the surface fines over time.When both processes operate simultaneously a slope can be weatheringdominated upslope (where runoff and therefore erosion is low) and armouring-dominated downslope.In all cases, for a constant gradient slope the armour coarsens downslope (i.e. as drainage area increases) as a result of a balance between erosion and weathering.Thus even for weatheringdominated slopes the surface grading catena is dependent on armouring through the balance between weathering and armouring (Cohen et al., 2009).They also observed that for many slopes the surface initially armours but, after some period of time (space and rate dependent), weathering begins to dominate and the grading of the soil surface subsequently fines.Depending on the relative rates of armouring and weathering the final equilibrium grading of the slope may be finer or coarser than the initial conditions but in all cases the surface coarsened with increasing area and slope.Subsequently, mARM3D was developed by Cohen et al. (2010) to incorporate soil profile evolution by using several soil profile layers and a semi-infinite bedrock layer into the mARM framework.They used exponential and humped exponential depth-dependent weathering functions (soil production functions) to quantify the weathering characteristics of the soil profile and the bedrock.They concluded that although the soil depth and the subsurface soil profiles are dependent on the depth-dependent weathering function, their effect on the spatial organisation of the grading of the soil surface was minimal.Their simulations showed that the area-slope-d 50 relationship was still present at the soil surface even with different depth-dependent weathering functions.
These results were in good agreement with the results of the ARMOUR model used by Sharmeen and Willgoose (2006).The work of both Sharmeen and Cohen used process parameters calibrated to observed field erosion (Willgoose and Riley, 1998) and laboratory weathering data (Wells et al., 2006(Wells et al., , 2008) ) for a site at Ranger Uranium Mine.Thus their conclusions only apply to the site at Ranger.
The aim of this paper is to present a new model (SSS-PAM) that extends this previous work and generalises the conclusions using a sensitivity analysis of its process parameters.In this way we test the robustness of the Cohen's areaslope-d 50 relationship under different conditions (e.g.different climates and soil production functions).Here we present (1) the extensions in SSSPAM, (2) calibration and validation of SSSPAM, and (3) exploration of the spatial and temporal patterns of soil grading and weathering and armouring processes.The model discussed here is the soilscape component of a coupled soil-landscape evolution model and this paper aims to better understand the behaviour of this soilscape model before examining the more complex coupled soillandform system.

The SSSPAM model
SSSPAM is a state-space matrix model simulating temporal and spatial variation of the grading of the soil profile through depth over a landscape and extends the approach of the mARM model (Cohen et al., 2009) and mARM3D (Cohen et al., 2010).It uses matrix equations to represent physical processes acting upon the soil grading through the soil profile.SSSPAM uses the interaction between a number of layers to simulate soil grading evolution (Fig. 1).These layers are the following: (1) a water layer flowing over the ground which moves soil particles laterally, (2) a surface soil layer from which the water entrains soil particles and which produces an armour over the soil below, (3) several soil layers representing the soil profile, and (4) a semi-infinite nonweathering bedrock/saprolite layer underlying the soil.In SSSPAM two processes are modelled: erosion due to overland flow, and weathering within the profile.The armouring module consists of three components.
The grading of the surface (armour) layer changes over time because of three competing processes, (1) selective entrainment of finer fractions by erosion, (2) the resupply of material from the subsurface (that balances the erosion to ensure mass conservation in the armour layer) and (3) the breakdown of the particles within the armour due to physical weathering.The erosion rate of the armour layer is calculated from the flow shear stress.The entrainment of particles into surface flow at each time step from the armour layer is determined by the erosion transition matrix, which is constructed using Shield's shear stress threshold.The Shield's shear stress threshold determines the maximum particle size that can be entrained in the surface water flow.For particles smaller than the Shield's shear stress threshold a selective entrainment mechanism is used which was found to be a good fit to field data (Willgoose and Sharmeen, 2006).Resupply of particles to the armour layer from underneath is mass conservative.The rate of resupply equals the rate of erosion, so the armour's mass is constant.
The weathering module simulates the disintegration of particles in the armour and underlying soil profile layers.Weathering is also modelled with a transition matrix.It defines the change in the armour grading as a result of the fracturing of particles through the weathering mechanism.The "Body Fracture" mechanism (Fig. 2) splits the parent particle into a number of daughter particles.Wells et al. (2008) found that a body fracture model with two equal-volume daughter fragments best fitted his laboratory salt weathering experiments.This does not guarantee that this fragmentation mechanism is appropriate for other rock types not tested by Wells,  Wells, et al., 2008).and one of the cases studied in this paper is a generalisation of this equal volume fragmentation geometry.Weathering in this paper is mass conservative so that when larger particles break into smaller particles the cumulative mass of the soil grading remains constant.Thus we do not model dissolution.
The state vector g defines the soil grading at any specific time and in every layer.Entries g i in the state vector g are the proportion of the material in the grading size range i.The evolution (of the state vector) from one state to another state during a single time step is defined using a matrix equation.This matrix (called the transition matrix) describes the relationship between the states at two times and defines the change in the state during a time step where g t 1 and g t 2 are state vectors defining the soil grading at time t 1 and t 2 , R is the marginal transition matrix, I is the identity matrix, and t is the timestep (Cohen et al., 2009).
For multiple processes Eq. ( 1) can be applied sequentially for each process, using the R matrix appropriate for each of the processes.
Within each layer the equation for weathering follows Eq. ( 1) where W is the rate of weathering (which is depthdependent), and B is the non-dimensional weathering marginal transition matrix.Parameter W determines the rate of weathering while B determines the grading characteristics of the weathered particles.
For the armour layer the mass in the layer is kept constant so that as fines are preferentially removed by erosion, the mass removed is balanced by new material added from the layer below, and with the grading of the layer below.For each layer in the profile mass conservation is applied, and any net deficit in mass is (typically) made up from the layer below (i.e. by removing material in the layer below).The only exception to this rule is the case of deposition at the surface where material is pushed down.In this latter case the pushing down results from an excess of mass in the armour layer and this excess propagates down through the profile.

Constitutive relationships for erosion and armouring
The erosion rate (E) of the armour is calculated by a detachment-limited incision model, where e is the erodibility rate, q is discharge per unit width (m 3 s −1 m −1 ), S is slope, d 50 a is the median diameter of the material in the armour (m), α 1 , α 2 and β are exponents governing the erosion process.It is possible to derive exponents α 1 and α 2 from the shear stress dependent erosion physics (Willgoose et al., 1991b) or they can be calibrated to field data (e.g.Willgoose and Riley, 1998).In this paper for simplicity we will consider a one-dimensional hillslope with a unit width, constant gradient, and a 2 m maximum soil depth.The discharge was calculated by where r is the runoff excess generation (m 3 s −1 ) and x is the distance down the slope (m) from the slope apex to each node.
The implementation details of the erosion physics (e.g.how selective entrainment of fines is incorporated into the marginal transition matrix for erosion) are identical to that of Cohen et al. (2009) and will not be discussed here.The primary process of relevance here is that a size-selective entrainment of fine fractions of the soil grading by erosion is used and it follows the approach of Parker and Klingeman (1982) as calibrated by Willgoose and Sharmeen (2006).The result is that for surfaces that are being eroded the surface becomes coarser with time (and thus why we call the top layer the armour layer).

Constitutive relationships for weathering
The fracturing geometry determines the weathering transition matrix B. Each grading size class will lose some of its mass to smaller grading size classes as larger parent particles are transformed into smaller daughter particles.The daughter products can fall in one or more smaller grading classes depending on the size range of particles produced by the breakdown of the larger parent particles.The amount of material received by each smaller size class is a function of size distribution of the grading classes, fracture mechanism, and the size characteristics of the daughter particles.Wells et al. (2008) found that for his material (a mining waste product from Ranger Uranium Mine) a simple symmetric fracture model with two equal volume daughter products best fitted his experimental data.While the formulation of the weathering transition matrix in Cohen et al. (2009) allows a general fragmentation geometry, Cohen only used the symmetric fragmentation found experimentally by Wells.This paper will generalise these results and examine a broader range of fracture geometries.
To generalise the fracture geometries we will assume that a parent particle with a diameter d breaks into a single daughter particle with diameter d 1 and n−1 smaller daughters with diameter d 2 (the total number of daughters being n).For simplicity all the particles considered are assumed to be spherical.Mass conservation implies If the single larger daughter with diameter d 1 accounts for α fraction of the parent then By changing the α fraction value and the number of daughters n we are able to simulate various fracture geometries such as symmetric fragmentation, asymmetric fragmentation, and granular disintegration (Wells et al., 2008).
For instance α = 0.5, n = 2 represents symmetric fragmentation with two daughter particles, α = 0.99, n = 11 represents a fracture mechanism resembling granular disintegration where a large daughter retains 99 % of the parent particle volume and 10 smaller daughters have 1 % of the parent volume collectively.The construction of the weathering transition matrix then follows the methodology outlined in Fig. 1 in Cohen et al. (2009).

Soil profile development through depth-dependent weathering
The weathering module of SSSPAM consists of two components.They are (1) the weathering geometry for the grading of the daughter particles discussed above, and (2) the weathering rate for the different soil layers which determines the rate at which the parent material is weathered.The weathering rate of each soil layer typically (though not always) depends on the depth below the soil surface.
To characterize the weathering rate with soil depth, depthdependent weathering functions are used.In their mARM3D model Cohen et al. (2010) used two depth-dependent weathering functions (Fig. 3), (1) exponential decline (called exponential) (Humphreys and Wilkinson, 2007) and (2) humped exponential decline (called humped) (Ahnert, 1977;Minasny and McBratney, 2006).For the exponential, the weathering rate declines exponentially with depth.The rationale underpinning the exponential function is that the surface soil layer is subjected to the high rates of weathering because it is closer to the surface where wetting and drying, and temperature fluctuations are greatest.The humped function has the maximum weathering rate at a finite depth below the surface instead of being at the surface itself and then declines exponentially below that depth.The rationale for the humped function is evidence that the weathering is highest at the water table surface which leads to a humped function.
We also examined another depth-dependent weathering function we call the dynamic reversed exponential function.In this function the highest weathering rate is located at the soil-bedrock/saprolite interface and exponentially decreases upwards toward the surface and downwards into the underlying bedrock.The soil-bedrock interface is defined as that layer above which the porosity increases abruptly, reflecting the transformation from bedrock/saprolite to soil (Anderson and Anderson, 2010).Unlike the exponential and humped functions the depth of the peak weathering rate in the dynamic reversed exponential function moves up and down with the ups and downs of the soil-bedrock interface.At the soil-bedrock interface the bedrock material is transformed from bedrock to soil.The bedrock has a higher potential for chemical weathering than the soil above the soil-bedrock interface that has been subjected to chemical weathering.The function declines below the soil-bedrock interface because the reduced porosity of the bedrock inhibits water flow.Although we do not model chemical weathering in this paper, we believe that the dynamic reversed exponential function can be used to conceptualise chemical weathering.
The three depth-dependent weathering functions are graphically represented by Fig. 3.The exponential function is (Cohen et al., 2010) where w h is the weathering rate at the soil layer at a depth of h (m) below the surface and δ 1 is the depth scaling factor (here δ 1 = 1.738).
The humped function used is (Minasny and McBratney, 2006) where P 0 and P a are the maximum weathering rate and the steady state weathering rate respectively, δ 2 and δ 3 are constants used to characterise the shape of the function, and M is the maximum weathering rate at the hump which is used to normalise the function.Values we used here were P 0 = 0.25, P a = 0.02, δ 2 = 4, δ 3 = 6, and M = 0.04.The dynamic reversed exponential function is where H is the depth (m) to the soil bedrock interface from the surface which is calculated from the soil grading distribution at each iteration during the simulation, λ is a constant which determines the function value at the asymptote, δ 4 and δ 5 are constants used to characterise the rate of decline with depth of the function.We used λ = 0.98, δ 4 = 3, δ 5 = 10.The non-zero weathering below the bedrock-soil interface in Eq. ( 10) represents a slower rate of chemical weathering within the bedrock due to its lower porosity and hydraulic conductivity.In general δ 5 > δ 4 .
The weathering rate of each layer is determined by modifying the base weathering rate W 0 (Eq.2) and the depthdependent weathering function used, f (h).The weathering rate of a soil layer at a depth of h from surface W h is given by 3 Data used in this study Four soil particle size distribution data sets were used as input data for SSSPAM simulations.Two particle size distribution data sets were collected from the Ranger Uranium Mine (Northern Territory, Australia) spoil site (Willgoose and Riley, 1998;Sharmeen and Willgoose, 2007;Cohen et al., 2009;Coulthard et al., 2012).The third and fourth gradings were created from the previous two gradings to simulate the subsurface bedrock conditions.The naming convention used here is "a" for the actual grading data set and "b" for the  1).
-Ranger1a: this grading distribution was first used by Willgoose and Riley (1998) for their landform evolution modelling experiments.This soil grading was subsequently used by Sharmeen and Willgoose (2007) and Cohen et al. (2009) for their armouring and weathering simulations.This grading distribution consists of stony metamorphic rocks of medium to coarse size produced by mechanical weathering breakdown, has a median diameter of about 3.5 mm, and has a maximum diameter of 19 mm.
-Ranger2a: the second grading distribution was used by Coulthard et al. (2012) in their soil erosion modelling experiments and has a maximum diameter of 200 mm.The Coulthard data set includes a coarse fraction which is not included in Ranger1a, has a median diameter of 40 mm, and has a maximum diameter of 200 mm.Nominally Gradings 1a and 2a are for the same site but the gradings are not identical in the overlapping part of the grading below 19 mm.
-Ranger1b and Ranger2b: these grading data sets were created using the particle distribution classes of Ranger1a and Ranger2a to represent the underlying bedrock for each of the grading distributions mentioned above.To represent the bedrock for these data sets 100 % of the material was assumed to be in the largest diameter class for each grading classes (19 mm for the 1b and 200 mm for 2b).
We divided our planar hillslope into nodes with 4 m spacing downslope and the armouring and weathering was simulated at these nodes.The soil profile at each node was www.earth-surf-dynam.net/4/607/2016/Earth Surf.Dynam., 4, 607-625, 2016 represented by 21 layers representing the armour layer and 20 subsurface layers.Initially the armour layer was set to either Ranger1a or Ranger2a grading data set (depending on the type of simulation) and all the subsurface layers were set to the corresponding bedrock layer (for Ranger1a surface grading Ranger1b was set as the bedrock grading for all other subsurface layers).For brevity henceforth simulations run with the "Ranger1 data set" used the Ranger1a grading for the initial surface layer and Ranger1b as the subsurface grading unless otherwise stated.Likewise "Ranger2 data set" means, Ranger2a for the initial surface and Ranger2b for the subsurface).We have used 30 years of measured pluviograph data (Willgoose and Riley, 1998) to calculate discharge.The 30 years of runoff was repeated to create a 100-year data set as was done in our earlier work (Sharmeen and Willgoose, 2006;Cohen et al., 2009).

SSSPAM calibration
To provide a realistic nominal parameter set around which parameters could be varied in the parametric study, SSSPAM was calibrated to mARM3D, which in turn had been calibrated to ARMOUR1D (Willgoose and Sharmeen, 2006) and we know ARMOUR1D corresponded well with field data.
Figure 5 shows a comparison between contour plots generated by mARM3D and SSSPAM using identical initial conditions (Ranger1 data set) and model parameters.The figure shows that mARM3D and SSSPAM produce similar d 50 values, though SSSPAM is very slightly coarser.The slight differences between the two contour plots result from the improved numerics of SSSPAM and an improved implementation of the matrix methodology in SSSPAM.We are thus confident that SSSPAM and mARM3D are comparable.The parameter values used for SSSPAM are α 1 = 1.0, α 2 = 1.2, β = 1.0, m = 4, e = 2.5 × 10 −8 and n = 0.1.et al. (2009, 2010) found a strong log-log linear relationship between contributing area, slope and the d 50 of the armour soil grading.They quantified the relationship between soil grading, local topographic gradient and drainage area by

Cohen
where A is the contributing area to the point of interest, S is the slope at the point of interest, d 50 is the 50th percentile Earth Surf.Dynam., 4, 607-625, 2016 www.earth-surf-dynam.net/4/607/2016/Cohen used only one parent material grading and one parameter set for his analyses.To explore the generality of Eq. ( 12), we have examined the behaviour of the contour plots with changes to (1) weathering parameters, (2) grading of the parent material, (3) process and climate parameters, and (4) armouring mechanisms.We also examined a broader range of area-slope combinations that would typically occur in nature (since we are interested in man-made landforms which may have far from natural geomorphology), and which Cohen examined.For the initial conditions, unless otherwise indicated, in each simulation the "a" grading was used for the initial surface layer and the corresponding "b" bedrock grading for all the initial subsurface layers (e.g.Ranger1a for the surface and Ranger1b for the subsurface).To ensure that the hillslopes had reached equilibrium, the model simulated 100 000 years with output every 200 years.Equilibrium was assessed to occur when the grading of all nodes on the hillslope stopped changing, typically well before 100 000 years.
Figure 4 shows a time series d 50 evolution of all the nodes with lowest slope gradient (2.1 %).It shows that equilibrium is reached well before 100 000 years.Hillslopes with higher gradients reached equilibrium even faster.

Interpretation of the grading contour plots
Before discussing the parametric study and its myriad of contour plots, Fig. 5 shows how the contour plots can be used to estimate soil properties for any hillslope type.Five profiles are illustrated: Curve 1: This is a hillslope where the slope is increasing down the hillslope so is concave down in profile and looks like a rounded hilltop.The d 50 increases down the hillslope (i.e.increasing area, moving from left to right in Fig. 5).All our contour plots increase from left to right and from bottom to top, so in general concave hillslopes will always coarsen downslope.
Curve 2: This hillslope has constant slope downslope and, as for slope 1, will always coarsen downslope.
Curve 3: This hillslope has slopes that are decreasing downslope and is concave up.Importantly the gradient of the line in Fig. 5 is less than the gradient of the contours so the hillslope coarsens downslope.
Curve 4: This hillslope is similar to 3 except that the rate of decrease of slope downstream is more severe (i.e.concavity is greater) so the gradient of the line in Fig. 5 is steeper than the gradient of the contours.This hillslope fines downstream.
Curve 5: This hillslope is a classic catena profile with a rounded hilltop and a concave profile downstream of the hilltop.By tracking this hillslope downstream the surface grading (or surface d 50 ) will initially coarsen.As it transitions to concave up it will continue to coarsen until the rate of reduction of the hillslope slope is severe enough that is starts to fine downstream.Whether this latter region of fining occurs will depend on the concavity of the hillslope and whether it is strong enough relative to the gradient of the soil contours in Fig. 5.
Note that the erosion model in SSSPAM is an incision model dependent on upstream area and slope.With this model the planar shape and slopes of the catchment upstream of the point are irrelevant, so while we derived Fig. 5 for a planar hillslope it is equally valid for a natural twodimensional catchment with flow divergence and convergence.Thus it should be clear that the spatial distribution of soils, and any questions of downslope fining or coarsening of those soils, must depend on the interaction between the pedogenesis processes that produce the soils (and thus drive the area-slope dependence of soil grading) and landform evolution processes that generate those profiles (and the area-slope relations for those slopes).Ultimately deeper understanding of these links will only come from a coupled landscapesoilscape evolution model, but in this paper we confine ourselves to better understanding of the soilscape processes and the area-slope dependence of grading.The coupled model will be discussed in a subsequent paper.

Parametric Study of SSSPAM
All the nominal parameters used in the parametric study are presented in Table 2.In order to fully explore the area-sloped 50 relationship a parametric study was carried out using SSSPAM.The area-slope-diameter relationship was derived by evolving the soil on a number of one-dimensional, constant width, planar hillslopes, each with a different slope, with evolution continuing until the soil reached equilibrium.A contour plot was then created where the soil grading metric (usually the median diameter, d 50 ) was contoured for a range of slopes and area.Because of the planar slope, only erosion occurs, no deposition.Erosion is a function of local discharge, slope and soil surface grading as indicated in Eq. ( 3), and is assumed to be detachment-limited.Detachment-limitation means that the upstream sediment loads do not impact on erosion rates.Hillslope elevations are not evolved (i.e.no landform evolution occurs) which is equivalent to assuming that the soil evolves more rapidly than the hillslope so that the soils equilibrate quickly to any landform changes.

Changing surface and subsurface gradings and weathering rate
Figure 6 shows the equilibrium contour plots generated for the two grading data sets and with different weathering rates.
The equilibrium d 50 decreases with increasing weathering rate.Higher weathering rates break down the larger particles more rapidly.The equilibrium d 50 values were the same even  if the initial surface grading was changed.For example, using the Ranger1a or Ranger2b grading data for the surface but with Ranger2b for the bedrock yielded identical equilibrium d 50 results.As weathering broke down the surface layer and it was eroded it was replaced by the weathered bedrock material, which was identical when the same subsurface grading and weathering mechanism was used.Finally a coarser subsurface grading led to a coarser armour.
These trends with weathering rate are consistent with Cohen et al. (2010) where the log-log linear area-slope-d 50 relationship was observed regardless of the weathering rate.Moreover the contour lines in Fig. 6 all have the same slope.This implies that although the magnitude of the coarseness of the equilibrium armour depends on the underlying soil grading and weathering mechanism, the slope of the contours is independent of the subsurface grading and weathering process.This result demonstrates that the area-slope-d 50 relationship is robust against changes in the grading of the source material, and the only change is in the absolute grading, not the grading trend with area and slope.

Changing the runoff rate
Erosion is a function of the discharge, and the discharge depends on the climate and rainfall.The effect of changing the runoff is shown in Fig. 7. To simulate a more arid climate the runoff generation parameter in Eq. ( 4) was halved.Figure 7 shows that a reduced discharge produced a finer armour.While not shown, higher discharge rates produced coarser armour.For lower discharges (1) the Shield's Stress threshold decreases thus allowing smaller particles to be retained in the armour layer, and (2) the rate of erosion decreases while the weathering rate remains constant so that weathering (i.e.fining) becomes more dominant.Both of these processes work in tandem to produce finer armour.This conclusion is qualitatively consistent with Cohen et al. (2013), where they applied natural climate variability over several ice-age cycles and observed switching between fining and coarsening of the soil surface depending on the relative dominance of erosion and weathering at different stages in the climate cycle.

Changing the erosion discharge and slope exponents
The influence of the exponents on area and slope in the erosion equation (Eq.3), α 1 and α 2 , is shown in Fig. 8.These contour plots used the Ranger1a surface grading for the surface grading and Ranger1b bedrock grading for the initial subsurface layers.Figure 8 shows that although the d 50 values changed with different α 1 and α 2 values, the slope of the contours only changed when α 1 /α 2 was changed.To investigate the generality of this conclusion, contours were then plotted for different α 1 /α 2 .The slope of the contours was strongly correlated with α 1 /α 2 .The slope of the contours increased for higher α 1 /α 2 ratios.Similar results were obtained for the Ranger2 data set.The α 1 /α 2 ratio not only influences the slope of the contour lines but also influences the equilibrium d 50 values.For low α 1 /α 2 , the equilibrium d 50 values at the hillslope nodes were coarser than for high α 1 /α 2 .These relationships allow us to generalise the area-sloped 50 relationship www.earth-surf-dynam.net/4/607/2016/Earth Surf.Dynam., 4, 607-625, 2016 where δ, γ and ε are exponents on contributing area, slope and d 50 respectively, and c is a constant, and where the ratio δ/γ is a function of the erosion dependence on area and slope.
Figure 9 shows that δ/γ was strongly correlated with the model α 1 /α 2 even though there was no correlation with the individual parameters (i.e.α 1 with δ, or α 2 with γ ).In the regression analysis the parameter ε was assumed to be 1 in order to calculate δ and γ constants.This assumption does not affect the δ/γ ratio.This result was independent of the subsurface grading.

Changing the erodibility and selectivity exponent β and e
This section examines the effect of changing erosion equation parameters, (1) the d 50 exponent β (Eq. 3) which relates the erosion rate to median sediment diameter, and (2) the erodibility rate e.The slope of the contours was independent of these parameters.The parameters β and e influence (1) the absolute value of d 50 , and (2) the spacing of the contours.These impact on the value of c in Eq. ( 13).For higher β, the equilibrium d 50 was coarser than for low β values.Increasing the erodibility factor e yields similar results.

Different weathering fragmentation geometries
To study different weathering mechanisms we used a fragmentation geometry (Fig. 2) that has two parameters, n and α (Eqs.5-7).The simulations in the previous sections used symmetric fragmentation with n = 2 and α = 0.5 (i.e.where a parent particle breaks down to two equal volume daughter particles).Here we examine four other geometries, (1) symmetric fragmentation with multiple daughter products (n = 5, α = 0.2; i.e. the parent breaks into five equal daughters each having 20 % of the volume of the parent), (2) moderately asymmetric (n = 2, α = 0.75; the parent breaks into two daughters, with 75 and 25 % of the parent volume), (3) granular disintegration (n = 11, α = 0.9; the parent breaks into 11 daughters, one with 90 % of the parent volume and the other 10 daughters each have 1 % of the parent volume), and (4) as for Geometry 3 but with the large daughter having 99 % of the parent particle volume (n = 11, α = 0.99).Figure 10 shows results using the Ranger1 data set.The corresponding symmetric results are in Fig. 6.Symmetric fragmentation with five equal daughter particles (Geometry 1) leads to the finest equilibrium contour plot but the contours are otherwise unchanged.The granular disintegration geometries produced coarser results with the coarsest armour from Geometry 4. We conclude that when fragmentation produces a number of symmetric daughters the equilibrium grading of a hillslope is finest.Finally the slope of the contours did not change for different fragmentation geometries.

Effect of initial conditions
The simulations in the sections above used the same grading for the initial surface and the subsurface.To explore the initial conditions we changed the initial surface and subsurface data sets.The equilibrium grading contour plot generated using Ranger2a surface grading and Ranger1b subsurface gradings was identical to the equilibrium grading contour plot generated using Ranger1a surface and Ranger1b subsurface grading.Likewise the equilibrium grading contour plot generated using Ranger1a surface and Ranger2b subsurface gradings was identical to the equilibrium grading contour plot generated using Ranger2a surface and Ranger2b for subsurface gradings.The results were slightly different for different subsurface gradings.These results also show that, as expected, there was no effect of the initial conditions on the equilibrium grading.Though not shown the influence of the initial grading is only felt during the dynamic phase of the simulation before the armour reaches equilibrium.

Generalising beyond median grain size
The results above have focussed on d 50 as a measure of soil grading.However, the model can provide any particle percentile or statistic of interest.Figure 11 shows area-slope results for d 10 (i.e. 10 % by mass is smaller than this diameter).
It shows that the general trends observed in the d 50 contour plots (Fig. 6b2) are also evident in d 10 .Though not shown, similar results were found for d 90 .The slope of the contours is independent of diameter but as expected the d 10 and d 90 values are ranked d 10 < d 50 < d 90 .We conclude that the areaslope-diameter relationship we have observed in our simulations is robust across the grading profile.

Influence of the depth-dependent weathering functions
In this section we consider the three different depthdependent weathering functions (Fig. 3, Eqs.8 to 10) for the weathering rate in the subsurface soil layers.All the simulations in the previous sections used the exponential function (Eq.8).Figures 5 and 12 show that the contour plots for all weathering functions are very similar.However, as slope and area are increased the humped function produces a more rapidly coarsening armour.Overall the reversed exponential produces the coarsest armour.For the reversed exponential after an initially high weathering rate at the surface, the weathering rate reduces rapidly as the soil-bedrock interface moves deeper into the soil profile.This low near surface weathering decreases the rate of fining of the armour and dramatically reduces the erosion.This reduction in erosion rate prevents weathered fine particles from reaching the surface.We also analysed the subsurface soil profile.Figure 13 shows the d 50 through the soil profile for our onedimensional hillslope of length 32 m, divided in to eight nodes at 4 m intervals, and with 10 % slope, and Ranger1 data set.The bedrock layers are those layers near the base of the profile with the d 50 = 19 mm.The exponential and humped functions produce similar soil profiles except that the humped function produces a shallower soil and a coarser armour compared with the exponential.In contrast, the reversed exponential produces a markedly different soil profile.It produces very coarse armour, a soil thickness beyond the modelled 2000 mm limit, and a more uniform soil grading through the profile.This latter result is because the weathering is greatest at the bedrock-soil interface so most of the soil grading change is focussed at the base of the profile and relatively less occurs within the profile.
A final question is whether the area-slope-grading relationship occurs only in the armour or exists throughout the profile using the exponential weathering function.We gener-ated area-slope-d 50 contours for four different depths within the profile extending down to the base of the soil profile (Fig. 14).The slope of the contours is the same for all depths and hence we believe that the area-slope-grading log-log linear relationship is exhibited for the entire soil profile, with the only change being the coarseness of the soil (which reflects the maturity of weathering of the soils) at any particular depth.This result is intriguing because while the armouring from erosion occurs at the surface it has an impact throughout the profile, it is not simply a property of the near surface layer directly impacted by erosion.Thus the act of soil profile generation, which is solely driven by the depth-dependent weath-

Discussion
Here we have used a new pedogenesis model, SSSPAM, to analyse the equilibrium soil grading and spatial organisation of soil profiles.This model extends the mARM3D model of Cohen et al. (2010) and improves the numerics.Our results have generalised previous studies (Cohen et al., 2009(Cohen et al., , 2010(Cohen et al., , 2013) ) that have found a log-log linear relationship between d 50 , contributing area and slope.Using a broader range of environmental conditions, we have found that log-log linear relationship for grading is robust against changes in environment and underlying geology and for hillslopes where the dominant processes are surface fluvial erosion and inprofile weathering.The main factors influencing the quantitative form of the relationship are the area and slope dependency of the erosion equation, and the relative rates of the weathering and erosion processes.Coarsening of the downslope nodes was observed in all the simulations.
Our parametric study has demonstrated the versatility of our model for studying the influence of different process parameters and the dynamics of hillslope evolution.Our d 10 and d 90 contour plots show that the area-slope-diameter relationship is not only true for d 50 but is also true for other aspects of the particle size grading of the soil.This strengthens our confidence in the generality of the area-slope-diameter relationship.This relationship provides us with a methodology to predict the characteristics of soil grading on a hillslope as a function of geomorphology.It also allows us to interpolate between field measurements.Furthermore, our parametric study showed how parameters of the armouring component affect the area-slope-diameter relationship.Particularly interesting was that the ratio of the erosion exponents (α 1 /α 2 ) changes the slope of the contours.This observation also hints at the importance of topographic and process characteristics in soil evolution and hillslope catena and how these topographical units may be used for predictive soil mapping and inference of erosion process.
Previous work (e.g.Willgoose, et al., 1991b;Tucker and Whipple, 2002) has shown that topography is also a function of α 1 /α 2 and this suggests a strong underlying process link between the spatial distribution of topography and the spatial distribution of soil grading that goes beyond the concept of soil catena.The soil catena concept says that systematic changes occur in soils as a function of their position on the hillslope.Our results suggest that the same processes that influence the equilibrium distribution of topography (e.g. the erosion process that determines α 1 /α 2 ) also influence the equilibrium distribution of soils.Thus while a soil catena presumes a causal link from topography, we postulate a causal link for both topography and soils from erosion processes.
Using our model we were able to explore the soil profile characteristics and how the soil profile will change depending on the weathering characteristics of the bedrock material.Another important insight is that the area-slope-d 50 relationship is present in all the subsurface layers as well as the surface armour.
In this paper SSSPAM did not model transport-limited erosion.The implication is that the eroded sediment from nodes upslope did not impact the erosion on the downslope nodes.We also did not model an interaction between grading and the infiltration of water so no coupled behaviour with hydrology was modelled.In this paper we have only considered erosion from overland fluid flow and physical weather- ing mechanisms to predict the equilibrium soil distribution of hillslopes.There is a need to explicitly incorporate chemical and biological weathering (Green et al., 2006;Lin, 2011;Riebe et al., 2004;Roering et al., 2002;Vanwalleghem et al., 2013).Another important aspect needed is accounting for deposition of sediments so that we can model alluvial soils which requires a transport-limited erosion model.A future task is to incorporate a soils model like SSSPAM into a landform evolution model such as SIBERIA (Willgoose et al., 1991a).This would allow the modelling of the interaction between the pedogenesis process in this paper with hillslope transport processes such as creep and bioturbation.If soils evolve rapidly then it may be possible to use the equilibrium grading results from this paper as the soilscape model, on the basis that the soil evolves fast enough to always be at, or near, equilibrium with the evolving landform.If soils evolve slowly then it may be necessary to fully couple the soils and landform evolution models.This is a subtle, and not fully re-solved, question of relative response times of the soils and the landforms (Willgoose et al., 2012).

Conclusions
The most important insight from this paper is that the areaslope-grading relationship observed from an earlier generation soil profile pedogenesis model by previous authors (Cohen et al., 2009(Cohen et al., , 2010) ) is general and robust across a range of climate and geologic conditions.Despite the wide range of parameters we used in our simulations, we always observed the log-log linear area-slope-diameter relationship in our simulations although the soil coarseness depended on the parameters used.In addition, contour plots of d 10 and d 90 indicated that the area-slope-diameter relationship is valid throughout the soil grading range, not just for d 50 .It was also true for depths below the surface.The parametric study conducted on the area-slope-diameter relationship demonstrated how this relationship would change with changes in the pedogenic processes.We found that the ratio of the erosion exponents on discharge and slope, α 1 /α 2 , changes the angle of the contours in the log-log contour plots (Fig. 7).This has application in the field of digital soil mapping where easily measurable topographical properties can be used to predict the characteristics of soil properties.Importantly, the contributing area and the slope data can be easily derived from a digital elevation model, which can be produced using remote sensing and GIS techniques.Coupling SSSPAM with a GIS system can potentially improve the field of digital soil mapping by providing a physical basis to existing empirical methods and potentially streamlining existing resource inten-sive and time-consuming soil mapping techniques as, for example, in the current initiatives in global digital soil mapping (Sanchez et al., 2009).
The simple physical processes currently implemented in SSSPAM also enables it to model the evolution of hillslope soil grading.A subsequent paper will focus on the dynamics of the soil profile evolution process.Although we used only armouring and weathering as soil forming factors in this study, other processes such as chemical weathering or biological influence on soil formation can also be included in our state-space matrix modelling framework (e.g.Willgoose, 2017).With its high computational efficiency and ability to incorporate various processes in to its modelling framework, www.earth-surf-dynam.net/4/607/2016/Earth Surf.Dynam., 4, 607-625, 2016

Figure 3 .
Figure 3. Graphical representation of all the depth-dependent weathering functions used in SSSPAM.

Figure 5 .
Figure 5. Log-log Area-Slope-d 50 contour plots generated using the Ranger1 data set.(a) mARM3D (Cohen et al., 2009), (b) SSS-PAM.The dotted lines in (b) are hypothetical long profiles down a drainage line showing how the contour figure can be used to generate soil properties down a drainage line.See the text for more detail.

Figure 6 .
Figure 6.Equilibrium contour plots of d 50 values (interpolated from 48 data values, the diamonds) simulated by SSSPAM for different surface and subsurface grading data and different weathering rates (Top to Bottom: 0.1, 1.0, 10.0).(Left Column) Ranger1 data set, (Right Column) Ranger2 data set.

Figure 7 .
Figure 7. Equilibrium contour plots of d 50 generated using the Ranger1 data set with identical model parameters as used in Fig. 6a2 except changing the runoff rate, half the nominal runoff rate.

Figure 11 .
Figure 11.Equilibrium contour plots of d 10 generated using Ranger1 data set with identical model parameters as Fig. 6a2 (where the d 50 results are presented).

Figure 12 .
Figure 12.Equilibrium contour plots of d 50 generated using Ranger1 data set with identical model parameters as Fig. 6a2 except changing the depth-dependent weathering function to (a) humped, (b) dynamic reversed exponential.

Figure 13 .
Figure 13.Equilibrium soil profile d 50 generated using the Ranger1 data set with a one-dimensional hillslope with 10 % slope and 32 m length using (a) exponential, (b) humped, (c) reversed exponential weathering functions.

Table 1 .
Size distribution of soil gradings used for SSSPAM4D simulation.

Table 2 .
Parameters used in the simulations generate Fig. 6a2.