Spatially Explicit Soil Compaction Risk Assessment of Arable Soils at Regional Scale: The SaSCiA-Model

Soil compaction caused by field traffic is one of the main threats to agricultural landscapes. Compacted soils have a reduced hydraulic conductivity, lower plant growth and increased surface runoff resulting in numerous environmental issues such as increased nutrient leaching and flood risk. Mitigating soil compaction, therefore, is a major goal for a sustainable agriculture and environmental protection. To prevent undesirable effects of field traffic, it is essential to know where and when soil compaction may occur. This study developed a model for soil compaction risk assessment of arable soils at regional scale. A combination of (i) soil, weather, crop type and machinery information; (ii) a soil moisture model and (iii) soil compaction models forms the SaSCiA-model (Spatially explicit Soil Compaction risk Assessment). The SaSCiA-model computes daily maps of soil compaction risk and associated area statistics for varying depths at actual field conditions and for entire regions. Applications with open access data in two different study areas in northern Germany demonstrated the model’s applicability. Soil compaction risks strongly varied in space and time throughout the year. SaSCiA allows a detailed spatio-temporal analysis of soil compaction risk at the regional scale, which exceed those of currently available models. Applying SaSCiA may support farmers, stakeholders and consultants in making decision for a more sustainable agriculture.


Introduction
Soil compaction represents one of the main threats to soils worldwide [1]; it causes increased surface runoff, soil erosion and nutrient leaching, while infiltration rate, plant growth, root growth and biological activity decrease [2][3][4]. Hence, soil compaction affects soil functionality, agricultural productivity, food security, flood risk and nutrient input to water bodies [5][6][7].
Soil compaction occurs when the exposure to field traffic exceeds soil strength [2]. It is commonly separated into topsoil and subsoil compaction. While primary tillage may reverse topsoil compaction, subsoil compaction persists in the long term [8,9]. Restoring compacted subsoil by e.g., deep loosening is cost and time-consuming; resulting unstable soil conditions prohibit wheeling with heavy machinery after loosening [10]. Preventing soil compaction, in particular subsoil compaction, is therefore the best way to preserve soil functionality.
A first step to prevent soil compaction is the knowledge about where and when soil compaction may occur. A spatial evaluation of soil compaction risk of arable soils is challenging. Highly dynamic natural processes such as the variation in soil moisture strongly influence soil compaction Figure 1 provides a schematic overview of the SaSCiA-model. Input data include weather information, soil data, crop information, machinery setups and field traffic days. Using this input data, the SaSCiA-model generates daily maps of soil compaction risk estimates for entire regions by coupling various models. In detail, it involves spatially explicit modelling of soil water content at daily intervals using the MONICA-model [34], -calculation of actual soil strength [11,15,35,36], -estimation of soil stress [32,37].

The SaSCiA-Model-Material and Methods
The model is implemented in "R" (version 3.3.1, [38]). Data and models employed are freely available to enable widespread application by interested users. The script including a test dataset is available on request from corresponding author. In case of adaptations or modifications, users can easily replace individual modules or data sets.
Sustainability 2018, 10, x FOR PEER REVIEW 3 of 29 models [15,32,35]. Considering the dynamic changes of soil properties depending on the present crop types and growing stage, we hypothesize that the SaSCiA-model supports a more realistic estimation of soil compaction risk at regional scale compared to currently available models. The term "regional scale" as it is used in the presented study means a spatial extent of landscape between tens to thousands of square kilometers. To demonstrate the applicability of the model, SaSCiA was applied at two study areas with different soil and weather conditions in northern Germany.
Applying the SaSCiA-model enables the production of daily maps of soil compaction risk, identification of areas exposed to soil compaction and estimation of optimal field traffic days. These regional maps are essential to understand where and to what extent soil compaction may occur. Furthermore, these maps may help to prevent soil compaction by e.g., showing the effects of adapted management strategies. Figure 1 provides a schematic overview of the SaSCiA-model. Input data include weather information, soil data, crop information, machinery setups and field traffic days. Using this input data, the SaSCiA-model generates daily maps of soil compaction risk estimates for entire regions by coupling various models. In detail, it involves spatially explicit modelling of soil water content at daily intervals using the MONICA-model [34], -calculation of actual soil strength [11,15,35,36], -estimation of soil stress [32,37].

The SaSCiA-Model-Material and Methods
The model is implemented in "R" (version 3.3.1, [38]). Data and models employed are freely available to enable widespread application by interested users. The script including a test dataset is available on request from corresponding author. In case of adaptations or modifications, users can easily replace individual modules or data sets.

Input Data and Data Preparation
The required input data for the model are weather data, soil data, present crop type, crop information (fertilizer applications; irrigation; dates of sowing, harvest and primary tillage), machinery information and days of field traffic. Soil data and present crop types must be provided as spatially explicit (raster) data. All other input data have to be stored as tables (CSV).

Weather Information
The weather information is necessary to calculate the actual soil moisture with the MONICA-model. The model requires daily values of temperature (average, minimum and maximum), wind speed, precipitation, relative humidity and sun duration (Table A1, Appendix A).

Soil Data
The soil data are used to calculate soil moisture and soil strength. A minimum set of soil parameters includes soil texture class, dry bulk density, gravel content, organic carbon content, soil aggregate structure, air capacity, field capacity, available field capacity and horizon depths. The soil data is required as soil profiles (soil data from soil surface to a depth of 2 m). Additionally, each soil profile must be identifiable by a unique soil identifier ("soil-ID"). The soil profile information has to be stored in a table (Table A2, Appendix A). Moreover, the spatial distribution of all soil profile/soil identifier must be provided as raster data.

Present Crop Types
Knowing the present crop types for the year of soil compaction risk assessment is essential. The cultivated crop type affects the soil moisture of a certain soil and determines the field traffic activities (e.g., times of sowing, fertilizing, spraying, harvesting, used machinery). The crop type information must also be provided as raster data. A reclassification of the crop types by unique crop identifiers ("crop-ID") converts the raster into a structure readable for the SaSCiA-model. Table A6 (Appendix A) lists the assigned crop-IDs.

Crop Rotation
To consider the effects of previously cultivated crop types on the water balance of the soil, a three-year crop rotation is assumed. The present crop type and the crop types of the preceding two years compose the crop rotation. Intersecting the spatial distribution of crop types from the three consecutive years results in a raster with all possible combinations of crop rotations; a unique identifier ("rotation-ID") clearly distinguishes all combinations. Figure A1 (Appendix A) demonstrates this principle exemplarily. Joining the "rotation-ID" and "crop-ID"-raster results in a "crop-rotation-ID"-raster, which is used to spatially allocate the crop-information in the SaSCiA-Model.
If the crop type of a previous year is unknown, a representative crop rotation has to be assumed, e.g., winter wheat, winter barley and sugar beet.
In addition to the spatial allocation of the individual crop-rotations, further crop information such as tillage, sowing and harvest dates and tillage depth needs to be added (Table A3, Appendix A).

Fertilizer and Irrigation Information
Information about nitrogen application and irrigation can be included to stimulate plant growth. This information must be stored in individual tables (Table A4, Appendix A). Input parameters include the "crop-rotation-ID", amount of nitrogen, kind of fertilizer, the application date and incorporation for fertilizer application. For irrigation, the "crop-rotation-ID", amount of irrigation water and date of irrigation are necessary. Fertilizer and irrigation information is assigned to the crop types based on the "crop-rotation-ID". The machinery used and machinery characteristic are crucial to calculate the mechanical stress imposed on the soil. Therefore, wheel load, tire inflation pressure and field traffic days need to be parameterized for each field. Due to a huge variety of agricultural machinery, the spatial allocation of machinery used on the different fields of a region is mostly unknown.
To enable both, a spatially and a crop type dependent differentiation, the user has to define typical machinery setups used in the study area for each crop type, based on e.g., literature review and manufacturer's information. In general, the same machinery setup will be applied for each crop type (e.g., maize, sugar beet) within the entire study area.
A similar assumption is required for field traffic activities. The time of field traffic is highly variable at regional scale. To achieve a region-specific approximation, reasonable periods of field traffic activities for each crop type have to be predetermined (if the real dates are unknown). Typical field operations, e.g., primary tillage, sowing and harvest form categories of periods. Each category contains a date for the beginning and ending of potential field traffic activity. The wheel loads and tire inflation pressures are related to each period and crop type. One value for wheel load and tire inflation pressure can be considered per period and crop type. Table A5 (Appendix A) shows an example of the table structure for the machinery and field traffic days.

Soil Moisture Modelling (MONICA)
Soil compaction risk mainly depends on the soil moisture. Wet soils have a low load bearing capacity and a high susceptibility to soil compaction [39,40]. Increasing desiccation reduces the susceptibility to soil compaction. The actual soil moisture therefore is essential to evaluate field traffic activities on soil functionality. Soil moisture varies strongly, temporally, spatially and with depth. Thus, a continuous update of soil moisture in space and time is necessary to increase the reliability of soil compaction risk assessment. At the regional scale, soil moisture models may overcome the limited field measurements to provide such an update. A variety of models exist to determine soil moisture; they often are included in crop-growth models (e.g., HYDRUS, DAISY). For the SaSCiA-model, the crop-model MONICA (version 2.0, [34]) was selected. MONICA computes daily soil moisture content at 10 cm depth intervals from top to a depth of 2 m. Several studies demonstrated the usability of the model (e.g., [41,42]). MONICA is applicable to freely available data and transferable to regional scale [41]. Nendel et al. [34] and Nendel and Specka [43] give a detailed description of MONICA.
Soil data, crop type and rotation data, fertilizer data and weather information form the input tables for soil moisture modelling with MONICA in the SaSCiA-model. Matching soil and crop tables results in all possible combinations between soil profiles and crop types. Each combination receives a unique identifier ("soil-crop-rotation-ID") and guarantees a spatial allocation at the end of calculation. The SaSCiA-model automatically converts all input-tables into the MONICA-specific table format. Afterwards, the MONICA model is applied to each "soil-crop-rotation-ID". The resulting tables include the calculated soil moisture for each soil and crop combination, forming the basis for subsequent soil strength calculation.

Soil Strength Modelling
Soil strength describes the resistance of a soil against an applied pressure [2]. Soil strength depends on several soil properties, e.g., soil texture, organic carbon content and soil structure. Wet soils have a reduced soil strength compared to dry soils; therefore, the soil strength changes continuously depending on soil moisture. Calculating soil strength in the SaSCiA-model follows an approach published by Rücknagel et al. [32], which considers varying soil moisture and therefore enables a soil compaction risk assessment for dynamic field conditions. Furthermore, the approach is transferable to the regional scale by combining spatially explicit soil and crop type data and soil moisture modelling. To calculate the soil moisture dependent soil strength, the soil strength at field capacity (water content at a suction rate of −6 kPa or pF 1.8) is necessary. Soil strength at field capacity may be measured directly, calculated from measured dry bulk density and aggregate density [44] or may be derived by pedotransfer functions [15].
At regional scale, measured values of soil strength and aggregate density are usually unavailable, necessitating the use of pedotransfer functions. Horn and Fleige [15] provide common pedotransfer functions (see introduction) to derive soil strength at field capacity. The structure of the pedotransfer functions enables an application with freely available soil data and the transfer to the regional scale. Lebert [45] and D'Or and Destain [28], however, demonstrated that the pedotransfer functions provided by Horn and Fleige [15] tend to overestimate soil strength of certain texture classes (e.g., loam).
Lebert [45] recommended the use of pedotransfer functions by Horn and Fleige [15] for sandy soils and DIN V 19688 [35] for loamy, silty and clay soils. The pedotransfer functions of DIN V 19688 [35] are based on the same data and assumptions as the approach by Horn and Fleige [15], but calculate more reliable values [45]. The used pedotransfer functions disregard organic soils. These soils are therefore excluded from further analyses, i.e., SaSCiA calculates the soil compaction risk only for mineral soils, not for organic soils. Table 1 lists the equations (Equations (1)-(4)) for calculating the soil strength at field capacity for various soil texture classes as implemented in the SaSCiA-model.
The cohesion (Coh) and angle of internal friction (Phi) are derived from the given soil texture class and the soil structure according to Horn and Fleige [15] and Lebert [45]. The calculated values represent the soil strength at field capacity (σ1.8) for each depth of the individual soil types in the study area.
Subsequently, soil strength at field capacity is joined with the actual soil moisture values calculated by the MONICA-model. Afterwards, SaSCiA-computes the moisture-dependent soil strength at daily intervals, using the pedotransfer functions provided by Rücknagel et al. [11] (Table 2, Equations (5)- (9)).
The SaSCiA-model calculates the total soil strength for all "soil-crop-rotation"-combinations in the study area and all defined periods/days of field traffic. The results are daily information of actual soil strength depending on actual soil moisture calculated. The actual soil moisture is calculated by the MONICA-model using the crop information (crop type and crop rotation), the soil data and weather information.

Soil Stress Modelling
Soil stress is applied by field traffic and depends, among other factors, on wheel load and tire inflation pressure. Several approaches exist to calculate the soil stress and stress distribution in the soil (e.g., [18,29,48]). These approaches, however, often need various input parameters (e.g., exact tire descriptions), which are unavailable at regional scale. Koolen et al. [37] published a generalized method, which only requires wheel load (WL; in kg), tire inflation pressure (TIP; in bar), concentration factor (ν) and the depth (Depth; in cm) to calculate the soil stress at a certain depth (σ Depth ; in kPa) (Equation (14)): This approach is transferable to the regional scale and therefore was integrated in the SaSCiA-model. Rücknagel et al. [32] developed an equation to calculate the concentration factor (ν) based on soil strength at field capacity and actual soil moisture in percent of field capacity (%FC) (Equation (15)): The concentration factor determines the stress distribution in the soil. Söhne [49] developed this approach and classified the concentration factor into 4 (for dry soils), 5 (for moist soils) and 6 (for wet soils). Equation (15) determines the concentration factor for any soil moisture. The SaSCiA-model contains both soil stress equations (Equations (14) and (15)). The daily soil moisture calculated by the MONICA-model and the calculated soil strength at field capacity (pF 1.8) is used to determine the concentration factor. The present crop type, provided field traffic days and the date of calculation define the wheel load and tire inflation pressure used for soil stress calculation. When the present crop type is sugar beet, for instance, a wheel load of 11.200 kg and tire inflation pressure of 2.2 bar is used for stress calculation at harvest time (e.g., 30 September 2016). At the same time, but with crop type maize, a wheel load of 5.100 kg and tire inflation pressure of 2.2 bar are selected. Thus, SaSCiA conducts the calculation of soil stress for each single field in the study area for each day; on days without traffic activity, no stress is applied to the soil.

Soil Compaction Risk Evaluation
Field traffic affects various soil functions negatively when the applied soil stress exceeds the soil strength. Subtracting the applied soil stress from soil strength (Equation (16)) results in the so called "Soil Compaction Index" (SCI; [32]), which is used in the SaSCiA-model: SCI values (in log(kPa)) smaller than or equal to 0 indicate that expected field traffic will not affect soil functionality (no soil compaction risk). Values greater than 0 point to an increasing soil compaction risk due to plastic soil deformation, going along with degrading soil functionality [32,33]. Table 4 lists the classification of the SCI; each SCI-class represents a soil compaction risk class, ranging from "no risk" to "extremely high" [32,33].
As the final result, SaSCiA computes the "Soil Compaction Index" (Equation (16)) ranked into different classes of compaction risk (Table 4) for each soil-crop combination on a daily basis. The model outcomes are spatially explicit daily maps of the soil compaction risk and associated tables summarizing the area statistics throughout the entire region. Table 4. Classification of "Soil Compaction Index" (SCI, in log(kPa)) [32,33].

SCI Soil Compaction Risk
Very high >0. 4 Extremely high

Study Areas
The SaSCiA-model was applied to two study areas in northern Germany. The first one, named "Adenstedt" with a total area of 336 km 2 , is located in the southern part of Lower Saxony, the Lower Saxon Loess Hill Country. The soil parent material ranges from shallow layers of clayey and sandy weathering residuals at the hilltops, deeply weathered loess along the hill slopes and loamy deposits in the valleys, forming a wide variety of soil types. Predominating soil types under arable use are Luvisols and stagnic Luvisols, which are highly susceptible to soil compaction. Typical crop rotations in this region include sugar beet, (Beta vulgaris L.), winter wheat (Triticum aestivum L.) and winter barley (Hordeum vulgare L.) with an increasing share of silage maize (Zea mays L.) in the last decades.
The second study area, named "Kummerow" with a total area of 2500 km 2 , is located in the center of Mecklenburg-Vorpommern. Kummerow is part of the Northern German Plain. The parent material is glacial till, pervaded by numerous rivers and lakes. Typical soil types are Luvisol, Cambisol, Fluvisol, Gleysol and Histosol. Most common crop types are winter wheat and rapeseed (Brassica napus L.).

Weather Information
In Germany, the German Weather Service (DWD) provides weather data on a web-based platform (WebWerdis, [50]). All necessary input data (minimum, maximum and average temperature, wind speed, precipitation, relative humidity, sun duration) are available free of charge. The data are updated daily, allowing the calculation of soil moisture and consequently the soil compaction risk at almost real-time field conditions. The weather stations Teterow and Sukow-Levitzow were chosen for the Kummerow area, weather station Alfeld for the Adenstedt area.

Soil Information
A soil map at a scale of 1:200,000 (BUEK 200, [51]) was used for both study areas. The Federal Institute for Geosciences and Natural Resources (BGR) provides the map free of charge. Soil maps with higher resolution are only available on payment. The BUEK 200 consists of a shape file including the geometries of the individual soil types and related tables with the elementary soil attributes (e.g., soil texture). Based on the soil texture class, the dry bulk density and content of organic matter the more complex physical soil properties such as air capacity, field capacity, available field capacity and the wilting point were derived using the pedotransfer functions provided by Wessolek et al. [52].
They provide classified tables (34 soil texture classes, five dry bulk density classes, four organic matter classes), which enable the derivation of the target values. An additional script was created to automatically calculate the necessary soil properties.

Present Crop Types and Crop Rotation
Freely available data of present and spatially explicit data on cultivated crop types are rare. Remote sensing data are a valuable source of present and former land use/land cover and therefore on cropland and crop types mapping (e.g., [53][54][55]). The 'Semi-Automated Classification' plug-in [56] for QGIS [57] provides a freely available implementation of classical supervised image classification approaches. These approaches are able to learn the spectral signatures of target classes from user-defined training pixels and categorize image pixels according to statistical similarities in their reflectance behavior [58]. To this end, a certain amount of field mapping data on crop types is required [59]. Crop type mappings in the field were conducted in Kummerow (2015) and in Adenstedt (2016) during the summer months from July to September. Field mapping data covered the following target classes: cereals (e.g., winter wheat, rye, barley), maize, winter rapeseed, sugar beets, and grassland. Google Earth imagery supported the identification of less dynamic land cover classes such as water bodies, evergreen/deciduous forest and sealed areas. The field mapping data were divided randomly into a validation and training data set. Recent studies demonstrated the suitability of multispectral Landsat 8 and Sentinel-2 datasets for crop type mapping with a spatial resolution between 10 and 30 m (e.g., [60][61][62]). The United States Geological Survey (USGS) and the European Space Agency (ESA) provide free of charge Landsat 8 [63] and Sentinel-2 [64] data on a regular basis (every 5-16 days).
For Kummerow, four Landsat 8 scenes were available during the vegetation period in 2015. Atmospherically corrected Landsat 8 Surface Reflectance products provided by USGS [65] were delivered geo-referenced on Transverse Mercator UTM Zone 33N (WGS84). The data showed a highly accurate co-registration (RMSE < 8 m).
For Adenstedt, Sentinel-2A data form the database for the vegetation period 2016; due to cloud coverage, however, three images are available (processing baseline: 02.04). Sentinel-2A products of processing baseline 02.04 showed highly accurate co-registration (RMSE < 11 m, Transverse Mercator UTM 32 N WGS 84, [66]). Atmospheric correction was conducted using Sen2Cor (version 2.1.2; [67]), which is implemented in the SNAP toolbox (version 5.0; [68]). Table 5 lists the acquisition dates for both study areas. Using the field mappings and satellite data, a supervised classification (minimum distance; distance threshold: 0.2) was conducted for each study area with the plug-in 'Semi-Automated Classification' [56] developed for the open source software QGIS (version 2.18.7; [57]). To this end, atmospherically corrected satellite data from the different dates were merged into one dataset and then classified. Thus, class separability benefited from phenological characteristics (multi-temporal classification). Figure 2a shows the classification for Adenstedt, Figure 2b for Kummerow.
Kummerow), which were relevant for the SaSCiA-model, performed well with class accuracies >80%. The spectral similarity of winter wheat, rye and barley hampered their differentiation; they therefore were aggregated as cereals.   Accuracy assessment followed the suggestions by Foody [59] using the independent field mapping data for validation. Table 6 lists the percentage share of each class and calculated accuracy measures separately for each study area. Both classifications exceeded a widely reported orientation value of overall accuracy, i.e., 85% [59] indicating sufficient classification accuracy of the entire data set. In particular, the classes cereals, winter rapeseed, maize (Kummerow) and sugar beets (except Kummerow), which were relevant for the SaSCiA-model, performed well with class accuracies >80%. The spectral similarity of winter wheat, rye and barley hampered their differentiation; they therefore were aggregated as cereals.
Both study areas showed a high share of arable land, i.e., 49.8% at Kummerow and 44.9% at Adenstedt. Cereals predominated in both study areas. In Kummerow, winter rapeseed had the second highest share, while sugar beets were second in Adenstedt. In both study areas, unclassified pixels appeared as entire fields, which probably resulted from spectral signatures (crop types) uncovered by the field mappings. PA = producer's accuracy is the share of all correctly classified pixels and the sum of all validation pixels of a class; UA = user's accuracy is the share of all correctly classified pixels and the sum of all classified pixels of a class; overall accuracy is the share of all correctly classified pixels and the sum of all validation pixels; the kappa coefficient additionally includes erroneously classified pixels and considers agreement by chance; * due to rare cultivation in Kummerow area all mapped sugar beet fields were used for classifier training impeding an independent validation.
Both study areas lack information on crop types mapping in the preceding two years impeding a remote sensing based derivation of crop rotations (e.g., [69]). Assuming common crop rotations based on crop types in 2015 (Kummerow) and 2016 (Adenstedt) allowed considering effects of preceding crop types on soil moisture calculation. For instance, when the present crop type was sugar beet, the crop types for former years were winter wheat, resulting in a crop rotation winter wheat, winter wheat and sugar beet.

Fertilizer and Irrigation Information
At a regional scale, fertilizer types and applied amount differ; therefore, the most common fertilizers and averages of applied amounts for the regions were chosen based on the information provided by Achilles et al. [70] and KTBL [71].
As no information about irrigation was available, the irrigation was set to zero for both study areas. Table 7 shows machinery setup used in this study, in particular wheel load and tire inflation pressure. The machinery selection is based on literature review (e.g., [33]) and manufacturer information (e.g., Grimme for sugar beet harvester). The maximum possible wheel load was assumed for soil compaction risk assessment. For instance, the wheel load of fully loaded combine harvester was used and not the empty or half-filled weight. Since explicit field traffic days were unavailable, general information provided by Achilles et al. [70] and KTBL [71] were used to define reasonable mean periods for field traffic activity (Table A5, Appendix A).

Results and Discussion
The SaSCiA-model as presented in the Material and Methods section was applied to the data sets listed in the Section 3 ("The SaSCiA-model application"). The following section contains selected results of both study areas, which are discussed to demonstrate the spatio-temporal advantages and limitations of the developed model. Figure 3 illustrates the spatial distribution of soil compaction risk in three different depths at Adenstedt on 7 August 2016 as modelled with SaSCiA. Table 8 summarizes the area percentage for each compaction risk class and each of the four crop types. The selected depths are 20 cm (topsoil), 35 cm (soil directly beneath the maximum tillage depth) and 50 cm (subsoil).

Spatial Distribution of Soil Compaction Risk
On 7 August 2016, field traffic potentially affected 80.13% (Table 8) of the total arable area at Adenstedt. The areas covered by maize (10.36%) and rapeseed (9.51%) were not endangered through soil compaction since no field traffic activity took place at this time.
In 20 cm depth, 73.52% showed "high" to "extremely high" risk. Only 3.99% were classified as "no risk" during field traffic activity. The prevailing red colors in Figure 3a clearly demonstrate the high share of areas classified as "high" to "extremely high". Bluish areas ("no risk" to "medium" risk) were distributed over the entire study area without showing any spatial pattern.
In 35 cm depth, the area proportions changed; only 2.30% of the total area, which related entirely to "cereals", was classified as "high" or "extremely high". The percentage area classified as "no risk" during field traffic increased up to 21.12%. Figure 3b highlights the shift with predominating bluish tones, especially light blue. Two regions with higher soil compaction risk are prominent, i.e., a region with "very high" soil compaction risk in the southwestern part and an area with "high" soil compaction risk along a north-to-south corridor in the East.
In 50 cm depth, most of the arable area (78.90%) showed "no risk"; only 1.12% was classified as "high". Dark blue colors therefore predominate in Figure 3c except in a region in the southwest with "high" and some fields in the northwest with "low" soil compaction risk.
The spatial detail of the soil compaction risk maps depends on the spatial resolution of the soil and crop type data used for soil moisture modelling. The spatial resolution of crop type mapping was 20 m (Adenstedt) and 30 m (Kummerow). As field sizes are often several hectares, the spatial resolution of the crop type information was sufficient for the soil compaction risk assessment. The soil map used in the study area, however, was only available at a scale of 1:200,000. Hence, the spatial distribution of soil properties is highly aggregated and small-scale variations in soil properties remain disregarded. Using a more detailed soil map within SaSCiA-model may enable a more differentiated soil compaction risk assessment. was 20 m (Adenstedt) and 30 m (Kummerow). As field sizes are often several hectares, the spatial resolution of the crop type information was sufficient for the soil compaction risk assessment. The soil map used in the study area, however, was only available at a scale of 1:200,000. Hence, the spatial distribution of soil properties is highly aggregated and small-scale variations in soil properties remain disregarded. Using a more detailed soil map within SaSCiA-model may enable a more differentiated soil compaction risk assessment.

Temporal Variation of Soil Compaction Risk
The spatial distribution of soil compaction risk for individual days enables the analysis of the temporal variation for an entire year. Temporal changes of soil compaction risk depend on weather conditions, related crop growth and varying wheel loads. Figure 4 illustrates an example for different depths (20, 35 and 50 cm) at Adenstedt. The area percentage of the individual soil compaction risk classes (ordinate) is highlighted for each day in 2016 (abscissa).
The percentage area affected by field traffic varied between 0% and 100% for arable land throughout the year. No field traffic occurred during winter (January, second half of November and December), in the first half of March and from mid-June to mid-July leading to 0% affected area. During the second half of September and second half of October, field traffic affected the entire arable area, i.e., 100%. During the remaining months, the percentage varied between 9.8% and 85%. Figure 4 depicts strong leaps of the percentage of soil compaction risk area between certain dates, e.g., 9.5% on 31 July 2016 and 80% on 1 August 2016. These sharp transitions marked the edges of field traffic periods defined in SaSCiA prior to the model run and depend on the particular crop type. This explains, why the area percentage of soil compaction risk at a given date remained the same for all depths (Figure 4a-c), while the share of the soil compaction risk classes varied at the different depths.
For the topsoil (20 cm), the SaSCiA-model revealed the highest soil compaction risk with a high proportion of "high" and "extremely high" risk classes throughout the year. During spring, almost each field traffic activity resulted in an "extremely high" soil compaction risk. In November, each field traffic activity resulted in "high" and "very high" soil compaction risk. In spring and autumn, the high percentage of "high" risk area resulted from high soil moisture. In northern Germany, most precipitation occurs from November to March, while plant transpiration is very low. The wet soil therefore has a low soil strength and is unable to withstand the applied stresses, even of low wheel loads. In summer, the amount of precipitation is low, while plant transpiration rate is high. Soil strength increases with decreasing soil moisture, resulting in a reduced soil compaction risk. Accordingly, the percentage of soil compaction risk classes "low" and "no risk" increased at the beginning of August. Nevertheless, Figure 4 shows that a high percentage of fields remained highly exposed to soil compaction in summer and early autumn. At this time, the harvest activities led to highest wheel loads (5000-10,000 kg), causing soil pressures that exceeded the soil strength even of dry soils.
Below the tillage depth (35 cm) and in the subsoil (50 cm), the soil compaction risk was low for the entire cropping season. At a depth of 35 cm, "no risk" areas predominated. From February to the end of June, SaSCiA revealed a "high" or "extremely high" soil compaction risk only for 3% of the arable land. "Low" to "medium" classes shared 23%. The majority of the trafficked area, however, showed "no risk". During the first half of August, only 3% were classified as "high" and "very high", 55% were classified as "low" to "medium". From the second half of August until the end of the year, all field traffic activities resulted in "no risk" of soil compaction. At a depth of 50 cm, the soil compaction risk decreased continuously compared to the depths at 35 cm. The class "extremely high" never occurred; classes "high" and "very high" had a maximum of 2%. Soil compaction risk classes "low" and "medium" resulted in a maximum of 7%. The majority of field traffic activities were classified as "no risk". The low soil compaction risk in autumn contradicts the expectation that maize and sugar beet harvest cause the highest soil compaction as reported by e.g., Nevens and Reheul [72], Peth et al. [13] and Destain et al. [40]. In 2016, precipitation sum from July to September was 50% lower (93.2 mm) compared to long-term measurements (184.7 mm) [50]. The low precipitation led to dry soils especially in the subsoil. The dry soil compensated most of the applied soil pressure in the topsoil, resulting in a less affected subsoil.
To demonstrate the influence of weather conditions, the soil compaction risk assessment for Adenstedt was repeated with the weather data from 2014. In 2014, precipitation sum was 240.6 mm [50] from July to September representing a year with wet soil conditions during harvest. All other input data remained the same for this scenario. Figure 5 shows the results of the wet scenario at a depth of 35 cm: a clear increase in soil compaction risk for harvest times (16 July to 15 August 2016 and 16 September to 15 October2016) emphasizes the effects of high wheel loads during wet soil conditions compared with dry soil conditions in Figure 4. The results of the wet weather scenario are in agreement with the findings of the former mentioned studies [13,40,72]: high wheel loads during maize and sugar beet harvest increase soil compaction risk. and 16 September to 15 October2016) emphasizes the effects of high wheel loads during wet soil conditions compared with dry soil conditions in Figure 4. The results of the wet weather scenario are in agreement with the findings of the former mentioned studies [13,40,72]: high wheel loads during maize and sugar beet harvest increase soil compaction risk.

Spatio-Temporal Variation for Single Crop Types
The SaSCiA-model enables analyses of varying soil compaction risk for single crop types. Selecting only one crop type (e.g., sugar beet) enables the analysis of varying soil properties on soil compaction risk. Figure 6 shows an example of the temporal dynamics in soil compaction risk at a depth of 20 cm for areas cultivated with cereals at Kummerow in 2015. Since cereals was the focused crop type and SaSCiA considers similar field traffic for one crop type, the area affected by field traffic activity was either 0% (no field traffic) or 100% (field traffic on all cereal cultivated areas). Six periods had no field traffic activity. During the year, however, all soil compaction risk classes occurred, whereas spring and autumn showed the highest soil compaction risk. Changes in soil compaction risk classes became apparent within a traffic period. As an example, Figure 7 shows the changes in soil compaction risk during the cereal-harvest period (1 August 2015-15 August 2015) in the Kummerow area.

Spatio-Temporal Variation for Single Crop Types
The SaSCiA-model enables analyses of varying soil compaction risk for single crop types. Selecting only one crop type (e.g., sugar beet) enables the analysis of varying soil properties on soil compaction risk. Figure 6 shows an example of the temporal dynamics in soil compaction risk at a depth of 20 cm for areas cultivated with cereals at Kummerow in 2015.

Spatio-Temporal Variation for Single Crop Types
The SaSCiA-model enables analyses of varying soil compaction risk for single crop types. Selecting only one crop type (e.g., sugar beet) enables the analysis of varying soil properties on soil compaction risk. Figure 6 shows an example of the temporal dynamics in soil compaction risk at a depth of 20 cm for areas cultivated with cereals at Kummerow in 2015. Since cereals was the focused crop type and SaSCiA considers similar field traffic for one crop type, the area affected by field traffic activity was either 0% (no field traffic) or 100% (field traffic on all cereal cultivated areas). Six periods had no field traffic activity. During the year, however, all soil compaction risk classes occurred, whereas spring and autumn showed the highest soil compaction risk. Changes in soil compaction risk classes became apparent within a traffic period. As an example, Figure 7 shows the changes in soil compaction risk during the cereal-harvest period (1 August 2015-15 August 2015) in the Kummerow area. Since cereals was the focused crop type and SaSCiA considers similar field traffic for one crop type, the area affected by field traffic activity was either 0% (no field traffic) or 100% (field traffic on all cereal cultivated areas). Six periods had no field traffic activity. During the year, however, all soil compaction risk classes occurred, whereas spring and autumn showed the highest soil compaction risk. Changes in soil compaction risk classes became apparent within a traffic period. As an example, Figure 7 shows the changes in soil compaction risk during the cereal-harvest period (1 August 2015-15 August 2015) in the Kummerow area. At the beginning of the harvest period (Figure 7a; 1 August 2015), soil compaction risk was heterogeneously distributed, whereas classes "high" to "extremely high" (81.52%) predominated. Only 3.67% of the area showed "no risk" and 14.81% "low" to "medium" soil compaction risk. Rainy days in July (precipitation sum of 72.4 mm; DWD) caused an increased soil moisture with decreased soil strength. In the following days, however, the desiccating soil led to a decreasing soil compaction risk in the entire region.
The analysis of the cereal-harvest in the Kummerow area ( Figure 7) demonstrated how weather conditions may influence the spatial distribution of soil compaction risk on short term. Furthermore, the different soil compaction risk classes for a particular day highlighted the influence of varying soil properties; as crop type, weather information and machinery setup were the same for the entire region, the spatio-temporal variations in soil compaction risk result from different soil conditions.

Advances and Limitations of the SaSCiA-Model
The applicability of the SaSCiA-model has been demonstrated for two study areas in Germany with different spatial extent and varying land cover, pedological characteristics and weather conditions. The study highlights that the SaSCiA-model generates daily maps of soil compaction risk, which can support a sustainable field management and field traffic strategies to maintain various environmental functions of soils. The potential applications of the SaSCiA-model are manifold: At the beginning of the harvest period (Figure 7a; 1 August 2015), soil compaction risk was heterogeneously distributed, whereas classes "high" to "extremely high" (81.52%) predominated. Only 3.67% of the area showed "no risk" and 14.81% "low" to "medium" soil compaction risk. Rainy days in July (precipitation sum of 72.4 mm; DWD) caused an increased soil moisture with decreased soil strength. In the following days, however, the desiccating soil led to a decreasing soil compaction risk in the entire region.
The analysis of the cereal-harvest in the Kummerow area ( Figure 7) demonstrated how weather conditions may influence the spatial distribution of soil compaction risk on short term. Furthermore, the different soil compaction risk classes for a particular day highlighted the influence of varying soil properties; as crop type, weather information and machinery setup were the same for the entire region, the spatio-temporal variations in soil compaction risk result from different soil conditions.

Advances and Limitations of the SaSCiA-Model
The applicability of the SaSCiA-model has been demonstrated for two study areas in Germany with different spatial extent and varying land cover, pedological characteristics and weather conditions. The study highlights that the SaSCiA-model generates daily maps of soil compaction risk, which can support a sustainable field management and field traffic strategies to maintain various environmental functions of soils. The potential applications of the SaSCiA-model are manifold: (i) Real-time soil compaction risk assessment The SaSCiA-model enables a real-time soil compaction risk assessment, if sufficient data (especially weather and crop type data) is available. Currently available spatial soil compaction models show only the potential soil compaction risk, which is of limited importance for decision-making under real field conditions. The spatio-temporal high-resolution maps produced by SaSCiA show areas with varying (time and space) soil compaction risk. SaSCiA, therefore, can be potentially used as communication and visualization tool. Thus, the up-to-date soil compaction risk maps may support farmers, stakeholders and consultants in making decision for a more sustainable agriculture.
(ii) Retrospective soil compaction risk assessment The SaSCiA-model is also applicable for retrospective soil compaction risk assessments, as demonstrated in the case studies. As long as crop type information is available (e.g., derived from field mappings or from remote sensing data), the SaSCiA-model can be used to obtain retrospective daily soil compaction risk maps for the last year, decade or even longer (the Landsat archive provides data back to the 1980s). Additional information of former crop types derived by farmers or public institutions may also contribute to generate typical crop rotations for study areas and, therefore, long-term analyses of soil compaction risk.

(iii) Hot-spot detection
The continuous application of the model for several years may identify areas with high soil compaction risk year by year. This hotspot identification may help to prevent further soil compaction by means of adapted field cultivation. Since soil compaction reduces the infiltration rate and increases runoff, detected hot-spot areas may be subject to further investigations, e.g., regarding risk assessment of soil water erosion. Accordingly, Alaoui et al. [7] recently published a review in which they advise the use of remote sensing and new methods of soil compaction mapping for flood analysis and flood prevention. The SaSCiA-model may contribute to this recommendation.
(iv) Deriving days of trafficability Apart from the assessment of soil compaction risk, the SaSCiA-model may be used to calculate the maximum wheel load day-by-day until no risk (or low or medium risk) occurs. Such analyses, for instance, may provide recommendations to farmers on the maximum payload during harvesting.
(v) Scenario calculations The study shows that the SaSCiA-model can also be used for scenario calculation. Incorporating modified weather data further enables the analysis of climate change on soil compaction risk. This may help to develop strategies for future machinery setup, field traffic behavior or crop type selection. Effects of crop types on soil compaction risk may be investigated by changing present crop types, e.g., winter wheat instead of maize. This may be of particular interest in regions with changing cultivation practices as observed in Germany; here, the amount of area cultivated with maize extremely increased in the last decade [73,74]. Since maize is one of the crop types causing highest soil compaction risk due to harvest in autumn under wet soil conditions, applying SaSCiA (present crop type (maize) versus an assumed crop type (e.g., winter wheat)) demonstrates changing soil compaction risks under different crop type cultivation.
The forecasting ability of the SaSCiA-model, however, depends on the quality and availability of input data and the assumptions necessary for a regional analysis. In addition, some limitations of the SaSCiA-model must be taken into account: The crop type is a key factor for soil compaction risk assessment. Without spatially explicit crop type information, SaSCiA is unable to calculate the soil moisture; related field traffic activities and wheel loads therefore cannot be assigned. To calculate real-time or near real-time soil compaction risk For the two study areas, Sentinel-2A and Landsat 8 satellite data accompanied by field mappings served as input for the crop type mappings. Nevertheless, satellite remote sensing data allow a retrospective assessment; it may also help to identify crop type for the present years. Separation of crop types using remote sensing further depends on their spectral variability, i.e., the spectral signatures must differ significantly at the acquisition day of the satellite. Furthermore, the spectral variability within a crop type has to be lower that the spectral variability between crop types.
During summer months, matured rapeseed may have a similar spectral behavior as cereals; during spring, rapeseed is in full bloom and differs significantly from cereals. After sowing in spring, both maize and sugar beets depict a similar spectral behavior as bare soil; consequently, they are spectrally easier to distinguish in late summer. After maize harvest, it is difficult to distinguish winter wheat as a subsequent crop and yellow mustard as catch crop, or weed growth. In general, images acquired directly after harvest or after sowing tend to impede an accurate crop type classification.
Moreover, the spectral/radiometric characteristics of the satellite sensors determine whether and which crop types can be identified [75,76]. Using satellite remote sensing for real-time crop monitoring therefore may be limited. Land use classifications using remote sensing data in spring or early summer, however, enable an assessment for the following harvest period. Using a series of land cover classifications for several consecutive years, typical crop rotations may also be derived (e.g., [69]). Applying this crop rotation scheme allows a prediction of the following crop type, which then enables a soil compaction risk assessment for the subsequent year.
(ii) Assumptions for input-data In SaSCiA, the spatial pattern of soil compaction risk results from spatially varying soil properties, crop types and derived soil moisture. Weather information, machinery information and field traffic days are assumed as constant values within the entire region.
For wheel load and tire inflation pressure, one value for each crop type and type of field traffic (e.g., sowing, spraying and harvest) can be specified. If more than one piece of machinery works in the field at the same time, e.g., during maize harvest, only one wheel load can be used for soil compaction risk assessment. In the case studies, the maximum wheel load, i.e., the worst case, was chosen. This may result in an overestimation of soil compaction risk, as wheel load changes dynamically; during sugar beet harvest, for instance, wheel load varies between the empty gross weight of 30 t and a fully loaded weight of 60 t. Some studies therefore focus on changes in wheel load at field scale (e.g., [29,77]). At regional scale, however, the wheel load affecting each raster cell remains unknown. Each part of a field may potentially be wheeled with the maximum load; hence, each raster cell may be affected with a maximum soil compaction risk. In SaSCiA, however, the wheel load and tire inflation pressure is user defined; it therefore depends on the user, which may also select the minimum or the average wheel load.
Due to a lack of information on the distribution of field traffic, the SaSCiA-model neglects rollover frequency. The number of wheel passages, however, affects soil functionality. Each wheel passage leads to a decrease of soil functions (e.g., [78]). Highest field traffic intensity occurs in the headlands (turning areas of the field) and the tramlines [77,79]. The position of tramlines and headlands is highly variable, depending e.g., on field geometry and working width, which aggravates considering its effects at a regional scale. Since each raster cell could potentially be trafficked, any raster cell may be susceptible to soil compaction. Some studies therefore focus on the prediction of field traffic distribution based on field geometry and used machinery setup (e.g., [80]). If these models are further developed, integration into the SaSCiA-model may overcome the disadvantage of missing rollover frequency.
The exact days of field traffic activity is another parameter, which is difficult to consider at a regional scale. The period for e.g., maize sowing lies within a range of weeks and the exact date of sowing depends on the farmer's decision. The SaSCiA-model therefore integrates periods for each field traffic operation to consider varying field traffic days. The real field traffic activity will most likely occur within the specified period resulting in a reliable soil compaction risk assessment. If available, exact field traffic days can be used.
Another important point is the assumption of static characteristic of further soil properties, for instance dry bulk density and soil structure. The highly dynamic change of soil moisture is part of the SaSCiA-model. Further soil properties such as dry bulk density may also vary during the season as a result of drying-wetting, settling after tillage or other processes. The same applies for the effects of root growth after sowing. Roots may change the soil structure and increase the soil strength during the season. A recently published study reviewed possible impacts of vegetation on soil strength and trafficability [81]. These changes are not considered automatically in the calculation of soil strength. However, effects of root/plant growth are considered in the MONICA-model. Thus, root/plant growth is considered in soil moisture calculation and, therefore, in the calculation of moisture dependent soil compaction risk. Apart from seasonal changes of soil properties, field traffic will affect soil properties. For instance, a traffic activity at the beginning of the year (e.g., fertilizer application) may increase the dry bulk density of the trafficked soil. The change in dry bulk density affects the subsequent calculation of soil strength and plant/root development and, therefore, the calculation of soil moisture.
As the SaSCiA-model uses the soil properties given in the soil table and calculates the soil compaction risk based on these properties, changes of soil properties as dry bulk density are not considered automatically. A manual change of soil properties is (until now) the only way to consider seasonal, tillage or traffic induced changes in certain soil properties.

(iii) Model validation
Additional work is required to provide a validation of the SaSCiA-model. Model accuracy, of course, depends on the quality of the input data. Soil information, for instance, is an important input for the SaSCiA-model, but is often only available in a low resolution, as was the situation in the case studies. Thus, the soil information is spatially aggregated and the depiction of the soil heterogeneity therefore is limited. This is, however, rather a problem of data availability than of the model itself. The individual components/models of the SaSCiA-model are well established and proven: the ability to calculate the soil compaction risk for two soil moisture contents (pF 1.8 and 2.5) using the approach by Horn and Fleige 2013 [15] at field, regional or higher scale was demonstrated by several studies [23,[26][27][28][29]. The MONICA-model was applied in studies to derive soil moisture and winter wheat yield at regional scale [34,41], which indicates its applicability in the SaSCiA-model. Götze et al. [33] and Jacobs et al. [82] successfully applied the approach from Rücknagel et al. [32] to evaluate the soil compaction risk at individual farms.
It is therefore necessary to verify the predicted soil compaction risk resulting from the interaction of the single model components. This requires further fieldwork and soil sampling conducted before and after any field traffic activity, as described by Batey [3]. The maps resulting from SaSCiA-modelling show the heterogeneous soil compaction risk at high spatial and temporal scales, thus enabling the detection of areas for appropriate soil sampling and model validation.

Conclusions
The SaSCiA-model enables a region-wide soil compaction risk assessment at a high spatial and temporal resolution. It combines well-established approaches and models [15,32,34,35] with spatially differentiated soil and crop information. Furthermore, the model integrates varying wheel loads, tire inflation pressures and field traffic days. The results are daily maps of soil compaction risk for entire regions in a spatial resolution exceeding already existing approaches. The applicability of the SaSCiA-model has been demonstrated for two study areas; soil compaction risk for both study areas has been calculated for entire years.
By using freely available data and open source software, a broad community of practicing experts, stakeholders and consultants involved in soil protection may apply the SaSCiA-model. Even if the model has some limitations as discussed, it currently is the only spatial approach concerning spatial and temporal changes in soil moisture, plant growth and wheel loads. Potential applications of the model are manifold and may range from retrospective, through current to potentially future soil compaction risk assessments.
Author Contributions: M.K. designed the study in collaboration with all co-authors. M.K. prepared the figures and wrote the majority of the text with contributions by K.D. M.K. developed the model and conducted data analyses; K.D. prepared and analyzed the remote sensing data. N.O. and R.D. contributed in reviewing and finalizing the manuscript.

Acknowledgments:
The authors thank the anonymous reviewers for their commitment and very helpful comments on improving the manuscript. We are grateful to James F. Petersen and Kilian Etter for the linguistic editing. The Federal Ministry of Education and Research (BMBF) supported this study within the framework of the BonaRes-initiative (Grant No.: 031A563C). We acknowledge financial support for publications cost by Land Schleswig-Holstein within the funding programme Open Access Publikationsfonds.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Structure and example of input weather information (date system: yyyy/mm/dd).

Date
T_min (     Definition of the Crop-Rotation-ID for the SaSCiA-Model Figure A1 demonstrates the principle of the creation of "crop-rotation-ID" exemplary. The "crop-rotation-ID" consists of six numbers and contains two information codes. The first three numbers represent the spatial distribution of a crop rotation ("rotation-ID"). The last three numbers represent the present crop type in the investigated year ("crop-ID").
A crop rotation is a sequence of crop types year by year at the same field/raster cell. The left-hand side of Figure A1 demonstrates the creation of the "rotation-ID". Each field/raster cell in each year exhibit a crop type (e.g., maize) and is classified by a three-digit number (e.g., 110) based on the classification listed in Table A6. Fields/raster cells with the same sequence of crop types for the three years resulting in a "rotation-ID". This is the case for the upper and lower examples in Figure A1. The field borders remain unchanged for the period, but crop types vary (101-110-101 and 112-101-114), resulting in two different rotation-IDs (101,000, 104,000). Looking at the example in the middle, the field geometries and crop types are the same for the first two-years. In the third year, the field is subdivided in two different crop types (110 and 112), resulting in two rotation-IDs (102000, 103000). The "rotation-ID" itself is result from a sequential numbering, starting at 101000.
The present crop in the investigated year (year 3) is appended to the rotation-ID. The result is the "crop-rotation-ID", containing the information of the previous crop type, the present crop type and their spatial distribution (right had side, Figure A1).