Use of Hydrological Models to Predict Risk for Rutting in Logging Operations

Abstract

Using hydrological models with a high temporal resolution to predict risk for rutting may be a possible method to improve planning of forwarder trails or to schedule logging operations in sites with low bearing capacity to periods when soil moisture content is at a minimum. We have studied whether descriptions of rut variations, collected in 27 logging sites, can be improved by using hydrological data, modeled by Swedish HYdrological Prediction for Environment (S-HYPE). Other explanatory variables, such as field-surveyed data and spatial data, were also used to describe rut variations within and across logging sites. The results indicated that inclusion of S-HYPE data led to only marginal improvement in explaining the observed variations of the ruts in terms of both “rut depths” within the logging sites and “proportion of forwarder trails with ruts” across the logging sites. However, application of S-HYPE data for adapting depth-to-water (DTW) maps to temporal changes of soil moisture content may be a way to develop more dynamic soil moisture maps for forestry applications.

1. Introduction

Increasing sustainable forest production to substitute non-renewable energy-intensive material and services is an effective way to reduce human-induced greenhouse gas emissions [1,2]. However, increased forest production in Sweden implies more intensive logging operations, utilizing heavy machinery all year round, to deliver saw logs, pulpwood, and other products with a continuous flow from forest to industry. Data from official Swedish statistics shows an increase in annual fellings, from 49.3 to 85.3 million cubic meters, between 1956 and 2016 in terms of five-year averages [3].

Fully mechanized cut-to-length (CTL) logging operations with heavy machines may cause soil compaction or rutting in forest soils [4]. Soil compaction implies reduced pore volume and pore connectivity between the soil particles, and a resultant negative impact on air infiltration and water drainage in soil layers [5,6]. Increased bulk density in compacted soil may increase root penetration resistance and reduce tree and seedling growth in extreme conditions [7,8,9]. Oxygen deficiency and increased rates of water retention in compacted layers may also increase mobilization of total mercury/methyl mercury from the forest soil to nearby water streams, posing a threat to aquatic organisms [10]. Rutting occurs when machines apply compression or shear forces to sensitive soils, which may cause displacement of the soil to the sides or to the middle of the tracks [11]. The extent of soil compaction and rutting depends on factors such as axle load, ground pressure, number of machine passages, and ground bearing capacity [4,12,13]. In turn, ground bearing capacity is determined mainly by soil texture and moisture content [12,14], and thereby varies spatially and temporally in forest landscapes. It is also influenced by factors such as the armoring effect of roots and aboveground biomass [15]. Soils with fine-grained texture and high moisture content are more prone to developing deep ruts during logging operations [16]. Soil moisture content has been shown to be an important factor for predicting soil resistance to penetration and the associated risk for rut formation in conjunction with logging operations in fine-grained soils [17].

Various technical and planning measures are used in Swedish forestry to reduce the risk for rutting. Important planning measures include scheduling of logging to avoid sensitive sites in wet conditions and trafficability prediction models to decide where to drive and where strip trails need to be reinforced with slash [18,19]. Sites with low bearing capacity can be scheduled for harvesting in dry periods when soil moisture content and, consequently, wet areas are at minimum, or in frozen periods. Sites with stronger bearing capacity may be planned for harvesting during thaw or rainy conditions. Forest companies in Sweden usually sort their sites to plan harvesting seasons based on rough estimations of bearing capacity. However, soils with the best bearing conditions constitute only 5.2% of forest land in Sweden [20], while 45.2% of forest land comprises soil in the second-highest bearing class. Sites in the second-highest soil class can be harvested during most of the year if operational adjustment for soil moisture and thawing is taken into account [21]. Planning tools are required for more precise prediction of risk for rutting, both for scheduling logging operation and for planning machine operating trails during logging operations in the second-highest soil class.

This is especially important, since indicators of climate change in Nordic regions with implications for forestry operations include a reduced number of winter days with snow cover, from 120 days in 1950 to 100 days in 2020 [22], as well as higher average winter temperatures, from approximately −4 °C in 1900 compared to approximately −2 °C in 2020 in Sweden [23]. The indicators are calculated as ten-year averages for the whole country, thereby covering variations from south to north. In Finland, projections of probable climate scenarios indicate that the duration of winter periods with suitable bearing capacity in peatlands, e.g., with 20 cm frost depth or 40 cm snow cover, will shorten by the end of the 21st century [24].

In the past decade, topography-based soil moisture models, such as depth-to-water (DTW) maps, have been used as trafficability prediction models by most Swedish forest companies [19]. The DTW maps provided by the Swedish Forest Agency are available as spatial data layers covering the entire country. DTW maps estimate soil wetness by calculating the least elevation differences between the land’s surface and the nearest open water surface, e.g., flow channels or lakes, based on digital elevation models [25]. DTW maps do not include precipitation and therefore provide temporally static descriptions of wet areas. Changing the threshold value for upstream areas, required for initiating modeled flow channels, has been used to adapt wet areas in DTW maps to seasonal and temporal changes [26,27,28]. However, attaining a proper threshold value, for required area to initiate flow channels or separate wet/dry areas, for different seasons is an iterative procedure requiring field verification for practical applications [27,29]. Reeves et al. [30] developed a predictive geospatial model for identifying areas more susceptible to soil disturbances in harvesting sites, based on data regarding topography, land cover, and harvesting season. The harvesting season (winter/non-winter) was quite coarsely classified, lacking the temporal resolution required for more detailed scheduling of harvest sites. Although these kinds of trafficability prediction models have been effective in reducing the extent of severe rut damage [31], their precision might be affected by the static nature of these models. Using hydrological models to incorporate temporal variation of soil condition in trafficability models could be a way to improve predicting the risk for rutting.

Jones and Arp [32] used observed and hydrologically modeled soil moisture content at different soil depths to show how soil resistance, measured as cone penetration depth, and daily soil moisture fluctuations influence rutting. The results were used to develop vehicle-specific maps of estimated rut depth using seasonally adjusted DTW maps in three Canadian study sites. Jones and Arp [33] further contributed to the modeling of spatiotemporal soil moisture and cone penetrability variability, using hydrologically predicted soil moisture, elevation, soil particle size, soil bulk density, and land cover. Salmivaara et al. [34] used hydrologically modeled soil moisture content in combination with empirical data, such as logging transportation mass and vehicle rolling resistance, to provide a framework for producing rut-depth maps across a Finnish study site. However, the general applicability of these frameworks at high spatial resolutions needs to be examined by studying the suggested models under different terrain conditions.

By studying various logging sites, being cut at different periods, and spread throughout Sweden, we have tried to cover both temporal and spatial variability of soil condition in relation to rut formations. The main objective of this study was to evaluate whether risks for rutting can be predicted based on hydrologically modeled data within and across logging sites. The hydrological model used in this study is called HYPE (HYdrological Prediction for Environment) and is normally used for daily predictions of soil condition. The version of HYPE adjusted for Swedish conditions is called S-HYPE. The following research questions were specifically addressed:

• Will inclusion of S-HYPE-modeled data improve existing descriptions of rut depth within logging sites?

• Will inclusion of S-HYPE-modeled data improve existing descriptions of proportion of forwarder trails with ruts across logging sites?

Improved prediction of risk for rutting will give possibilities to improve planning of forwarder trails and scheduling logging sites by considering temporal variation of soil trafficability.

2. Materials and Methods

An empirical database containing data regarding ruts and forwarder trails [35] was used in this study. The database contained data regarding 35 logging sites, of which 27 sites (341 ha) had the complete data required for this study. The logging sites were well distributed in terms of operation season, size, and location in Sweden (Figure 1) and mainly comprised Norwegian spruce (Picea abies (L.) Karst.) and Scots pine (Pinus sylvestris L.). The sites were harvested during 2014 and 2015, using cut-to-length (CTL) mechanized systems. Each system involved a single-grip harvester, weighing approximately 20 Mg, and a forwarder with a laden weight of approximately 40–45 Mg. Detailed technical specification of machines used on these logging sites were missing in the database, since the data had been collected to address other research questions than those studied here.

2.1. Input Data

Input data used in this study can be divided into three main categories (Table 1).

Field-surveyed data: Information about forwarder trails and visible ruts along the trails was procured in the field using the ArcGIS collector application. Ruts, defined as wheel-tracks with exposed mineral soil having a minimum depth of 10 cm at the deepest part and minimum length of 1 m, were surveyed as GPS-positioned point observations. Tracks on peat soil did not require exposed mineral soil to be counted as ruts. To simplify the field work, rut depth and length were estimated and collected in predefined classes. Depth classes were defined as 10–20 cm, 21–50 cm, and >51 cm, and length classes as 1–5 m, 6–10 m, 11–20 m, and >20 m.

Data about forwarder trails were also collected in the field, describing the estimated number of machine passages including both loaded and unloaded forwarder passages. This information was used to define forwarder trail type in classes: strip forwarder trails with 1–5 passages, base forwarder trails with 6–10 passages, and main base forwarder trails with >10 passages. Any ground protection measures, such as slash or temporary timber bridges over wet areas, were also recorded in the field for both ruts and the forwarder trails. These data were later reclassed into a binary categorical attribute, ground protection or no ground protection, in the inferential analysis. A total of 5021 GPS-positioned rut-points and 234 km of forwarder trails were used as field-surveyed data across the logging sites. A more detailed description of the sites, forwarder trails, and ruts is presented in Table A1.

Spatial data: Digital maps were consulted to collect spatial data about soil type, elevation, and estimated soil wetness at the studied logging sites. Soil type information was extracted from Quaternary soil maps of the Geological Survey of Sweden (SGU). Elevation and estimated soil wetness, modeled in DTW maps, were extracted from raster layers with 2 m resolutions available at the Swedish Mapping, Cadastral and Land Registration Authority and the Swedish Forest Agency, respectively. DTW maps were reclassed into binary classes, where areas with DTW ≤ 1 m were considered wet, and areas with DTW > 1 m as dry.

Hydrologically modeled data, (S-HYPE): The HYPE model has been developed by the Swedish Meteorological and Hydrological Institute (SMHI) to simulate flow channels and circulation of nutrients such as nitrogen and phosphorous [36]. The Swedish version of HYPE, S-HYPE [37], has been used as a platform for hydrological forecasting and warning services in Sweden since 2013. S-HYPE subdivides the Swedish landscape into 37,000 sub-basins, at the time of study, with areas of 700–1000 hectares. The sub-basins are, in turn, classified into 65 combinations of soil type and land use. The soil type data used in S-HYPE are extracted from soil maps of the Swedish Geological Survey (SGU), and are further simplified into classes of peatland, clay, till, thin soil or uncovered bedrocks, and fluvioglacial sediments. Using these data together with altitude and climatological data, the S-HYPE provides various hydrologically modeled variables, including soil moisture, ground water level, soil temperature, and frost depth, which were utilized in this study.

S-HYPE variables were provided as time series, covering a time span from 1981 to 2015. Hydrological data corresponding to logging periods for each of the sites were derived from S-HYPE time series. Wherever study units, e.g., ruts or logging sites, overlay more than one soil type or sub-basin, S-HYPE variables were weighted accordingly through the inferential analysis (Table A1). The mathematical operations were performed using ArcGIS 10.7 and MS Access software. A detailed presentation of all S-HYPE variables (a total of 30 variables) used in the study is provided in Table A2.

2.2. Statistical Inference

The sampled ruts in logging sites and sub-basin areas were analyzed by inferring variations of “rut depth” and “proportion of forwarder trails with ruts” using two regressor subsets:

• Field and spatial data only.

• Field, spatial, and S-HYPE data.

The subsets 1 and 2 differ mainly with regard to the regressors that reflect soil conditions. The soil moisture content estimated by DTW indices was only included in the first subset, whereas the S-HYPE variables were used to estimate soil conditions in the second subset of regressors.

2.2.1. Rut Depth Variation

The ruts used in this study had been originally surveyed for other purposes than the research questions addressed here and were therefore sampled at any distance (with an average of 9.5 m) between sampling positions. To avoid autocorrelative redundancy in sampled ruts, the dataset was resampled, and a minimum distance of 35 m was assigned between the sample points. This autocorrelative threshold distance was identified with Kriging technique using an anisotropic semi-variogram [38] and may be considered as a generally valid minimum sample distance in the type of boreal forest environments covered in the present study. The resampled dataset contained 1756 field-collected rut samples and was supplemented with map-collected sample points along the existing forwarder trails, where no rut damage had been identified. The latter type of sample points was called null ruts and they were extracted retroactively on top of digital logging maps. The addition of null ruts facilitated inferences regarding differences across disturbed and undisturbed forwarder trails and increased the sample size to 2063 sample points.

The soil moisture condition was described by DTW indices, extracted at each sample position, in regressor subset 1. Weighted-average S-HYPE variables at buffer zones around each sample position described the soil condition, using regressor subset 2. The sample points were matched with S-HYPE data by calculating the associated weighted averages of the S-HYPE time series variables, using information regarding logging periods, sub-basin areas, and soil types within buffer zones around each sample position. The buffer zones were made using a radius value (r) = rut length/2. The dependent variable was defined as the number of ruts within the predefined rut-depth classes (Section 2.1) and was assumed to follow a Poisson distribution [39]. The regressor variables shared by the two subsets included trail type, ground protection, and elevation.

Partial least squares (PLS) analysis was used [40] to estimate the true dimensionality of the regressor matrix, and to identify the associated regressor base. PLS was also used to create rank-ordered lists of relatively independent regressor variables. The top-ranked PLS regressors were used to fit a generalized mixed log-linked model [38] describing the observed variation of rut depths. Using a nested model design made it possible to account for the covariance structures assumed to stratify the dataset due to rut measurements within logging sites across sub-basins. In the mixed model design, logging site identities were exclusively used as random regressors when using field and spatial data as only regressors (subset 1). The logging site identities were nested within random sub-basin identities when S-HYPE variables were included in the rut depth prediction model (subset 2). More explicitly, the reason for the mixed effect design is that logging sites were harvested in different periods, by different operational teams that may have used different routines in their logging operations. Ruts within logging sites may therefore have properties not shared by ruts in other sites. The nested design was used because inclusion of S-HYPE variables required considering that individual sub-basins shared individual S-HYPE model outputs, presumably not shared by logging sites across different sub-basins.

The two subsets were compared in terms of efficiency in explaining the observed variation of the response variable, e.g., rut depth within logging site scale.

2.2.2. Proportion of Forwarder Trails with Ruts Variation

Potential regressor variables were compiled to infer the proportion of forwarder trails with ruts across logging sites. The proportion of forwarder trails with ruts was assessed by dividing the total rut lengths from all sample points by the corresponding length of forwarder trail within the respective classes of trail types and ground-protective measures. The proportion of wet areas within logging sites was also assessed using DTW maps, to reflect the soil moisture condition, when using the regressor subset 1. Logging sites with peat as dominant soil type (two sites) were excluded from the analysis due to large differences in bearing capacity between mineral soil and peat soil.

Three new variables were derived from S-HYPE data to assess the relative values of precipitation (Rel-prec), groundwater level (Rel-gwat), and soil moisture at root zone (Rel-srfd) to describe the soil condition in regressor subset 2. The relative values were calculated by dividing the average values of the associated S-HYPE variables during the logging operation (3–10 days) by averages of these variables for the reference period 1 January 2000 to 31 December 2014. The S-HYPE variables were calculated as weighted averages, using the proportion of soil type per logging site area as a reference prior to the assessment of relative values.

Using the method of parameterization described above, the dependent response variable, proportion of forwarder trails with ruts, was now a real number within the interval [0, 1]. The variation of this response variable was inferred with a logit-linked generalized mixed linear model [38]. The two regressor subsets (1 and 2) were again tested and compared with regard to their efficiency in explaining the variation of proportion of forwarder trails with ruts, over logging site scale. All explanatory regressors apart from logging site identities were considered as fixed effects. Logging site identities were considered as random variables in the evaluation of subset 1 and were (again) nested within random sub-basin identities in the evaluation of regressor subset 2.

Analyses of spatial data, including kriging analysis, were performed in ArcGIS 10.7. The statistical inferences (PLS analysis) were performed in Dell Statistica 13.0, and the generalized mixed linear model analysis using GENMOD procedure in SAS^® statistical software (Version 9.4, SAS Institute Inc., Cary, NC, USA), ArcGIS 10.7 (Esri, Redlands, CA, USA), MS Excel, and MS Access (Microsoft, Redmond, WA, USA) were used for database management.

3. Results

3.1. Rut Depth Variation

With PLS applied to the resampled data, using only spatial and field data as regressors (subset 1), a list of relatively independent rank-ordered variables explaining the observed rut depth variation was created. The high-ranked PLS regressors identified and used as input variables to the generalized linear models were trail type, ground protection, elevation, and DTW soil moisture. Trail type, ground protection, and elevation were identified as significant regressors (p < 0.05), explaining 18.8% of rut depth variation observed within logging sites according to the PLS analysis.

When S-HYPE variables were used to describe the soil condition in regressor subset 2, estimated second layer soil moisture (Sml2) was added to the top significant variables (p < 0.05) identified by PLS. Accordingly, a linear combination of trail type, ground protection, elevation, and Sml2 explained 19.3% of the rut-depth variation observed within logging sites (Table 2).

3.2. Proportion of Forwarder Trails with Ruts Variation

The proportion of forwarder trails with ruts was inferred (with a generalized mixed linear model) across logging sites, using the first subset regressor including trail type, ground protection, soil type, and proportion of wet areas in logging sites, e.g., proportion of site areas with DTW ≤ 1 m. The two (fixed) factors trail type and ground protection were once again identified as significant regressors (p < 0.05), explaining 33.1% of the observed variation in proportion of forwarder trails across the sites. The inference across logging sites was repeated with second regressor subset, including the relative S-HYPE variables Rel-srfd, Rel-prec, and Rel-gwat, in the regressor matrix instead of soil moisture estimation by DTW. The relative soil moistures (Rel_srfd) and (Rel_gwat) were added to the list of significant (p < 0.05) variables. A linear combination of relatively independent field, spatial, and S-HYPE regressors explained 35.4% of the response variation observed across logging sites (Table 3).

4. Discussion

An existing database on ruts and forwarder trail over 27 boreal logging sites throughout Sweden was used to analyze the hydrological S-HYPE variables within and across logging sites to evaluate their potential to predict risk for rutting. The PLS results of assessing the observed variation of ruts with different regressor subsets—(1) field and spatial data, and (2) a combination of field, spatial, and S-HYPE data—indicated that the S-HYPE variables could contribute to additional explanatory power, but at very small magnitudes.

When S-hype variables were added to the regressor basis, the explanatory power of the models improved from 18.8% to 19.3% for rut depth, and from 33.1% to 35.4% for proportion of trails with ruts inferences. The improvement was slightly more pronounced when the proportion of trails with ruts was inferred, which indicated that S-HYPE variables can perform better at larger scale, such as logging site areas compared to finer scales and rut depth positions within logging sites. The S-HYPE variable sml2, soil moisture content at the second soil layer, and relative values of S-Hype Rel_srfd and Rel_gwat were found to be significant factors (p < 0.05) in our study. The effect is logical, since both the rut depths and lengths, i.e., proportion of forwarder trail with ruts, increase at higher soil moisture/ground water levels. However, due to the material used, we cannot say that this is an effect of the weather and not of some other difference between the sites. This needs to be investigated further using a controlled experimental design, and not survey material, as used here.

Trail type and ground protection were identified as significant factors in both subsets, explaining rut depth/proportion of forwarder trails with ruts. Higher number of machine passages means in practice the passage of higher accumulated load over trail segments, which encourages rut development. The significance of trail type in explaining the ruts variations is in agreement with results reported by Marra et al. [41] and Eliasson [4]. The use of slash for ground protection distributes the machine loads over larger areas and can thereby reduce the risk for rutting when the measures are properly applied [42]. The field data of ruts and forwarder trail used in this study had been collected previously to evaluate how application of DTW maps could contribute to minimized rut formation. Based on this essential decision support material, logging planners and drivers were able to avoid wet parts of the terrain, or to strengthen them prior to passage to minimize rut formation. The logging operation method indicated the effectiveness of ground protection in minimizing the risk for rutting, but this has faded out the possible effect of soil moisture estimations by DTW maps for the same objective. Soil moisture estimations by DTW maps could not explain the rut depth variations within logging sites, nor the proportion of forwarder trail with ruts. The logging operation method also implied that the choice of forwarder trail, and hence the location of ruts, is far from random, which needs to be kept in mind when interpreting the results. Elevation values, extracted from DEM layers, were more effective for evaluation of rut variations at finer scale, e.g., rut depths within logging sites, but were not found effective for predicting the proportion of trails with ruts when elevation was aggregated over the whole area of logging sites.

S-HYPE variables, when added to field and spatial data, improved description of the observed variations of the ruts in terms of both “rut depths” within the logging sites and “proportion of forwarder trail with ruts” across the logging sites. The improvements were, however, quite marginal, which makes it difficult to extract more general conclusions. The low effect of S-HYPE variables may also be because they are partially derived from the field and spatial variables already in the model, and they therefore introduce covariance to the regressor matrix rather than add unique non-redundant information. We suggest that the main reason for the limited effect of S-HYPE variables is that the S-HYPE model is designed for operating at relatively large landscape scales, so it cannot reproduce the local effects that regulate rut formation. This might depend on uncertainties incorporated in the S-HYPE model. Uncertainties are intrinsic parts of hydrological modeling and may have a number of causes: data uncertainty, model parameter uncertainty, and model structure uncertainty [43].

Data uncertainty: The S-HYPE model is built on different classes of land use (e.g., forest/agriculture land), elevation, and soil type. The land use classification is quite coarse and does not account for spatial variation in standing trees or the understory vegetation type. Soil information is extracted from Quaternary soil maps of the Geological Survey of Sweden (SGU), whose quality and resolution vary across the nation, and may affect the outcomes of the hydrological modeled data.

Model parameter uncertainty: The parameters in S-HYPE are primarily linked to soil type and land use. There are, however, also local deviations of key parameters, called super parameters [44], by which the accuracy in streamflow is improved considerably. The focus of the parameter calibration is streamflow. The model simulates other internal variables, e.g., lake levels, groundwater levels, and snow depth, with a varying degree of accuracy. It can be calibrated to simulate internal variables with a high level of accuracy [36,45]; but without local adaption, the results are quite uncertain at the fine scale used here. The variation in time is usually more reliable than the absolute values, for instance, for groundwater levels. In turn, this means that the detailed information required within the scales of typical logging sites is lacking. In a study by Tyystjärvi [46], the spatial and temporal variation of soil moisture was estimated using three different hydrological models, JSBACH [47], SpaFHy [48], and Ecohydrotools [49], where modeled and observed soil moisture was compared across a study area in northwest Finland. The study indicated that all the models had difficulties in modeling small-scale spatial variation, particularly at the driest part of the terrain, but performed better when simulating the temporal variation of soil moisture over longer time spans. More spatially detailed input data were recommended to achieve accurate estimations.

Model structure uncertainty: The S-HYPE model is mainly developed to estimate the flow and turnover of water and nutrients, and for monitoring the quality of water resources. Detailed information, such as land cover and soil type, is upscaled by the S-HYPE model without affecting its main application areas, whereas the estimation of such variables at smaller landscape scales may suffer as a result. Furthermore, the location of the different soil types is not resolved within the S-HYPE model sub-basin units, with areas of typically 700–1000 ha. This implies that S-HYPE variables lack the spatial resolution required for estimating the local ground-bearing capacity, based on soil condition estimations within the logging sites or capturing the variation of bearing capacity across the logging sites.

The strength of the S-HYPE model lies in describing how hydrological landscape characteristics may vary with time across a specific area, rather than capturing the differences among parts of different sub-basin areas over the same time span. The dataset used in the evaluation consisted of information from different points in space rather than different points in time over the same area. The predictive power of the hydrological model is likely to be advice on when, rather than where, to avoid logging operations. This strength could be applied to improve DTW maps by adjusting threshold values concerning flow initiation areas (FIA), required for mapping flow channels, or to decide proper DTW limit to distinguish wet/dry areas. A 1 ha FIA threshold value and a DTW limit of ≤ 1 m to separate wet/dry areas are currently used in DTW maps available at the Swedish Forest Agency [19]. Higher values of relative soil moisture (Rel_srfd) or relative ground water level (Rel_gwat), compared to averaged conditions over longer periods, can indicate wetter soil moisture conditions. This information can be used to apply a lower FIA or higher wet/dry DTW threshold limit and thereby develop more dynamic soil moisture maps for forestry applications. However, new studies are required to evaluate exactly how to use S-HYPE variables for such purposes.

5. Conclusions

To predict risk for rutting, using a hydrological model applicable for planning of logging operations, the small-scale spatial variations of topography, soil type, and land cover need to be represented in the model. Further improvement of the S-HYPE model may be needed to adapt it to the requirements of identifying logging sites and forwarder trail with adequate bearing capacity. However, the use of S-HYPE data to adapt the DTW maps to temporal changes in soil moisture may be a possible alternative for developing more dynamic soil moisture maps for forestry applications.

Author Contributions

Conceptualization: J.S. and L.E.; methodology, S.M. and T.T.; software, S.M. and T.T.; validation, L.E., J.S. and T.T.; formal analysis, S.M. and T.T.; investigation, S.M.; resources, J.S.; data curation, S.M.; writing—original draft preparation, S.M.; writing—review and editing, all authors; visualization, S.M.; supervision: L.E. and J.S.; project administration, G.L.; data acquisition. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

We would like to thank Tomas Nordfjell, Carola Häggström, Ph.D. (Swedish University of Agricultural Sciences, SLU), and Linnea Hansson, Ph.D. (The Swedish Forestry Research Institute, Skogforsk) for their valuable comments on the manuscript.

Conflicts of Interest

Authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figure and Tables

Figure 1. (Left) Surveyed logging sites in Sweden. (Right) Distribution of ruts per hectare logging site area with the season of logging site operations indicated. Ruts were defined as machine tracks with a minimum depth of 10 cm and a minimum length of 1 m.

Table 1

Summary of input data used in this study: field-surveyed, spatial, and S-HYPE-data.

Input-Data	Contained Data	Data Category	Data Type	Reference
Rut length and depth	Length classes:1–5 m, 6–10 m, 11–20 m, >20 m.Depth classes:10–20 cm, 21–50 cm, >51 cm	Field surveyed	Shape file
Slash protection	Yes/No	Field surveyed
Trail type	Three classes: strip trails = 1–5 passes, base trails = 6–10 passes, main base trails > 10 passes.	Field surveyed	Shape file
Soil type	Soil type in top 50 cm of the soil according to Quaternary maps	Spatial data	Shape file of varying resolution: 1:25–100,000 and 1:750,000	Geological Survey of Sweden, https://www.sgu.se/en/ (accessed on 15 August 2020).
Soil moisture	Estimations by depth-to-water (DTW) index, converted to two classes, wet (DTW < 1 m) and dry (DTW > 1 m)	Spatial data	Raster (2 m)	The Swedish Forest Agency, https://www.skogsstyrelsen.se/ (accessed on 1 April 2020).
Elevation	Digital elevation model created from high-resolution laser scanning of Sweden	Spatial data	Raster (2 m)	The Swedish Mapping, Cadastral and Land Registration Authority, https://www.lantmateriet.se/ (accessed on 1 April 2020).
S-HYPE	Modeled hydrological variables, see Table A2	Hydrologically modeled data	Text datasheets per soil type and sub-basin area	https://www.smhi.se/forskning/forskningsenheter/hydrologisk-forskning/hype-1.557 (accessed on 12 August 2017).

Table 2

Description of regressor subsets and models efficiency, R², for inferring variation of rut depth within logging sites. Regressors marked with (*) were identified as significant at p < 0.05 level.

Regressor Subset	Fixed Regressors				Random Regressors	R2
(1) Field and spatial data	Trail type *	Ground protection *	DTM (elevation) *	DTW (soil moisture)	Logging sites	18.8%
(2) Field, spatial, and S-HYPE data	Trail type *	Ground protection *	DTM (elevation) *	Sml2 (soil moisture at second soil layer by S-HYPE) *	Logging sites nested within sub-basin areas	19.3%

Table 3

Description of regressor subsets and models efficiency, R², for inferring variation of forwarder trails with ruts across logging sites. Regressors marked with (*) were identified as significant at p < 0.05 level.

Regressor Subset	Fixed Regressors						Random Regressors	R2
(1) Field and spatial data	Forwarder trail type *	Ground protection *	Soil type	Proportion of wet area in logging sites, % (areas with DTW ≤ 1 m)			Logging sites	33.1%
(2) Field, spatial, and S-HYPE data	Forwarder trail type *	Ground protection *	Soil type	Relative precipitation (Rel_prec)	Relative soil moisture * (Rel_srfd)	Relative ground water level (Rel_gwat)	Logging sites nested with sub-basin areas	35.4%

Appendix A

Table A1

Description of logging sites, showing number of sub-basin areas in S-HYPE, total area (ha), number of ruts per hectare, total length of forwarder trail, total length of ruts, forwarder trail types, any ground protection of the forwarder trail, and soil type distribution according to SGU Quaternary maps.

Logging Site ID	Number of Sub-Basins Sustaining Sites	Area (ha)	Number of Ruts per ha	Total Length of Ruts (m)	Total Length of Forwarder Trails (m)	% Forwarder Trails Protected with Slash	Forwarder Trail Type Distribution			Soil Type Distribution
							% Main Base Forwarder Trails	% Base Forwarder Trails	% Strip Forwarder Trails	% Clay	% Glaciofluvial Sediment	% Till	% Peat	% Thin Soil and Bare Bedrocks
1	1	1.78	2.8	23	1760	65.94	1.89	60.07	38.04	-	-	100.00	-	-
2	1	2.33	34.3	303	1752	59.94	-	52.04	47.96	-	-	100.00	-	-
3	1	1.38	4.34	20	1388	84.41	-	60.43	39.57	-	-	85.52	14.48	-
4	1	6.74	10.69	240	7467	74.11	5.63	61.37	33.00	-	-	100.00	-	-
5	1	7.24	16.86	315	5406	86.68	0.00	67.04	32.96	-	-	100.00	-	-
6	2	28.87	36.34	3895	14,291	81.82	1.66	59.28	39.06	-	-	100.00	-	-
7	1	13.89	42.12	2160	6855	68.05	0.00	56.83	43.17	-	-	100.00	-	-
8	1	14.99	20.75	840	9406	80.20	2.46	63.12	34.42	-	-	44.82	51.26	3.93
9	1	16.1	14.16	740	12,529	84.28	0.00	64.77	35.23	-	-	98.25	1.75	-
10	1	8.92	2.69	75	7692	65.39	3.76	64.49	31.75	-	-	69.42	0.92	29.66
11	1	11.39	11.85	475	7432	63.48	0.00	51.86	48.14	-	-	99.63	0.37	-
12	1	48.56	13.43	2579	35,183	40.18	1.70	39.32	58.98	-	-	99.49	0.51	-
13	1	17.67	2.77	148	13,634	67.59	0.00	26.60	73.40	-	-	98.10	1.90	-
14	1	20.53	35.8	2235	14,805	71.61	1.83	46.69	51.48	-	-	100.00	-	-
15	1	3.06	15.38	163	1806	67.30	-	53.71	46.29	-	90.35	9.65	-	-
16	2	9.75	1.74	76	7493	85.67	3.30	58.37	38.33	-	-	2.71	3.68	93.61
17	1	19.74	22.95	1620	15,096	76.65	0.49	74.03	25.47	-	-	27.80	72.20	-
18	1	4.37	2.29	30	2922	80.04	3.07	58.22	38.71	-	-	100.00	-	-
19	1	7.36	1.36	50	4697	91.98	-	58.44	41.56	-	-	100.00	-	-
20	1	21.19	5.71	463	16,272	36.91	2.00	74.39	23.60	-	-	99.59	0.41	-
21	1	7.49	20.96	500	6137	76.48	2.04	64.85	33.11	-	-	98.39	1.61	-
22	2	5.03	42.72	1128	2739	40.27	12.74	13.48	73.78	8.70	-	78.39	12.91	-
23	1	15.95	45.15	2258	10,908	37.67	1.59	35.58	62.83	51.27	-	31.82	1.98	14.92
24	1	9	10	510	5480	12.06	0.00	27.74	72.26	4.34	-	65.42	0.72	29.52
25	1	15.92	40.77	3317	8600	51.74	2.45	24.20	73.34	-	-	79.94	4.53	15.53
26	1	12.1	47.6	5635	6944	38.13	1.01	18.46	80.53	-	-	82.74	1.22	16.03
27	1	9.66	38.71	2680	5633	56.59	1.48	15.62	82.90	-	-	97.64	1.28	1.08

Table A2

List of S-HYPE variables used in the study to evaluate the possibility of describing the risk of rutting during logging operations. Three soil layers (upper, middle, and lower), with varying thickness for each soil type, are used in the model.

Variable Name	Scale	Definition
temp	°C	air temperature, provided in Tobs.txt/Tobs_nnn.txt
ctmp	°C	corrected air temperature
snow	mm	snow water equivalent
sdep	cm	snow depth
soim	mm	computed soil moisture (including standing water)
som2	mm	soil water of upper two soil layers (including standing water)
sml1	mm	soil moisture upper soil layer (not including standing water)
sml2	mm	soil moisture second soil layer
sml3	mm	soil moisture third soil layer
smrz	mm	soil moisture root zone (upper two soil layers) (not including standing water)
sm13	mm	soil moisture all soil layers (not including standing water)
stsw	mm	standing soil water
srff	-	soil moisture root zone (upper two soil layers) (not including standing water) as fraction of wcfc volume
smfd	-	soil moisture (not including standing water) as fraction of soil depth
srfd	-	soil moisture root zone (upper two soil layers) (not including standing water) as fraction of root depth
smfp	-	soil moisture (not including standing water) as fraction of pore volume
srfp	-	soil moisture root zone (upper two soil layers) (not including standing water) as fraction of pore volume
smdf	mm	soil moisture deficit to field capacity of upper two soil layers
gwat	m	groundwater level
sfst	cm	frost depth
stmp	°C	soil temperature
stm1	°C	upper soil layer temperature
stm2	°C	middle soil layer temperature
stm3	°C	lowest soil layer temperature
cout	m3/s	simulated outflow from lake/subcatchment
prec	mm/[period]	precipitation as provided in Pobs.txt
cprc	mm/[period]	corrected precipitation
crun	mm/[period]	calculated local runoff from land area. Note that this is not the same as the flow to the local stream if floodplains are used.
cros	mm/[period]	simulated surface runoff. Note that this is not the same as the flow to the local stream if floodplains are used
psim	mm/[period]	precipitation including water that will be removed as “interception losses”
temp	°C	air temperature, provided in Tobs.txt/Tobs_nnn.txt
ctmp	°C	corrected air temperature
snow	mm	snow water equivalent
sdep	cm	snow depth
soim	mm	computed soil moisture (including standing water)
som2	mm	soil water of upper two soil layers (including standing water)

References

1. Lundmark, T.; Bergh, J.; Hofer, P.; Lundström, A.; Nordin, A.; Poudel, B.; Sathre, R.; Taverna, R.; Werner, F. Potential Roles of Swedish Forestry in the Context of Climate Change Mitigation. Forests; 2014; 5, pp. 557-578. [DOI: https://dx.doi.org/10.3390/f5040557]

2. Högberg, P.; Ceder, L.A.; Astrup, R.; Binkley, D.; Bright, R.; Egnell, L.D.G.; Filipchuk, A.; Genet, H.; Ilintsev, A.; Kurz, W.A. et al. Sustainable Boreal Forest Management—Challenges and Opportunities for Climate Change Mitigation; Swedish Forest Agency: Sweden, 2021; 60.

3. Anon. Riksskogstaxeringen_Officiell Statistik Om Den Svenska Skogen. Available online: https://skogsstatistik.slu.se/pxweb/sv/OvrStat/OvrStat__Avverkning/AVV_arlig_avverkning_landsdelar_tab.px/ (accessed on 15 November 2020).

4. Eliasson, L. Effects of Forwarder Tyre Pressure on Rut Formation and Soil Compaction. Silva Fenn.; 2005; 39, pp. 549-557. [DOI: https://dx.doi.org/10.14214/sf.366]

5. Hansson, L.J.; Koestel, J.; Ring, E.; Gärdenäs, A.I. Impacts of off-road traffic on soil physical properties of forest clear-cuts: X-ray and laboratory analysis. Scand. J. For. Res.; 2017; 33, pp. 166-177. [DOI: https://dx.doi.org/10.1080/02827581.2017.1339121]

6. Cambi, M.; Certini, G.; Neri, F.; Marchi, E. The impact of heavy traffic on forest soils: A review. For. Ecol. Manag.; 2015; 338, pp. 124-138. [DOI: https://dx.doi.org/10.1016/j.foreco.2014.11.022]

7. Kozlowski, T.T. Soil Compaction and Growth of Woody Plants. Scand. J. For. Res.; 1999; 14, pp. 596-619. [DOI: https://dx.doi.org/10.1080/02827589908540825]

8. Wästerlund, I. Growth reduction of trees near strip roads resulting from soil compaction and damaged roots—A literature survey. Sver. Skogsvårdsförbunds Tidskr.; 1983; 81, pp. 97-109.

9. Mariotti, B.; Hoshika, Y.; Cambi, M.; Marra, E.; Feng, Z.; Paoletti, E.; Marchi, E. Vehicle-induced compaction of forest soil affects plant morphological and physiological attributes: A meta-analysis. For. Ecol. Manag.; 2020; 462, 118004. [DOI: https://dx.doi.org/10.1016/j.foreco.2020.118004]

10. Eklöf, K.; Lidskog, R.; Bishop, K. Managing Swedish forestry’s impact on mercury in fish: Defining the impact and mitigation measures. Ambio; 2016; 45, (Suppl. 2), pp. 163-174. [DOI: https://dx.doi.org/10.1007/s13280-015-0752-7]

11. Horn, R.; Vossbrink, J.; Peth, S.; Becker, S. Impact of modern forest vehicles on soil physical properties. For. Ecol. Manag.; 2007; 248, pp. 56-63. [DOI: https://dx.doi.org/10.1016/j.foreco.2007.02.037]

12. Håkansson, I. Packning av åkermark vid maskindrift. Omfattning-Effekter-Motåtgärder. Machinary induced compaction of arable soils- Incidence, consequence and counter measures. Soil Science; Sveriges lantbruks universitet (SLU): Uppsala, Sweden, 2000; 123.

13. Sakai, H.; Nordfjell, T.; Suadicani, K.; Talbot, B.; Bøllehuus, E. Soil Compaction on Forest Soils from Different Kinds of Tires and Tracks and Possiblity of Accurate Estimate.pdf. Croat. J. For. Eng.; 2008; 29, pp. 15-27.

14. Toivio, J.; Helmisaari, H.-S.; Palviainen, M.; Lindeman, H.; Ala-Ilomäki, J.; Sirén, M.; Uusitalo, J. Impacts of timber forwarding on physical properties of forest soils in southern Finland. For. Ecol. Manag.; 2017; 405, pp. 22-30. [DOI: https://dx.doi.org/10.1016/j.foreco.2017.09.022]

15. Uusitalo, J.; Ala-Ilomäki, J. The significance of above-ground biomass, moisture content and mechanical properties of peat layer on the bearing capacity of ditched pine bogs. Silva Fenn.; 2013; 47, [DOI: https://dx.doi.org/10.14214/sf.993]

16. Sirén, M.; Salmivaara, A.; Ala-Ilomäki, J.; Launiainen, S.; Lindeman, H.; Uusitalo, J.; Sutinen, R.; Hänninen, P. Predicting forwarder rut formation on fine-grained mineral soils. Scand. J. For. Res.; 2019; 34, pp. 145-154. [DOI: https://dx.doi.org/10.1080/02827581.2018.1562567]

17. Uusitalo, J.; Ala-Ilomäki, J.; Lindeman, H.; Toivio, J.; Siren, M. Modelling soil moisture—Soil strength relationship of fine-grained upland forest soils. Silva Fenn.; 2019; 53, [DOI: https://dx.doi.org/10.14214/sf.10050]

18. Eliasson, L.; Wästerlund, I. Effects of slash reinforcement of strip roads on rutting and soil compaction on a moist fine-grained soil. For. Ecol. Manag.; 2007; 252, pp. 118-123. [DOI: https://dx.doi.org/10.1016/j.foreco.2007.06.037]

19. Hoffmann, S.; Schönauer, M.; Heppelmann, J.; Asikainen, A.; Cacot, E.; Eberhard, B.; Hasenauer, H.; Ivanovs, J.; Jaeger, D.; Lazdins, A. et al. Trafficability Prediction Using Depth-to-Water Maps: The Status of Application in Northern and Central European Forestry. Curr. For. Rep.; 2022; 8, pp. 55-71. [DOI: https://dx.doi.org/10.1007/s40725-021-00153-8]

20. Svensson, S.A.; Braide, A. Tekniska Skogsdata [Träd, Bestånd, Mark, Ståndort; Inkl. Trädfunktioner]; Sveriges Lantbruksuniverisitet: Umeå, Sweden, 1987.

21. Berg, S. Terrängtypsschemat för Skogsarbete; Skogforsk: Gävle, Sweden, 1995; ISBN 91-7614-035-0

22. Anon. Klimatindikator—Antal Dagar Med Snötäcke. Available online: https://www.smhi.se/klimat/klimatet-da-och-nu/klimatindikatorer/klimatindikator-antal-dagar-med-snotacke-1.91081 (accessed on 1 December 2020).

23. Anon. Klimatindikator_Temperatur. Available online: https://www.smhi.se/klimat/klimatet-da-och-nu/klimatindikatorer/klimatindikator-temperatur-1.2430 (accessed on 1 December 2020).

24. Lehtonen, I.; Venäläinen, A.; Kämäräinen, M.; Asikainen, A.; Laitila, J.; Anttila, P.; Peltola, H. Projected decrease in wintertime bearing capacity on different forest and soil types in Finland under a warming climate. Hydrol. Earth Syst. Sc.; 2019; 23, pp. 1611-1631. [DOI: https://dx.doi.org/10.5194/hess-23-1611-2019]

25. Murphy, P.N.C.; Ogilvie, J.; Castunguay, M.; Connors, T.; Meng, F.R.; Arp, P.A. DEM-derived flow channel and wet area mapping: A new tool for forest operations planning. Proceedings of the Sustainable Forest Management Network, Fourth International Conference; Edmonton, AB, Canada, 20–22 June 2006.

26. Ågren, A.; Lidberg, W.; Ring, E. Mapping Temporal Dynamics in a Forest Stream Network—Implications for Riparian Forest Management. Forests; 2015; 6, pp. 2982-3001. [DOI: https://dx.doi.org/10.3390/f6092982]

27. White, B.; Ogilvie, J.; Campbell, D.M.H.M.H.; Hiltz, D.; Gauthier, B.; Chisholm, H.K.H.; Wen, H.K.; Murphy, P.N.C.N.C.; Arp, P.A.A. Using the Cartographic Depth-to-Water Index to Locate Small Streams and Associated Wet Areas across Landscapes. Can. Water Resour. J. Rev. Can. Des Ressour. Hydr.; 2013; 37, pp. 333-347. [DOI: https://dx.doi.org/10.4296/cwrj2011-909]

28. Murphy, P.N.C.; Ogilvie, J.; Arp, P. Topographic modelling of soil moisture conditions: A comparison and verification of two models. Eur. J. Soil Sci.; 2009; 60, pp. 94-109. [DOI: https://dx.doi.org/10.1111/j.1365-2389.2008.01094.x]

29. Mohtashami, S.; Eliasson, L.; Hansson, L.; Willén, E.; Thierfelder, T.; Nordfjell, T. Evaluating the effect of DEM resolution on performance of cartographic depth-to-water maps, for planning logging operations. Int. J. Appl. Earth Obs. Geoinf.; 2022; 108, 102728. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102728]

30. Reeves, D.A.; Reeves, M.C.; Abbott, A.M.; Page-Dumroese, D.S.; Coleman, M.D. A detrimental soil disturbance prediction model for ground-based timber harvesting. Can. J. For. Res.; 2012; 42, pp. 821-830. [DOI: https://dx.doi.org/10.1139/x2012-034]

31. Kankare, V.; Luoma, V.; Saarinen, N.; Peuhkurinen, J.; Holopainen, M.; Vastaranta, M. Assessing feasibility of the forest trafficability map for avoiding rutting—A case study. Silva Fenn.; 2019; 53, [DOI: https://dx.doi.org/10.14214/sf.10197]

32. Jones, M.-F.; Arp, P.A. Relating Cone Penetration and Rutting Resistance to Variations in Forest Soil Properties and Daily Moisture Fluctuations. Open J. Soil Sci.; 2017; 7, pp. 149-171. [DOI: https://dx.doi.org/10.4236/ojss.2017.77012]

33. Jones, M.-F.; Arp, P. Analyzing and Projecting Soil Moisture and Cone Penetrability Variations in Forest Soils. Open J. For.; 2019; 9, pp. 109-142. [DOI: https://dx.doi.org/10.4236/ojf.2019.92005]

34. Salmivaara, A.; Launiainen, S.; Perttunen, J.; Nevalainen, P.; Pohjankukka, J.; Ala-Ilomäki, J.; Sirén, M.; Laurén, A.; Tuominen, S.; Uusitalo, J. et al. Towards dynamic forest trafficability prediction using open spatial data, hydrological modelling and sensor technology. For. Int. J. For. Res.; 2020; 93, pp. 662-674. [DOI: https://dx.doi.org/10.1093/forestry/cpaa010]

35. Friberg, G.; Bergkvist, I. Så Påverkar Arbetsrutiner Och Markfuktighetskartor Körskador I Skogsbruket [How Operational Procedures and Depth-To-Water Maps Can Reduce Damage on Soil And Water And Rutting in the Swedish Forestry]; 904-2016 Forestry Research Institute of Sweden: Uppsala, Sweden, 2016; pp. 1-36.

36. Lindström, G.; Pers, C.; Rosberg, J.; Strömqvist, J.; Arheimer, B. Development and testing of the HYPE (Hydrological Predictions for the Environment) water quality model for different spatial scales. Hydrol. Res.; 2010; 41, pp. 295-319. [DOI: https://dx.doi.org/10.2166/nh.2010.007]

37. Strömqvist, J.; Arheimer, B.; Dahné, J.; Donnelly, C.; Lindström, G. Water and nutrient predictions in ungauged basins: Set-up and evaluation of a model at the national scale. Hydrol. Sci. J.; 2012; 57, pp. 229-247. [DOI: https://dx.doi.org/10.1080/02626667.2011.637497]

38. Cressie, N. Statistics for Spatial Data; 2nd ed. John Wiley & Sons, Incorporated: Hoboken, NJ, USA, 1993.

39. McCullagh, P.; Nelder, J.A. Generallized Linear Models; 2nd ed. CRC Press LLC: Boca Raton, FL, USA, 1989; Volume 37, 532.

40. Haenlein, M.; Kaplan, A.M. A Beginner’s Guide to Partial Least Squares Analysis. Underst. Stat.; 2004; 3, pp. 283-297. [DOI: https://dx.doi.org/10.1207/s15328031us0304_4]

41. Marra, E.; Cambi, M.; Fernandez-Lacruz, R.; Giannetti, F.; Marchi, E.; Nordfjell, T. Photogrammetric estimation of wheel rut dimensions and soil compaction after increasing numbers of forwarder passes. Scand. J. For. Res.; 2018; 33, pp. 613-620. [DOI: https://dx.doi.org/10.1080/02827581.2018.1427789]

42. Labelle, E.R.; Poltorak, B.J.; Jaeger, D. The role of brush mats in mitigating machine-induced soil disturbances: An assessment using absolute and relative soil bulk density and penetration resistance. Can. J. For. Res.; 2019; 49, pp. 164-178. [DOI: https://dx.doi.org/10.1139/cjfr-2018-0324]

43. Montanari, A. Uncertainty of hydrological predictions. Treatise on Water Science; University of Bologna: Bologna, Italy, 2011; Volume 2, pp. 459-478.

44. Lindström, G. Lake water levels for calibration of the S-HYPE model. Hydrol. Res.; 2016; 47, pp. 672-682. [DOI: https://dx.doi.org/10.2166/nh.2016.019]

45. Arheimer, B.; Dahné, J.; Lindström, G.; Marklund, L.; Strömqvist, J. Multi-variable evaluation of an integrated model system covering Sweden (S-HYPE). Proceedings of the Conceptual and Modelling Studies of Integrated Groundwater, Surface Water, and Ecological Systems; Melbourne, Australia, 28 June–7 July 2011; pp. 145-150.

46. Tyystjärvi, V. Soil Moisture in Process Based Modelling; University of Helsinki: Helsinki, Finland, 2019.

47. Roeckner, E.; Bäuml, G.; Bonaventura, L.; Brokopf, R.; Giorgetta, M.E.M.; Hagemann, S.; Kirchner, I.; Kornblueh, L.; Manzini, E.; Rhodin, A. et al. The Atmospheric General Circulation Model Echman5, Part 1, Model Description; Max Planc Institute for Meteorology: Hamburg, Germany, 2003; 349.

48. Launiainen, S.; Guan, M.; Salmivaara, A.; Kieloaho, A.-J. Modeling boreal forest evapotranspiration and water balance at stand and catchment scales: A spatial approach. Hydrol. Earth Syst. Sc.; 2019; 23, pp. 3457-3480. [DOI: https://dx.doi.org/10.5194/hess-23-3457-2019]

49. Maclean, I.M.D. Ecohydrotools. Available online: https://github.com/ilyamaclean/ecohydrotools (accessed on 11 August 2020).