Articles | Volume 20, issue 15
Research article
14 Aug 2023
Research article |  | 14 Aug 2023

Mammalian bioturbation amplifies rates of both hillslope sediment erosion and accumulation along the Chilean climate gradient

Paulina Grigusova, Annegret Larsen, Roland Brandl, Camilo del Río, Nina Farwig, Diana Kraus, Leandro Paulino, Patricio Pliscoff, and Jörg Bendix

Animal burrowing activity affects soil texture, bulk density, soil water content, and redistribution of nutrients. All of these parameters in turn influence sediment redistribution, which shapes the earth's surface. Hence it is important to include bioturbation into hillslope sediment transport models. However, the inclusion of burrowing animals into hillslope-wide models has thus far been limited and has largely omitted vertebrate bioturbators, which can be major agents of bioturbation, especially in drier areas.

Here, we included vertebrate bioturbator burrows into a semi-empirical Morgan–Morgan–Finney soil erosion model to allow a general approach to the assessment of the impacts of bioturbation on sediment redistribution within four sites along the Chilean climate gradient. For this, we predicted the distribution of burrows by applying machine learning techniques in combination with remotely sensed data in the hillslope catchment. Then, we adjusted the spatial model parameters at predicted burrow locations based on field and laboratory measurements. We validated the model using field sediment fences. We estimated the impact of bioturbator burrows on surface processes. Lastly, we analyzed how the impact of bioturbation on sediment redistribution depends on the burrow structure, climate, topography, and adjacent vegetation.

Including bioturbation greatly increased model performance and demonstrates the overall importance of vertebrate bioturbators in enhancing both sediment erosion and accumulation along hillslopes, though this impact is clearly staggered according to climatic conditions. Burrowing vertebrates increased sediment accumulation by 137.8 % ± 16.4 % in the arid zone (3.53 kg ha−1 yr−1 vs. 48.79 kg ha−1 yr−1), sediment erosion by 6.5 % ± 0.7 % in the semi-arid zone (129.16 kg ha−1 yr−1 vs. 122.05 kg ha−1 yr−1), and sediment erosion by 15.6 % ± 0.3 % in the Mediterranean zone (4602.69 kg ha−1 yr−1 vs. 3980.96 kg ha−1 yr−1). Bioturbating animals seem to play only a negligible role in the humid zone. Within all climate zones, bioturbation did not uniformly increase erosion or accumulation within the whole hillslope catchment. This depended on adjusting environmental parameters. Bioturbation increased erosion with increasing slope, sink connectivity, and topography ruggedness and decreasing vegetation cover and soil wetness. Bioturbation increased sediment accumulation with increasing surface roughness, soil wetness, and vegetation cover.

1 Introduction

Bioturbation was shown to shape the land surface (Hazelhoff et al., 1981; Istanbulluoglu, 2005; Taylor et al., 2019; Tucker and Hancock, 2010; Whitesides and Butler, 2016; Wilkinson et al., 2009; Corenblit et al., 2021) by influencing surface microtopography (Reichman and Seabloom, 2002; Kinlaw and Grasmueck, 2012; Debruyn and Conacher, 1994) and soil properties such as soil porosity, permeability, and infiltration (Reichman and Seabloom, 2002; Yair, 1995; Hancock and Lowry, 2021; Ridd, 1996; Hall et al., 1999; Coombes, 2016; Larsen et al., 2021). Cumulatively, these modifications lead to changes in sediment redistribution (Gabet et al., 2003; Nkem et al., 2000; Wilkinson et al., 2009) and hence have the potential to affect surface topography and nutrient redistribution on large spatial and temporal scales. To quantify these effects, the shared role of climate, landscape characteristics, and burrowing dynamics on sediment redistribution needs to be understood.

On a local scale, currently used field methods to monitor sediment redistribution under real-life conditions are mainly erosion pins, splash boards, and rainfall simulators (Imeson and Kwaad, 1976; Wei et al., 2007; Le Hir et al., 2007; G. Li et al., 2019; T. C. Li et al., 2019; T. Li et al., 2018; Voiculescu et al., 2019; Chen et al., 2021; Übernickel et al., 2021a). The monitoring of box experiments yields a high spatiotemporal resolution and can also be linked to mathematical equations, such as random walks (Boudreau, 1986; Wheatcroft et al., 1990), stochastic differential equations (Boudreau, 1989; Milstead et al., 2007), finite-difference mass balancing (Soetaert et al., 1996; François et al., 1997), and Markov chain theory (Jumars et al., 1981; Foster, 1985; Trauth, 1998; Shull, 2001) to describe sediment redistribution.

Previously used methods have, however, several limitations when studying bioturbation. Field measurements likely lead to an underestimation of sediment fluxes, as they are one-time or seasonal measurements and thus do not capture the continuous excavation of the sediment by the animal (Grigusova et al., 2022) at a high temporal resolution. Box experiments and the mathematical equations derived from them describe bioturbation as an isolated process and ignore adjacent environmental parameters (such as climate or vegetation). However, the field measurements showed both positive (Hazelhoff et al., 1981; Black and Montgomery, 1991; Chen et al., 2021) and negative impact of bioturbation on erosion (Imeson and Kwaad, 1976; Hakonson, 1999). Also, previous field-based studies observed an increased bioturbation activity with higher (Milstead et al., 2007; Meserve, 1981; Tews et al., 2004; Wu et al., 2021; Ferro and Barquez, 2009) and lower vegetation cover (Simonetti, 1989; S. Zhang et al., 2020; Q. Zhang et al., 2019; Qin et al., 2021). Furthermore, soil mixing rates are not homogenous throughout the year; they depend on the animal phenological cycles (Eccard and Herde, 2013; Jimenez et al., 1992; Katzman et al., 2018; Malizia, 1998; Morgan and Duzant, 2008; Monteverde and Piudo, 2011; Gray et al., 2020; Yu et al., 2017).

Another approach offers raster-based soil erosion and landscape evolution models which integrate co-dependencies between bioturbation-relevant environmental parameters (Black and Montgomery, 1991; Meysman et al., 2003; Yoo et al., 2005; Schiffers et al., 2011). The most common soil erosion models are empirical (Wischmeier and Smith, 1978; Williams, 1975; Renard et al., 1991), process-based (Morgan et al., 1998; Roo et al., 1996; Nearing et al., 1989; Beasley et al., 1980), and semi-empirical models, the latter of which are a combination of both (Morgan et al., 1984; Beven and Kirkby, 1979).

Process-based models are based on a mechanistic understanding of the underlying physical, chemical, and biological processes that govern the behavior of the system being studied. They must be parameterized for each site; however, these models explicitly represent the governing equations and simulate the system's behavior by numerically solving these equations. Process-based models are generally considered to be more realistic and accurate than empirical models because they capture the fundamental processes that drive the system's behavior. However, process-based models can be computationally expensive, require more data and knowledge of system properties, and may require complex numerical algorithms (Morgan et al., 1998; Roo et al., 1996; Nearing et al., 1989; Beasley et al., 1980).

Within empirical models, on the other hand, the physical equations are completely replaced by empirically determined equations which only hold for the specific area they are derived for. These models are generally simpler, are less computationally expensive, and require more data and knowledge of system properties than process-based models. However, empirical models also tend to be less accurate than process-based models, particularly when applying beyond the range of data used to fit the model. In contrast to physical-based models, empirical models may not be applicable to new or different conditions, as they are based on observed relationships and do not capture the underlying processes that govern system behavior (Wischmeier and Smith, 1978; Williams, 1975; Renard et al., 1991).

Semi-empirical models combine the advantages of the both model types (Morgan et al., 1984; Morgan, 2001; Morgan and Duzant, 2008; Devia et al., 2015; Lilhare et al., 2015). Most landscape models have not yet implemented the impacts of bioturbators on water and sediment fluxes (Brosens et al., 2020; Anderson et al., 2019; Braun et al., 2016; Cohen et al., 2010, 2015; Carretier et al., 2014; Welivitiya et al., 2019). There are numerous models describing benthic soil mixing (Francois et al., 1997, 2002; Kadko and Heath 1984; Croix et al., 2002), biodiffusion caused by all invertebrate bioturbators (Meysman et al., 2005; Rakotomalala et al., 2015; Morris et al., 2006), and vertical soil mixing and lateral sediment redistribution caused by single invertebrate species (Orvain et al., 2006; Román-Sánchez et al., 2019; Orvain, 2003, 2005; Sanford, 2008). However, there are also models which described the impact of bioturbation on sediment redistribution by vertebrate animal species, such as the impact of pocket gophers on non-linear hillslope diffusion (Gabet, 2000) or on the creation of Mima mounds (Gabet et al., 2014). Several models include soil vertical mixing caused by bioturbation and its effect on landscape evolution on a millennial scale. This rather large spatiotemporal scale, however, means an omission of the natural variability in burrow sizes and densities, climate zones, and seasonality. In these models, soil erosion increases proportionally with increasing bioturbation, vertical soil mixing rates are uniform, and bioturbation is positively linked with vegetation cover (Temme and Vanwalleghem, 2016; Vanwalleghem et al., 2013; Yoo and Mudd, 2008; Pelletier et al., 2013). None of the previous studies included vertebrate bioturbator burrows of various sizes and spatial distribution by adjusting the soil properties and topography into a raster-based area-wide soil erosion model. This approach would enable us to understand the impact of all vertebrate bioturbators by considering the spatial distribution and variable impacts of bioturbator burrows on sediment redistribution. For this, bioturbation has to be included into erosion models at a spatial resolution which allows the imitation of the surface processes occurring within and near the burrow and at a temporal resolution which captures the animal daily burrowing behavior.

A suitable model which can be extended to include continuous bioturbating activity is the semi-empirical Morgan–Morgan–Finney soil erosion model (Morgan et al., 1984; Morgan, 2001). This model was successfully tested in several climate zones and land use types, such as Mediterranean sites (Jong et al., 1999); rainfed agrosystems, fields, and pastures (López-Vicente et al., 2008); the East African Highlands (Vigiak et al., 2005); and humid forests (Vieira et al., 2014). One of the recently developed improvements of this model is the daily Morgan–Morgan–Finney model (DMMF), which introduces subsurface flow and vegetation structures (type, size, height, root depth) and enables modeling at a high spatial (0.5 m) and temporal (daily) resolution (Choi et al., 2017). These improvements yield the potential to integrate the bioturbation into the model, as the burrowing activity is not constant and depends on vegetation structure (Tews et al., 2004; Ferro and Barquez, 2009).

In this study, we include vertebrate bioturbator burrows into a semi-empirical soil erosion model (DMMF) at a daily temporal and 0.5 m spatial resolution. For this, we predict the distribution of burrows by applying machine learning techniques in combination with using remotely sensed data as predictors. Then, we adjust soil properties, topography, and vegetation properties at predicted burrow locations based on field and laboratory measurements. We validate the model using field sediment fences. We run the model for a time period of 6 years, once with and without burrow adjustments. We estimate the impact of bioturbator burrows on sediment redistribution (including accumulation, erosion, and excavation) and surface runoff within four sites along the Chilean climate gradient. Lastly, we analyze how the impact of bioturbation on sediment redistribution depends on the burrow structure, climate, topography, and adjacent vegetation. Our study shows the importance of including bioturbation into erosion modeling and describes the interplay between bioturbation, environmental parameters such vegetation and topography, and sediment redistribution.

2 Study area

Our study was performed along a climate and ecological gradient in Chile (Übernickel et al., 2021b), comprising four study sites in the Chilean Coastal Cordillera: Pan de Azúcar (PdA) National Park (NP), Santa Gracia (SG), La Campana (LC) NP, and Nahuelbuta (NA) NP (Fig. 1). PdA NP is located in the arid zone in a fog-laden environment in the southern part of the Atacama Desert, with almost no rainfall. The vegetation cover is less than 5 % and dominated by small desert shrubs, several types of cacti, and biocrusts (Lehnert et al., 2018). SG is a natural reserve located in the semi-arid zone near La Serena, which is dominated by goat grazing. The vegetation consists of shrubs and cacti, covering up to 40 % of the study area. LC NP is part of the Mediterranean-type climate zone in the Valparaíso Region and is also affected by cattle. The study site is dominated by an evergreen sclerophyllous forest with endemic palms. The canopy reaches a height of up to 9 m, and the understory consists of deciduous shrubs and herbs. NA is located in the humid–temperate zone and characterized by a dense evergreen Araucaria forest comprising broadleaved trees with heights of up to 14 m. The ground is covered by bamboo, shrubs, and herbs (Bernhard et al., 2018; Oeser et al., 2018). The most common bioturbating vertebrate animal species recorded within these sites are carnivores of the family Canidae (Lycalopex culpaeus, Lycalopex griseus) as well as rodents of the families Abrocomidae (Abrocoma bennetti), Chinchillidae (Lagidium viscacia), Cricetidae (Abrothrix andinus, Phyllotis xanthopygus, Phyllotis limatus, Phyllotis darwini), and Octodontidae (Cerqueira, 1985; Jimenez et al., 1992; Übernickel et al., 2021a).

Figure 1Study area and study sites. Black lines outline the hillslope catchments. Along the blue lines, the in situ data (mound locations, soil samples, vegetation mapping) were collected. (a) Position of the study sites along the climate gradient. PdA – Pan de Azúcar, SG – Santa Gracia, LC – La Campana, NA – Nahuelbuta. Positions of plots in (b) PdA, (c) SG, (d) LC, and (e) NA. The background image is an RGB composite calculated from WorldView-2 satellite imagery. Images were obtained with a single license from GAF AG. Scale bar is the same for (b), (c), (d), and (e).

3 Methodology

We combined semi-empirical soil erosion modeling with in situ measurements, remote sensing data, and machine learning methods (Fig. 2). Along eight hillslope catchments within four climate zones, we mapped locations of burrows, estimated the vegetation cover, and extracted soil samples. We analyzed the soil samples in the laboratory. Then we used remote sensing datasets and machine learning to upscale burrow distribution, vegetation cover, and soil properties into the hillslope catchments. The hillslope catchment-wide predictions, the topographical information retrieved from lidar data (Kügler et al., 2022), and the climate information retrieved from climate stations were the input parameters for our soil erosion model. We ran the model with and without bioturbation. We included the bioturbation into the model by adjusting the input parameters at the predicted burrow locations. We also included continuous burrowing activity and soil mixing (Grigusova et al., 2021), the seasonality (Kraus et al., 2022), and the animal phenological cycle as found in Jimenez et al. (1992). The models were validated using self-constructed sediment traps. We studied the modeled surface runoff and sediment redistribution. Lastly, we analyzed whether and how the impact of bioturbation on sediment redistribution depends on environmental parameters (topography, landscape connectivity, and vegetation).

Figure 2Flowchart of our study. Green indicates in situ input data; blue indicates remote sensing input data. Red indicates model parametrization. Yellow indicates model output and analysis. Gray indicates model validation.


3.1 In situ data

The study setup consisted of eight hillslope catchments: one north-facing and one south-facing hillslope catchment per study site. We defined a line with a width of 1 m from the top to the base of each hillslope catchment (see blue line, Fig. 1). We subdivided the track into tiles of 1 m2. We saved the GPS information of each tile.

Within each tile of the line, we mapped burrow presence, land cover, and the extracted soil samples. A burrow consisted of an entrance and a mound (Fig. 3a). Each 1 m2 tile with a burrow was described as a presence data point and tiles without a burrow as absence data points. We noted the size of the burrow, vegetation cover, and land cover types (bare soil, herbs, shrubs, trees) within the tile. We extracted 162 soil samples from soil without a mound at a depth of 10 cm. Additionally, we took a photo of the surface every second tile along the track.

To validate the model output, we set up sediment traps (Fig. 3b), with six traps per site, two of which were located at the hillslope catchment base and four were located on two random positions within the hillslope catchment. The sediment traps consisted of geotextile and wooden poles and had a length of 2–5 m. A total of 1.5 m of geotextile was laid horizontally down at the surface, and 1 m of geotextile was vertically attached to wooden poles to enable the collection of sediment (Fig. 3b).

The sediment accumulated within the traps was collected after 1 year, and its mass (cm3) and dry weight (kg) were estimated.

Climate information was retrieved from climate stations located adjacent to the hillslope catchments, which provide climate data in 5 min intervals (Übernickel et al., 2021). To force the model on an hourly basis, hourly air temperature, precipitation total and intensity, wind speed, wind direction, and humidity were calculated for the study period from 1 April 2016 to 1 December 2021. Evapotranspiration was estimated with the Penman–Monteith equation (Penman, 1948).

Figure 3In situ constructions. (a) Example of a burrow consisting of burrow entrance and mound. (b) Fence construction used for the collection of eroded sediment to validate the model. Both photos by Paulina Grigusova.


3.2 Estimation of soil properties

We estimated several soil properties from the soil samples and photos collected in situ (Grigusova et al., 2022). We estimated the rock coverage on the surface and debris from the photos taken every second tile. For this, the photos were firstly classified into five classes. The classification was unsupervised using k means (Fig. A1). Then we calculated the ratio of pixels classified as skeleton and/or debris to the overall number of all pixels to determine the proportion of both parameters in percent.

In the lab, we estimated soil water content, bulk density, soil particle density, soil texture (sand; silt; clay; coarse, middle, and fine sand; coarse, middle, and fine silt), soil skeleton, organic matter, and organic carbon.

Gravimetric soil water content (%) (GSWC) described the mass of water within the soil sample and was estimated as in Eq. (1):

(1) GSWC = S m - S d S d × 100 ,

where Sm (g) is the mass of moist soil measured directly after the extraction and Sd (g) is the mass of soil dried at 105 C for at least 24 h. Bulk density (g cm−3) (BD) was calculated as follows:

(2) BD = S d S v ,

where Sv (cm−3) is the volume of the sample. Soil particle density (g cm−3) (SPD) was calculated as in Eq. (3):

(3) SPD = dm S v ,

where dm (g) is the dry mass of soil particles excluding pores.

Particle size distribution (%) of clay (< 0.002 mm); coarse, middle, and fine silt (0.002 to 0.02 mm); and coarse, middle, and fine sand (0.02 to 2 mm) was estimated according to Durner et al. (2017). Soil skeleton was estimated as the ratio of particles with a diameter above 2 mm. Ratio of organic matter (OM) was estimated as in Eq. (4):

(4) OM = 1 - S c S d ,

where Sc is the weight (g) of the sample dried at 500 C for 16 h.

We used pedotransfer functions to determine porosity, saturated soil moisture, hydraulic conductivity, water content at field capacity, and permanent wilting point. Pore ratio (θs) was estimated from bulk and particle density as in Eq. (5):

(5) θ s = BD SPD .

Saturated water content (g g−1) (Ws) was estimated as in Eq. (6):

(6) W s = θ s p w BD ,

where pw (g cm−3) is the density of water, which is set to be 1 g cm−3 (Pollacco, 2008).

Hydraulic conductivity Ks (m s−1) was estimated as in Eq. (7):

(7) K s = 1.15741 × 0.0000001 × exp ( x ) ,

where x for sandy soil is

(8) x = 9.5 - 1.471 × BD × BD - 0.688 × OM + 0.0369 × OM × OM - 0.332 × CS ,

and x for loamy and clayey soils is

(9) x = - 43.1 + 64.8 × BD - 22.21 × BD × BD + 7.02 × OM - 0.1562 × OM × OM + 0.985 × ln OM - 0.01332 × C × OM - 4.71 × BD × CS ,

where C is percentage of clay and CS is percentage of clay and silt (Wösten, 1997). To estimate water content at field capacity (%) (FC) and permanent wilting point (PWP), we applied functions by Tomasella et al. (2000) as these were developed for South American soils:


where Si is the percentage of silt.

3.3 Processing of remote sensing data

The digital elevation models (DEMs) were calculated from the lidar data (Kügler et al., 2022; Horn, 1981) at a resolution of 0.5 m. Slope was calculated according to Horn (1981). Manning's surface roughness coefficient was estimated following Li and Zhang (2001). The topographic position index (TPI) and the topographic ruggedness index (TRI) were calculated according to Wilson et al. (2007). To calculate the TPI, the average elevation of pixels within a range specified by the user needs to be subtracted from the elevation of the central pixel. Positive values represent hills, while negative values represent valleys. The TRI adds together the elevation differences between a grid cell and its eight neighbors. It measures the relative level of topography irregularity: the higher the value, the more irregular the topography. Plan and profile curvature were determined after Zevenbergen and Thorne (1987). Connectivity indices, sinks, wetness index, flow direction, flow path, catchment slope, and catchment were calculated in SAGA GIS.

Single license stereo WorldView-2 images with a resolution of 0.5 m were retrieved from GAF AG Munich GmbH. The topographic correction of WorldView-2 images was done using the lidar data, solar elevation angle, solar zenith angle, and azimuth angle according to Goslee (2012). The digital surface models (DSMs) were calculated from the stereo images. Additionally, we extracted single bands and calculated the normalized difference vegetation index (NDVI).

3.4 The erosion model

3.4.1 Daily Morgan–Morgan–Finney model

The DMMF model is a combined soil erosion model used to estimate surface runoff and sediment flux on a field scale on a daily basis. Spatially, the DMMF model represents an area as several interconnected elements (e.g., pixels) of uniform topography, soil characteristics, land cover type, and vegetation structure. Through coupling, the model operates with flow direction algorithms: each element receives water and sediments from upslope elements and delivers the generated surface runoff and eroded soils to downslope elements. On a temporal scale, the model estimates surface runoff and sediment flux of each element on a daily basis. The model input parameters include climate, topography, soil properties, and land cover information (Choi et al., 2017). Data pre-processing, modeling, and analysis (see Fig. 2) were done in the R statistical environment. The raster data were cropped to the size of the hillslope catchments (Fig. 1). Input parameters are listed in Table 1 and plotted in Fig. A2.

During the model simulation, water and sediment are transferred from pixels located at higher elevations to pixels situated at lower elevations. This occurs in two stages: the first stage is the hydrological phase where the model calculates surface runoff, which happens when the amount of surface water input exceeds the water-holding capacity. The amount of surface runoff is computed by taking the infiltration capacity of the surface, the volume of surface water input, and the fraction of the impervious area of a pixel into account. Infiltration capacity represents the maximum amount of surface water that can penetrate the subsurface layer. It is determined by the percentage of the impervious area and the available pore space.

The second stage is the sediment phase, where the model estimates the sediment budget for each particle size class, based on the surface conditions. The model calculates the detachment and deposition of sediments in a step-by-step process. The sources of sediments are detached particles from the pixel itself due to rainfall and surface runoff and delivered soil particles from higher-elevation pixels. The detachment of soil particles by rainfall occurs when raindrops hit the ground with enough energy to detach soil particles from the surface. Rainfall has different impacts on areas with and without canopy cover, as canopy cover changes the kinetic energy of raindrops.

The quantity of soil particles detached by raindrops is calculated based on the soil particle detachability, the percentage of each particle size class, the bare soil surface area, and the kinetic energy of effective rainfall. The quantity of detached soil particles by surface runoff is calculated based on the soil particle detachability, the amount of runoff, the slope angle of the pixel, and the proportion of the bare surface area. The third source of sediment is from higher-elevation pixels and is averaged by the surface area of the pixel.

Once sediments are delivered to the surface runoff, a portion of the suspended sediments settle to the bottom due to gravitational force. To calculate this settling, the model requires the flow velocity of the runoff and the settling velocity of each particle size class, which are influenced by the flow depth, slope angle of the pixel, and Manning's roughness coefficient (Choi et al., 2019).

3.4.2 Estimation of spatial parameters

For spatial parameterization of the DMMF model, we predicted land cover, soil properties, and burrow distribution onto the hillslope catchments using machine learning techniques. We used the approach of Meyer et al. (2018). The most important predictors were selected by forward feature selection. The quality of the random forest (RF) models was assessed by leave-location-out cross-validation. We trained the model stepwise, using in situ data collected from seven of the hillslope catchments and validated the model using in situ data from the remaining hillslope catchment (Meyer et al., 2018). The prediction was done at 0.5 m spatial resolution. We used the WorldView-2 layers obtained with a single license with GAF, NDVI, DEM, DSM, slope, and roughness as predictors. The PAN-sharpening of the WV-2 layers was done by GAF. The accuracy of the classifications was estimated by dividing the number of correctly classified pixels to the number of all pixels.

For the area-wide prediction of burrow locations across the hillslope catchments, we used the burrow presence and absence data (Sect. 3.1) as the response data within the RF models. The accuracy was 0.82 for PdA, 0.77 for SG, 0.75 for LC, and 0.85 for NA. The prediction of soil properties was done using soil properties estimated along the track line (see Sect. 3.1) as response data within the RF models. All of the models reached a high accuracy (see Table A1).

To obtain land cover classification, we used, as the response within the RF models, the land cover measured in situ. The classes were soil without rocks, rocks, biocrusts, grass/herbs, shrubs, and trees. Predictor values for each class were extracted from at least 100 polygons per site and class. The accuracy of the RF models was 0.71 for PdA, 0.81 for SG, 0.83 for LC, and 0.75 for NA.

The vegetation height measured in plots was averaged for each class per site. All pixels classified as the respective class were assigned the same vegetation height information. Vegetation density was estimated per hillslope catchment as the number of vegetation individuals per square meter. Vegetation diversity was calculated with the Shannon index (Shannon, 1948). The interception area was the area not covered by vegetation (herbs, shrubs, or trees).

3.4.3 Inclusion of bioturbation

In the grid cells with predicted burrow locations, we adapted the values of input parameters to include bioturbation. The adaptations varied with climate zone and burrow size. The size, geometric structure, and excavation rates of burrowing animals were previously estimated at a high spatial and temporal resolution (Grigusova et al., 2022). Based on these results, we firstly adjusted the microtopography. We modified the layer depth to represent burrow entrance and elevation to represent animal mound. Mounds were always located downslope of burrow entrances in the direction of flow.

Secondly, we adjusted the soil properties. The soil properties texture and organic carbon were estimated from soil extracted from mounds in Kraus et al. (2022). In this study we additionally estimated bulk density, initial water content, soil skeleton, porosity, saturated water content, available water capacity, and water content at field capacity from the same dataset (see Sect. 3.2). We calculated the median value of each property for the samples extracted from mounds and for the samples extracted from soil without mounds. Then, we estimated the change in percent between these two values. This was then used to adjust the soil property for each pixel including a mound.

Thirdly, modeled mound pixels had to be cleared from ground vegetation cover. For this, we removed ground vegetation cover from pixels with burrow locations and decreased ground vegetation cover, height, diameter, and number of ground vegetation individuals from adjacent pixels as measured in situ. Then, the number of rocks and amount of debris were set as estimated from soil samples (Sect. 3.2)

Animal activity has been found to be highly variable throughout the year (Grigusova et al., 2022; Kraus et al., 2022). The density of burrows does not stay stable throughout the year but increases or decreases depending on the season and climate zone. We therefore artificially removed or added burrows into the hillslope catchments during the particular seasons. For this, we adapted the density of soil, the topography, and vegetation cover accordingly. We created a 3D model of the burrow structure and adjusted subsurface soil properties and properties of soil excavated to the surface: the removed vegetation within the pixel with a predicted burrow and decreased adjacent vegetation cover.

Lastly, we also included the vertical movement of sediment particles from deeper soil layers to the surface depending on climate. Animals were found to reconstruct their burrows after each rainfall event (Grigusova et al., 2022). Corresponding with these findings, we increased the entrance depth and mound height by 30 % after each rainfall event, which represents the averaged value found in the previous study (Grigusova et al., 2022).

For the validation, we ran the model for the time periods between the installation of sediment fences and the collections of sediment. We compared the mass and weight of modeled and collected sediment and estimated R2 and RMSE. To test the importance of the inclusion of individual bioturbation parameters into the model, we ran the model under four conditions: (i) no burrows, (ii) solely entrances, (iii) solely mounds, and (iv) entire burrows (entrances and mounds).

Table 1Model input layers and respective changes to layer values at the predicted burrow locations. Ground vegetation was removed from the respective pixels, while tree canopy was not changed. The values were estimated as described in Sect. 3.5.2. Using the adjusted values, we calculated evapotranspiration using the Penman–Monteith equation, surface roughness from the elevation layer, hydraulic conductivity, water content at field capacity, and saturated water content using pedotransfer functions.

The ns means not significant.

Download Print Version | Download XLSX

3.5 DMMF model sensitivity test

We conducted a sensitivity test to identify those input parameters which significantly influence the model output. For this, we first estimated the mean value of each input parameter. Then, we created an artificial hillslope catchment of 100 m × 100 m. To start the test, each pixel received the mean value of each parameter. We ran the model for one rainfall event. Then, stepwise, we changed the single input parameter values from their minimum to their maximum values while not adjusting any other parameters. To quantify the significance of the input variations, we conducted a t test (Table A2). For this, we compared the amount of redistributed sediment of each model run to the first model run.

3.6 Impact of burrows on surface processes

We estimated burrow density as the ratio of pixels with predicted burrows to all pixels. Additionally, we calculated the ratio of pixels which are part of a burrow aggregation to all pixels which include a burrow. Burrow aggregation describes at least four neighboring pixels with predicted burrows. We calculated the amount of excavated sediment as a sum of burrow density and the burrow excavation rate as estimated in Grigusova et al. (2022).

To estimate the impact of burrows on sediment redistribution and surface runoff, we ran the DMMF model for the time period from 1 April 2016 to 31 December 2021 for all hillslope catchments. We ran the model (i) with no burrows and (ii) with entire burrows. We estimated (i) sediment redistribution (accumulation minus erosion) and (ii) surface runoff. We analyzed the redistribution and runoff on the plot (1 m2) and hillslope catchment (1 ha) scale.

Lastly, to analyze under which biotic and abiotic environmental parameters (topography, vegetation cover) the bioturbation enhances sediment erosion or accumulation, we set up a generalized additive model (GAM) (Wood, 2006). For this, we first subtracted the output of the model with no burrows from the output of the model with entire burrows. Within each pixel, two processes are happening simultaneously: a certain amount of sediment erodes, and a certain amount of sediment accumulates. To estimate the sediment redistribution for each pixel of each model run, we estimate which of these processes dominated. Positive pixel values thus mean that bioturbation enhanced sediment accumulation; negative pixel values mean that bioturbation enhanced sediment erosion. We tested the following environmental parameters: mound density, vegetation cover, elevation, slope, aspect, TRI, TPI, curvature and connectivity, and wetness index. The model performance was evaluated by the percentage of explained data variance. We analyzed the impact of environmental parameters within 1 m and within 10 m from the burrows.

4 Results

4.1 Model sensitivity test and accuracy

Parameters which significantly influenced the model output were precipitation, slope, vegetation cover, surface roughness, silt content, and water content (Table A2). There was a correlation between some of the spatial model parameters (Fig. A10), especially between the initial and saturated water content, between water content and vegetation cover, and between clay content and field capacity. However, a high correlation between spatial parameters does not mean that these parameters impact the sediment redistribution in a similar way.

We quantified the model performance by comparing the modeled and measured sediment redistribution. The performance varied depending on the burrow inclusion (Figs. 4 and 5). The performance of the model without any bioturbation was lower (R2= 0.73, RMSE = 1.50, MSE = 2.27), as when burrow entrances (R2= 0.81, RMSE = 1.34, MSE = 1.16) or mounds (R2= 0.83, RMSE = 1.10, MSE = 1.22) were included. The model had the highest performance when entire burrows were included (R2= 0.85, RMSE = 1.01, MSE = 1.01). However, as the scatterplots showed, the model performance seemed to be determined strongly by one measurement (Fig. 5). For this reason, we calculated the metrics without this measurement (Fig. A2). The model without any burrows (R2= 0.17, RMSE = 1.18, MSE = 1.39) in this case performed much lower than models with burrows. The model performance increased when burrow entrances (R2= 0.48, RMSE = 0.61, MSE = 0.78) or mounds (R2= 0.51, RMSE = 0.75, MSE = 0.57) were included. The model with whole burrows reached the highest performance (R2= 0.71, RMSE = 0.63, MSE = 0.39). When we compare the modeled redistribution to the sediment redistribution estimated using time-of-flight cameras in Grigusova et al., (2022), the differences appear to be minor (R2= 0.62, RMSE = 0.12, MSE = 0.35).

Figure 4R2 and RMSE of the Morgan–Morgan–Finney soil erosion model. For dataset A, we compared the amount of sediment collected in all sediment fences with the modeled eroded sediment (see Fig. A3). For dataset B, we removed one measurement, as the R2 seemed to be defined by this measurement (see Fig. A4). For scenario A, we did not include any burrows into the model. For scenario B, we included burrow entrances, and for scenario C, we included mounds. For scenario D, we included whole burrows into the model. The adjustments made to include entrances, mounds, and burrows into the model are described in Sect. 3.5.2.


Figure 5Measured and modeled redistributed sediment without an outlier. (a) Model without bioturbation. (b) Model with entrances. (c) Model with mounds. (d) Model with burrows.


4.2 Model output: surface runoff and sediment redistribution

Hillslope catchment-wide sediment redistribution (1 ha resolution) was the highest in humid NA, followed by Mediterranean LC, semi-arid SG, and arid PdA (Figs. 6a, b, 8). In NA, LC, and SG, the erosion processes dominated, while in PdA, more sediment accumulated than eroded. The impact of burrows on sediment redistribution was significant in arid PdA, semi-arid SG, and Mediterranean LC. Burrows increased sediment redistribution by 137.8 % ± 16.4 % in arid PdA (3.53 kg ha−1 yr−1 vs. 48.79 kg ha−1 yr−1), by 6.5 % ± 0.7 % in semi-arid SG (129.16 kg ha−1 yr−1 vs. 122.05 kg ha−1 yr−1), and by 15.6 % ± 0.3 % in Mediterranean LC (4602.69 kg ha−1 yr−1 vs. 3980.96 kg ha−1 yr−1). Overall, bioturbation increased sediment accumulation in the arid zone (as the magnitude of the sediment excavation by the animals exceeded sediment erosion which occurs during rainfall events), but it increased sediment erosion in semi-arid and Mediterranean climate (where animal burrowing activity and rainfall are present). The largest impact was found under Mediterranean conditions. We found no significant effect on redistribution in the humid zone (Fig. 7). However, impact of bioturbation varied throughout the hillslope catchment (Figs. 7, 8, and 9).

Surface runoff was the highest in humid NA, followed by Mediterranean LC, arid PdA, and semi-arid SG (Fig. 6c). The impact of burrows on surface runoff was significant in all climate zones. Burrows increased surface runoff in PdA by 34 %, in SG by 40%, and in LC by 4.1 %, but it decreased surface runoff by 5.9 % in NA. Hillslope catchment-wide maps are shown in Figs. A6–A8.

Figure 6Summary of model outputs across the climate gradient. PdA is arid Pan de Azúcar. SG is semi-arid Santa Gracia. LC is Mediterranean La Campana. NA is humid Nahuelbuta. Graphs (a) and (b) show the modeled sediment redistribution. Positive values indicate sediment accumulation; negative values indicate sediment erosion. In (a) sediment redistribution is shown on a pixel scale (in kg m−2 yr−1), while in (b) sediment redistribution is shown on the hillslope catchment scale (in kg ha−1 yr−1). The impact of bioturbation on sediment redistribution was estimated by a t test and was significant in three sites: PdA***, SG**, and LC***. Bioturbation increased sediment redistribution by 137.8 % in PdA, by 6.5 % in SG, and by 15.6 % in LC. For hillslope catchment-wide maps see Figs. A6–A8. Graph (c) represents the modeled surface runoff on the hillslope catchment scale (in m3 ha−1 yr−1). The impact of bioturbation on surface runoff was estimated by a t test and was significant at all sites. Bioturbation increased surface runoff in PdA by 34 %, in SG by 40 %, and in LC by 4.1 %, but it decreased surface runoff by 5.9 % in NA. For hillslope catchment-wide maps, see Fig. A6.


Figure 7Comparison of the model outputs with and without bioturbation of each pixel (0.5 m) in all study sites. The x axis shows the output of the model with bioturbation, and the y axis shows the model output without bioturbation. PdA is arid Pan de Azúcar. SG is semi-arid Santa Gracia. LC is Mediterranean La Campana. NA is humid Nahuelbuta. Points represent single pixel values; lines show linear regressions for the sites. The lower the R, the higher the impact of burrows is on sediment redistribution at the resolution of 0.5 m. The dashed black line symbolizes a perfect correlation – along this line the bioturbation would have no effect on sediment redistribution. Bioturbation leads to more accumulation if the regression line representing results from a particular climate zone is steeper than the perfect correlation line. Bioturbation leads to more erosion if the regression line representing results from a particular climate zone is flatter than the perfect correlation line. Bioturbation increases sediment accumulation in arid PdA (through the high burrowing rate, more sediment is accumulated on the surface than eroded during rainfall events). Bioturbation increases sediment erosion in semi-arid SG and Mediterranean LC. In absolute terms, the highest impact on sediment redistribution is in the Mediterranean climate zone. The lowest impact is in the humid zone.


Figure 8Hillslope catchment-wide predicted sediment redistribution. Colors indicate sediment redistribution. Gray indicates the hill shading calculated from lidar data. (a) Pan de Azúcar, (b) Santa Gracia, (c) La Campana, and (d) Nahuelbuta.


4.3 Role of continuous burrowing activity on sediment redistribution

We included transport of the sediment to the surface by animal excavation into the model. The density of burrows was the highest in the arid PdA, then Mediterranean LC, and then semi-arid SG, and it was the lowest in humid NA. Burrows were mostly distributed within groups of several burrows in Mediterranean LC and semi-arid SG, while they were more evenly distributed in the arid PdA and humid NA. The burrows were of the largest in Mediterranean LC, followed by arid PdA, semi-arid SG, and humid NA. Similarly, the highest volume of excavated sediment at the beginning of the modeling period was in Mediterranean LC and arid PdA. The volume of excavated sediment during the burrow reconstruction after rainfall events was the highest in humid NA, followed by Mediterranean LC, semi-arid SG, and arid PdA. The percentage of sediment excavated by the animal to sediment redistributed during rainfall events was 128 % in PdA, 24 % in SG, 33.5 % in LC, and 5.6 % in NA.

Table 2Impact of animal bioturbation activity on overall sediment redistribution on various scales. The bioturbation activity was estimated using time-of-flight-based cameras in Grigusova et al. (2022). This study showed that animals reconstruct their burrows after each rainfall event. During this process, 10 % of the overall sediment burrow volume is relocated from within the burrow to the surface. We integrated this process into our model and calculated the percentage of newly excavated sediment by the animals to the amount of sediment which was redistributed during rainfall for the period of 1 year.

Download Print Version | Download XLSX

4.4 Role of adjacent environment

We subtracted the output of the model with included burrows from the output of the model without burrows (Fig. A8). Although the burrows on average enhanced sediment erosion on the hillslope catchment scale, the high-resolution maps unveiled that burrows enhance sediment erosion within some pixels, while they rather increased sediment accumulation within others.

The amount of data variance explained by the GAM models (see Sect. 3.6.) differed between models (Table A3). Models estimating the impact of environmental parameters on sediment redistribution within 1 m from the burrows explained 3.84 % of the variance in PdA, 37.1 % in SG, 46 % in LC, and 42. % in NA. Models estimating the impact of environmental parameters on sediment redistribution within 10 m from the burrows explained 1.99 % of the variance in PdA, 12.8 % in SG, 52 % in LC, and 72.9 % in NA. The parameters selected for SG were slope, roughness, curvature, TRI, and NDVI. Parameters selected for LC were elevation, slope, NDVI, sinks, and roughness. Parameters selected for NA were elevation, slope, aspect, TRI, sinks, and roughness (Fig. 10).

Bioturbation strongly increased sediment redistribution (erosion and accumulation) at high values of elevation, slope, surface roughness, TRI, sinks, and topographic wetness index; at the middle values of elevation and aspect; and at low values of profile curvature and NDVI. From these parameters, bioturbation increased sediment erosion at high and middle values of elevation; at high values of slope, sinks, and TRI; and at low values of profile curvature. Bioturbation increased sediment accumulation at high values of surface roughness and topographic wetness index and at low values of NDVI (Figs. A3–A8).

Bioturbation somewhat enhanced sediment erosion at medium values of surface roughness, NDVI, and sinks and at low values of topographic wetness index. Bioturbation somewhat increased sediment accumulation at low values of slope and TRI, at low and medium values of elevation, and at high values of profile curvature.

Figure 9Hillslope catchment-wide impact of bioturbation on sediment redistribution. Color indicates the impact. Positive values indicate that bioturbation enhanced sediment accumulation; negative values indicate that bioturbation enhanced sediment erosion. Gray indicates the hill shading calculated from lidar data. (a) Pan de Azúcar, (b) Santa Gracia, (c) La Campana, and (d) Nahuelbuta.


Figure 10This figure is a conceptual summary of the detailed results from Figs. A3–A8. Bioturbation increases erosion or accumulation depending on the values of environmental parameters. The dependencies are the same for all climate zones. The figure is the conceptual summary for all climate zones; therefore, there are no values stated on the x and y axes. The x axis shows whether bioturbation increases erosion or accumulation. The y axes are environmental parameters. Line thicknesses indicate the magnitude of impact. Please note that bioturbation has no impact on sediment redistribution in regions with low sink connectivity and topographic ruggedness. The relationship between the values of environmental parameters and the impact of bioturbation is not linear: bioturbation can have the same impact on sediment redistribution at high or low values of an environmental parameter but a contrasting impact at middle values of this parameter (as in this case for elevation, slope, or surface roughness).


5 Discussion

5.1 The inclusion of bioturbation increases model performance

Overall, our DMMF model including bioturbation performed much better than the model without bioturbation. The DMMF model without bioturbation performed worse (RMSE of 1.18 kg ha−1 yr−1 and R2 of 0.17) than the model with bioturbation (RMSE was 0.63 kg ha−1 yr−1 and R2 was 0.71).

We hence argue that the higher accuracy of our model can be explained with the inclusion of bioturbation. This is confirmed by the fact that our model run without bioturbation performed similarly to previously run models without bioturbation: in earlier studies, the accuracy of the MMF model reached an RMSE between 4.9 and 8.2 kg ha−1 yr−1, with an estimated R2 of between 0.21 and 0.57 (Jong et al., 1999; Vigiak et al., 2005; López-Vicente et al., 2008; Vieira et al., 2014; Choi et al., 2017). However, we acknowledge that previous studies were all conducted in more temperate climate zones. To be able to compare our results with previous studies, we calculated the model performance considering solely the Mediterranean and humid climate zone, which are more similar in climate to the more temperate locations of previous studies. The performance of the model was still high (R2= 0.72, RMSE = 0.45 kg ha−1 yr−1), confirming the conclusion that bioturbation increased model performance.

We compared the modeled impact of bioturbation on sediment redistribution with the impact of bioturbation estimated in previous studies. In the humid zone, our model predicted an erosion up to 3.5 kg m−2 yr−1. This estimation is in line with erosion rates established by in situ measurements in other studies conducted in a more humid climate zone (between 1.5 and 3.7 kg m−2 yr−1) (Black and Montgomery, 1991; Yoo and Mudd, 2008; Yoo et al., 2005; Rutin, 1996). This also confirms the reliability of our approach. Previous authors estimated the impacts using rainfall simulators, erosion pins, or splash boards. The measurements were conducted for a time period between 3 months and 3 years, and the sites were revisited for each estimation. We do not compare our results with studies which previously applied models to estimate impacts of bioturbation, as, to our knowledge, none of the previous studies integrated vertebrate burrow structures into a soil erosion model and ran the model on a daily basis.

5.2 The relevance of bioturbation for sediment redistribution depends on the environmental context

On the hillslope catchment scale (1 ha), our study finds that bioturbation increases erosion in semi-arid and Mediterranean zone, increases accumulation in the arid zone, and has no impact within the humid zone (Fig. 6b). In contrast, bioturbation increases both erosion and accumulation on the plot scale (1 m2) (Fig. 6a). On this scale, in the arid and semi-arid zone, sediment erosion and accumulation were predicted to be about equal (erosion and accumulation both up to 0.1 kg m−2 yr−1 in the arid zone and erosion and accumulation both up to 0.2 kg m−2 yr−1 in the semi-arid zone; see Fig. 6a). Bioturbation marginally increased erosion and decreased accumulation in the semi-arid zone but reduced accumulation 2-fold in the arid zone. In contrast, in the Mediterranean and humid zone, erosion was predicted to be almost double when compared to accumulation (predicted erosion up to 2.5 kg m−2 yr−1 and accumulation up to 1.4 kg m−2 yr−1). Inclusion of bioturbation increased erosion up to 3 kg m−2 yr-1, and accumulation up to 1.6 kg m−2 yr−1 in the Mediterranean zone, while it had no significant effect in the humid zone. We argue that sediment redistribution due to bioturbation is heavily influenced by mesotopographic structures which determine the flow path of surface runoff and influence the infiltration processes. Due to this, the erosion and accumulation on the plots scale are more heavily impacted by bioturbation with increasing surface runoff.

Our study found an increase of erosion in the semi-arid and Mediterranean climate zone to be between 6.5 % and 15.6 % due to bioturbation. Previous studies found that even a small increase of erosion has significant impacts on the whole hillslope catchment. A 10 % increase in erosion rates over a 10-year period can lead to significant changes in the landscape, including, for example, a 20 %–30 % reduction in soil thickness and an increase in sediment transport in nearby rivers (Kuhn, 2016).

According to our analysis, bioturbation increases erosion or accumulation of sediment mostly based on an interplay between topographic structures elevation, slope, and TRI (Fig. 10). Over all research sites, this study found that bioturbation leads to an increase in surface erosion in areas where erosional processes dominate (upper and/or steeper slopes) and tends to increase sediment accumulation in areas where sediment is naturally deposited (e.g., lower slopes or shallow depressions; Fig. 10). This finding is based on the fact that erosion in general is positively affected by slope and negatively by surface roughness and vegetation (Rodríguez-Caballero et al., 2012; Wang et al., 2013; Kirols et al., 2015). Additionally, the redistribution of sediment is largely affected by topographic meso-/macroforms, such as rills or cliffs. These can be quantified by topographic ruggedness index (TRI), which describes the amount of elevation drop between adjusting cells of DEM (Wilson et al., 2007). At high values of this index, we would therefore expect high erosion rates, due to concentrated runoff within the connected rills or undisturbed flow of runoff from the cliffs downslope.

Our data show that one burrow provides up to 0.43 m3 of additional loose sediment at the surface (Table 2), while the surface roughness increases up to 200 % (Grigusova et al., 2022). When including burrows into the model, with slope values from 0 to 5, the presence of burrows had no impact on sediment redistribution. From 5 onwards it increased sediment erosion proportionally to the slope of the hillside (an increased erosion from 0.4 g ha−1 yr−1 in the semi-arid zone to 150 kg ha−1 yr−1 in the Mediterranean zone; Figs. A3–A6). Similarly, at locations with elevation drops ranging from 0 m to 0.2 m (lower TRI values), the presence of burrows had no impact. However, at locations with elevation drops of 0.2 to 0.5 m (higher TRI values), bioturbation increases sediment erosion by 1.5 kg ha−1 yr−1 (Figs. A3–A8). Lastly, bioturbation proportionally increased accumulation when the surface roughness values were above 0.5 (an increased accumulation from 0.2 g ha−1 yr−1 in semi-arid zone to 5000 kg ha−1 yr−1 in the Mediterranean zone; Figs. A3–A6).

We conclude that, in locations with slope values over 5 or at locations with sudden drops in elevation (high TRI) and connected rills, more sediment is eroding than accumulating. Here, additional surface sediments generated by bioturbators provide more source material for erosion, and thus bioturbation increases sediment erosion at these locations (Figs. 10 and 11). In contrast, at locations with a slope below 5, where processes are dominantly controlled by surface roughness, sediment accumulation caused by bioturbation increases proportionally when the surface roughness has a value above 0.5. This is likely because burrows through their above-ground structures heavily increase surface roughness (Grigusova et al., 2022), and hence the presence of bioturbating animals leads to an increase in sediment accumulation.

Additionally, we hypothesize that it is not only the additional availability of sediment on the surface and the topography of the vicinity which controls the contribution of bioturbation to sediment surface flux, but also the spatial distribution of animal burrows. We interpret that, in locations with high burrow aggregation, surface flow might be redirected and centralized around the aggregates and thus increase sediment erosion in the areas adjacent burrow aggregates (Fig. 11). This mechanism could explain why bioturbation promotes sediment erosion especially in the Mediterranean zone, where burrows are more aggregated. The relative role of burrow aggregation should be studied in detail and included in future studies.

Figure 11Context dependency of sediment redistribution. (a) Pan de Azúcar, (b) Santa Gracia, (c) La Campana, and (d) Nahuelbuta. Brown arrows indicate the direction and magnitude of overall sediment redistribution within each climate zone. Blue arrows indicate the direction of flow (runoff vs. infiltration). Semicircles indicate the distribution and size of the burrows. The dashed line indicates the median value of each parameter for the first four parameters.


6 Conclusions

Our study found that the inclusion of vertebrate bioturbators' burrows into a soil erosion model significantly increases its reliability. Vertebrate bioturbators increase sediment accumulation in the arid climate zone, increase sediment erosion in the semi-arid and Mediterranean zone, and have no impact on sediment redistribution in the humid zone. Our study furthermore shows that the impact of bioturbation heavily depends on adjacent environmental parameters. The burrows increase sediment erosion at high and low values of elevation; at high values of slope, sink connectivity, and topography ruggedness; and at low values of vegetation cover. The burrows increase accumulation at high values of surface roughness and soil wetness. This means that, overall, on geological timescales, as burrowing animals increase both erosion in steeper zones and accumulation in areas with gentler slopes and higher roughness, hillslope relief should become more quickly equalized and overall more flat. This tendency is most pronounced in the Mediterranean zone with high burrow density and excavation rates, as well as comparably high precipitation rates.

Appendix A

Table A1R2 and RMSE of random forest models trained for the prediction of soil properties needed for model parametrization. RMSE is root mean square error.

Download Print Version | Download XLSX

Table A2Model sensitivity analysis. For the analysis, the minimum, maximum, and mean values of each parameter were calculated. The model was run for a hillslope catchment of 1 km2 with homogenous mean parameters. Then, the minimum and maximum values of each parameter were tested. Each parameter was changed stepwise to its minimum or maximum value while the remaining parameters stayed homogenous. The significance of the parameter was estimated by a t test conducted between the erosion estimated by the model with homogenous mean parameters and the erosion estimated by the model with varying minimum and maximum parameter values. Only significant parameters are shown.

Download Print Version

Table A3Summary of GAM models. We analyzed the impact of parameters within a 1 and 10 m from burrows. The asterisks indicate p values of the selected parameters. **** p<0.001. *** p<0.01. ** p<0.05. * p<0.1. One GAM model was run per parameter. Only results for models with an explained variance above 5 % are shown.

Download Print Version | Download XLSX

Figure A1An example of the unsupervised k-means classification of the surface photo from La Campana. Original photo was taken by Paulina Grigusova. The collection of in situ data is explained in Sect. 3.1. and the estimation of soil properties in Sect. 3.2. The image was classified into five classes using unsupervised k-means classification; the land cover was then assigned manually. In some cases, like in this case for rocks, multiple k-means classes stand for the same land cover. These were then unified to the class “rocks”.


Table A4Review of studies which integrated any kind of bioturbation into models. Previous models integrated either benthic, invertebrate, or single species of vertebrate bioturbators. Models applied described either the vertical soil mixing or long-term landscape evolution models. None of the previous studies included vertebrate burrows of bioturbators into an erosion model which would be capable of capturing the daily redistribution processes.

Download Print Version | Download XLSX

Figure A2Measured and modeled redistributed sediment for different scenarios. (a) Model without bioturbation. (b) Model with entrances. (c) Model with mounds. (d) Model with burrows.


Figure A3Environmental parameters influencing impact of bioturbation on sediment redistribution in Santa Gracia within 1 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A4Environmental parameters influencing impact of bioturbation on sediment redistribution in Santa Gracia within 10 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A5Environmental parameters influencing impact of bioturbation on sediment redistribution in La Campana within 1 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A6Environmental parameters influencing impact of bioturbation on sediment redistribution in La Campana within 10 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A7Environmental parameters influencing impact of bioturbation on sediment redistribution in Nahuelbuta within 1 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A8Environmental parameters influencing impact of bioturbation on sediment redistribution in Nahuelbuta within 10 m from the burrows. Positive values indicate bioturbation enhances sediment accumulation at the respective parameter values; negative values indicate bioturbation enhances sediment erosion at the respective parameter values.


Figure A9Burrow aggregation concentrates the runoff and increases erosion. An example of a north-facing hillside in Mediterranean La Campana for the time period of 1 year. (a) Sediment erosion as estimated by the model without bioturbation. (b) Sediment erosion as estimated by the model with bioturbation. (c) Sediment erosion as estimated by the model with bioturbation with predicted burrow locations. (d) Surface runoff as estimated by the model without bioturbation. (e) Surface runoff as estimated by the model with bioturbation. (f) Surface runoff as estimated by the model including bioturbation and predicted burrow locations. Black indicates that at least one burrow was located within this pixel. Four neighboring pixels which contain a burrow form a burrow aggregation.


Figure A10Correlation matrix between the model input parameters.


Author contributions

PG set up the model, analyzed the data, and wrote the manuscript draft; PG and AL performed the measurements; AL, JB, NF, RB, DK, PP, LP, and CdR reviewed and edited the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Earth surface shaping by biota (Esurf/BG/SOIL/ESD/ESSD inter-journal SI)”. It is not associated with a conference.


We thank CONAF for the kind support provided during our field campaign. We thank our two reviewers for their valuable comments that helped to significantly improve the manuscript.

Financial support

This study was funded by the German Research Foundation, DFG (grant nos. BE1780/52-1, LA3521/1-1, FA 925/12-1, and BR 1293-18-1) and is part of the DFG Priority Programme SPP 1803: EarthShape: Earth Surface Shaping by Biota, sub-project “Effects of bioturbation on rates of vertical and horizontal sediment and nutrient fluxes”.

Review statement

This paper was edited by Todd A. Ehlers and reviewed by Emmanuel Gabet and one anonymous referee.


Anderson, R. S., Rajaram, H., and Anderson, S. P.: Climate driven coevolution of weathering profiles and hillslope topography generates dramatic differences in critical zone architecture, Hydrol. Process., 33, 4–19,, 2019. 

Beasley, D. B., Huggins, L. F., and Monke, E. J.: ANSWERS: A Model for Watershed Planning, T. ASAE, 23, 938–944,, 1980. 

Bernhard, N., Moskwa, L.-M., Schmidt, K., Oeser, R. A., Aburto, F., Bader, M. Y., Baumann, K., Blanckenburg, F. von, Boy, J., van den Brink, L., Brucker, E., Büdel, B., Canessa, R., Dippold, M. A., Ehlers, T. A., Fuentes, J. P., Godoy, R., Jung, P., Karsten, U., Köster, M., Kuzyakov, Y., Leinweber, P., Neidhardt, H., Matus, F., Mueller, C. W., Oelmann, Y., Oses, R., Osses, P., Paulino, L., Samolov, E., Schaller, M., Schmid, M., Spielvogel, S., Spohn, M., Stock, S., Stroncik, N., Tielbörger, K., Übernickel, K., Scholten, T., Seguel, O., Wagner, D., and Kühn, P.: Pedogenic and microbial interrelations to regional climate and local topography: New insights from a climate gradient (arid to humid) along the Coastal Cordillera of Chile, CATENA, 170, 335–355,, 2018. 

Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant, Hydrol. Sci. Bull., 24, 43–69,, 1979. 

Black, T. A. and Montgomery, D. R.: Sediment transport by burrowing mammals, Marin County, California, Earth Surf. Proc. Land., 16, 163–172,, 1991. 

Boudreau, B. P.: Mathematics of tracer mixing in sediments; I, Spatially-dependent, diffusive mixing, Am. J. Sci., 286, 161–198,, 1986. 

Boudreau, B. P.: The diffusion and telegraph equations in diagenetic modelling, Geochim. Cosmochim. Ac., 53, 1857–1866,, 1989. 

Braun, J., Mercier, J., Guillocheau, F., and Robin, C.: A simple model for regolith formation by chemical weathering, J. Geophys. Res.-Earth, 121, 2140–2171,, 2016. 

Brosens, L., Campforts, B., Robinet, J., Vanacker, V., Opfergelt, S., Ameijeiras-Mariño, Y., Minella, J. P. G., and Govers, G.: Slope Gradient Controls Soil Thickness and Chemical Weathering in Subtropical Brazil: Understanding Rates and Timescales of Regional Soilscape Evolution Through a Combination of Field Data and Modeling, J. Geophys. Res.-Earth, 125, 5,, 2020. 

Carretier, S., Goddéris, Y., Delannoy, T., and Rouby, D.: Mean bedrock-to-saprolite conversion and erosion rates during mountain growth and decline, Geomorphology, 209, 39–52,, 2014. 

Cerqueira, R.: The Distribution of Didelphis in South America (Polyprotodontia, Didelphidae), J. Biogeogr., 12, 135–145,, 1985. 

Chen, M., Ma, L., Shao, M.'a., Wei, X., Jia, Y., Sun, S., Zhang, Q., Li, T., Yang, X., and Gan, M.: Chinese zokor (Myospalax fontanierii) excavating activities lessen runoff but facilitate soil erosion – A simulation experiment, CATENA, 202, 105248,, 2021. 

Choi, K., Arnhold, S., Huwe, B., and Reineking, B.: Daily Based Morgan–Morgan–Finney (DMMF) Model: A Spatially Distributed Conceptual Soil Erosion Model to Simulate Complex Soil Surface Configurations, Water, 9, 278,, 2017. 

Cohen, S., Willgoose, G., and Hancock, G.: The mARM3D spatially distributed soil evolution model: Three-dimensional model framework and analysis of hillslope and landform responses, J. Geophys. Res., 115, 191,, 2010. 

Cohen, S., Willgoose, G., Svoray, T., Hancock, G., and Sela, S.: The effects of sediment transport, weathering, and aeolian mechanisms on soil evolution, J. Geophys. Res.-Earth, 120, 260–274,, 2015. 

Coombes, M. A.: Biogeomorphology: diverse, integrative and useful, Earth Surf. Proc. Land., 41, 2296–2300,, 2016. 

Corenblit, D., Corbara, B., and Steiger, J.: Biogeomorphological eco-evolutionary feedback between life and geomorphology: a theoretical framework using fossorial mammals, Naturwissenschaften, 108, 55, 2021. 

Debruyn, L. A. L. and Conacher, A. J.: The bioturbation activity of ants in agricultural and naturally vegetated habitats in semiarid environments, Soil Res., 32, 555–570,, 1994. 

Devia, G. K., Ganasri, B. P., and Dwarakish, G. S.: A Review on Hydrological Models, Aquat. Pr., 4, 1001–1007,, 2015. 

Durner, W., Iden, S. C., and von Unold, G.: The integral suspension pressure method (ISP) for precise particle-size analysis by gravitational sedimentation, Water Resour. Res., 53, 33–48,, 2017. 

Eccard, J. A. and Herde, A.: Seasonal variation in the behaviour of a short-lived rodent, BMC Ecol., 13, 43,, 2013. 

Ferro, L. I. and Barquez, R. M.: Species Richness of Nonvolant Small Mammals Along Elevational Gradients in Northwestern Argentina, Biotropica, 41, 759–767,, 2009. 

Foster, D. W.: BIOTURB: A FORTRAN program to simulate the effects of bioturbation on the vertical distribution of sediment, Comput. Geosci., 11, 39–54,, 1985. 

François, F., Poggiale, J.-C., Durbec, J.-P., and Stora, G.: A New Approach for the Modelling of Sediment Reworking Induced by a Macrobenthic Community, Acta Biotheor., 45, 295–319,, 1997. 

Gabet, E. J.: Gopher bioturbation: field evidence for non-linear hillslope diffusion, Earth Surf. Proc. Land., 25, 1419–1428, 2000. 

Gabet, E. J., Reichman, O. J., and Seabloom, E. W.: The Effects of Bioturbation on Soil Processes and Sediment Transport, Annu. Rev. Earth Pl. Sci., 31, 249–273,, 2003. 

Gabet, E. J., Perron, J. T., and Johnson, D. L.: Biotic origin for Mima mounds supported by numerical modeling, Geomorphology, 206, 58–66,, 2014. 

Goslee, S. C.: Topographic Corrections of Satellite Data for Regional Monitoring, Photogr. Remote Sens., 78, 973–981,, 2012. 

Gray, H. J., Keen-Zebert, A., Furbish, D. J., Tucker, G. E., and Mahan, S. A.: Depth-dependent soil mixing persists across climate zones, P. Natl. Acad. Sci. USA, 117, 8750–8756,, 2020. 

Grigusova, P.: Soil properties along Chilean climate gradient, [data set],,, last access: 3 August 2023a. 

Grigusova, P.: Modelled sediment redistribution along climate gradient, [data set],,;jsessionid=22F870744C71E3DAB58C6201A5026656?citid=521, last access: 3 August 2023b. 

Grigusova, P.: Model sediment redistribution caused by bioturbating animals, [code],, last access: 3 August 2023c. 

Grigusova, P., Larsen, A., Achilles, S., Klug, A., Fischer, R., Kraus, D., Übernickel, K., Paulino, L., Pliscoff, P., Brandl, R., Farwig, N., and Bendix, J.: Area-Wide Prediction of Vertebrate and Invertebrate Hole Density and Depth across a Climate Gradient in Chile Based on UAV and Machine Learning, Drones, 5, 86,, 2021. 

Grigusova, P., Larsen, A., Achilles, S., Brandl, R., del Río, C., Farwig, N., Kraus, D., Paulino, L., Pliscoff, P., Übernickel, K., and Bendix, J.: Higher sediment redistribution rates related to burrowing animals than previously assumed as revealed by time-of-flight-based monitoring, Earth Surf. Dynam., 10, 1273–1301,, 2022. 

Hakonson, T. E.: The Effects of Pocket Gopher Burrowing on Water Balance and Erosion from Landfill Covers, J. Environ. Qual., 28, 659–665,, 1999. 

Hall, K., Boelhouwers, J., and Driscoll, K.: Animals as Erosion Agents in the Alpine Zone: Some Data and Observations from Canada, Lesotho, and Tibet, Arct. Antarct. Alp. Res., 31, 436–446,, 1999. 

Hancock, G. and Lowry, J.: Quantifying the influence of rainfall, vegetation and animals on soil erosion and hillslope connectivity in the monsoonal tropics of northern Australia, Earth Surf. Proc. Land., 46, 2110–2123,, 2021. 

Hazelhoff, L., van Hoof, P., Imeson, A. C., and Kwaad, F. J. P. M.: The exposure of forest soil to erosion by earthworms, Earth Surf. Proc. Land., 6, 235–250,, 1981. 

Horn, B. K. P.: Hill shading and the reflectance map, Proc. IEEE, 69, 14–47,, 1981. 

Imeson, A. C. and Kwaad, F. J. P. M.: Some Effects of Burrowing Animals on Slope Processes in the Luxembourg Ardennes, Geografiska Annaler: Series A, Phys. Geogr., 58, 317–328,, 1976. 

Istanbulluoglu, E.: Vegetation-modulated landscape evolution: Effects of vegetation on landscape processes, drainage density, and topography, J. Geophys. Res., 110, 11,, 2005. 

Jimenez, J. E., Feinsinger, P., and Jaksi, F. M.: Spatiotemporal Patterns of an Irruption and Decline of Small Mammals in Northcentral Chile, J. Mammal., 73, 356–364,, 1992. 

Jong, S. M. de, Paracchini, M. L., Bertolo, F., Folving, S., Megier, J., and de Roo, A. P. J.: Regional assessment of soil erosion using the distributed model SEMMED and remotely sensed data, CATENA, 37, 291–308,, 1999. 

Jumars, P. A., Nowell, A. R. M., and Self, R. F. L.: A simple model of flow – Sediment – Organism interaction, Mar. Geol., 42, 155–172,, 1981. 

Kadko, D. and Heath, G. R.: Models of depth-dependent bioturbation at MANOP Site H in the eastern equatorial Pacific, J. Geophys. Res., 89, 6567,, 1984. 

Katzman, E. A., Zaytseva, E. A., Feoktistova, N. Y., Tovpinetz, N. N., Bogomolov, P. L., Potashnikova, E. V., and Surov, A. V.: Seasonal Changes in Burrowing of the Common Hamster (Cricetus cricetus L., 1758) (Rodentia: Cricetidae) in the City, Povolzhskiy Journal of Ecology, 17, 251–258,, 2018. 

Kinlaw, A. and Grasmueck, M.: Evidence for and geomorphologic consequences of a reptilian ecosystem engineer: The burrowing cascade initiated by the Gopher Tortoise, Geomorphology, 157/158, 108–121,, 2012. 

Kirols, H. S., Kevorkov, D., Uihlein, A., and Medraj, M.: The effect of initial surface roughness on water droplet erosion behaviour, Wear, 342/343, 198–209,, 2015. 

Kraus, D., Brandl, R., Achilles, S., Bendix, J., Grigusova, P., Larsen, A., Pliscoff, P., Übernickel, K., and Farwig, N.: Vegetation and vertebrate abundance as drivers of bioturbation patterns along a climate gradient, Plos One, 17, e0264408,, 2022. 

Kuhn, W.: Bioerosion auf Rhyolith im westlichen Mainzer Becken, Mainzer geowissenschaftliche Mitteilungen, 44, 63–72,, 2016. 

Kügler, M., Hoffmann, T. O., Beer, A. R., Übernickel, K., Ehlers, T. A., Scherler, D., and Eichel, J.: (LiDAR) 3D Point Clouds and Topographic Data from the Chilean Coastal Cordillera,, 2022. 

La Croix, A. D., Gingras, M. K., Dashtgard, S. E., and Pemberton, S. G.: Computer modeling bioturbation: The creation of porous and permeable fluid-flow pathways, Bulletin, 96, 545–556,, 2012. 

Larsen, A., Nardin, W., Lageweg, W. I., and Bätz, N.: Biogeomorphology, quo vadis? On processes, time, and space in biogeomorphology, Earth Surf. Proc. Land., 46, 12–23,, 2021. 

Le Hir, P., Monbet, Y., and Orvain, F.: Sediment erodability in sediment transport modelling: Can we account for biota effects?, Cont. Shelf Res., 27, 1116–1142,, 2007. 

Lehnert, L. W., Thies, B., Trachte, K., Achilles, S., Osses, P., Baumann, K., Schmidt, J., Samolov, E., Jung, P., Leinweber, P., Karsten, U., Büdel, B., and Bendix, J.: A Case Study on Fog/Low Stratus Occurrence at Las Lomitas, Atacama Desert (Chile) as a Water Source for Biological Soil Crusts, Aerosol Air Qual. Res., 18, 254–269,, 2018. 

Li, G., Li, X., Li, J., Chen, W., Zhu, H., Zhao, J., and Hu, X.: Influences of Plateau Zokor Burrowing on Soil Erosion and Nutrient Loss in Alpine Meadows in the Yellow River Source Zone of West China, Water, 11, 2258,, 2019. 

Li, T., Shao, M.'a., Jia, Y., Jia, X., and Huang, L.: Small-scale observation on the effects of the burrowing activities of mole crickets on soil erosion and hydrologic processes, Agriculture, Ecosyst. Environ., 261, 136–143,, 2018. 

Li, T. C., Shao, M. A., Jia, Y. H., Jia, X. X., Huang, L. M., and Gan, M.: Small-scale observation on the effects of burrowing activities of ants on soil hydraulic processes, Eur. J. Soil Sci., 70, 236–244,, 2019. 

Li, Z. and Zhang, J.: Calculation of Field Manning's Roughness Coefficient, Agr. Water Manag., 49, 153–161,, 2001. 

Lilhare, R., Garg, V., and Nikam, B. R.: Application of GIS-Coupled Modified MMF Model to Estimate Sediment Yield on a Watershed Scale, J. Hydrol. Eng., 20, 1443–1459,, 2015. 

López-Vicente, M., Navas, A., and Machín, J.: Modelling soil detachment rates in rainfed agrosystems in the south-central Pyrenees, Agr. Water Manag., 95, 1079–1089,, 2008. 

Malizia, A. I.: Population dynamics of the fossorial rodent Ctenomys talarum (Rodentia: Octodontidae), J. Zool., 244, 545–551,, 1998. 

Meserve, P. L.: Trophic Relationships among Small Mammals in a Chilean Semiarid Thorn Scrub Community, J. Mammal., 62, 304–314,, 1981. 

Meyer, H., Reudenbach, C., Hengl, T., Katurji, M., and Nauss, T.: Improving performance of spatio-temporal machine learning models using forward feature selection and target-oriented validation, Environ. Model. Softw., 101, 1–9,, 2018. 

Meysman, F. J. R., Boudreau, B. P., and Middelburg, J. J.: Relations between local, nonlocal, discrete and continuous models of bioturbation, J. Mar. Res., 61, 391–410,, 2003. 

Meysman, F. J. R., Boudreau, B. P., und Middelburg, J. J.: Modeling reactive transport in sediments subject to bioturbation and compaction, Geochim. Cosmochim. Ac., 69, 3601–3617,, 2005. 

Milstead, W. B., Meserve, P. L., Campanella, A., Previtali, M. A., Kelt, D. A., and Gutiérrez, J. R.: Spatial Ecology of Small Mammals in North-central Chile: Role of Precipitation and Refuges, J. Mammal., 88, 1532–1538,, 2007. 

Monteverde, M. J. and Piudo, L.: Activity Patterns of the Culpeo Fox (Lycalopex Culpaeus Magellanica ) in a Non-Hunting Area of Northwestern Patagonia, Argentina, Mamm. Study, 36, 119–125,, 2011. 

Morgan, R. P. C. and Duzant, J. H.: Modified MMF (Morgan–Morgan–Finney) model for evaluating effects of crops and vegetation cover on soil erosion, Earth Surf. Proc. Land., 33, 90–106,, 2008. 

Morgan, R. P. C.: A simple approach to soil loss prediction: a revised Morgan–Morgan–Finney model, CATENA, 44, 305–322,, 2001. 

Morgan, R. P. C., Morgan, D. D. V., and Finney, H. J.: A predictive model for the assessment of soil erosion risk, J. Agr. Eng. Res., 30, 245–253,, 1984. 

Morgan, R. P. C., Quinton, J. N., Smith, R. E., Govers, G., Poesen, J. W. A., Auerswald, K., Chisci, G., Torri, D., and Styczen, M. E.: The European Soil Erosion Model (EUROSEM): a dynamic approach for predicting sediment transport from fields and small catchments, Earth Surf. Proc. Land., 23, 527–544,<527:AID-ESP868>3.0.CO;2-5, 1998. 

Morris, J. E., Hampson, G. J., und Johnson, H. D.: A sequence stratigraphic model for an intensely bioturbated shallow-marine sandstone: the Bridport Sand Formation, Wessex Basin, UK, Sedimentology, 53, 1229–1263,, 2006. 

Nearing, M. A., Foster, G. R., Lane, L. J., and Finkner, S. C.: A Process-Based Soil Erosion Model for USDA-Water Erosion Prediction Project Technology, T. ASAE, 32, 1587–1593,, 1989. 

Nkem, J. N., Lobry de Bruyn, L. A., Grant, C. D., and Hulugalle, N. R.: The impact of ant bioturbation and foraging activities on adjacent soil properties, Pedobiologia, 44, 609–621,, 2000. 

Oeser, R. A., Stroncik, N., Moskwa, L.-M., Bernhard, N., Schaller, M., Canessa, R., van den Brink, L., Köster, M., Brucker, E., Stock, S., Fuentes, J. P., Godoy, R., Matus, F. J., Oses Pedraza, R., Osses McIntyre, P., Paulino, L., Seguel, O., Bader, M. Y., Boy, J., Dippold, M. A., Ehlers, T. A., Kühn, P., Kuzyakov, Y., Leinweber, P., Scholten, T., Spielvogel, S., Spohn, M., Übernickel, K., Tielbörger, K., Wagner, D., and Blanckenburg, F. von: Chemistry and microbiology of the Critical Zone along a steep climate and vegetation gradient in the Chilean Coastal Cordillera, CATENA, 170, 183–203,, 2018. 

Orvain, F.: A model of sediment transport under the influence of surface bioturbation: generalisation to the facultative suspension-feeder Scrobicularia plana, Mar. Ecol. Prog. Ser., 286, 43–56,, 2005. 

Orvain, F., Le Hir, P., und Sauriau, P.-G.: A model of fluff layer erosion and subsequent bed erosion in the presence of the bioturbator, Hydrobia ulvae, J. Mar. Res., 61, 821–849,, 2003. 

Orvain, F., Sauriau, P.-G., Bacher, C., and Prineau, M.: The influence of sediment cohesiveness on bioturbation effects due to Hydrobia ulvae on the initial erosion of intertidal sediments: A study combining flume and model approaches, J. Sea Res., 55, 54–73,, 2006. 

Pelletier, J. D., Barron-Gafford, G. A., Breshears, D. D., Brooks, P. D., Chorover, J., Durcik, M., Harman, C. J., Huxman, T. E., Lohse, K. A., Lybrand, R., Meixner, T., McIntosh, J. C., Papuga, S. A., Rasmussen, C., Schaap, M., Swetnam, T. L., and Troch, P. A.: Coevolution of nonlinear trends in vegetation, soils, and topography with elevation and slope aspect: A case study in the sky islands of southern Arizona, J. Geophys. Res.-Earth, 118, 741–758,, 2013. 

Penman, H.: Natural evaporation from open water, hare soil and grass, Proc. Roy. Soc. Lond. Ser. A, 193, 120–145,, 1948. 

Pollacco, J. A. P.: A generally applicable pedotransfer function that estimates field capacity and permanent wilting point from soil texture and bulk density, Can. J. Soil. Sci., 88, 761–774,, 2008. 

Qin, Y., Yi, S., Ding, Y., Qin, Y., Zhang, W., Sun, Y., Hou, X., Yu, H., Meng, B., Zhang, H., Chen, J., and Wang, Z.: Effects of plateau pikas' foraging and burrowing activities on vegetation biomass and soil organic carbon of alpine grasslands, Plant Soil, 458, 201–216,, 2021. 

Rakotomalala, C., Grangeré, K., Ubertini, M., Forêt, M. und Orvain, F.: Modelling the effect of Cerastoderma edule bioturbation on microphytobenthos resuspension towards the planktonic food web of estuarine ecosystem, Ecol. Modell., 316, 155–167,, 2015. 

Reichman, O. J. and Seabloom, E. W.: The role of pocket gophers as subterranean ecosystem engineers, Trend. Ecol. Evol., 17, 44–49,, 2002. 

Renard, K., Foster, G., Weesies, G., and Porter, J.: RUSLE: The Revised Universal Soil Loss Equation, J. Soil Water Conserv., 46, 30–33, 1991. 

Ridd, P. V.: Flow Through Animal Burrows in Mangrove Creeks, Estuar. Coast. Shelf Sci., 43, 617–625,, 1996. 

Rodríguez-Caballero, E., Cantón, Y., Chamizo, S., Afana, A., and Solé-Benet, A.: Effects of biological soil crusts on surface roughness and implications for runoff and erosion, Geomorphology, 145/146, 81–89,, 2012. 

Román-Sánchez, A., Reimann, T., Wallinga, J., and Vanwalleghem, T.: Bioturbation and erosion rates along the soil-hillslope conveyor belt, part 1: Insights from single-grain feldspar luminescence, Earth Surf. Proc. Land., 44, 2051–2065,, 2019. 

Roo, A. P. J. De, Wesseling, C. G., and Ritsema, C. J.: lisem: a single-event physically based hydrological and soil erosion model for drainage basins, I: theory, input and output, Hydrol. Process., 10, 1107–1117,<1107:AID-HYP415>3.0.CO;2-4, 1996. 

Rutin, J.: The burrowing activity of scorpions (Scorpio maurus palmatus) and their potential contribution to the erosion of Hamra soils in Karkur, central Israel, Geomorphology, 15, 159–168,, 1996. 

Sanford, L. P.: Modeling a dynamically varying mixed sediment bed with erosion, deposition, bioturbation, consolidation, and armoring, Comput. Geosci., 34, 1263–1283,, 2008. 

Schiffers, K., Teal, L. R., Travis, J. M. J., and Solan, M.: An open source simulation model for soil and sediment bioturbation, Plos One, 6, e28028,, 2011. 

Shannon, C. E.: A Mathematical Theory of Communication, Bell Syst. Tech. J., 27, 379–423,, 1948. 

Shull, D. H.: Transition-matrix model of bioturbation and radionuclide diagenesis, Limnol. Oceanogr., 46, 905–916,, 2001. 

Simonetti, J. A.: Microhabitat Use by Small Mammals in Central Chile, Oikos, 56, 309–318,, 1989. 

Soetaert, K., Herman, P. M. J., Middelburg, J. J., Heip, C., deStigter, H. S., van Weering, T. C. E., Epping, E., and Helder, W.: Modeling 210Pb-derived mixing activity in ocean margin sediments: Diffusive versus nonlocal mixing, J. Mar. Res., 54, 1207–1227,, 1996. 

Taylor, A. R., Lenoir, L., Vegerfors, B., and Persson, T.: Ant and Earthworm Bioturbation in Cold-Temperate Ecosystems, Ecosystems, 22, 981–994,, 2019. 

Temme, A. J. A. M. and Vanwalleghem, T.: LORICA – A new model for linking landscape and soil profile evolution: Development and sensitivity analysis, Comput. Geosci., 90, 131–143,, 2016. 

Tews, J., Brose, U., Grimm, V., Tielbörger, K., Wichmann, M. C., Schwager, M., and Jeltsch, F.: Animal species diversity driven by habitat heterogeneity/diversity: the importance of keystone structures, J. Biogeogr., 31, 79–92,, 2004. 

Tomasella, J., Hodnett, M. G., and Rossato, L.: Pedotransfer Functions for the Estimation of Soil Water Retention in Brazilian Soils, Soil Sci. Soc. Am. J., 64, 327–338,, 2000. 

Trauth, M. H.: TURBO: a dynamic-probabilistic simulation to study the effects of bioturbation on paleoceanographic time series, Comput. Geosci., 24, 433–441,, 1998. 

Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50,, 2010. 

Übernickel, K., Pizarro-Araya, J., Bhagavathula, S., Paulino, L., and Ehlers, T. A.: Reviews and syntheses: Composition and characteristics of burrowing animals along a climate and ecological gradient, Chile, Biogeosciences, 18, 5573–5594,, 2021a. 

Übernickel, K., Ehlers, T. A., Paulino, L., and Fuentes Espoz, J.-P.: Time series of meteorological stations on an elevational gradient in National Park La Campana, Chile, GFZ Data Services,, 2021b. 

Vanwalleghem, T., Stockmann, U., Minasny, B., and McBratney, A. B.: A quantitative model for integrating landscape evolution and soil formation, J. Geophys. Res.-Earth, 118, 331–347,, 2013. 

Vieira, D. C. S., Prats, S. A., Nunes, J. P., Shakesby, R. A., Coelho, C. O. A., and Keizer, J. J.: Modelling runoff and erosion, and their mitigation, in burned Portuguese forest using the revised Morgan–Morgan–Finney model, Forest Ecol. Manag., 314, 150–165,, 2014. 

Vigiak, O., Okoba, B. O., Sterk, G., and Groenenberg, S.: Modelling catchment-scale erosion patterns in the East African Highlands, Earth Surf. Proc. Land., 30, 183–196,, 2005. 

Voiculescu, M., Ianăş, A.-N., and Germain, D.: Exploring the impact of snow vole (Chionomys nivalis) burrowing activity in the Făgăra? Mountains, Southern Carpathians (Romania): Geomorphic characteristics and sediment budget, CATENA, 181, 104070,, 2019. 

Wang, B., Zheng, F., Römkens, M. J.M., and Darboux, F.: Soil erodibility for water erosion: A perspective and Chinese experiences, Geomorphology, 187, 1–10,, 2013. 

Wei, X., Li, S., Yang, P., and Cheng, H.: Soil erosion and vegetation succession in alpine Kobresia steppe meadow caused by plateau pika – A case study of Nagqu County, Tibet, Chin. Geograph. Sc., 17, 75–81,, 2007. 

Welivitiya, W. D. D. P., Willgoose, G. R., and Hancock, G. R.: A coupled soilscape–landform evolution model: model formulation and initial results, Earth Surf. Dynam., 7, 591–607,, 2019. 

Wheatcroft, R. A., Jumars, P. A., Smith, C. R., and Nowell, A. R. M.: A mechanistic view of the particulate biodiffusion coefficient: Step lengths, rest periods and transport directions, J. Mar. Res., 48, 177–207,, 1990. 

Whitesides, C. J. and Butler, D. R.: Bioturbation by gophers and marmots and its effects on conifer germination, Earth Surf. Proc. Land., 41, 2269–2281,, 2016. 

Wilkinson, M. T., Richards, P. J., and Humphreys, G. S.: Breaking ground: Pedological, geological, and ecological implications of soil bioturbation, Earth-Sci. Rev., 97, 257–272,, 2009. 

Williams, J. R. (Ed.): Sediment-yield prediction with Universal Equation using runoff energy factor. In Present and prospective technology for predicting sediment yield and sources: Proceedings of the Sediment-Yield Workshop, ARS-S-40, United States Department of Agriculture (USDA), New Orleans, USA, 1975. 

Wilson, M. F. J., O'Connell, B., Brown, C., Guinan, J. C., and Grehan, A. J.: Multiscale Terrain Analysis of Multibeam Bathymetry Data for Habitat Mapping on the Continental Slope, Mar. Geod., 30, 3–35,, 2007. 

Wischmeier, W. and Smith, D. D.: Predicting rainfall erosion losses – A guide to conservation planning, Agriculture Handbook, US Department of Agriculture, 1–58, 1978.  

Wood, S. N.: Generalized Additive Models, Chapman and Hall/CRC, ISBN 9781315370279, 2006. 

Wösten, J. H. M. (Ed.): Soil Quality for Crop Production and Ecosystem Health, Developments in Soil Science, Elsevier, ISBN: 9780080541402, 1997. 

Wu, C., Wu, H., Liu, D., Han, G., Zhao, P., and Kang, Y.: Crab bioturbation significantly alters sediment microbial composition and function in an intertidal marsh, Estuar. Coast. Shelf Sci., 249, 107116,, 2021. 

Yair, A.: Short and long term effects of bioturbation on soil erosion, water resources and soil development in an arid environment, Geomorphology, 13, 87–99,, 1995. 

Yoo, K. and Mudd, S. M.: Toward process-based modeling of geochemical soil formation across diverse landforms: A new mathematical framework, Geoderma, 146, 248–260,, 2008. 

Yoo, K., Amundson, R., Heimsath, A. M., and Dietrich, W. E.: Process-based model linking pocket gopher (Thomomys bottae) activity to sediment transport and soil thickness, J. Geophys. Res., 33, 917,, 2005. 

Yu, C., Zhang, J., Pang, X. P., Wang, Q., Zhou, Y. P., and Guo, Z. G.: Soil disturbance and disturbance intensity: Response of soil nutrient concentrations of alpine meadow to plateau pika bioturbation in the Qinghai-Tibetan Plateau, China, Geoderma, 307, 98–106,, 2017. 

Zevenbergen, L. W. and Thorne, C. R.: Quantitative analysis of land surface topography, Earth Surf. Proc. Land., 12, 47–56,, 1987. 

Zhang, Q., Li, J., Hu, G., and Zhang, Z.: Bioturbation potential of a macrofaunal community in Bohai Bay, northern China, Mar. Pollut. Bull., 140, 281–286,, 2019. 

Zhang, S., Fang, X., Zhang, J., Yin, F., Zhang, H., Wu, L., and Kitazawa, D.: The Effect of Bioturbation Activity of the Ark Clam Scapharca subcrenata on the Fluxes of Nutrient Exchange at the Sediment-Water Interface, J. Ocean Univ. China, 19, 232–240,, 2020. 

Short summary
In our study, we included bioturbation into a soil erosion model and ran the model for several years under two conditions: with and without bioturbation. We validated the model using several sediment fences in the field. We estimated the modeled sediment redistribution and surface runoff and the impact of bioturbation on these along a climate gradient. Lastly, we identified environmental parameters determining the positive or negative impact of bioturbation on sediment redistribution.
Final-revised paper