Next Article in Journal
The Role of Actors in Social Innovation in Rural Areas
Next Article in Special Issue
How Eco-Efficiency Is the Forestry Ecological Restoration Program? The Case of the Sloping Land Conversion Program in the Loess Plateau, China
Previous Article in Journal
Evaluating the Impact of Forest Tenure Reform on Farmers’ Investment in Public Welfare Forest Areas: A Case Study of Gansu Province, China
Previous Article in Special Issue
Tourist Agroforestry Landscape from the Perception of Local Communities: A Case Study of Rwenzori, Uganda
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Depth-to-Water Maps to Identify Soil Areas That Are Potentially Sensitive to Logging Disturbance: Initial Evaluations in the Mediterranean Forest Context

1
Consiglio per la Ricerca in Agricoltura e l’Analisi dell’Economia Agraria (CREA), Centro di Ricerca Ingegneria e Trasformazioni Agroalimentari, Via della Pascolare 16, 00015 Monterotondo, Italy
2
Department of Agricultural and Forest Sciences, University of Tuscia, Via San Camillo de Lellis, 01100 Viterbo, Italy
*
Author to whom correspondence should be addressed.
Submission received: 30 March 2022 / Revised: 26 April 2022 / Accepted: 5 May 2022 / Published: 9 May 2022
(This article belongs to the Special Issue Advances in Sustainable Forest Management)

Abstract

:
Scientific research on reduced-impact logging has been addressed to develop effective approaches and methodologies to limit soil disturbance caused by forest operations. In recent years, the development of soil trafficability maps based on soil wetness indices is the approach that has been extensively used in the context of the Boreal forests. In particular, the depth-to-water (DTW) index has been identified as an interesting solution for the identification of areas particularly sensitive to soil disturbance. This study aimed to evaluate the cost-benefit factor of DTW maps for the identification of soil-sensitive areas in the Mediterranean context. In particular, a DTW map was developed for two oak coppice areas located in Italy and harvested over a period of 2–4 years with different mechanisation levels. Soil surveys concerning soil moisture, physico-chemical properties (bulk density, penetration resistance, shear resistance, organic matter), and biological properties (soil microarthropods community measure via soil biological quality (QBS-ar) index) were carried out in these forests, checking for significant differences between the zones at DTW index ≤1 (which should be more sensitive to soil disturbance) and the other areas of the forest soil. The results obtained revealed the efficiency of a DTW index in potential areas at a higher level of soil moisture. On the other hand, the values of soil physico-chemical properties in the areas at a DTW index ≤1 did not differ significantly from the ones in other zones. However, the values of the QBS-ar index in areas with a low DTW index were significantly lower than the ones in zones at the DTW index >1. Therefore, the obtained findings reveal that the DTW index is a reliable tool to identify and predict which areas are more prone to impact soil biological properties.

Graphical Abstract

1. Introduction

The implementation of sustainable forest management (SFM) is one of the main goals specified in the New European Forest Strategies for 2030 [1]. One of the most important aspects to be taken into account to achieve SFM consists of implementing sustainable forest operations (SFOs). The concept of SFOs involves the development of logging activities in a way that ensures economic efficiency without compromising the health of the environment and the safety of the operators [2]. One of the main environmental targets of the SFOs is to limit disturbance to the forest soil caused by logging activities [3,4,5]. Soil compaction caused by machinery traffic can indeed lead to hydrological concerns such as increased runoff and erosion [6,7]. Furthermore, machinery-induced soil compaction can negatively influence several features of natural regeneration, for instance, decreasing root length and the overall biomass of the seedlings [8].
Therefore, scientific research in the sector of reduced impact logging has been focused on the development of possible solutions to limit soil disturbance related to forest operations [3,6,9]. The suggested best practices converge on the importance of carrying out forest harvesting in optimal climatic conditions or better with the soil in a situation of adequate bearing capacity, which is strictly related to soil moisture. In this framework, the application of the precision forest harvesting (PFH) approach has been showing interesting results [10]. In particular, the application of geographic information system (GIS) and global navigation satellite system (GNSS) are valuable instruments for forest harvesting planning [11,12,13,14]. Picchio et al. demonstrated, for instance, that an optimised strip road pattern developed via GIS can decrease by more than half the soil area affected by the passage of the machineries [15].
One of the most recent and interesting applications of GIS to tackle the issue of reducing soil disturbances caused by logging consists of the development of trafficability maps [16,17,18]. These maps are generally developed starting from a digital terrain model (DTM) and indicate areas that are expected to be particularly sensitive to soil compaction and therefore should be avoided by the operator who is driving the machine [19,20]. For the development of trafficability maps, sensitive areas are identified as the ones with high moisture content [21,22]. In fact, wet soils are more prone to soil compaction [17,23].
Therefore, trafficability maps are generally developed based on soil wetness indices, and among these, one of the most applied in the last years has been the depth-to-water (DTW) index. The depth-to-water (DTW) concept was developed and tested by Murphy et al. [24].
DTW maps have been applied and tested in a Boreal climate for the prediction of perennial streams, wet or water-saturated areas, and sensitive areas for ground-based trafficking [25,26,27,28]. These instruments are well-established in the Boreal climates and are assumed to be good predictors of soil moisture conditions [27,29].
More recently, DTW was applied to predict soil strength and areas prone to higher machinery-induced compaction in European temperate forests; however, the results were not fully satisfactory [29,30].
It is worth mentioning that DTW maps have never been applied in the context of Mediterranean forests. Instead, these instruments could be particularly useful in this zone. Coppicing interventions, which are one of the most common forest management in Italy [31,32], are carried out during the winter when soil moisture is higher. These have a subsequent lower bearing capacity, considering that soils are generally not frozen during the winter months. Therefore, the possibility of developing trafficability maps starting from the DTW index is very interesting for implementing SFOs in Mediterranean forests.
In this framework, the present study aimed to evaluate the efficiency of DTW maps in the prediction of soil-sensitive areas in oak coppice forests. A preliminary application of a DTW model was superimposed in a GIS environment on the skid trail network identified after logging activities. Then, a comparison of soil physico-chemical and biological properties (microarthropods community measured via QBS-ar index) among skid trails located in areas at low DTW index (≤1) and skid trails in areas at higher DTW value (>1) was carried out. The authors thus checked the effectiveness of the DTW index for the further development of soil trafficability maps by testing the following research hypotheses: (i) the DTW index is a clear descriptor of more sensitive soil areas; (ii) the magnitude of machinery-induced disturbance in areas at low DTW value is higher than the one in zones with high DTW value.
There are two new aspects to this study. This represents the first attempt to apply DTW maps in the context of Mediterranean forests. Moreover, this study represents a critical assessment of the efficiency of DTW maps in the prediction of areas more prone to soil disturbance while taking into account the biological aspect as well.

2. Materials and Methods

2.1. Study Area

The study area is located in Castel Giorgio, a municipality in the region of Umbria in Italy (coordinates WGS84UTM33T 249968 E; 4730069 N). It consists of two different parcels belonging to a private owner. Both parcels are turkey oak (Quercus cerris L.) coppices with standards (Figure 1). The intervention consisted of both sub-compartments in a coppicing intervention at a growth age of 19 years, releasing 95 standards per hectare.
The two parcels (Figure 1) were harvested in two harvesting seasons, 2017 and 2019. One parcel (the one reported in red in Figure 1) was harvested with a high mechanisation level via the whole-tree harvesting system (WTS), felling operation was motor-manual with a chainsaw, and extraction was performed with a grapple skidder. The intervention surface was 16.59 ha. The other parcel (in purple in Figure 1) was harvested with a medium mechanisation level via tree length system (TLS), felling operations were done by motor-manual with chainsaw, and extraction was performed with a forestry-fitted farm tractor equipped with a winch. The intervention surface was 20.24 ha.
In both parcels, soil is volcanic with a depth ranging from 10 cm to 100 cm. According to USDA classification, it can be classified as “Typical Haplustepts fine, mixed, mesic”. It consists of brown soil with an acidic sub-acidic reaction. Soil texture is the same in both parcels, and the soil can be classified as silty loam. Topographic features are also very similar in the two parcels, with limited roughness and low slope (prevalent slope 15–20%).
The overall surface of areas at a low DTW index accounted for 13% in the parcel harvested with a high mechanisation level and 11% in the one harvested with a medium level of mechanisation. Of the overall surface at a low DTW index, 6% was directly affected by machinery passage in the parcel harvested with a high mechanisation level, while this percentage was 7% in the ones harvested with a medium mechanisation level (Table 1).

2.2. DTW Map Development

The DTW index indicates the least elevation difference between surface flow channels and nearby landscape areas [33].
DTW values are zero in the zones consisting of surface flow channels. Moving upwards from flow channels in the landscape, DTW values increase, indicating decreased soil moisture away from surface waters. To develop DTW maps, therefore, it is necessary to define two thresholds: (1) flow initiation area (FIA), i.e., a catchment area required to form flow channels; (2) a DTW threshold value to consider soil as wet [34], which is generally set up as 1 [25,26]. The main advantage of DTW in comparison to other wetness indices is the function of simulating different soil moisture scenarios by changing the FIA. Indeed, small FIAs (i.e., 0.25 ha) represent high soil moisture at a site, whereas a larger FIA (i.e., 4 ha) would represent dryer conditions [17].
The DTW map was calculated from a 1 m DTM derived from LiDAR data and provided by the Italian National Geodatabase. All the procedures were implemented using the open-access software Quantum Gis 3.16, following the methodology proposed by Schönauer et al. [29].
The DTM was first processed by removing local depressions [35]. A flow direction tool was used to identify the preferential water path from each cell of the DTM. Flow channels were subsequently extracted from the DTM with a deterministic 8 (D8) flow accumulation tool [36], using an FIA value of 0.25 ha. The 0.25 ha FIA value was established to simulate the soil moisture conditions during the period in which extraction operations were carried out. This is done during the months of December–March, which is the period of time when the condition of the soil is at the highest content of moisture in the Mediterranean context.
Subsequently, the output of the flow accumulation tool was used as a starting point to calculate the shortest vertical difference between each surface grid cell and adjacent flow paths with the application of a lowest-cost function. The output raster is in the DTW index map. This map was reclassified into a binary map, assigning value “1” to all the cells with DTW index ≤1 and “0” to the other cells of the raster.
The reclassified map was uploaded on a handheld Trimble Juno GNSS receiver for the field surveys.

2.3. Experimental Design

The experimental design consisted of comparing the magnitude of disturbance to the forest soil, which occurred in areas predicted as sensitive to soil compaction by the DTW algorithm, as well as in areas that were not considered as sensitive. In both the investigated parcels, disturbed soil in the skid trails was identified within DTW-sensitive areas and non-sensitive ones. To avoid the influence of variability due to different factors, such as the number of machine passages on the skid trail, the sampling was carried out in skid trails crossing both DTW sensitive and non-sensitive areas. Thus, it can be assumed that the number of passages in those skid trails was the same in the two areas and that the presence of differences in compaction is related to different moisture conditions as predicted by the DTW algorithm. The acronyms of the different experimental treatments and the related explanation are reported in Table 2.
In each treatment, soil sampling was carried out, collecting information on the soil physico-chemical and biological parameters described in the following sub-sections. Field surveys were carried out at the beginning of March 2022 during the same period of the year in which logging activities took place. Thus, it is possible to assume that the soil moisture conditions at the moment of the field survey were comparable to the ones during logging interventions.

2.4. Soil Physico-Chemical Properties

Penetration resistance (MPa) and shear resistance (Mg m−2) were assessed with a handheld specific instrument in the first 3–10 cm of soil. The values for the two variables referred to the soil water-holding capacity according to Saxton et al. [37]. Eighteen measurements of both penetration and shear resistance were carried out in each experimental treatment for a total of 90 samplings per feature.
Soil bulk density (g cm−3) was calculated by sampling the soil via a dedicated instrument (30 soil samples in each experimental treatment). Soil samples were then put in plastic bags and sealed for shipping to the laboratory. The samples were then weighed after oven drying at 105 °C to constant weight (dry weight). The ratio between the dry weight and the volume of the cylinder of the corer (100 cm3) gives the measure of bulk density [38]. The same samples were used to assess soil moisture.
The percentage of soil organic matter content was performed by collecting 30 soil samples in each experimental treatment. The same corer as for bulk density sampling was applied. Collected samples were put into plastic bags and taken to the laboratory.
The organic matter content was determined by incineration in a mitten at 400 °C for 4 h [38].

2.5. QBS-ar Index (Soil Biological Quality Based on Microarthropods)

The impact on the biological component of soil was evaluated by applying the QBS-ar Index. This is a qualitative index that assesses the complexity of the microarthropod community in the soil. It is based on the assumption that a higher soil quality implies a higher number of microarthropod groups which are well adapted to the soil environment.
Therefore, soil microarthropods are grouped into different biological forms taking into account their morphological adaptation to soil habitats. Each form presents a score called EMI (ecomorphological index), which ranges from 1 to 20 in proportion to the degree of adaptation [38]. The QBS-ar Index is the sum of the EMI of all the groups within the collected soil sample. To assess the QBS-ar index, 30 soil cores 100 cm2 and 10 cm deep were collected with a specific corer in each treatment. Microarthropods were subsequently extracted from the soil with a Berlese–Tüllgren funnel. The specimens were further collected in a preserving solution (75% ethyl alcohol and 25% glycerol by volume) and identified at different taxonomic levels (class for Myriapoda and order for Insecta, Chelicerata, and Crustacea) via a stereo microscope.

2.6. Statistical Analysis

After a preliminary check of data normality with Shapiro-Wilk test [39] and homoscedasticity via the Levene test [40], the presence of statistically significant differences among the mean values of treatments was investigated via unpaired samples T-test (for soil moisture) [41] and one-way ANOVA (for bulk density, penetration resistance, shear resistance and organic matter) [42] applying HSD Tukey test as post hoc [43]. Data with non-normal distribution, or those which presented heteroscedasticity, were analysed using the non-parametric ANOVA Kruskal-Wallis test [44], applying the Duncan test [45] as a post hoc. Furthermore, non-metric multidimensional scaling (NMDS) was applied to show in a clear way the differences in the average soil parameters for the different experimental treatments. This technique is usually applied to collapse information from multiple dimensions into just a few so that they can be visualised and interpreted. The similarity measure applied was “Gower” in the form of centred axes and convex hulls as a geometric approximation. Statistical analysis was performed with Statsoft Statistica 7.0 (Statsoft, Tulsa, OK, USA) [46] and PAST software [47].

3. Results

The values of soil moisture in the experimental treatments (Table 3) showed the presence of statistically significant differences between areas with low DTW index and other areas in both parcels. In particular, in both sub-compartments, soil moisture in the areas at DTW index ≤1 resulted in being significantly higher.
Concerning soil bulk density, no statistically significant differences were detected among the mean values of the treatments (Figure 2). Therefore, soil compaction in the areas at DTW index ≤1 was not significantly different from the one obtained in the areas at higher scores of DTW index in both of the investigated forest parcels. Interestingly, there is no difference in comparison to the control area. This suggests that from three to four years after the intervention, the soil bulk density in oak coppices is fully recovered to the theoretical pre-harvesting values represented by the control zone.
Regarding penetration resistance and shear resistance, no difference was detected for both variables in the HM parcel. There are no differences between the values in HM-DTW and HM. A slight difference is, however, shown by the values of shear resistance between MM-DTW and MM, with statistically higher values in the zones with DTW index ≤1. On the other hand, no difference between MM-DTW and MM was revealed for penetration resistance. Focusing on the comparison with the control area, values of penetration resistance are fully recovered to the pre-intervention situation in HM, while there are still significant differences between the values in MM parcel and the control area. Values of shear resistance seem to be recovered to the pre-intervention situation in both parcels, except for MM-DTW (Figure 3).
Values of soil organic matter are not fully recovered to pre-intervention values in both parcels, as the control is still significantly higher (Figure 4). No difference in the mean values of organic matter was detected for both the investigated parcels between the zones at DTW index ≤1 and the ones at DTW index >1. Interestingly, in all the treatments affected by harvesting, the values of organic matter are statistically different from the control area. This demonstrates that this parameter was not recovered to pre-intervention values after a period of 3–4 years.
Interestingly, although physico-chemical variables did not show any differences between the zones at DTW index ≤1 and the ones at higher DTW index values, the biological impact seems statistically higher in the zones identified as sensitive to soil disturbance by the DTW algorithm. In fact, in both investigated parcels, the QBS-ar values in zones at DTW index ≤1 are significantly lower than the values obtained in areas at DTW index >1 (Figure 5).
Non-metric multidimensional scaling (nMDs) tests produced a two-dimensional ranking (Figure 6), which provided a significantly greater reduction in statistical stresses (i.e., the disagreement between the 2D figure obtained with nMDS and the predicted values from the regression) than expected by chance (α = 0.05). When considering bulk density (BD), penetration resistance (PR), shear resistance (SR), organic matter (OM), and QBS-ar index, the two axes explained 97.5% of the overall variance (coordinate 1 > 96%; coordinate 2 > 1%). The five variables showed the maximum correlation with the ordination axis 1. The soil condition along axis 2 was dominated mainly by OM and BD (Figure 6). The nMDS for the four disturbed scenarios showed a clear difference among the DTW areas as compared to the others, mostly in MM parcel. A clear difference is shown between HM and MM, with a greater level of disturbance in the latter case.

4. Discussions

This study was the first attempt to evaluate the performance of DTW maps in the Mediterranean context from this perspective. The obtained results are rather controversial but can be useful to address future research.
The DTW index has been confirmed to be a good predictor of soil moisture conditions also in the Mediterranean context. The obtained findings showed that soil moisture in zones at low DTW index had statistically significant higher moisture than the ones located in zones at DTW index >1. Therefore, from these preliminary findings, the first research hypothesis seems to be confirmed. These findings are in line with previous literature. Indeed DTW maps are considered reliable predictors of soil moisture conditions in the Boreal context [26]. In current literature accuracy of prediction is, in fact, reported to be in the range between 51% and 92% [19,24,25,26].
However, DTW maps did not identify soil sites particularly prone to soil compaction. The values of soil bulk density, penetration resistance, and shear resistance did not differ significantly on skid trails located in zones at a low DTW index in comparison to the values obtained in areas at higher DTW values. Therefore, soil compaction in zones at DTW index ≤ 1 was not higher than in other areas not identified as sensitive to soil disturbance, leading to the rejection of the second research hypothesis when dealing with soil physico-chemical properties. Additionally, in this case, the findings are in agreement with previous literature on the topic. A sufficient level of reliability in the prediction of soil bearing capacity via DTW maps has not been achieved yet. Indeed, recent studies failed to delineate a relationship between DTW index values and soil strength measures with different terramechanical indices, including penetration resistance [29,48]. Moreover, previous research from Ågren et al. revealed a questionable relationship between soil strength, measured via cone index, and DTW index [27].
However, in contrast with inferential statistics, the results of non-parametric statistics (nMDS) showed that soil impacts in the zones at DTW index ≤ 1 are higher than in the other areas, suggesting the need to carry out further investigation on this topic.
On the other hand, DTW maps have been shown to predict zones in which the biological impact is higher. Specifically, it can be speculated from the obtained results that in the areas at a low DTW index, the same amount of compaction led to higher disturbance to the microarthropod community of the soil. In both the investigated parcels, the values of the QBS-ar index were significantly lower in the areas at a low DTW index (HM-DTW and MM-DTW) than in the zones at DTW values > 1 (HM and MM). There may be more than one explanation for these findings, and further study is necessary. Assuming that DTW maps effectively predict zones at higher soil moisture, as confirmed by the inferential statistic applied to the obtained soil moisture data, the difference in the QBS-ar index in the soil affected by machinery passage in the areas at DTW index ≤ 1 in comparison to the other areas of the forest parcel can be related to (1) the higher moisture of the areas at a low DTW index triggers higher disturbance to soil microarthropods even if the level of compaction is the same as in other areas of the forest soil; (2) in the areas at a low DTW index there is a lower soil microarthropod biodiversity apart from the influence of logging activities.
Literature seems to indicate rather that the first explanation is probably closer to being correct. Indeed, soil moisture conditions have been identified as a major driver of the soil microarthropod community, with higher soil moisture leading to increased biodiversity of soil fauna [49,50,51]. Additionally, in the case of waterlogged conditions, the soil microarthropod community showed the higher variability in the presence of higher moisture content in the soil [52], but further investigation is needed to address the questions raised by the results of the present study. In each case, the obtained results can have important implications for SFM applications. The efficiency and reliability of DTW maps to identify specific areas as particularly prone to disturbance of the soil microarthropod community should be extensively applied to avoid machinery passages in these zones, limiting the impacts on soil biological properties and enhancing, therefore, the forest operations sustainability.
Indeed, uploading and displaying GIS-developed DTW maps on the onboard computers of modern forest machineries, integrated into mobile GIS applications, can provide the machine driver with site-specific information to select the extraction route, which combines consideration for both benefits for soil conservation and operational efficiency [22]. Moreover, dedicated GIS tools and methodologies can be applied to identify an optimal skid trails or strip roads network, calculated by excluding areas at DTW index ≤ 1 from the passage of the machinery [15]. The obtained results suggest another interesting finding, which is indirect because this evaluation was beyond the scope of the research, but which is however evident. This consists of the low recovery time of soil physical and biological properties after coppicing intervention in Mediterranean conditions. After a three-to-four-year period from logging, soil bulk density, penetration resistance, shear resistance, and the QBS-ar index returned almost entirely to the pre-intervention values (CON). This confirms the high resilience to soil disturbances related to logging of this kind of stand [53]. Differently, soil organic matter was not re-established. Furthermore, recovery in HM parcel seems to be stronger, confirming that the application of machineries specifically developed for forest activities, such as skidders and forwarders, can decrease the magnitude of soil disturbance [9,22].
Being this a preliminary study in the Mediterranean context, further research has to be carried out, repeating the study in different study areas with different soils and forest typologies.
In summary, however, although not able to predict areas particularly sensitive to soil compaction, DTW maps were revealed to be effective in identifying areas with higher soil moisture and, most of all, areas prone to be greatly disturbed from the biological point of view.

5. Conclusions

This study represents the first trial of applying DTW maps to identify soil areas that could be particularly prone to machinery-induced compaction in a Mediterranean context.
The obtained findings partially confirmed the research hypotheses. Indeed, DTW maps performed well in predicting soil areas with higher moisture content. On the other hand, soil physico-chemical properties in the areas identified as particularly sensitive to compaction according to DTW index did not show significant differences with other zones.
However, the most important finding of the study was that the disturbance to the soil microarthropod community in areas at a low DTW index was statistically higher than in the other zones of the parcel. More studies are needed to confirm these findings, which are potentially very important for forest management.

Author Contributions

Conceptualisation, F.L., R.V. and R.P.; methodology, F.L., R.V., D.T. and R.P.; formal analysis, F.L. and R.P.; writing—original draft preparation, F.L., R.V. and D.T.; writing—review and editing, F.L., R.V. and R.P.; supervision, R.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding authors.

Acknowledgments

The authors gratefully acknowledge the Forest Consortium of Amiata, Arcidosso, Grosseto, Italy. This research was in part carried out within the framework of the MIUR (Italian Ministry for Education, University and Research) initiative “Departments of Excellence” (Law 232/2016), WP3, which financed the Department of Agriculture and Forest Science at the University of Tuscia. This research was in part carried out within the project PLANNING CHANGE (CAMBIO PIANO), PSR 2014–2020 Tuscany Region, met. Leader. Meas. 16.2-Gal F.A.R. Maremma.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. European Union. New European Forest Strategy for 2030. Available online: https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:52021DC0572 (accessed on 15 August 2021).
  2. Marchi, E.; Chung, W.; Visser, R.; Abbas, D.; Nordfjell, T.; Mederski, P.S.; McEwan, A.; Brink, M.; Laschi, A. Sustainable Forest Operations (SFO): A new paradigm in a changing world and climate. Sci. Total Environ. 2018, 634, 1385–1397. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Picchio, R.; Mederski, P.S.; Tavankar, F. How and How Much, Do Harvesting Activities Affect Forest Soil, Regeneration and Stands? Curr. For. Rep. 2020, 6, 115–128. [Google Scholar] [CrossRef] [Green Version]
  4. Ampoorter, E.; Van Nevel, L.; De Vos, B.; Hermy, M.; Verheyen, K. Assessing the effects of initial soil characteristics, machine mass and traffic intensity on forest soil compaction. For. Ecol. Manag. 2010, 260, 1664–1676. [Google Scholar] [CrossRef] [Green Version]
  5. Ampoorter, E.; Goris, R.; Cornelis, W.M.; Verheyen, K. Impact of mechanized logging on compaction status of sandy forest soils. For. Ecol. Manag. 2007, 241, 162–174. [Google Scholar] [CrossRef]
  6. Picchio, R.; Jourgholami, M.; Zenner, E.K. Effects of Forest Harvesting on Water and Sediment Yields: A Review Toward Better Mitigation and Rehabilitation Strategies. Curr. For. Rep. 2021, 7, 214–219. [Google Scholar] [CrossRef]
  7. Jourgholami, M.; Khoramizadeh, A.; Venanzi, R.; Latterini, F.; Tavankar, F.; Picchio, R. Evaluation of Leaf Litter Mulching and Incorporation on Skid Trails for the Recovery of Soil Physico-Chemical and Biological Properties of Mixed Broadleaved Forests. Land 2021, 10, 625. [Google Scholar] [CrossRef]
  8. Tavankar, F.; Picchio, R.; Nikooy, M.; Jourgholami, M.; Naghdi, R.; Latterini, F.; Venanzi, R. Soil natural recovery process and Fagus orientalis lipsky seedling growth after timber extraction by wheeled skidder. Land 2021, 10, 113. [Google Scholar] [CrossRef]
  9. Labelle, E.R.; Hansson, L.; Högbom, L.; Jourgholami, M.; Laschi, A. Strategies to Mitigate the Effects of Soil Physical Disturbances Caused by Forest Machinery: A Comprehensive Review. Curr. For. Rep. 2022, 8, 20–37. [Google Scholar] [CrossRef]
  10. Keefe, R.F.; Wempe, A.M.; Becker, R.M.; Zimbelman, E.G.; Nagler, E.S.; Gilbert, S.L.; Caudill, C.C. Positioning Methods and the Use of Location and Activity Data in Forests. Forests 2019, 10, 458. [Google Scholar] [CrossRef] [Green Version]
  11. Bont, L.G.; Fraefel, M.; Frutig, F.; Holm, S.; Ginzler, C.; Fischer, C. Improving forest management by implementing best suitable timber harvesting methods. J. Environ. Manag. 2022, 302, 114099. [Google Scholar] [CrossRef]
  12. Phelps, K.; Hiesl, P.; Hagan, D.; Hotaling Hagan, A. The Harvest Operability Index (HOI): A Decision Support Tool for Mechanized Timber Harvesting in Mountainous Terrain. Forests 2021, 12, 1307. [Google Scholar] [CrossRef]
  13. Palander, T.; Kärhä, K. Utilization of Image, LiDAR and Gamma-Ray Information to Improve Environmental Sustainability of Cut-to-Length Wood Harvesting Operations in Peatlands: A Management Systems Perspective. ISPRS Int. J. Geo-Inf. 2021, 10, 273. [Google Scholar] [CrossRef]
  14. Marčeta, D.; Petković, V.; Ljubojević, D.; Potočnik, I. Harvesting system suitability as decision support in selection cutting forest management in northwest Bosnia and Herzegovina. Croat. J. For. Eng. 2020, 41, 251–265. [Google Scholar] [CrossRef]
  15. Picchio, R.; Latterini, F.; Mederski, P.S.; Tocci, D.; Venanzi, R.; Stefanoni, W.; Pari, L. Applications of GIS-Based Software to Improve the Sustainability of a Forwarding Operation in Central Italy. Sustainability 2020, 12, 5716. [Google Scholar] [CrossRef]
  16. Flisberg, P.; Rönnqvist, M.; Willén, E.; Frisk, M.; Friberg, G. Spatial optimization of ground-based primary extraction routes using the bestway decision support system. Can. J. For. Res. 2021, 51, 675–691. [Google Scholar] [CrossRef]
  17. Jones, M.-F.; Arp, P. Soil Trafficability Forecasting. Open J. For. 2019, 9, 296–322. [Google Scholar] [CrossRef] [Green Version]
  18. 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, 10197. [Google Scholar] [CrossRef]
  19. Lidberg, W.; Nilsson, M.; Ågren, A. Using machine learning to generate high-resolution wet area maps for planning forest management: A study in a boreal forest landscape. Ambio 2020, 49, 475–486. [Google Scholar] [CrossRef] [Green Version]
  20. 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. Forestry 2021, 93, 662–674. [Google Scholar] [CrossRef]
  21. Campbell, D.M.H.; White, B.; Arp, P.A. Modeling and mapping Soil resistance to penetration and rutting using LiDAR-derived digital elevation data. J. Soil Water Conserv. 2013, 68, 460–473. [Google Scholar] [CrossRef] [Green Version]
  22. 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, 55–71. [Google Scholar] [CrossRef]
  23. Tavankar, F.; Picchio, R.; Nikooy, M.; Jourgholami, M.; Latterini, F.; Venanzi, R. Effect of soil moisture on soil compaction during skidding operations in poplar plantation. Int. J. For. Eng. 2021, 32, 128–139. [Google Scholar] [CrossRef]
  24. Murphy, P.N.C.; Ogilvie, J.; Connor, K.; Arp, P.A. Mapping wetlands: A comparison of two different approaches for New Brunswick, Canada. Wetlands 2007, 27, 846–854. [Google Scholar] [CrossRef]
  25. Murphy, P.N.C.; Ogilvie, J.; Meng, F.R.; White, B.; Bhatti, J.S.; Arp, P.A. Modelling and mapping topographic variations in forest soils at high resolution: A case study. Ecol. Model. 2011, 222, 2314–2332. [Google Scholar] [CrossRef]
  26. Ågren, A.M.; Lidberg, W.; Strömgren, M.; Ogilvie, J.; Arp, P.A. Evaluating digital terrain indices for soil wetness mapping—A Swedish case study. Hydrol. Earth Syst. Sci. 2014, 18, 3623–3634. [Google Scholar] [CrossRef] [Green Version]
  27. Ågren, A.M.; Lidberg, W.; Ring, E. Mapping temporal dynamics in a forest stream network-implications for riparian forest management. Forests 2015, 6, 2982–3001. [Google Scholar] [CrossRef] [Green Version]
  28. Mohtashami, S.; Eliasson, L.; Jansson, G.; Sonesson, J. Influence of soil type, cartographic depth-to-water, road reinforcement and traffic intensity on rut formation in logging operations: A survey study in Sweden. Silva Fenn. 2017, 51, 2018. [Google Scholar] [CrossRef] [Green Version]
  29. Schönauer, M.; Väätäinen, K.; Prinz, R.; Lindeman, H.; Pszenny, D.; Jansen, M.; Maack, J.; Talbot, B.; Astrup, R.; Jaeger, D. Spatio-temporal prediction of soil moisture and soil strength by depth-to-water maps. Int. J. Appl. Earth Obs. Geoinf. 2021, 105, 102614. [Google Scholar] [CrossRef]
  30. Schönauer, M.; Prinz, R.; Väätäinen, K.; Astrup, R.; Pszenny, D.; Lindeman, H.; Jaeger, D. Spatio-temporal prediction of soil moisture using soil maps, topographic indices and SMAP retrievals. Int. J. Appl. Earth Obs. Geoinf. 2022, 108, 102730. [Google Scholar] [CrossRef]
  31. Di Marzio, N. An overview of forest cover and management in Italy. Nov. Meh. Sumar. 2020, 41, 63–71. [Google Scholar] [CrossRef]
  32. Camponi, L.; Cardelli, V.; Cocco, S.; Serrani, D.; Salvucci, A.; Cutini, A.; Agnelli, A.; Fabbio, G.; Bertini, G.; Roggero, P.P. Effect of coppice conversion into high forest on soil organic C and nutrients stock in a Turkey oak (Quercus cerris L.) forest in Italy. J. Environ. Manag. 2022, 312, 114935. [Google Scholar] [CrossRef] [PubMed]
  33. Murphy, P.N.C.; Ogilvie, J.; Meng, F.-R.; Arp, P. Stream network modelling using lidar and photogrammetric digital elevation models: A comparison and field verification. Hydrol. Process. 2008, 22, 1747–1754. [Google Scholar] [CrossRef]
  34. 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. [Google Scholar] [CrossRef]
  35. Tarboron, D.G. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resour. Res. 1997, 33, 309–319. [Google Scholar] [CrossRef] [Green Version]
  36. O’Callaghan, J.F.; Mark, D.M. The extraction of drainage networks from digital elevation data. Comput. Vision Graph. Image Process. 1984, 28, 323–344. [Google Scholar] [CrossRef]
  37. Saxton, K.E.; Rawls, W.; Romberger, J.S.; Papendick, R.I. Estimating generalized soil-water characteristics from texture. Soil Sci. Soc. Am. J. 1986, 50, 1031–1036. [Google Scholar] [CrossRef]
  38. Venanzi, R.; Picchio, R.; Piovesan, G. Silvicultural and logging impact on soil characteristics in Chestnut (Castanea sativa Mill.) Mediterranean coppice. Ecol. Eng. 2016, 92, 82–89. [Google Scholar] [CrossRef]
  39. Shapiro, S.S.; Wilk, M.B. An analysis of variance test for normality (complete samples). Biometrika 1965, 52, 591–611. [Google Scholar] [CrossRef]
  40. Glass, G.V. Testing homogeneity of variances. Am. Educ. Res. J. 1966, 3, 187–190. [Google Scholar] [CrossRef]
  41. Pfanzagl, J.; Sheynin, O. Studies in the history of probability and statistics XLIV A forerunner of the t-distribution. Biometrika 1996, 83, 891–898. [Google Scholar] [CrossRef]
  42. Fisher, R.A. The Correlation Between Relatives on the Supposition of Mendelian Inheritance. Philos. Trans. R. Soc. Edinb. 1918, 52, 399–433. [Google Scholar] [CrossRef] [Green Version]
  43. Tukey, J.W. Comparing individual means in the analysis of variance. Biometrics 1949, 5, 99–114. [Google Scholar] [CrossRef] [PubMed]
  44. Kruskal, W.H.; Wallis, W.A. Use of Ranks in One-Criterion Variance Analysis. J. Am. Stat. Assoc. 1952, 47, 583–621. [Google Scholar] [CrossRef]
  45. Duncan, D.B. Multiple range and multiple F tests. Biometrics 1955, 11, 1–42. [Google Scholar] [CrossRef]
  46. STATISTICA, Version 7.0. Data Analysis Software System. StatSoft, Inc.: Tulsa, OK, USA, 2007. Available online: www.statsoft.com(accessed on 10 April 2021).
  47. Hammer, Ø.; Harper, D.A.T.; Ryan, P.D. PAST: Paleontological statistics software package for education and data analysis. Palaeontol. Electron. 2001, 4, 9. [Google Scholar]
  48. Schönauer, M.; Hoffmann, S.; Maack, J.; Jansen, M.; Jaeger, D. Comparison of Selected Terramechanical Test Procedures and Cartographic Indices to Predict Rutting Caused by Machine Traffic during a Cut-to-Length Thinning Operation. Forests 2021, 12, 113. [Google Scholar] [CrossRef]
  49. Sterzyńska, M.; Pižl, V.; Tajovský, K.; Stelmaszczyk, M.; Okruszko, T. Soil Fauna of Peat-Forming Wetlands in a Natural River Floodplain. Wetlands 2015, 35, 815–829. [Google Scholar] [CrossRef]
  50. Bhagawati, S.; Bhattacharyya, B.; Medhi, B.K.; Bhattacharjee, S.; Mishra, H. Diversity and density of Collembola as influenced by soil physico-chemical properties in fallow land ecosystem of Assam, India. J. Environ. Biol. 2020, 41, 1626–1631. [Google Scholar] [CrossRef]
  51. Aupic-Samain, A.; Baldy, V.; Delcourt, N.; Krogh, P.H.; Gauquelin, T.; Fernandez, C.; Santonja, M. Water availability rather than temperature control soil fauna community structure and prey–predator interactions. Funct. Ecol. 2021, 35, 1550–1559. [Google Scholar] [CrossRef]
  52. Jakšová, P.; Ľuptáčik, P.; Miklisová, D.; Horváthová, F.; Hlavatá, H. Oribatida (Acari) communities in arable soils formed under waterlogged conditions: The influence of a soil moisture gradient. Biologia 2020, 75, 243–257. [Google Scholar] [CrossRef]
  53. Venanzi, R.; Picchio, R.; Grigolato, S.; Latterini, F. Soil and forest regeneration after different extraction methods in coppice forests. For. Ecol. Manag. 2019, 454, 117666. [Google Scholar] [CrossRef]
Figure 1. Google satellite image of the study area (date 7 July 2019). The boundaries of the two parcels are reported in purple (MM—medium mechanisation level) and red (HM—high mechanisation level). Control area is reported in green. Blu zones represent the areas in which DTW index is ≤1.
Figure 1. Google satellite image of the study area (date 7 July 2019). The boundaries of the two parcels are reported in purple (MM—medium mechanisation level) and red (HM—high mechanisation level). Control area is reported in green. Blu zones represent the areas in which DTW index is ≤1.
Land 11 00709 g001
Figure 2. Bulk density for the different treatments. Obtained p-value from ANOVA was p = 0.37 (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels).
Figure 2. Bulk density for the different treatments. Obtained p-value from ANOVA was p = 0.37 (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels).
Land 11 00709 g002
Figure 3. Penetration resistance (A) and shear resistance (B) for the different treatments. Obtained p-values from ANOVA were p = 0.00041 for penetration resistance and p = 0.00012 for shear resistance. Different letters indicate homogeneous groups (p < 0.05) according to HSD Tukey test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Figure 3. Penetration resistance (A) and shear resistance (B) for the different treatments. Obtained p-values from ANOVA were p = 0.00041 for penetration resistance and p = 0.00012 for shear resistance. Different letters indicate homogeneous groups (p < 0.05) according to HSD Tukey test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Land 11 00709 g003
Figure 4. Organic matter for the different treatments. Obtained p-value from ANOVA was p = 0.00001. Different letters indicate homogeneous groups (p < 0.05) according to HSD Tukey test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Figure 4. Organic matter for the different treatments. Obtained p-value from ANOVA was p = 0.00001. Different letters indicate homogeneous groups (p < 0.05) according to HSD Tukey test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Land 11 00709 g004
Figure 5. QBS-ar index for the different treatments. Obtained p-value from Kruskal–Wallis test was p = 0.0013. Different letters indicate homogeneous groups (p < 0.05) according to Duncan test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Figure 5. QBS-ar index for the different treatments. Obtained p-value from Kruskal–Wallis test was p = 0.0013. Different letters indicate homogeneous groups (p < 0.05) according to Duncan test (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels). Different lowercase letters indicate different homogeneous groups at p < 0.05 according HSD Tukey Test.
Land 11 00709 g005
Figure 6. Non-metric multidimensional scaling (nMDS) analysis of the main indexes and indicators of the disturbed soil and the control (BD—bulk density; PR—penetration resistance; SR—shear resistance; OM—organic matter; QBS-ar—soil biological quality based on microarthropods soil biological quality) (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels).
Figure 6. Non-metric multidimensional scaling (nMDS) analysis of the main indexes and indicators of the disturbed soil and the control (BD—bulk density; PR—penetration resistance; SR—shear resistance; OM—organic matter; QBS-ar—soil biological quality based on microarthropods soil biological quality) (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level; CON—control area consisting of unharvested oak coppice in the last 20 years located close to the two parcels).
Land 11 00709 g006
Table 1. Details of areas at DTW index ≤1 in the two investigated parcels.
Table 1. Details of areas at DTW index ≤1 in the two investigated parcels.
ParcelOverall Surface (ha)Surface at DTW Index ≤1 (ha)Surface at DTW Index ≤1 Affected by Machinery Passage (ha)
High mechanisation (HM)16.592.090.13
Medium mechanisation (MM)20.242.260.16
Table 2. Acronyms for the experimental treatments.
Table 2. Acronyms for the experimental treatments.
AcronymTreatment
HM-DTWAreas identified as sensitive by DTW algorithm (DTW index ≤1) in the forest parcels harvested with high mechanisation level
HMAreas not identified as sensitive by DTW algorithm (DTW index >1) in the forest parcels harvested with high mechanisation level
MM-DTWAreas identified as sensitive by DTW algorithm (DTW index ≤1) in the forest parcels harvested with medium mechanisation level
MMAreas not identified as sensitive by DTW algorithm (DTW index >1) in the forest parcels harvested with medium mechanisation level
CONControl area consisting of unharvested oak coppice in the last 20 years located close to the two parcels
Table 3. Soil moisture values referred to the different conditions analysed and results of the independent T-test analysis (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level).
Table 3. Soil moisture values referred to the different conditions analysed and results of the independent T-test analysis (HM-DTW—areas identified as sensitive by DTW algorithm and harvested with high mechanisation level; HM—areas not identified as sensitive by DTW algorithm and harvested with high mechanisation level; MM-DTW—areas identified as sensitive by DTW algorithm and harvested with medium mechanisation level; MM—areas not identified as sensitive by DTW algorithm and harvested with medium mechanisation level).
TreatmentMoisture %
(AVG ± SD)
t-Valuep-Value
HM-DTW39.6 ± 2.26.566<0.001
HM35.7 ± 1.9
MM-DTW40.0 ± 3.97.930<0.001
MM32.9 ± 1.9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Latterini, F.; Venanzi, R.; Tocci, D.; Picchio, R. Depth-to-Water Maps to Identify Soil Areas That Are Potentially Sensitive to Logging Disturbance: Initial Evaluations in the Mediterranean Forest Context. Land 2022, 11, 709. https://0-doi-org.brum.beds.ac.uk/10.3390/land11050709

AMA Style

Latterini F, Venanzi R, Tocci D, Picchio R. Depth-to-Water Maps to Identify Soil Areas That Are Potentially Sensitive to Logging Disturbance: Initial Evaluations in the Mediterranean Forest Context. Land. 2022; 11(5):709. https://0-doi-org.brum.beds.ac.uk/10.3390/land11050709

Chicago/Turabian Style

Latterini, Francesco, Rachele Venanzi, Damiano Tocci, and Rodolfo Picchio. 2022. "Depth-to-Water Maps to Identify Soil Areas That Are Potentially Sensitive to Logging Disturbance: Initial Evaluations in the Mediterranean Forest Context" Land 11, no. 5: 709. https://0-doi-org.brum.beds.ac.uk/10.3390/land11050709

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop