The current and future role of biota in soil-landscape evolution models

Biota are major drivers of geomorphological development. Vegetation and soil fauna act as ecosystem engineers, changing the environment through physical structures and individual activities such as litter layering, tree uprooting


Introduction
Soil, landscape and biota are closely interrelated in complex ways. For a better understanding of landscape evolution, we need to consider soil formation. Soil influences the landscape in direct and indirect ways, e.g., through soil texture differences leading to the variation of erosion processes and vice versa (Minasny et al., 2015). Biota also play an important role, as vegetation is the primary source of soil organic matter (Gabet et al., 2003;Lidicker, 2000). Plants indirectly convert solar energy into geomorphologically acting forces and modify land from the micro-to global scale. The appearance and evolution of plants are considered to be a major component of the geomorphic history of the Earth (Corenblit and Steiger, 2009). Also, in modern ecological theory, bioturbation by soil animals is considered a typical example of "ecosystem engineering" (Meysman et al., 2006). From a geomorphological perspective, biotic processes can impact rock weathering, soil formation and erosion, slope stability, and hydrology over short time scales. For example, uprooted trees (such as tree windthrows) can create glades with a particular meso-or micro-climate and lead to the start of a new succession of the ground vegetation, due to disappearance of the tree cover increasing light density and rainfall reaching the soil surface (Langohr, 1993). Over geological times, biota affect climate and climatic conditions, and dictate the mechanisms and rates of erosion that control topographic evolution (Dietrich and Perron, 2006). The importance of biota has been highlighted by the US National Research Council (NRC). Two of the NRC's nine "Grand challenges" in earth surface science are related to biota. Its fifth challenge is related to transport laws that govern the evolution of the Earth's surface and the sixth challenge is related to the question of how ecosystems and landscapes coevolve. The fifth challenge states that a mechanistic understanding of links among climate, hydrology, geology, biota, topography, rates of erosion and deposition is a fundamental goal of earth surface process research. The sixth challenge highlights the efforts in exploring life-landscape interactions and outlines a need existing for achieving a greater mechanistic understanding of these interactions (NRC, 2010). To achieve such understanding, biota-soil-landscape interactions need to be accounted for in mathematical earth surface models.
Mathematical models are often used to simulate the dynamics involved in soil and landscape development at a given spatial and temporal scale (Tucker and Hancock, 2010;van der Beek, 2013). Most mathematical earth surface models apply the same set of formulae to many locations (cells, or nodes), from which landscape-level behavior emerges. This emergence at higher scale levels is a particular strength of such models because emergent behavior cannot easily be studied in other ways. The way in which mathematical earth surface models are used varies widely. On one end of the spectrum, sandbox-type models are used as tools to visualize and freely explore the long-term implications of process and parameter choices. On the other end, models are calibrated and validated for specific landscapes, with the aim to correctly simulate a set of measured properties or process rates in those landscapes (Temme et al., 2017). Once the correct or acceptable correspondence between model-predicted and observed properties or rates is achieved, such models can be used to infer optimal parameter values (van der Meij et al., 2016) or optimal boundary conditions (Barnhart et al., 2020) or to predict future development (Willgoose, 2005).
In geomorphology, soil profile development models such as SoilGen2 (Keyvanshokouhi et al., 2016) and landscape evolution models (LEMs) have been developed and applied separately for a long time. However, they both have their limitations in presenting complete interactions between soil and landscape (Vanwalleghem et al., 2013). Soil-landscape evolution models (SLEMs) such as MILESD (Vanwalleghem et al., 2013) and Lorica (Temme and Vanwalleghem, 2016) overcome these limitations to some extent by integrating soil and landscape processes. These SLEMs succeed in simulating the development of soil properties and at the same time keeping track of surface and sediments changes, and therefore have the potential to explore the complicated soil-landscape system (Temme and Vanwalleghem, 2016;Vanwalleghem et al., 2013). Previous valuable review papers give readers insights into mathematical aspects of soil and landscape processes modeling. For instance, Willgoose (2018) gave a brief review of existing SLEMs as examples of coupled soilscape and landscape evolution models. Minasny et al. (2015) reviewed data and models on soil weathering and transport and explored the possibility of integrating pedogenesis with landscape evolution. Tucker and Hancock (2010) provided a systematic overview of common process ingredients of LEMs and SLEMs. However, the current role of biota in SLEMs is still not reviewed systematically, and knowledge gaps and future developments are not clear. Our objective is to complement existing reviews of SLEMs to help develop biota-soillandscape evolution models by highlighting neglected roles of biota in geomorphic and pedologic processes and giving suggestions for improvement. In this contribution, we will address the following research questions: • What are the roles that plants and animals play in the development of soils and landscapes? • Which of these roles are currently simulated in SLEMs?
• What is needed to fill the gaps?
We do so by reviewing the roles of biota in landscape and soil processes first and then discussing current challenges and future improvement. We realize that other processes like microbial activities and chemical processes may be also related, however physical disturbance of soil animals and vegetation is considered as a key driver of soil production and transport (Wilkinson et al., 2009). To limit the scope, we focus on vegetation and soil animal bioturbation and some chemical properties of vegetation leaf litter in non-organic terrestrial environments in this study.

Scoping review
In this review, we collected literature about SLEMs by a scoping review, which is used to obtain comprehensive coverage of available literature and identify research gaps in the existing literature (Arksey and O'Malley, 2005). Based on the methodological framework proposed by Arksey and O'Malley (2005), we included the following stages ( Fig. 1) in the scoping review: (1) Identification: this stage included determining the scoping review aims, formulating searching strategies, and identifying relevant studies. This scoping review was aimed to find current studies about SLEMs in an as comprehensive as possible scope and identify knowledge gaps from the identified studies. With these aims, we formulated research strategies including advanced searching in comprehensive electronic databases, and manual searching of relevant citations from references lists of articles, related journals and conferences. The search query applied in databases included the following terms: soil profile, regolith, soil layer, soil depth, soil evolution, formation, weathering, pedogenesis, landscape, hillslope, landform, topography, model, simulation. The search query was modified each time to meet the requirement of each database. The initial search was implemented in November 2019, and updated searches were conducted in June 2020 and December 2020. All identified citations were imported into the bibliographic manager software Zotero 5.0 (Corporation for Digital Scholarship, Virginia, USA) and duplicate items were removed manually with further duplicates removed when found later in the process. (2) Screening: with the search strategy, we picked up many irrelevant papers at the first "Identification" stage. The first-level screening, consisting of an only title and abstract relevance screening, was conducted at this stage to exclude articles that did not meet the minimum inclusion criteria. We exported citations from Zotero to the web-based systematic review software Rayyan QCRI (Qatar Computing Research Institute, Doha, Qatar) for the title and abstract relevance screening. Within Rayyan QCRI, we removed irrelevant citations and duplicate citations which were not found in Zotero. (3) Eligibility: to obtain a more precise scope of SLEMs and separate SLEMs from LEMs and soil models, a second-level screening, fulltext screening was conducted. Studies were eligible for inclusion if they discussed the development or application of SLEMs. (4) Included: papers meeting the eligibility criteria were included at this stage and charted and analyzed in later sections. Papers that were excluded at this selection stage but being helpful for this review were also referenced in other sections. No distinction was made between studies using named models such as Lorica (Temme and Vanwalleghem, 2016) on the one hand, and studies using ad-hoc unnamed models such as the numerical model developed by Anderson et al. (2019).

Effects of biota on soil-landscape processes
Biota are an important geomorphological factor. The concept of ecosystem engineering was proposed by Jones et al. (1994) and ecosystem engineers were defined as "organisms that directly or indirectly modify the availability of resources to other species, by causing physical state changes in biotic or abiotic materials. In so doing, they modify, maintain and create habitats". Ecosystem engineering caused by vegetation and animals can impact hydrology, soil production, soil mixing and sediment transport, which may greatly change soil properties and shape the land surface (de Bruyn and Conacher, 1994;Gabet et al., 2003;Gray et al., 2020;Hazelhoff et al., 1981;Istanbulluoglu and Bras, 2005;Taylor et al., 2019;Tucker and Hancock, 2010;Whitesides and Butler, 2016;Wilkinson et al., 2009; Fig. 2A). As this field grows, there is an increasing body of studies revealing the long-term effect of ecosystem engineering on landscapes (Wright and Jones, 2006). From a biogeomorphological perspective, there are complex biotic-biotic and biotic-abiotic feedbacks between ecosystem engineers and the changes they cause in the abiotic environment (Corenblit et al., 2011). These complex feedbacks often lead to assemblages of specific species and formation of species-specific landforms at finer scales ( Fig. 2B), resulting in a spatially heterogenous landscape pattern (Field, 2006;Jungerius and van Zon, 1982;Kooijman and Cammeraat, 2010;Walsh and Voigt, 1977). Therefore, understanding the geomorphological effects of ecosystem engineering and biogeomorphological feedbacks is important.

Hydrology
Biota can influence soil hydrology through rainfall interception and infiltration enhancement. For vegetation-covered land, a proportion of rainfall will be shortly stored by the canopy or the litter layer in the forest floor. The intercepted water may evaporate or pass to the ground through stemflow and root flow afterwards, depending on the vegetation characteristics. Combining canopy and forest floor interception yields a total interception flux that can amount to 37% of the rainfall, or close to 50% of the total evaporation (Tsiko et al., 2012). As the first process in the chain from rainfall to runoff, interception may importantly influence subsequent hydrological processes including infiltration and runoff generation by changing the rate of rainfall arriving at the surface (Tsiko et al., 2012;Wang et al., 2020). Soil moisture is also impacted by water uptake by plants driven by transpiration, and plant transpiration is one of the central water-cycle components of the terrestrial hydrosphere, which causes the freshwater flux of 4.5 × 10 4 to 6.2 × 10 4 km 3 yr − 1 at the global scale (Evaristo et al., 2015). Rainfall reaching the ground may infiltrate into the soil, and the infiltration can also be enhanced by the litter cover and macropores left by dead plant roots and animals. The litter cover can protect the soil surface from crust formation and water ponding, which will increase soil infiltration (Meeuwig, 1970;Xin et al., 2016). Then enhanced infiltration water fluxes move deeper into the soil, following preferential paths and conduits, such as macropores left by dead plant roots and burrowing animals (Devitt and Smith, 2002). Tree uprooting also plays an important role in soil hydrology. Uprooted trees can affect hydrology by decreasing interception of rainfall and transpiration due to the disappearance of tree cover (Langohr, 1993). Furthermore, pits and mounds left by uprooted trees can influence surface hydrological connectivity and infiltration (Phillips et al., 2017;Shouse and Phillips, 2016;van der Meij et al., 2020). Soil animals, such as earthworms also influence infiltration through activities such as burrowing, ingesting, and digesting. Ingestion and digestion products from earthworms evolve into the soil as crumbs, and crumbs and burrows together increase infiltration. Bouché and Al-Addan (1997) for instance found that the infiltration rate is correlated to earthworm biomass, burrow length, surface and volume.

Soil production
Biota are an important factor in bedrock weathering. Disturbances caused by the plant (tree) roots are considered a primary process in soil development and bedrock erosion (Pawlik et al., 2016;Phillips and Marion, 2006). Roots can enter intact, unweathered bedrock through micro-spaces created by microbial communities and fungi in the rhizosphere and mycorrhizosphere. Once roots have entered bedrock, the biochemical weathering of bedrock will be enhanced by moisture fluxes along with roots (Pawlik et al., 2016). As roots penetrate and grow in the cracks or weaker portions of the bedrock, their radial pressures can be sufficient to split up weak bedrock and subsequently increase bedrock weathering (Gabet et al., 2003). If trees are uprooted, root mass with rocks and soils will be torn out of the ground forming pit-mound topography which causes bedrock disruption and erosion (Gabet and Mudd, 2010). Vegetation litter also plays an important role in weathering. Its role in the translocation of silica and other elements from the parent rock to the upper soil layers and the formation of humic acids is well established, all of which are important in controlling weathering rates within the soil and parent rock (Walsh and Voigt, 1977). Soil animals also influence soil production when they displace weak rocks in search of food or shelter (Wilkinson et al., 2009). The importance of animals in soil production was first studied by Darwin (1840), who noted that disturbance caused by earthworms can affect subsoils up to 2.5 m depth and speculated that intestinal acid from earthworms played a role in the mineral dissolution. Darwin also calculated earthworms mounding rate and estimated that the downslope soil flux resulting from the replacement of their casts by rain can be up to 0.224 cm 2 /yr (Wilkinson et al., 2009). Most invertebrates are considered to influence soil formation processes by increasing mineral weathering and humus formation (Blouin et al., 2013;Lavelle, 1988), while vertebrates such as gophers and mountain beavers are however thought less important in soil production as most of them act no deeper than 30 cm (Gabet et al., 2003). Most vertebrates are considered significant geomorphic workers in resurfacing the land through mounding and burrowing (Winchell et al., 2016).

Soil mixing
According to Schaetzl et al. (1989), plants can mix soil by several processes: "(1) root expansion during growth; (2) decay and infilling of former root channels; (3) settling of the soil due to water extraction by roots; (4) agitation of the plant during storms, which promotes root movement; and (5) uprooting". The turbation caused by tree windthrows can even reach to 1.8-2.0 m depth, and it needs a sufficiently long time for the complete decay of the dead roots and trucks and stabilizing the soil (Langohr, 1993). Additionally, sediment accumulated on litter layers by splash processes can be transported and deposited, leading to the mixing of the plant and soil material (van Zon, 1980). Invertebrates move through the soil and push particles when seeking sources and nutrients in the soil, and during these processes, soil and organic matter may be displaced to the surface, which is known as mounding, or in any directions within soil below the surface, which is defined as mixing (Wilkinson et al., 2009). The mounding rate and mixing rate can be up to 22-28 t hm − 2 yr − 1 and 730-1100 t hm − 2 yr − 1 respectively in tropical regions (Wilkinson et al., 2009), while they might be much lower in temperate regions (Gobat et al., 2004). This significantly impacts soil development and soil texture (Phillips, 2007). For example, de Bruyn and Conacher (1994) found that in ant-disturbed plots, soil profiles at 2 cm and 26 cm have a similar texture, which was attributed to size-selective soil (clay) displacement caused by ants.

Sediment transport
Vegetation plays a vital role in sediment transport. Vegetation roots can increase slope stability which is generally expressed as additive cohesion (Istanbulluoglu and Bras, 2005;Wu, 1976Wu, , 1984. The uprooting of trees is often considered as a dramatic example of floralturbation, which significantly influences soil surface processes and soil stability (Schaetzl et al., 1989). Vegetation cover can reduce the flow velocity, enhancing soil resistance to erosion, and disperse effective shear stress, reducing the detachment potential of water (Istanbulluoglu and Bras, 2005;Prosser and Dietrich, 1995). On a given slope, areas with dense vegetation cover have lower soil denudation rates compared with areas with low cover of vegetation (Acosta et al., 2015). In vegetation sparse steeplands and drylands, vegetation can work as sediment capacitors through their physical structures such as canopies and stems (Lamb et al., 2011). In steep bedrock landscapes, vegetation "dams" can trap sediment and build mounds upslope (Lamb et al., 2013(Lamb et al., , 2011, and in semiarid environments shrub canopy protection leads to differential rain splash, which causes a flux of particles from outside the radius of the shrub towards it, and mounding of soil materials including nutrients beneath shrubs, forming a "resources island" (Furbish et al., 2009;Worman and Furbish, 2019). Besides the living part of vegetation, leaf litter on the floor also plays a substantial role in erosional processes. The litter layer can intercept rain, neutralizing its kinetic energy and splash erosion, and slow down runoff, reducing its carrying capacity (Blanco Sepúlveda and Aguilar Carrillo, 2015;Elwell and Stocking, 1976). Field observations confirm that a thick intact litter layer is an important soil protection agent, and it contributes more significantly to reducing runoff erosion than undergrowth vegetation in some circumstances (Brandt, 1988;Liu et al., 2017;Park et al., 2010). Besides covering the soil surface, litter can be incorporated into top-layer soil through burrowing and mixing caused by wind, water and soil animals. This litter-soil layer can reduce the effective erosive energy and the soil detachment capacity of overland flow, decreasing soil erosion (Wang et al., 2019). Burning of vegetational cover caused by forest fires can significantly increase yields of runoff and sediment until its recovery (Inbar et al., 1998;Wittenberg and Inbar, 2009). Animals are generally considered to be major agents of soil transport. A substantial number of mounds, burrows, and casts generated by animals are not protected by vegetation cover and are easily eroded by rain splash or overland flow. Also, soil creep can be influenced by biota. In some circumstances, soil creep is dominated by bioturbation such as animal burrowing combined with slope wash (Heimsath et al., 2002).

Biogeomorphological feedbacks
As we reviewed in the above sections, vegetation litter plays a substantial role in hydrology and geomorphology and litter quantity and quality are the most important variables of litter functioning. The reduction of soil erosion from the litter layer is positively related to coverage of the litter layer, and when the coverage is beyond the threshold level (60-65%), the reduction in erosion is significant (Blanco Sepúlveda and Aguilar Carrillo, 2015). The cover and thickness of a litter layer are influenced by its decomposition rate, which is well related to litter quality. Litter rich in nutrient content or soluble substances with less acid-soluble and acid-insoluble substances, and lower secondary compounds can normally decompose faster (Martínez-Yrízar et al., 2007;McClaugherty et al., 1985). Easily decomposed litter will lead to a thinner Leaves-Fragmented-horizon, which results in less resistance to splash and overland erosion. Furthermore, leaf litter is a major source of soil organic matter and animal food, and litter quality plays a major role in the litter decomposition food web (Staaf, 1987). Litter quality together with soil properties (e.g., pH, bulk density) can determine the richness, abundance, and species composition of soil invertebrates. Some invertebrates, such as ants and earthworms are regarded as the most important bioturbators and they can significantly influence geomorphological processes through bioturbation (Bjørnlund and Christensen, 2005;Reich et al., 2005;Staaf, 1987;Szlavecz et al., 2018). Soil fauna are more sensitive to litter quality than soil attributes, such as pH and bulk density (Ilieva-Makulec et al., 2006;Moço et al., 2010). Moço et al. (2010) revealed that the richness of fauna was directly affected by lignin and polyphenol content in the litter positively and negatively respectively and the chemical components related to acidity, nutrition, and palatability are most decisive for the abundance and diversity of soil and litter fauna. Litter quality can serve as a useful environmental filtering tool, through which soil animals feeding on fresh soil organic matter and highly processed soil organic matter can be separated into forests with litter differing in decomposing rates (Szlavecz et al., 2018;Taylor et al., 2019). Soil animals, then, change the structure of litter layers and underlying soils in seeking food, mates and habitats (Szlavecz et al., 2018). During these processes, soil animals not only translocate, reorganize and loosen soil but also can create soil through biochemical and biomechanical processes (Darwin, 1840), which are considered as major drivers of soil and landscape development (Black et al., 1990;Ellison, 1946;Gabet, 2000;Gray et al., 2020). Under specific vegetation species with different litter quality, the level of faunal turbation is different. These specific vegetation-animals interactions can lead to spatially heterogeneous abiotic-biotic feedbacks within ecosystems for a long time, often resulting in characteristic biogeomorphological patterns (Baas and Nield, 2007;Corenblit et al., 2011;Fig. 2B).

A case study
This well-documented interaction between vegetation-soil, biota-soil formation-hydrology, and geomorphological development has been observed and documented in a Keuper forest of the Gutland region in Luxembourg (Cammeraat, 2002;Duijsings, 1987;Hazelhoff et al., 1981;Hendriks, 1993;Imeson and Jungerius, 1977;Kooijman and Cammeraat, 2010;van Hooff, 1983), where European hornbeam (Carpinus betulus L.) and European beech (Fagus sylvatica L.), deciduous broadleaf trees with different litter quality, are dominant tree species. Hornbeam has high litter quality and leads to high earthworm activity while beech has low litter quality and leads to low earthworm activity (Kooijman et al., 2019). After long-term interactions, there are obvious spatial differences in microtopography and soil profiles under beech and hornbeam patches (Fig. 3). Hornbeam plots are wet, with a shallow impermeable clayey soil layer (Bg horizon) and undergrowth vegetation, and the litter layer is present only in autumn and winter (Hazelhoff et al., 1981), while beech plots are dry, with an impermeable clayey soil layer at slightly greater depth, less vegetation and whole-year dense litter layers (Kooijman and Cammeraat, 2010). Litter quality plays an important role in making differences between hornbeam plots and beech plots (Kooijman et al., 2019). The topsoil under hornbeam trees is bare or covered by a thin litter layer for only a few months, which is sensitive to splash erosion. Furthermore, soil animals such as earthworms and moles create channels and pipes in the topsoil layer, strongly increasing pipe flow and soil erosion (Hazelhoff et al., 1981;van Hooff, 1983). Due to the reduced impact by the litter layer, higher soil moisture, and better soil conditions, hornbeam plots have more undergrowth vegetation. In beech plots, the soil is relatively dry because of the high water retention in the dense litter layers and the absence of clay brought by soil animals from deeper soils. The impermeable subsoils are deeper because animal bioturbation is very low and thick litter layers decrease splash and overland erosion from the topsoil. Meanwhile, leaching and throughflow are enhanced, resulting in a greater loss of clay particles from the top of the subsoil (Bg horizon) (Cammeraat, 2002). Undergrowth vegetation is very low as their growth is hampered by dense organic layers and relatively dry conditions. In this forest, differences in litter quality together with soil fauna, through interacting with hydrology, geomorphic and soil processes, lead to a spatially heterogeneous pattern of trees patches.

Overview of identified SLEMs
From the 35 papers included in the last stage of the scoping review, we identified 19 unique SLEMs. Though it is realized that some models were not classified as SLEMs by their authors at first, we present them as SLEMs here because they combine soil production or formation and landscape evolution processes. The extent to which SLEMs have presented the range of soil-landscape-biota processes can be expressed as the process coverage (Minasny et al., 2015). We selected 4 soil-related processes (weathering and soil production, clay translocation, underground water flow, soil organic matter dynamics) and 3 landscape processes (hillslope sediment transport, overland water erosion and deposition, particle size selectivity) which are considered key in pedogenesis and landscape evolution (Minasny et al., 2015;Tucker and Hancock, 2010) (Fig. 4). The soil layering system is also included in Fig. 4 as it represents soil profiles and may influence the universal application of the model (Temme and Vanwalleghem, 2016). Some soilrelated processes such as soil organic matter (carbon) dynamics are rarely included in pure-geomorphic LEMs, whereas most soil models do not consider geomorphic processes (Minasny et al., 2015). For biotic processes, we selected additional 6 processes reviewed in Section 3 (sediment transport influenced by vegetation, canopy interception and permeability reinforcement, biota influencing weathering, soil mixing  caused by animals, bioturbation-dominated colluvial transport, the geomorphic effect of litter and biotic interactions). Each model's coverage of these 14 processes was then assessed. The description of the presence of the model process is based on the most relevant papers we identified and functions of individual models that are unreferred by their papers or beyond our scope are not included here. Software packages such as Landlab presenting a framework to help create individual models and enable interaction studies in Earth surface dynamics (Barnhart et al., 2020;Hobley et al., 2017) are not included in this study.
In the 19 selected models, the most commonly simulated processes are weathering, overland water erosion and deposition, and hillslope sediment transport. These are fundamental processes in soil production and soil erosion, which together determine the thickness of the soil. Only a few models additionally include clay translocation, organic matter dynamics, and bioturbation (except the effect of vegetation on sediment transport), not only because of the different focuses of individual models but also due to the difficulty of measuring and hence confidently simulating these underground processes. The models also differ in simulating fixed or dynamic soil layers numbers and thickness. It seems that SSSPAM (Welivitiya et al., 2019), mARM3D (Cohen et al., 2010), mARM5D (Cohen et al., 2015), Lorica (Temme and Vanwalleghem, 2016) and HydroLorica (van der Meij et al., 2020) enable users to change the number of soil layers. SSSPAM, mARM3D and mARM5D simulate soil profiles with the constant user-defined layer numbers and thickness in time, while Lorica, HydroLorica include soil layers that vary in number, and vary in thickness down profile and with time, and other models use fixed soil layers numbers and changing soil thickness over the simulation time (Fig. 4). In general, it seems that having a variable number of layers and dynamic layer thickness, allows the best balance between recording pedological detail and minimizing computer memory requirements (by using fewer or thicker layers where soils are less variable) -but this balance comes with a computational cost as dynamic layer systems need to be evaluated and updated regularly during model runs.
It seems that some SLEMs use empirical or other functional relationships to describe the whole/partial biotic and soil-landscape processes (Hoosbeek and Bryant, 1992). For example, SSSPAM (Welivitiya et al., 2019) involved empirical relationships in characterizing erosion and deposition. Minasny and McBratney (2001) incorporated an empirical relationship between the rate of physical weathering and soil thickness. In the model developed by Pelletier et al. (2013), they linked the potential bedrock weathering with the effective energy and mass transfer using best-fit coefficients obtained from published data. In an ideal situation, a mechanistic understanding would describe and predict the behavior of a system in changing environment (Hoosbeek and Bryant, 1992) and more research is still needed to achieve more mechanistic understandings of the underlying soil-landscape process.
The models also differ in their coverage of processes. The models MILESD (Vanwalleghem et al., 2013) and Lorica (Temme and Vanwalleghem, 2016) have the higher process coverage, and could therefore be considered better from the perspective of integrating soil processes and landscape processes. The cost of this better ability to incorporate interactions is a larger number of model parameters that need to be measured or estimated. MILESD is one of the first models to highlight the link between soil production and landscape evolution. This model was devised to not only simulate soil horizon depth and soil properties but also to predict sediment export and its particle size distribution ( Vanwalleghem et al., 2013). However, it may have potential limitations, as only three layers in the model are simulated (A, B and C horizons) and the landscape evolution module is simplified (Minasny et al., 2015). Temme and Vanwalleghem (2016) developed SLEM Lorica, based on MILESD, and aspects of landscape evolution model LAPSUS (Schoorl et al., 2002). Lorica complements MILESD by introducing a dynamic soil layering system and new landscape processes. Lorica uses the stream power incision model, considering the spatial difference of topsoil layers and taking an approach in which sediment transport is limited by both detachment and transport. Erosion or deposition occurs depending on the comparison between streamflow transport capacity and the actual amount of sediment in transport at every transition between model cells. This approach (Davy and Lague, 2009) comes with a computational cost. Recently, the Lorica version HydroLorica was developed to include water flow and water availability over timescales of days rather than years as drivers of natural soil, landscape, and vegetation change (van der Meij et al., 2018). van der Meij et al. (2020) simulated diffusive fluxes caused by soil creep. In their simulation, the potential creep rate decreased exponentially with the soil depth, and the transport of sediment from a layer to lower neighbor layers was proportional to the surface slope and shared layer boundaries (van der Meij et al., 2020).

Biotic processes in SLEMs
As discussed earlier, plants and animals are important factors in soil and landscape evolution (Gabet et al., 2003). However, for a large fraction of LEMs (Tucker and Hancock, 2010) and SLEMs, abiotic soil and landscape processes are central issues, while biota are included as merely a minor parameter or processas indicated in Fig. 4. We acknowledge that the simplification of biotic processes is sometimes due to the intended scale of models. Specifically, some SLEMs and LEMs are intended to operate at orogen scale and therefore do not include detailed biotic processes. However, most SLEMs operate at watershed or hillslope scale, where biota impact geomorphology at patch and landscape scale (Dietrich and Perron, 2006), therefore biota deserve more inclusion. The mathematical descriptions of 6 biotic processes in individual SLEMs (Fig. 4) are presented in Table 1.

Vegetation influencing sediment transport
Root cohesion can influence hillslope stability and sediment transport. The Dietrich-et-al model (1995) output maps of relative slope stability with changing root cohesion and revealed how different vegetation altered the relative instability over the landscape. Brooks et al. (1995) simulated the effect of root cohesion of pine and grass on slope stability and failure, as shown in eq. 1 (Table 1). Other models consider vegetation cover as a variable of the sediment transport rate function, as a greater fraction of the flow shear stress is exerted on vegetation than on soils, reducing the detachment potential (Istanbulluoglu and Bras, 2005). MILESD combined both cover and roots of vegetation by introducing the biological activity index (Vanwalleghem et al., 2013) seeing eq. 2 (Table 1). Lorica assumed that erosion caused by overland flow has a negative exponential relationship with vegetation cover (Temme and Vanwalleghem, 2016) as shown in eq. 3 (Table 1). While these models succeed in incorporating vegetation in a simplified way, they may have potential limitations. The models do not consider other bioturbation disturbances such as tree throw, growth, and death of roots (Gabet and Mudd, 2010;Kirwan and Shugart, 2008), and vegetation is only assumed to reduce transport rates by protecting the surface from water flow. Hydrolorica recently included tree throw, which was assumed to form pit and mound topography and influence vertical and lateral soil transport and disrupt overland water flow by creating depressions (van der Meij et al., 2018). CAESAR-Lisflood, a coupled hydrodynamic LEM, linked vegetation types and dynamics with rates of erosion in the floodplain by considering the effect of different maturity time and lifecycles of trees and grass on the protection from erosion (Feeney et al., 2020). However, the effect of the litter layer in protecting the soil from erosion is still not considered. In some forest areas, the litter layer plays a more important role in reducing runoff erosion than undergrowth vegetation (Liu et al., 2017). Implementing such knowledge in models would also increase our understanding of such processes in the interplay between vegetation and geomorphological processes. Specifically, nonlocal process descriptions have the potential to incorporate the roles of biota in hillslope sediment transport in a more mechanistic manner These nonlocal sediment transport formulations include the probability distribution of particle travel distances, which incorporates all possibilities of individual particle motions (Doane et al., 2018). In principle, the parameters of such nonlocal sediment transport formulations could be tuned to different roles for different biological elements such as litter cover, which would help improve the detail involved with SLEMs. Brooks et al. (1995) tested the role of vegetation among other factors in influencing slope stability. They considered two scenarios, pine trees and grass cover, and assumed that pine trees cover provided larger canopy interception, soil cohesion, and permeability enhancement than grass cover (eqs. 4 and 5 in Table 1). Their results showed that vegetation's effect on slope stability is clear. In their study, water in the soil column was routed vertically cell by cell. However, soil water movement under vegetation is more complicated due to root effects. Studies show that when roots are dead and decomposed, they can leave behind macropores and water will move along these macropores rather than in a uniform pattern (Gabet et al., 2003). These preferential flow paths of water are heterogeneous and interconnected in a complex way (Germann and Beven, 2006) and are still difficult to quantify (Allaire et al., 2009).

Biota influencing weathering
The study conducted by Phillips et al. (2008) shows that rapid initial bedrock weathering rate is largely attributed to aggressive plants colonization, and favorable climate and soil moisture. However, due to the complexity of biological weathering, it is still not fully appreciated and modeled, which limits its inclusion in modeling geomorphic and soil development (Pawlik et al., 2016). Pelletier et al. (2013) tried to include  ( Brooks et al., 1995) (Pelletier et al., 2013) 4. Soil mixing caused by animals   (Pelletier et al., 2013) 6. Geomorphological effect of litter and biotic interactions Lacking partial biological effects on physical weathering in the model, by assuming that the potential weathering rate of bedrock increased exponentially with the effective energy and mass transfer (EEMT), the sum of energy from effective precipitation, and net primary production (eqs. 6,7,8 in Table 1). EEMT was considered as an optimal, climaticallybased variable for quantitative prediction of weathering rates, but the relation between EEMT and p 0 is of a functional type and more studies are still needed to achieve a mechanistic process description.

Soil mixing caused by animals
Animal bioturbation contributes greatly to the mixing of differentsize soil particles and soil carbon within and between soil layers. MILESD and Lorica considered that bioturbation causes soil mixing within each layer and leads to soil and carbon fluxes between layers (Temme and Vanwalleghem, 2016;Vanwalleghem et al., 2013). Fluxes were calculated as being proportional to biological activity and inversely proportional to the distance between layers (eqs. 9, 10 in Table 1). In MILESD, biological activity is expressed as a biological activity index varying according to depth in the profile, soil thickness, and soil carbon content. Due to the lack of experimental data, the selectivity of bioturbation transport is based on the assumption that coarse fragments cannot be moved by animals while other size classes are moved in the same proportion (Temme and Vanwalleghem, 2016;Vanwalleghem et al., 2013).The ability of mixing varies among different animal species, from worms to ants (Wilkinson et al., 2009;Winchell et al., 2016), it would be better if different groups of animals are considered.

Bioturbation-dominated colluvial transport
In the model frame proposed by Yoo and Mudd (2008), animal bioturbation is considered as the major cause of colluvial transport in a physically disturbed zone (PDZ) of the soil (eq. 11 in Table 1). Colluvial fluxes are calculated as a function of soil thickness and slope gradient on sloping grounds, while the gross colluvial fluxes are negligible on level grounds. Soil organic matter contents are usually higher in vegetated ecosystems than in bare soil, which may favor soil animal activity, and increase colluvial transport (Gabet, 2000). Based on this theory, Pelletier et al. (2013) assumed that creep includes abiotic and bioturbationdriven transport, and more vegetation cover can lead to more bioturbation. Using aboveground biomass as a measure of vegetation cover, they calculated bioturbation-driven transport by expecting that bioturbation increases linearly with aboveground biomass (eq. 12 in Table 1). Their results show that soil-vegetation-topography variables are closely correlated. This model provides a preliminary quantification of interaction between vegetation and bioturbation-dominated creep. In some cases, such as in tropical rainforests, high above-ground biomass may even be associated with low below-ground biological activity. Besides that, differences in litter quality can differentiate animal activity, while most SLEMs expected a uniform effect of vegetation on bioturbation. It deserves attention when working on the interactions between animals and vegetation.

Lack of leaf litter-related processes and vegetation-animal interactions
Comparing the role of biota in soil-landscape systems and biotic processes simulated, we found that the role of the litter layer and the effect of interactions between litter and animals are still lacking in current SLEMs (Fig. 4 and Table 1). Given that the litter layer is vital in protecting soil from water erosion, it deserves further efforts to incorporate its effects in SLEMs. Besides that, as we presented in previous sections, soil fauna and vegetation can act as important ecosystem engineers, driving both soil and geomorphological evolution. Simulating these complicated and specific biotic and biotic-abiotic feedbacks over different timescales are fundamental for linking geomorphology, ecology, and food web ecology, and understanding and predicting biogeomorphic systems development. However, vegetation and soil animals are not well connected in current models. Interactions between vegetation and animals are complicated and heterogenous, while the biotic effect is always considered to be uniform and spatially homogenous in current models. To better understand how biota and environment coevolve, more efforts are needed to complete existing models and present more biogeomorphological feedbacks between biota and the soil-landscape system.

Challenges
Through the scoping review, we identified 19 unique SLEMs, published from 1995 to 2020. Though these models succeeded in modeling landscape and soil processes, they gave less attention to a key factor in landscape and soils evolutionbiota (Fig. 2). It is not only because of the different focuses of these models but due to challenges existing in including biotic interactions. Developing a biota-soil-landscape model is complicated by the complexity and heterogeneity of biological activity, limitations of measurement methods, and lack of data.
Plants and animals are strongly influenced by environmental conditions. The influence of plants and animals on soil and landscape processes and interactions between each other vary in space and time (Cohen et al., 2010). This variation increases the difficulty of integrating biota with long-term geomorphological processes in modeling. Some vegetation-related parameters such as root cohesion may vary widely in time because plants are sensitive to land-use change, climate change, wildfires, floods and biological diseases . Protection of the soil by leaf litter also changes with litter cover. The thickness and presence of the litter layer vary with litter quality and seasons changing (Kooijman et al., 2019). For animals, processes are even more complicated, as their impacts are strongly species-dependent and difficult to predict. For instance, studies on the effect of pocket gophers show that their soil deposits and their burrows differ significantly from the background soil matrix. For instance, Reichman et al. (2001) found that in some cases, soil movement caused by pocket gophers was too complex to be predicted by physically-based models. Besides that, in many situations, a large variation in bioturbation exists when there are more than one type of bioturbators. Different soil animals such as ants, termites, and earthworms behave quite differently even at the species level, and they often vary in rates of mounding, mixing, burial, and distribution of deposited soil and nest presence (Cammeraat and Risch, 2008;de Bruyn and Conacher, 1994;Quintana et al., 2007;Taylor et al., 2019;Wilkinson et al., 2009). These species-specific activities can have unique effects on sediment transport processes and soil properties leading to spatial heterogeneity. Vegetation type and litter quality are important in determining the abundance and composition of bioturbators (Taylor et al., 2019). On the other hand, soil animal activity can also influence vegetation growth, for example, earthworm digesting fungi and bacteria can have a beneficial effect on root functioning (Ruiz et al., 2015). Given that biological activities and interactions are species-related, complicated, and heterogeneous, it is not surprising that these biological processes are difficult to generalize and include in the models.
Translating biotic processes into quantitative rates of soil movement is never easy. Timescales over which different animal activities rates are measured vary according to different measuring methods. In previous studies, to determine reliable estimates of soil turnover rates, long-term monitoring was necessary (de Bruyn and Conacher, 1994;Jungerius et al., 1989), which is difficult and cost-intensive. In addition, belowground biological activities such as the growth of roots, as well as soil mixing cannot be observed directly. Recently, luminescence (a lightsensitive property of mineral matter) measurements combined with existing theories have been applied as a long-timescale sediment tracer in soils to reveal the structure of soil mixing (Gray et al., 2020). Evidence showed that single-grain optically stimulated luminescence (OSL) dating can offer an efficient solution for quantifying pedogenic processes such as soil mixing, however, it is not widely applicable due to time cost and intensive labor needed (Román-Sánchez et al., 2019;Wilkinson et al., 2009). Román-Sánchez et al. (2019) provide evidence that a novel method, single-grain post-infrared infrared stimulated luminescence (post-IR IRSL) measurements on sand-sized grains of feldspar from the soil matrix, can provide direct information on bioturbation processes. And their results confirmed an assumption that bioturbation declines exponentially with soil depth below the surface. In their study, bioturbation rates are assumed to be homogeneous along the catena for simplification. However, it will be less realistic when there are more bioturbator species in the system, and bioturbation at the landscape scale will be more heterogeneous. Furthermore, biotic mixing processes (such as soil mixing by vegetation roots, animals) and abiotic mixing processes are not easily distinguished and further collection of soil mixing data is still needed (Gray et al., 2020).
Last but not least, limited available data and models of bioturbation also lead to a simplification of biotic processes in SLEMs. Darwin first realized that soil reworking by invertebrates can have a significant influence on soil formation and landscape development, however, his concept was picked up by geologists and geographers only in recent decades (Meysman et al., 2006). So, compared with research on physical and chemical processes in the Earth sciences, research on bioturbation is still relatively unexplored. Though soil turnover rates by animals and plants are better known than before, insights into controlling processes are still lacking because of the complex processes of bioturbation and their interactions. Modeling biological processes is challenging, more complicated biota-soil-landscape models are still lacking so far (Temme and Vanwalleghem, 2016;Vanwalleghem et al., 2013). There are some soil bioturbation models such as the SoilGen1-model, incorporating bioturbation with pedogenetic processes (Finke and Hutson, 2008), and OC-VGEN, testing the effect of bioturbation on soil organic carbon depth distribution (Keyvanshokouhi et al., 2019), and SLEMs may benefit from progress made by these soil models in the future. Although there is a growing body of literature about the leaf litter layer and interactions between litter and soil animals (e.g., Bjørnlund and Christensen, 2005;Hazelhoff et al., 1981;Reich et al., 2005;Staaf, 1987;Szlavecz et al., 2018), their mechanistic descriptions are still lacking.

Suggested improvements for future SLEMs
It is thus attractive to consider which future improvements in SLEMs would be needed to overcome the lack of biotic representation and investigate more involved biota-soil-landscape interactions. Approaching this question from the point of view of model architecture, we identified possible improvements in SLEM's vertical discretization (increasing vertical resolution) and state variables (increasing the number of soil properties that are simulated). We then defined five levels of SLEM complexity that would be made possible by such improvements that represent increased detail in biota-soil-landscape interactions (Fig. 5). It is recognized that model improvement is also possible within levels, by improving process formulations to bring them closer to biophysical reality. An example of such improvement would be the replacement of local diffusion formulations by non-local transport formulations that allow a more mechanistic exploration of the impacts of vegetation on hillslope sediment transport (see above).
The simplest structure (level 1) for SLEMs has one soil layer of a known thickness (level 1) and allows productive simulation of vegetation growth as a function of soil thickness and potentially landscape location (Collins et al., 2004). Improving on that, some current SLEMs discretize the soil as a set of multiple layers and simulate the change in mass across several grain size classes due to pedogenic and geomorphic processes. Based on a layer's texture, the first estimate of bulk density and soil hydraulic parameters can also be derived. Then, simplified soil hydrological functions such as the bucket-model approach (Guswa et al., 2002) can be used be simulate soil wetness in one temporal and three spatial dimensions (van der Meij et al., 2018). Compared to level 1, at level 2, plant-available soil water (PASW) can be simulated in substantially more detail, which will in most cases mean that simple rules relating PASW to landscape position and soil thickness are no longer neededthose relations would now emerge from simulations (Table 2). Equally thick, but differently textured or developed soils in similar landscape positions can thus support different vegetation intensity, different vegetation types and biotic activity, which in turn can impact pedogenic and geomorphic processes differentlyamong others through different rates or depths of bioturbation. We are not aware of currently available SLEMs that incorporate all aspects of level 2.
We define level 3 SLEMs as SLEMs that in addition to level 2 SLEMs; 1) incorporate a separate litter layer, 2) keep track of the amount of litter in that layer and the amount of soil organic matter in the mineral soil layers, and 3) allow vegetation, biotic activity and the various types of organic matter to influence soil structure and hence soil hydraulic properties. Although some existing SLEMs incorporate the organic matter aspects of level 3 (Rosenbloom et al., 2001;Temme and Vanwalleghem, 2016;Vanwalleghem et al., 2013), detailed representation of litter and organic matter so far is restricted to 1-dimensional soil models (Harden et al.,1999;Liu et al., 2003;van Oost et al., 2005;Yadav and Malanson, 2009;Yoo et al., 2006) and soil structure was recently Fig. 5. Overview of possible improvements in SLEM vertical discretization and state variables, with different levels of improvement identified with different levels of shading and numbers in pentagons. The inclusion of a property as a state variable means that that property will be changed over time and space during model runs, using equations or rules set in the model. included in a 1-dimensional simulation model for the first time (Meurer et al., 2020). The Meurer-et-al model (2020) productively relates organic matter storage to soil porosity, soil pore size distribution, bulk density and hence to layer thickness. It thus provides a great basis for the inclusion of structure also in SLEMs. If our understanding of the combined processes that create, maintain and destroy litter layers and soil structure would increase to a level that would allow such 3D model simulation, then vegetation growth and biotic activity could be much more closely coupled to soil development. Among others, this would allow process-based exploration of soil-biota co-evolution (Porder, 2014). The possible development and use of level 4 and 5 SLEMs appear currently limited not only by pedon-scale process knowledge but also by computing limitations. Level 4 SLEMs are defined as those that in addition to level 3 SLEMs have soil layers thin enough to incorporate Reynolds equation soil hydrology. At the time of writing, it seems possible to mathematically define such models (with SoilGen as inspiration, Keyvanshokouhi et al., 2016) but simulations over relevant timescales and reasonable resolutions are restricted to super-computing environments. Level 5 SLEMs, finally, would include state variables that express soil chemistry, including acidity and nutrient status. Although basic formulations for such SLEMs have been proposed (Kirkby, 1977, but without the detailed level 4 hydrology), no implementation is known to us. Vegetation formulations that can take soil nutrient status into account, do exist, for instance as crop growth models (Williams et al., 1989). As recognized in Fig. 5, not all aspects of a SLEM level need to be satisfied at the same time. Instead, we see the levels as aspirational goals that are worthy of partial adoption.

Conclusions
Bioturbation is considered important in soil development and landscape evolution. A mechanistic biota-soil-landscape model will increase our knowledge of the coevolution of ecosystem and landscape and help us predict and manage the environmental future.
• Among the identified 19 models, Lorica and MILESD have the highest coverage of key soil-landscape processes. Weathering, overland water transport, and hillslope sediment transport are commonly modeled by SLEMs, while clay translocation, organic matter dynamics and bioturbation are less included. • In current SLEMs, vegetation is assumed to influence overland water erosion and deposition, partial bedrock weathering, soil moisture regime and permeability. Soil animals are considered in soil mixing and colluvial transport. However, compared with geomorphic processes, biotic factors are still underappreciated and biotic processes are oversimplified. • None of the SLEMs incorporate the geomorphological role of litter and biogeomorphological feedbacks with modeling so far. Given the importance of litter quality and quantity, and interactions between vegetation and animals, more efforts are needed. • From the point of model architecture and existing challenges, we identified possible improvements in SLEM vertical discretization and state variables and defined five levels of SLEM complexity that allow an enhanced biota-soil-landscape model.

Funding
This work was supported by the China Scholarship Council (grant No. 201806390020).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.