The potential risk of combined effects of water and tillage erosion on the agricultural landscape in Czechia

ABSTRACT Tillage erosion is considered a major contributor to total soil erosion besides water erosion. However, tillage erosion has been neglected in Czechia. Therefore, our main goal was to analyse and compare the effect of tillage erosion with that of water erosion. The combined effect of both types of soil erosion was modelled at high spatial resolution. The Universal Soil Loss Equation approach was used to calculate water erosion, while a tillage model (diffusion approach), which takes into account the slope curvature and the tillage transport coefficient, was used to calculate tillage erosion. The constructed map showed that 48% of agricultural soils is threatened by total soil erosion at high rates. Furthermore, the area threatened by tillage erosion is almost 1.5 times larger than that endangered only by water erosion, and the mean contribution of tillage erosion to total soil erosion ranges between 20% and 30%.


Introduction
Soil erosion is generally considered to be the most problematic, widespread, and visible form of soil degradation in Europe (Boardman & Poesen, 2006;EEA and JRC, 2010;Grimm et al., 2002;Panagos et al., 2016;Stolte et al., 2016) and globally (FAO, 2011;FAO and ITPS, 2015). Erosion processes negatively influence a variety of soil properties, including changes in the thickness of the humic horizon and the soil organic carbon content and stock, which have major impacts on soil productivity. Furthermore, soil erosion has crucial effects on biogeochemical cycles due to intensive lateral fluxes (Van Oost et al., 2007;Quinton et al., 2010;Doetterl et al., 2016).
Traditionally, soil erosion has been related to water and wind processes that redistribute soil particles within the landscape. A potential negative effect of erosion and redistribution by tillage was rediscovered in the 1990s (Lindstrom et al., 1990;Govers et al., 1993;Lobb et al., 1995). Tillage erosion refers to the continuous downhill movement of soil particles caused by tillage operations. However, unlike water erosion, whose effect in a field is well recognised and the rates at which it can occur are widely known and available, tillage erosion remains poorly documented to date. Experimental studies conducted on individual fields or in small micro-catchments have shown the net spatial effect of the entire process of tillage erosion, a process controlled by the slope gradient (Lindstrom et al., 1990;1992), whereby its effects are predominant in hilly landscapes. The eroded parts are located on the convex landscape positions, whereas deposition takes place along the concave landscape positions (e.g. Govers et al., 1994;Hrabalíková et al., 2016;Lobb et al., 1995;Van Muysen et al., 2000). A strong difference between tillage and water erosion deals mainly with the spatial effects on the landscape in each case. The influence of water erosion is more visible on the sloping parts of the field, whereas evidence of tillage erosion is more visible on summits and shoulder slopes.
The rate of tillage erosion has dramatically increased over the last 70 years due to the widespread introduction of intensive, mechanised agriculture (Li et al., 2007). Van Oost et al. (2009) for example found that 45%-50% of the cropland area in Europe is threatened by tillage erosion. Many experimental and modelling studies have clearly shown that the rates of tillage erosion reach the same values as those for water erosion (Li et al., 2007;Lobb et al., 1995;Öttl et al., 2021;Wilken et al., 2017;2020). Thus, the advent of mechanised agriculture has caused tillage erosion to become an important component of total soil erosion, together with water and wind erosion. Despite the known facts verified in local studies, assessment of the influence of tillage erosion is frequently still a neglected process, especially in Czechia. Spatial information on the rate of water erosion is available mainly in the form of potential erosion risk, often calculated using empirical models of potential soil loss rates such as Universal Soil Loss Equation (USLE) or Revised USLE Panagos et al., 2015). Model calculations have also been used to estimate tillage erosion. Govers et al. (1994) introduced a mathematical formulaon which the most erosion models are builtfor simulating erosion due to tillage. During the last decades, a few models have been proposed to estimate erosion by tillage and have been successfully applied (e.g. De Alba, 2003;Li et al., 2007;Lindstrom et al., 2000;Vieira & Dabney, 2011).
Information about the rate of tillage erosion and the overall impact of erosion, together with water erosion, is available at a detailed level, either locally or on a gross spatial scale. However, detailed information on the spatial variability of these phenomena is very important for designing anti-erosion and decisionmaking strategies for soil protection. The main objectives of this study and of the constructed map were: i) to model the combined effects of water and tillage erosion to identify the most threatened areas in an agricultural landscape in Czechia at high spatial resolution, and ii) to analyse the importance of tillage erosion in relation to water erosion. One of the most common primary tillage implements is the mouldboard, whose operation settings were chosen to run the model. The main map is based on calculations using the concept proposed by Govers et al. (1994) and modified by Van Oost et al. (2009), which is the most widely used approach in tillage erosion modelling. We did not include the effect of wind erosion, which is less obvious because of the extensive spatial impact involved, whereby it is also very difficult to determine the rate of soil loss by wind erosion and its contribution to the total soil erosion rate. The final version of the map constructed was also provided as an official map certified by the Ministry of Agriculture of the Czech Republic as a suitable basis for the design and implementation of specific agricultural management strategies for a more effective soil protection in the most threatened areas at a national level.

Study area
The agricultural landscape of Czechia covers 4,16 mil. ha, approximately 42% of the total national territory. Most of this area is characterised by a rolling topography with mainly a hilly, highland character. Diverse geological structure resulted in the heterogeneous soil cover in Czechia. The landscape is predominantly covered by quaternary sediments (Chlupáč et al., 2002). 54% of the area is represented by colluvial material with Cambisols as a dominant soil unit. Lowlands and middle positions are covered by other significant aeolian sediments with main soil units Luvisols and Chernozems. Eluvial sediments occurring in river floodplains and terraces with the most presence of Fluvisols and Regosols Figure 1.
The rolling topogprahy together with the presence of soil units that are vulnerable to soil erosion (silty soils such as Haplic Luvisols, Albic Luvisols, and Chernozems), there is a potentially strong impact of water and tillage erosion ( Figure 2). Czechia is affected by a predominantly westerly atmospheric circulation with a temperate continental climate. The annual mean precipitation fluctuates within the 400-1500 mm range with peak rain in the summer months and in the mountainous area.

Modelling tillage erosion
The process of tillage erosion was calculated following the diffusion approach of Govers et al. (1994), in which the rate of tillage erosion is expressed by a single parameter, the tillage transport coefficient (k till ). The concept is based on the relationship between the mean displacement distance of soil particles and the slope gradient derived from field measurements in the early 1990s (Lindstrom et al., 1992;1990). Mathematical relationships for tillage erosion have been described previously (e.g. Govers et al., 1999;Govers et al., 1994;Li et al., 2009;Van Oost et al., 2006;Van Oost et al., 2000b), and a detailed review of the controlling factors involved is available . Tillage erosion or accumulation rate is described by the proportional factor k till and the rate of change of slope in the direction of operation. The relationship for the rate of soil mass translocation (Q s in kg·m −2 ) is expressed as follows: where D is the tillage depth (m), ρ b is the soil bulk density (kg·m −3 ), b is the empirically derived regression coefficient, and k till is the tillage transport coefficient (kg·m −1 ), h is the local elevation of the given point, and x is the direction of the movement. Total erosion or accumulation rate of the soil mass is then described using the one-dimensional mass conservation equation: where E is the rate of soil erosion or accumulation (kg·m −2 ·y −1 ), t is time . This means that the soil erosion rate can be described using the proportional factor k till and the rate of change of slope in the direction of operation which is using the slope curvature. It also follows from equation (3) that, while the soil loss from a particular point (Q s ) depends on the slope, the total mass balance (E) depends on the slope change. Thus, there is soil loss on convex slopes and accumulation on concave slopes. When some simplifying assumptions are applied, namely, the assumption of displacement of all soil particles in the direction of greatest slope, a linear relationship between translocation length and slope, tillage in the opposite direction in two subsequent operations, and spatially invariant ploughing depth and bulk density, equation (3) describes the one-dimensional movement of soil mass and can be subsequently used to simulate erosion by tillage in a two-dimensional space.
The tillage transport coefficient, a fundamental parameter, is difficult to obtain due to the time consuming and costly, extensive field measurements required. Very few empirical studies have been reported in Czechia (Hrabalíková et al., 2016;Novák & Hůla, 2018). Therefore, herein, we used the non-linear relationship developed by Van Oost et al. (2006), based on a regression analysis of available data of k till . This regression relationship is expressed as follows: where ρ b is the soil bulk density (kg·m −3 ), D is the soil tillage depth (m), V is the velocity of operation (km·h −1 ), and T represents the direction of operation (a value of 1 is assigned to contour tillage and a value of 2 is assigned to up and down slope tillage). The coefficients a, α, β, and γ are the regression coefficients for mouldboard tillage (the corresponding values are 0.97, 2.21, 0.57, and 0.67, respectively). Operation settings followed the most widely used values for mouldboard tillage. Soil depth and velocity of operation were 0,3 m and 8 km·h −1 (2.22 m·s −1 ), respectively. The simulation was run for a single annual mouldboard operation, considering the up and downslope directions of tillage. Secondary tillage operations, such as harrowing or discing and counter tillage were not considered. Owing to the lack of accurate data for the overall soil bulk density of agricultural soils, the values were derived based on the results of National basal soil monitoring (Kňákal, 2000). The monitoring showed overall topsoil average bulk density 1,39 g.cm −3 . Bulk density was evaluated for texture classes ranging from 1.17 g.cm −3 for clayey class to 1.45 g.cm −3 for clay loam (see Table 1). A 1:50,000 scale soil texture map based on data from a Systematic soil survey of agricultural soils in Czechoslovakia  was used as input for this conversion. The resulting k till values (Figure 3a) were confronted with empirical derived values on experimental plots (Hrabalíková et al., 2016;Novák & Hůla, 2018) and did not show any significant bias.
The second main input to the model of tillage erosion is the curvature of the slope, which is represented by the second derivative of the elevation of the terrain in Equation 3. The total curvature of the terrain is considered, which does not take into account the direction of tillage (data on tillage directions are not available nationwide). Curvature was calculated using the Curvature tool in ArcGIS software (Moore et al., 1991;Zevenbergen & Thorne, 1987). The resulting values ranged from −1 to 1, with negative values characterising convex slopes and positive values characterising concave slopes (Figure 3b). Elevation data were obtained from the Digital Terrain Model of the Czech Republic of the 4th generation (DMR 4G) with 5 m horizontal resolution and a total standard error of 0.3 m height in bare ground.

Modelling water erosion
Water erosion was calculated on the basis of the USLE (Wischmeier & Smith, 1965, 1978. This method is the most widely used to quantify the average annual soil loss (Mg·ha −1 ·year −1 ), and is based on the wellknown empirical equation: where R is the rainfall erosivity factor, defined as the product of rainfall kinetic energy and its highest 30-min intensity, K is the soil erodibility factor, expressing soil susceptibility to erosion, and L is the length factor, which expresses the influence of continual slope length on the extent of soil loss; S is the slope factor expressing the influence of slope gradient on the extent of soil loss, C is the cover-management factor expressing the influence of sowing methods and agro-technology, and P is the support-practice factor, which expresses the influence of erosion control measures. A more detailed description of the calculation of individual factors can be found in Novotný et al. (2016). In contrast to the aforementioned article, updates of the combined LS factor used in the DMR 4G elevation model and the regionalised R factor (Rožnovský et al., 2015) were included here. Values for C were determined according to climatic regions (Kadlec & Toman, 2002). Values of 0.005 and 0.8 were assigned to areas of permanent grassland and Hop gardens, respectively, while 0.45 was assigned to both orchards and vineyards.

Final map composition
Raster of potential annual soil loss by both tillage and water erosion were calculated with a 5 m spatial resolution (Main map). As the study focussed on potential annual soil loss, the entire agricultural landscape was considered. For permanent grassland, the annual soil loss was estimated at 0 Mg·ha −1 ·y −1 , as the area is not affected by tillage operations at this moment. For map readability at the final scale of 1:500,000 printed in A0 format, rasters were generalised on a mean aggregation strategy for individual fields from the database of the Land Parcel Identification System (LPIS -© Ministry of Agriculture of the Czech Republic). As soil accumulation values are not available for water erosion, only the rates of soil loss were considered in the assessment of tillage  The final map was prepared as a bivariate choropleth map that illustrates the combined effect of water and tillage erosion on agricultural land in Czechia and shows the places most prone to both the individual erosion factor and the combined erosion factors. Each raster was classified into four individual classes with equal distribution of pixel values based on histogram statistics (Figure 4). Analyses were implemented in an R environment (R Core Team, 2020). A classification based on interval limits into the combined legend was made in a raster calculator in QGIS 3.10.2. The colour scheme of the combined legend was based on a combination of two complementary colour schemes in which the values are represented by different colour saturation depending on the magnitude of soil loss. The combined legend was prepared in GIMP 2.10.8 by multiplying individual colour schemes in the 2D space. Subsequently, the colour for each class was read and applied to the individual pixels of the class in GIS.

Results and discussion
The calculated distribution of water and tillage erosion is presented in Figure 5. Markedly high rates (>4 Mg·ha −1 ·y −1 ) of tillage erosion were determined for 12% of the agricultural landscape in Czechia. However, when the calculated area was considered only for arable land, which is affected by tillage operations, the most endangered area reached 16%. This represents landscapes that are closely related to the rolling topography regions. The regions most affected by tillage erosion are located in south and northeast Moravia, and belong to the most fertile arable land in Czechia (Žížala et al., 2019). In contrast, very low erosion rates (<1 Mg·ha −1 ·y −1 ) were calculated for the major lowland regions with low relief intensity. The same applied to water erosion, consistently with previous findings reported for water erosion by .
The overall average tillage erosion rate for agricultural land was estimated at 1.5 and 2.4 Mg·ha −1 ·y −1 for arable land alone. As tillage erosion typically affected approximately half of the entire arable land, while the other half of the arable land was affected by deposition, average tillage erosion might correspond to a soil loss rate of 3.0 Mg ha −1 y −1 (4.8 Mg·ha −1 ·y −1 ) on the eroding parts of the agricultural landscape. Van Oost et al. (2009) estimated a mean tillage erosion rate of 4.1 Mg·ha −1 ·y −1 for Czechia. This discrepancy with our study might be mainly caused by the calculation of tillage erosion for permanent grasslands, which Van Oost et al. (2009) considered arable land where tillage operations are performed. Other potential effects on the results obtained in each study might derive from the input data. Van Oost et al. (2009) used the SRTM model with a resolution of 90 m, modelled in a 1 km grid with constant k till (500 kg·m −1 ·y −1 ) and used the CORINE land cover database to select land use classes. On the other hand, we used local and more accurate data; furthermore, our calculation was performed in a 5 m grid.
In terms of water erosion, based on the results of Cerdan et al. (2006), Van Oost calculated an average soil loss of 5.7 Mg·ha −1 ·y −1 in Czechia. In turn, Panagos et al. (2015) calculated an average soil loss rate due water erosion of 1.65 Mg·ha −1 ·y −1 for all lands and of 2.52 Mg·ha −1 ·y −1 for arable land in Czechia. Interestingly, according to our own calculations using detailed topography and land-use data incorporating information on parcels borders and objects interrupting the outflow, the average value was 4.8 Mg·ha −1 ·y −1 (7.1 Mg·ha −1 ·y −1 for arable land alone), i.e. more than double the value presented by Panagos et al. (2015). Because all mentioned values are based on model outputs, they only consider interrill and rill erosion processes but not gullies or rills in the concentrated drain paths. However, we believe that our estimates are more legitimate due to the use of local data, i.e. accurate topography data (5 m DEM versus 25 m DEM, along with the use of the local land-use information versus the use of the CORINE database), and the use the regionalised precipitation erosivity factor according to local measurements or the use of distributed values of the vegetation protection factor, which altogether logically lead to more valid results. Despite these improvements, however, based on experience of erosion monitoring (Žížala et al., 2016) and simulations at the local level, we know that model values can often be underestimated. This finding is in line with other authors criticising the approach of Panagos et al. (2015), e.g. Boardman (2016a, 2016b), Auerswald (2016), andFiener et al. (2020).
Total soil erosion rate was calculated as the sum of water and tillage erosion rates. Total soil erosion reached markedly high rates on 48% of the agricultural landscape (65% when only arable land was considered). The most endangered area is the Ždánický forest region in south Moravia. These results suggest that the spatial extent of the most threatened area is almost 1.5 times higher than the area vulnerable to the single effect of water erosion. Similar results were found by Van Oost et al. (2009) on the European scale. Tillage erosion rates reach similar values as water erosion rates on the sloping parts of the landscape, a finding that has already been demonstrated in many experimental and modelling studies. For the European agricultural landscape, Van Oost et al. (2009) estimated an average tillage erosion rate only 15% lower than the water erosion rate. Based on calculations at a resolution of 1 km, these authors found a tillage/water erosion ratio of 0.72 for Czechia. On the other hand, according to the study reported herein, the mean relative contribution of tillage  erosion to total soil erosion (sum of water and tillage erosion) was 28% ( Figure 6) and the resulting ratio of tillage/water erosion was 0.61. However, these values can only be considered as an indication of the method of calculation of both layers. Furthermore, it should be noted that the highest contribution of tillage erosion (more than 50%) was investigated in relatively flat regions with very low rates of water erosion and, concomitantly, slightly higher values of tillage erosion. Figure 7 shows rates of water and tillage erosion according to the spatial distribution of soil types. The distribution of soil erosion rates highly depends on the topographic position in which the individual soil types are located. It also shows dependency with the soil erodibility which is related topsoil characteristics such as grain size or humus content. High values are thus associated mainly with the occurrence of soil types determined by their location on slopes (Leptosols or erosion-related Regosols) or occurring on fine-grained substrates, e.g. on loess (Chernozems, Kastanozems, Luvisols).
Our results suggest that the impact of tillage erosion on the overall soil redistribution is very high. Nevertheless, the calculation of tillage erosion has a few limitations, and it is important to note several simplifications. Tillage erosion is primarily controlled by the erosivity of a specific tillage operation and landscape erodibility . As summarised by Van Oost et al. (2006), the most important factors influencing landscape erodibility include the settings of tillage implements (depth, speed, and direction of tillage), implement type (i. e., tillage erosivity and Figure 6. Spatial contribution of tillage erosion to total soil erosion for the agricultural landscape in Czechia.  soil properties (bulk density, texture, structure), and field topography (slope gradient and curvature). The model used in this study is invariable with respect to operation settings, and the final map was calculated only for the most widely used implement in Czechia, namely, the mouldboard with a single operation setting. Additionally, only one direction of tillage operation was considered, that is, down-and up-slope. Unfortunately, at present, suitable empirical data on the agricultural techniques used and directions of tillage operations are not available nationwide, whereby, we could not include secondary tillage operations either (e.g. harrowing) that would adequately increase the tillage erosion rate. Additionally, variation in the value for tillage transport coefficient for multiple tillage implements and operation settings has been demonstrated in many experimental studies (summarised in Van Oost et al., 2006). In our calculations, we used spatially distributed k till values that reflect differences in soil bulk density. The mean value of k till was estimated at 488 kg·m −1 ·y −1 , which is higher than in most previous experimental studies for mouldboard tillage (e.g. Lindstrom et al., 1992;Lobb et al., 1995Quine et al., 1994;Van Muysen et al., 2002). Similarly, Fiener et al. (2018) pointed the uncertainty involved in comparing different methods to determine the value for the tillage transport coefficient. They found a relatively wide range of soil translocation rates. Thus, a comparison of the mean soil translocation rate calculated by six different techniques revealed a considerable range among the magnitude of underestimation (−32.8 kg·m −1 ) and that of the overestimation (41.3 kg·m −1 ) for different techniques at the same slope position. Finally, tillage erosion values can be further affected by disregarding the effect of plot edges.

Conclusions
The final main map of the potential risk of the combined effects of water and tillage erosion showed the area of potential soil loss for individual fields. The combined visualisation showed the places where the potential risk is greatest, where a given type of erosion prevails, and where overall risk is minimal. Our results suggest that the impact of tillage erosion on overall soil redistribution is significantly high. We found that 12% of the agricultural land is affected by tillage erosion, with rates reaching as high as 4 Mg·ha −1 ·y −1 on 48% of agricultural land, when the combined effect of water and tillage erosion is considered. Thus, the most endangered area is almost 1.5 times higher than the area endangered only by the effects of water erosion. The relative contribution of tillage erosion on total soil loss was estimated at 20%-30%. These findings will surely be useful in planning anti-erosion and soil protection strategies in general and in the design and implementation of adequate anti-erosion measures. An innovation of the study reported herein, compared to the previously presented erosion risk maps, is the quantification and visualisation of potential hazards due to tillage erosion, largely neglected under the conditions of the Czechia to date, despite its significant influence. Furthermore, the results will surely prove useful to farmers as well as to the State administration to effectively design and establish protection policies for mitigating the impact of soil erosion on crop land use and productivity.
Software QGIS 3.10.2, Esri ArcGIS 10.4 (curvature calculation), environment for statistical computing R (version 3.5.1), and RStudio (Integrated development environment for Rversion 1.1.456) were used for calculations, analysis, and pre-processing. QGIS 3.10.2 was used for the construction of the main map. GNU Image Manipulation Program (GIMP, version 2.10.18) was used for preparing raster graphics used in the map and vector graphics editor Inkscape 1.0 was used for final graphic adjustments.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data Availability statement
The data used for the composition of the final map are available on request from the corresponding author. The data are not publicly available due to internal institutional restrictions, but it is possible to view them in the web-GIS application in the geoportal of the Research Institute for Soil and Water Conservation (https://geoportal.vumop.cz). The data used for map layer calculations are available from Czech Office for Surveying, Mapping and Cadastre (DMR 4G) and the Ministry of Agriculture of the Czech Republic (Land Parcel Identification System, data from Systematic soil survey of agricultural soils in Czechoslovakia). Restrictions apply to the availability of these data, which were used under license for this study.