Post-Harvest Evaluation of Soil Physical Properties and Natural Regeneration Growth in Steep-Slope Terrains

Abstract

Protection of forest soils during harvesting operations is necessary to reduce damage and accelerate recovery time. The current study aims to evaluate soil physical properties, natural regeneration, and its recovery process in treatments including slope gradient, traffic intensity and skid trail after long periods of time, after ceasing the timber harvesting operations. The most recent skidding operations within each 5 years recovery period were studied for a chronosequence of 20 years. Soil samples were taken in abandoned skid trails and data were recorded on naturally regenerated species and density. The results revealed that most soil disturbances occurred on the slopes >20%, as well as the highest levels of traffic intensity. Bulk density and penetration resistance were still higher than the control area, with a significant difference between them, while total porosity was partially recovered. Twenty years after the skidding operation, soil bulk density and penetration resistance were 13.2% and 23.7% higher than the control area, while total porosity was 9.78% lower the than value of the control area. Seedlings of 50–150 cm and >150 cm in height on skid trails had significantly lower density than those in the control. The number of seedlings per m² was less than the control area in all skid trails and for all height classes. The proportion of seedlings present in low traffic intensity was higher than in medium and high traffic intensities. The findings confirmed that full recovery rates are lengthy, and more time than 20 years is required to fully recover, especially with regards to penetration resistance.

1. Introduction

The use of heavy machinery in timber harvesting operations leads to increase soil disturbance and degradation [1,2,3]. However, profitable and sustainable harvesting operations require modern forestry mechanization. Machine traffic on forest soil leads to increased soil compaction, reduces soil porosity, and condenses the pore connectivity, which ultimately results in an increase in soil bulk density and penetration resistance [4,5,6,7]. Soil bulk density, soil strength and soil porosity are often considered as the most direct quantitative measurements of soil compaction [1,4,8,9,10]. Increases in bulk density or soil strength often coincide with decreasing total aeration, which adversely affects the air and water exchange and the microbial activity in topsoil layers [4,7]. For example, a small increase in soil bulk density and soil strength can cause a considerable reduction in soil porosity, particularly macroporosity volumes [11], and an increase in the proportion of micropore spaces [10]. When soil bulk density increases and the air pores are below 10%, the volume of soil pores reduces, which results in a decrease in the presence of microarthropod communities [12] and the intensity of microbial activity [13] in the soil. These changes have negative effects on the growth in height and diameter of regeneration located near the skid trail, and ultimately lead to a reduction in the fertility of the habitat over a long-term period [4,14,15]. In addition, Zenner et al. [16] reported that a partial recovery in the upper 10 cm layer occurred three years after harvesting operations, however, the density, height, and basal diameter of the early growth of quaking aspen (Populus tremuloides Michx.) were significantly reduced. Likewise, Whalley et al. [17] and Ampoorter et al. [2] found that the increase in penetration resistance can lead to a decrease in rooting growth and may restrict under soil strength to 3 MPa. Soil disturbance and its short-term effects have been widely studied across the world [7,18,19,20,21,22,23], but the information on the long-term impacts of soil disturbance and regeneration rates after this disturbance, especially in steep-slope regions, is scarce. In mountainous regions, due to the reduction of the speed of timber skidding machinery, and the increase in slope gradient, meaning the vibrations of surface soil take more time, this compaction can reach the soil beneath [24,25,26], which certainly affects the recovery time of disturbed soil in these regions. Recovery levels of soil physical properties vary depending on the degree of soil compaction, profile depth, soil type, vegetation and climate condition [27]. The regeneration process often begins with the formation of cracks, fissures and shear fractures, due to alternating shrinking and swelling induced by freeze/thaw and drying/wetting cycles [28,29]. These fragmentation processes are much more effective in the topsoil than in the subsoil [28,29]. Similarly, Ezzati et al. [11] reported that soil bulk density, total porosity and air permeability were significantly affected by soil compaction, compared with undisturbed sites, 20 years after the cessation of harvesting operations in mountainous forests. Meanwhile, Ebeling et al. [29] found that the recovery level of physical properties of soils in compacted forests was partially completed after 20 years at the sites, with high biological activity and high clay content. In some mountainous forests, Picchio et al. [30] showed that soil recovery after reduced impact logging activities can occur within 10 years, and this can be considered fairly fast. More frequently, the recovery of compacted forest soil in the absence of rehabilitation treatments is influenced by climatic conditions, thus the activities of roots and soil fauna are very slow and require a long time to fully recover. As a result, obtaining an accurate measure of soil disturbance for post-harvesting assessment of logging is necessary. Forest soil maintenance is a key factor in sustaining productive forests. On the other hand, the current soil condition is an important factor in the establishment of tree species. Investigating the reaction of different species in the establishment and distribution of regeneration in these conditions can play an important role in determining the nature and function of these species. The main objectives of this study were to determine (1) the effects of slopes and traffic on soil physical properties, (2) the process of change in and recovery of the soil’s physical properties after skidding operations, compared with undisturbed area, (3) the role of soil compaction on regrowth of forest stands after 20 years of harvesting operations in abandoned skid trails. The novelty and innovation of this study concerns the long-term investigation of important soil properties and their changes after skidding with the aim of accelerating and modifying recovery, in order to establish regeneration of skid trails.

2. Materials and Methods

2.1. Study Sites

The current research was conducted in Kheyrud Educational and Research Forest in the Hyrcanian forests in northern Iran. The research station lies between 51°33′12″ E and 51°39′56″ E longitude and 36°32′08″ N and 36°36′45.5″ N latitude. The altitude ranges from 790 m to 1135 m a.s.l. An average rainfall ranges from 1420 mm to 1530 mm, with the heaviest precipitation occurring in summer and fall. The average daily temperature ranges from a few degrees below 0 °C in winter to +25 °C during summer. This area is dominantly covered by natural forests containing native mixed deciduous tree species such as beech (Fagus orientalis Lipsky), hornbeam (Carpinus betulus L.), maple (Acer velutinum Boiss.), Cappadocian maple (Acer cappadocicum C.A.M.), Lime tree (Tilia Begonifolia Stev.), oak (Quercus castanifolia C.A.M.), elm (Ulmus glabra Huds.) and Caucasian alder (Alnus subcordata C.A. Mey.). The silvicultural management system is a mixed, uneven aged system with both single and group selective-cutting regimes, and the average timber volume is 329.25 m³ ha⁻¹. The skidding operations were carried out by a rubber-tired skidder TAF E655 equipped with a winch (Romanian rubber-tired skidder). This machine is an articulated four-wheel-drive vehicle with a weight of 6.8 tons (55% on the front and 45% on the rear axle) and engine power of 65 hp. The skidder was fitted with size 18.4–26 tires inflated to 659 kPa on both front and rear axles, and it had a ground clearance approximately 0.45 m with an overall width of 2.85 m. The soil texture of the study site ranges from silt loam to loamy and is classified as a brown forest soil (Alfisols) and well-drained. More general information of the experiments is given in

Table 1. Location and description of the study sites. Four skid trails were selected for study in nearby compartments, one in each of the time-since-most-recent-logging classes of 5 years, 10 years, 15 years, and 20 years.

Recovery Period (No. of Compartment)	Locations	Skid Trail Length (m)	Monthly Temperature (°C)	Elevation (m)
5 years (C. 118)	51°33′26.51″ E Long 36°32′17.24″ N Lat	850	14	850–900
10 years (C. 214)	51°34′31.24″ E Long 36°34′09.15″ N Lat	1100	15	920–975
15 years (C. 212)	51°34′20.41″ E Long 36°34′15.01″ N Lat	900	15	940–975
20 years (C. 217)	51°34′02.86″ E Long 36°34′24.12″ N Lat	1050	15	955–975

.

2.2. Sampling Strategy

Four abandoned skid trails were chosen with a prevalent downhill skidding direction. The trails ranged from 5 years, 10 years, 15 years and 20 years since their last forest harvesting. In order to have complete details of the experimental design of this retrospective study, it is important to specify that the exact number of passes and the age of the trail are unknown. The detailed information reported was found through official historical harvest reports and by interviews with local experts and forest operators as a proxy for determining the traffic intensity and the age of the trail. The motivation behind this research was to gain a retrospective view of the research conducted by Ezzati et al. [11] in another part of the Hyrcanian forests, which claimed that the recovery of soil physical properties can be lengthy. However, this precedent research did not investigate regeneration in response to soil compaction over time, which was one of the most interesting ideas developed in the current research. Three different skid trail types were selected, based on their distance from the log landings, the number of branches from the main trail, and the expert opinion of the forest manager. Each trail was divided into three traffic intensities (i.e., high (HST), medium (MST), and low (LST) skid trails). High traffic intensity (HST) originated from the log landing or main road and was exposed to a high level of traffic intensity. Medium traffic intensity (MST) branched off from an HST or primary skid trail and was exposed to a medium intensity of traffic. Low traffic intensity (LST) branched off from an MST or secondary skid trail and was exposed to the lowest intensity of harvesting traffic (e.g., the farther part of the skid trail network, originating from the interior of the forest). Each traffic level had to be further away on different slope gradients (i.e., 0–20% or gentle (G) and more than 21% steep (S) slopes of the same skid trail and four classes of age). Thus, research plots were located in three traffic intensity classes and two slope gradient classes, with three replications of each factor combination (n = 18 in each recovery period, Figure 1A). In each treatment (e.g., a combination of slope and traffic), six sample plots 4 × 10 m were measured. Plots 10 m long by 4 m wide were delineated prior to sampling with buffer zones between plots of at least 5 m to avoid interactions. Soil samples were taken at three locations in each plot—the left wheel track (LT), between the tracks (BT), and the right wheel track (RT)—along three randomized lines across the skid trail and perpendicular to the direction of travel, with a 2 m buffer zone between lines to avoid interactions (Figure 1B). Thus, 162 soil samples within 18 sample plots of 40 square meter size were taken in each recovery period for a grand total of 648 soil samples within 72 sample plots. Further, soil samples were taken in nearby undisturbed areas with no skidding impact that were at least 20–30 m away from the skid trail (the size of the average height of the dominant trees in the area) to avoid side effects (n = 54 in each recovery period). In order to study the trend and compare the regeneration in skid trails in undisturbed forests after skidding, for each skid trail and traffic intensity class, a 4 m² plot in slope of 0%–20% and a plot on a slope of >20% was designed to evaluate the regeneration density and height of all tree species combined to determine if there was a treatment effect from the same light intensity (intensity of light was determined with a luxury meter at height of 1 m). Along each plot, at a distance of 20–30 m from the skid trail in undisturbed forests, a plot of the same size was designated (Figure 1C). The treatment effect in this case would be evaluated as any difference in the seedling densities that could potentially be attributed to disturbed or compacted soil. All seedlings were counted, their heights measured, and then they were classified into three groups: less than 50 cm for small seedlings, 50–150 cm for medium seedlings and >150 cm and dbh <10 cm for large seedlings [31]. Due to the effect of seed trees on regeneration, the distance between plots and seedlings is considered to be the same in the skid trail and the undisturbed forests.

2.3. Measurements

Given that soil measurements were affected by soil moisture, a dry day was selected to make certain that the same conditions were created (according to soil moisture measurements) and all measurements were carried out on the same day. Soil sample cores were taken from the top mineral soil (from soil surface down to a depth of 10 cm) using a thin-walled steel cylinder, 100 mm long and 50 mm in diameter, driven horizontally into the soil by a hammer-driven device. After extracting the steel cylinder, soil cores were trimmed flush with the cylinder ends and extruded into a plastic bag for transport to the laboratory. Samples were weighed on the day they were collected, and again after oven drying at 105 °C, until constant mass was reached, to determine water content and bulk density. Eventually, soil bulk density was calculated as Equation (1):

$$\mathrm{BD}=\frac{\mathrm{WD}}{\mathrm{VC}}$$

(1)

where BD is the dry bulk density (g cm⁻³), WD is the weight of the dry soil (g), and VC is the volume of the cylinder (cm³).

Soil texture was also determined using the hydrometric method. The description of the soil size distribution in each trail, by traffic intensity, is presented in

Table 2. Soil texture classes at a depth of 0–100 mm for skid trails after different recovery lengths since the most recent harvest operation on skid trails with different traffic intensities. Un = Undisturbed area, HST = High traffic intensity, MST = Medium traffic intensity, LST = Low traffic intensity. The range of particle size was <0.002, 0.05–0.002 and 2–0.05 mm for clay, silt and sand, respectively.

Sample Site
Recovery Length	HST	MST	LST	Un
5 years	Clay	Clay loam	Clay loam	Clay loam
10 years	Clay	Silty clay	Silty clay loam	Clay loam
15 years	Loam	Loam	Loam	Silt loam
20 years	Clay loam	Silt loam	Silty clay	Loam

. Soil penetration resistance (PR) was measured using an analog hand-held soil penetrometer that was inserted vertically into the soil. Each sample point contained one PR measurement (i.e., the first 100 mm of soil depth). As PR measurements are highly influenced by soil moisture, a dry day was chosen to ensure soils were at dry conditions, and all measurements were performed on the same day to ensure that soil moisture conditions were relatively uniform within all sites. To calculate total porosity, first soil particle density was determined according to the ASTM D854-00 2000 standard, and then it computed as Equation (2):

$$\mathrm{TP}=1-\frac{\frac{{\mathrm{M}}_{\mathrm{s}}}{\mathrm{VC}}}{2.65}\times 100$$

(2)

where TP is the apparent total porosity (%), Ms is the mass of soil (g), 2.65 (g cm⁻³) is the soil particle density, and VC is the volume of the intact soil cores (196.25 cm³).

2.4. Statistical Analyses

We used a factorial experiment based on a complete block design that included a combination of three factors: recovery length (four classes), traffic intensity (four classes, including non-traffic control areas), and slope gradient (two classes); each combination was replicated three times. Generalized linear modeling (GLM, one-way analysis of variance) was applied to relate soil property responses with treatment. The analysis of variance, using a full factorial model with fixed factors, was used to determine the effects of traffic intensity and slope gradient for all soil variables. One-way ANOVA was also used to compare the density and height of seedlings in three traffic intensities and two slope classes. The normality and homogeneity of variance were verified with the Kolmogorov–Smirnov and the Levene tests (α = 0.05). Post hoc comparisons of the treatment group means were performed using Duncan’s multiple range tests with a 95% confidence level. Treatment effects were considered statistically significant when p ≤ 0.05. SPSS (release 17.0; SPSS, Chicago, IL, USA) statistical package was used for analyses. Non-metric multidimensional scaling (nMDS) was applied to compare the three soil param studied, related to traffic intensity and slope class. The results of the analyses were also presented using descriptive statistics.

3. Results

Results showed that the recovery period, traffic intensity and the slope of the skid trail have a significant effect on bulk density, penetration resistance and total porosity. The interaction between the age of the skid trail and traffic intensity was significant on bulk density and penetration resistance, while the interaction between the age of the skid trail and slope was not significant on the bulk density and total porosity at 0.01%, and, for penetration, resistance at 0.01% is significant (

Table 3. Two-way ANOVA results for the effects of traffic intensity (T), slope gradient (S) and skid trail ages or recovery period (A) on soil bulk density (BD), penetration resistance (PR), and total porosity (TP).

	p Value
Source of Variation	BD (g cm−3)	PR (MPa)	TP (%)
A	0.000 **	0.000 **	0.000 **
T	0.000 **	0.000 **	0.000 **
S	0.000 **	0.000 **	0.000 **
A × T	0.990 ns	0.000 **	0.997 ns
A × S	0.003 **	0.000 **	0.004 **

Note: ** significant difference at 0.01% level; ^ns not statistically significant.

).

3.1. The Effect of Traffic Intensity

In all skid trails, by increasing traffic intensity from low to high, soil bulk density and penetration resistance increased, but total porosity decreased (

Table 4. Means ± standard deviation of different physical soil properties on skid trails in different treatments. BD: soil bulk density; PR: penetration resistance; TP: total porosity. LST: skid trail with low traffic intensity; MST: skid trail with medium traffic intensity; HST: skid trail with high traffic intensity, Un: Undisturbed area.

Sample Site/Recovery Length	Traffic Intensity	Soil Physical Properties
		BD (g cm−3)	PR (MPa)	TP (%)
A (5 year)	LST	1.09 ± 0.01 b	2.94 ± 0.06 b	53.45 ± 0.77 b
	MST	1.11 ± 0.02 ab	3.22 ± 0.05 ab	52.66 ± 0.76 b
	HST	1.15 ± 0.01 a	3.74 ± 0.06 a	51.26 ± 0.76 b
	Un	0.91 ± 0.02 c	1.66 ± 0.04 c	61.39 ± 0.75 a
B (10 year)	LST	1.11 ± 0.01 a	2.77 ± 0.06 b	50.35 ± 0.76 b
	MST	1.13 ± 0.01 a	2.87 ± 0.05 b	49.46 ± 0.77 b
	HST	1.17 ± 0.02 a	3.01 ± 0.05 a	47.67 ± 0.76 bc
	Un	0.93 ± 0.01 b	1.51 ± 0.04 c	58.49 ± 0.77 b
C (15 year)	LST	1.01 ± 0.02 a	2.40 ± 0.05 b	53.66 ± 0.77 b
	MST	1.05 ± 0.01 a	2.86 ± 0.06 ab	51.76 ± 0.77 b
	HST	1.08 ± 0.02 a	3.21 ± 0.05 a	50.08 ± 0.76 b
	Un	0.89 ± 0.01 b	1.62 ± 0.04 c	59.05 ± 0.77 a
D (20 year)	LST	1.01 ± 0.01 a	1.89 ± 0.06 a	54.46 ± 0.75 a
	MST	1.04 ± 0.02 a	1.98 ± 0.05 a	53.11 ± 0.76 a
	HST	1.06 ± 0.01 a	1.99 ± 0.06 a	52.01 ± 0.77 ab
	Un	0.93 ± 0.01 b	1.46 ± 0.04 b	57.96 ± 0.77 a

Note: Different letters indicate statistically significant differences within columns for each traffic intensity separately (p < 0.01), based on analysis of variance (GLM).

). The highest and lowest amounts of soil bulk density were obtained from the A skid trail and D skid trail, respectively. In all skid trails, and three traffic intensities, the soil bulk density significantly differed with the control area (Table 4). The highest and lowest amounts of penetration resistance were obtained for the high and low traffic intensities. In all skid trails, penetration resistance has a significant difference in all three traffic intensities within the control area. In A and C skid trails, soil porosity has a significant difference between the control area and all three traffic intensity classes, while, in B and D skid trails, there was no significant difference for soil porosity between the control area and traffic intensity classes. The average soil porosity in different traffic intensities was not significantly different in the all four recovery periods. In B and D skid trails, there was no significant difference between traffic intensities, and between the control area and traffic intensities. In D skid trail (except for in high traffic areas), soil porosity is less than in the control area, but this difference is not significant. These three variables showed the maximum correlation with the ordination axes. The variables, bulk density, and total porosity values illustrate the soil scenario on the weighted scale of axis 1 (Figure 2). The impact arrangement along axis 2 is partially dominated by penetration resistance (Figure 2). The nMDS for the three soil impact scenarios (Figure 2) show a clear positive relationship between the traffic intensities and the impact levels.

In all skid trails, the percentage increase of bulk density and penetration resistance in high traffic intensities was higher than in low and medium traffic intensities (

Table 5. Soil properties’ recovery rate (change percent in different years after the skidding operations) in different traffic intensities.

Sample Site/Recovery Length	Traffic Intensity	Soil Physical Properties
		BD (g cm−3)	PR (MPa)	TP (%)
A (5 years)	LST	19.8	77.1	12.9
	MST	22	94	14.2
	HST	26.4	125.3	16.5
B (10 years)	LST	19.3	83.4	13.9
	MST	21.5	90	15.4
	HST	25.8	99.3	18.5
C (15 years)	LST	13.5	48.1	9.1
	MST	18	76.5	12.3
	HST	21.3	98.1	15.2
D (20 years)	LST	8.6	29.4	6
	MST	11.8	35.6	8.4
	HST	14	36.3	10.3

). The highest and lowest percentages of bulk density increase, in skidding A and D in high and low traffic intensities, were 26.4% and 8.6%, respectively. The highest and lowest percentage of increase in penetration resistance compared to the control area, on A and D skid trails in high and low traffic intensities, were 125.3% and 29.4%, respectively. In Table 5, total porosity reduction was higher in high traffic intensities than in low traffic. The highest and lowest percentages of porosity reduction in B and D skid trail in high and low traffic intensities were 18.5% and 6.0%, respectively (Table 5).

3.2. The Effect of Slope Gradients

Longitudinal slopes cause an increase in soil bulk density, especially on steep slopes, which was one of the factors affecting the soil recovery process. In general, in all skid trails (except D), the bulk density on the slope >20% was higher than on slopes of 0–20% (

Table 6. Means ± standard deviation of different physical soil properties on skid trails in different treatments. BD: soil bulk density; PR: penetration resistance; TP: total porosity. A, B, C, D, and Un represent, respectively, skid trails after a period of 5 years, 10 years, 15 years, 20 years and Undisturbed area.

Sample Site/Recovery Length	Slope (%)	Soil Physical Properties
		BD (g cm−3)	PR (MPa)	TP (%)
A (5 years)	0–20	1.1 ± 0.014 ab	2.915 ± 0.052 b	53.425 ± 0.63 b
	>20	1.145 ± 0.014 a	3.686 ± 0.041 a	51.498 ± 0.61 b
	Un	0.933 ± 0.018 c	1.667 ± 0.064 c	61.391 ± 0.77 a
B (10 years)	0–20	1.104 ± 0.02 a	2.368 ± 0.051 b	50.922 ± 0.64 b
	>20	1.183 ± 0.02 a	3.405 ± 0.052 a	47.402 ± 0.63 c
	Un	0.934 ± 0.018 b	1.511 ± 0.064 c	58.492 ± 0.78 b
C (15 years)	0–20	1.009 ± 0.02 b	2.696 ± 0.052 ab	53.723 ± 0.65 ab
	>20	1.091 ± 0.03 a	2.963 ± 0.043 a	49.951 ± 0.63 b
	Un	0.893 ± 0.018 b	1.626 ± 0.061 c	59.05 ± 0.77 a
D (20 years)	0–20	1.046 ± 0.031 a	1.965 ± 0.052 a	53.004 ± 0.63 a
	>20	1.035 ± 0.014 a	1.951 ± 0.043 a	53.994 ± 0.64 a
	Un	0.933 ± 0.018 b	1.463 ± 0.064 b	57.964 ± 0.77 a

Note: Different letters indicate statistically significant differences within columns for each slope class separately (p < 0.01), based on analysis of variance (GLM).

). The highest and lowest bulk densities were located in the B skid trail, with a slope of >20% and the C skid trail, with a slope of 0–20%. Table 6 shows that penetration resistance in all skid trails and for both slope classes was more than in the control area and has a significant difference with the control area. The highest and lowest penetration resistances were found in the A and D skid trails at slopes of >20%, with values 3.686 and 1.951 MPa, respectively (Table 6). Soil total porosity decreased with an increase in the slope, so that soil porosity decreased relative to the control area on all four skid trails and in two slope classes. In A and B skid trails, soil porosity in both slope classes has a significant difference with the control area. In the C skid trail and in the slope class, 0–20%, soil porosity has no significant difference with the control area, but there was a significant difference between the control area and skid trail in the slope class >20%. Soil porosity in the D skid trail, in both slope classes, is less than the control area, but with no significant difference. The variables, bulk density, and total porosity values illustrate the soil scenario on the weighted scale of axis 1 (Figure 3). The impact arrangement along axis 2 is dominated by penetration resistance (Figure 3). The nMDS for the three soil impact scenarios (Figure 3) show a slight positive relationship between the traffic intensities and the impact levels, however, the impacted areas showed a considerable closeness to the control areas. The three soil impact scenarios (Figure 3) did not show any relationship between the slope gradient and the impact level.

The highest and lowest bulk density increases in comparison with the control area were observed on the B and D skid trails with a slope of >20% (26.6% and 10.9%, respectively), which indicated that the effect of a steep slope increased bulk density, especially in the younger skid trail (

Table 7. Soil properties recovery rate (change percent in different years after the skidding operations) in two slope classes.

Sample Site/Recovery Length	Slope (%)	Soil Physical Properties
		BD (g cm−3)	PR (MPa)	TP (%)
A (5 years)	0–20	20.6	74.9	13
	>20	25.7	121.1	16.1
B (10 years)	0–20	19.5	56.7	12.9
	>20	26.7	125.3	19
C (15 years)	0–20	13	65.8	9
	>20	22.2	82.2	15.4
D (20 years)	0–20	12.1	34.3	8.5
	>20	10.9	33.3	7.9

). The results show that in all skid trails (except D) the penetration resistance was greater in the slope of >20% than slopes of 0–20% (Table 7). The highest and lowest increases in penetration resistance were seen on the A and D skid trails at a slope of >20%, with 121.1% and 33.3%, respectively. The decrease in soil porosity in comparison with the control area in both slope classes was lower in the D skid trail than the other skid trails (Table 7). The highest and lowest porosity decreases were seen in the A and D skid trails, in the slope class >20%, with 19% and 7.9%, respectively, compared with the control area.

3.3. The Effect of Recovery Period

We also analyze soil physical properties with respect to the number of years since the last harvesting operation (

Table 8. Changes in soil bulk density, penetration resistance and total porosity (mean ± standard error) in different years after the skidding operations. A, B, C, D, and Un represent, respectively, skid trails after a period of 5 years, 10 years, 15 years, 20 years and Undisturbed area.

Different Years After the Skidding Operations
Soil Physical Properties	Un	A	B	C	D
BD (g cm−3)	0.91 ± 0.01 c	1.11 ± 0.01 a	1.13 ± 0.01 a	1.04 ± 0.01 ab	1.03 ± 0.01 b
PR (MPa)	1.56 ± 0.03 c	3.21 ± 0.03 a	2.81 ± 0.03 b	2.76 ± 0.03 b	1.93 ± 0.03 bc
TP (%)	59.22 ± 0.38 a	52.93 ± 0.42 ab	49.65 ± 0.42 b	52.21 ± 0.42 ab	53.45 ± 0.42 a

Note: Different letters in a row indicate significant differences among intensities of soil physical properties (p < 0.05), based on Duncan’s multiple range tests.

and Figure 4). The rate of recovery was very slow for young skid trails (i.e., 5 and 10 years since the last harvesting operation), while partial recovery was observed for mid-aged and old trails, with 15 and 20 years since the last harvesting operation. Twenty years after logging operations (skid trails, i.e., 15 and 20 years), soil bulk density was significantly increased compared to control areas (Table 8). The highest and lowest penetration resistances were obtained in the younger skid trail (5 years), with 105.77% more than the control area, and older skid trail (i.e., 20 years), which was 23.72% greater than the control area. Although the soil porosities in different years after skidding operations were less than the value of the control area, there was no significant difference compared to the control area. Results show that soil properties were partially recovered over a 20-year period after skidding operations within different lengths of time, although these recovery values were still less than the values in the control (Undisturbed) area (Figure 4).

3.4. Density and Height of Seedling

The ANOVA showed that the number and heights of seedlings were significantly different according to the age of the skid trail, traffic and slope size, compared to the control area. However, seedling density was not significantly different in comparison to the wheel track in the control area, although the height of seedlings showed a significant difference (

Table 9. Results of ANOVA applied on density and height of seedlings at different treatments.

Source of Variation	Seedling Density		Seedling Height
	F Test	p Value	F Test	p Value
Age of skid trail × control	10.93 **	0.004	8.31 *	0.034
Traffic × control	11.46 **	0.001	6.35 **	0.007
Slope × control	8.96 **	0.002	7.16 *	0.027
Wheel track × control	6.70 ns	0.059	7.60 **	0.002

Note: ** significant difference at 0.01% level; * significant difference at 0.05% level; ^ns not statistically significant.

). The evaluation of seedling densities of the height classes of <50 cm revealed that there were no statistical differences between the control area and the skid trail D (Figure 5). However, there was a statistically significant difference in the seedling height class of 50–150 and >150 cm between the control and all skid trails. The greatest number per m² was in the seedling height class of <50 cm (Figure 5). In addition, the density of the seedlings was evaluated for three traffic and two slope classes in the skid trails (Figure 6). As the height classes increased from <50 cm to >150 cm, the proportion of seedlings in the low traffic intensities was higher than in the medium and high traffic intensities. In all three height classes, the proportion of seedlings occupying a slope of >20% was greater than the gradient of 0–20% (except for low traffic intensities).

4. Discussion

4.1. The Effect of Traffic Intensity

In this study, with the increase in traffic intensity, the values of BD were higher than those of the control treatments. Non-metric multidimensional scaling (nMDS) analysis, applied to the period “5 years”, produced a two-dimensional ordination that provided a significantly greater reduction in statistical stress than expected by chance (α = 0.05). When considering soil bulk density, penetration resistance, and total porosity, the two axes explained 88.9% of the overall variance. These findings are consistent with Zenner et al. [16], Ezzati et al. [11] and Picchio et al. [32]. The higher value of BD in high traffic intensities rather than low traffic intensities could be explained by the high soil moisture content of soil during harvesting operations, more timber from these locations [33], the thickness of litter layer, the mixing of organic matter with the soil and, ultimately, the low activity of soil organisms. BD has a significant difference in the traffic densities of high, medium and low compared with the control area; it was 14%, 11.6% and 8.6% higher than the control area and, in all traffic intensities, was less than the threshold of 15% [34]. PR in high traffic was higher than in medium and low traffic, which is consistent with the research by Zenner et al. [16] and Picchio et al. [32]. PR in all three classes of traffic intensity has a significant difference compared with the control area. The increase in BD at low and high traffic intensities was one of the reasons for this relationship. Where BD was constant, the PR values were still high and tended to increase. A longer time was needed for PR to fully recover. Ground-based skidding machines insert more pressure on surface soil and this pressure would lead to a decrease in the pore volume of the soil. This stress is much more effective at reducing volumes of macroporosity than microporosity. A reduction in the soil porosity is more associated with increasing soil bulk density. Unlike BD and PR, the TP, with a regular process in high traffic intensities, is less than medium and low traffic intensities. Thus, on every skid trail, high traffic has the least TP, which is consistent with the results of Rab [14]. Changes in TP for different traffic intensities indicate that the TP does not have a significant difference with the control area. However, the values were less than the control area after 20 years. Partial recovery of the TP in this study may be associated with the humid climatic conditions of the research area, restoring soil bulk density over a 20-years period, falling leaf litter, and the presence of soil organisms.

4.2. The Effect of Slope Gradients

The current study confirmed that the slope gradient was one of the critical factors affecting the soil recovery process, which could persist for a long period of time [11,14,29]. However, non-metric multidimensional scaling (nMDS) analysis, applied to the period “5 years”, for the three soil impact scenarios (Figure 2) did not show a clearly continuous positive relationship between the slope gradient and the impact level. Twenty years after timber harvesting, the values of BD were still higher in the steeper skid trail than the gentle slope. In the steeper terrains, skidders often slip repeatedly on the surface of the soil and do not remain within the skid trails, which may cause a disturbance of the topsoil layers [35]. Twenty years after the skidding operation, the BD at a slope of >20% and 0–20% was 10.9% and 12.1% higher than the control area, which was less than the threshold value of 15% [34]. The slow recovery rates of BD in steep trails may be associates with the instability of soil surface, low content of soil moisture and poor activity of soil organisms. By increasing the slope gradient of the skid trail, the PR increased, and its highest value occurred at the trail with slope gradient of >20%. PR in all skid trails and in both slope classes was significantly more than the value of the control area. There was a positive relationship between BD and PR during soil compaction. However, relative to the control, the increase in PR was always greater than that of BD. This increase in PR over BD was greatest in the steeper slope gradient and high traffic. The greater increase in PR in this present study, was possibly due to age-hardening [31]. According to Dexter et al. [36] age-hardening can occur from particle rearrangement or cementation, and possibly from both. After soil disturbance, the process of cementation occurs over time, as the broken bonds between particles renew [37]. De Moraes et al. [38] found that age-hardening caused the soil strength to increase in an untilled soil of similar clay content. So, this age-hardening process may have contributed to the greater increase in PR relative to BD [31]. TP in all skid trails was lower in steeper slopes compared to slopes of 0–20%, and there was a significant difference with the control area. TP after 20 years, at slopes of 0–20% and >20% is more than the minimum threshold of 10% [39], and there was no significant difference with the control area. In line with the current study, Von Wilpert and Schäffer [40], in the study of the ecological effects of the compaction and initial recovery of machinery in forest stands between 40 and 30 years old, concluded that at least 10 years are needed to recover soil porosity.

4.3. The Effect of Recovery Period

Soil’s physical properties, in different years after skidding operations, tend to recover, which confirms the results obtained by Ezzati et al. [11], DeArmond et al. [31] and Mohieddinne et al. [41]. Non-metric multidimensional scaling (nMDS) analysis, applied to the period “20 years”, produced a two-dimensional ordination that provided a significantly greater reduction in statistical stress than expected by chance (α = 0.05). When considering soil bulk density, penetration resistance, and total porosity, the two axes explained 87.4% of the overall variance. These three variables showed the maximum correlation with the ordination axes. Twenty years after the skidding operation, BD was 13.2% more than in the control area; this difference was significant, which indicated that BD was not fully recovered. The recovery process begins at the surface of the soil, and then gradually expands to the depth of the soil [40]. In this study, soil samples were taken from surface layers with depth of 10 cm), so BD was not recovered on the skid trail with 20 year old. However, with more recovery time, a full recovery of the soil’s physical properties can be achieved. Twenty years after skidding operations, PR is 23.7% more than in the control area; a significant difference that indicates the lack of recovery for PR. PR tended to decrease gradually from 3.21 MPa to 1.93 MPa for 5 and 20 year old skid trails. So, it is clear that more time would be needed for the recovery of PR than for BD. In the younger skid trails (e.g., A and B), lack of recovery of TP was significant as compared to the control areas. Shoulders and Terry [42] reported that changes in the soil porosity resulting from the preparation of infrastructure had persisted up to six years after harvesting operations. According to the results of this research, TP was 10.8% less than the control area, but their difference was not significant 20 years after skidding operations, which does not seem to be consistent with the findings of Tiarks et al. [43], Venanzi et al. [44] and Ebeling et al. [2]. McNabb and Startsev [45] concluded that significant changes were made in soil permeability, which attributes this increased penetration to the increase in the proportion of macroporosity.

4.4. Density and Height of Seedling

Analysis of the effects of treatments on regeneration showed a greater effect on seedling growth than seedling density. Seedling density showed recovery in all three height classes in different years after skidding, consistent with the other studies [31,46,47]. It appears that seedling establishment in the height class of <50 cm was recovered and no significant difference was observed with the control area in the undisturbed forest. However, the number of seedlings in this height class was less than in the control area. Improvement in the number of seedlings in this height class, and especially in the D skid trail, are likely due to the logging-induced changes of soil scarification [47], as the increasing the entrance of light radiation to the forest floor due to the openness of the cover canopy over as 20-years period after logging operations [48], which resulted in an increase the recruitment of pioneer species compared to the undisturbed forest [31]. In the case of the number of seedlings in the height classes of 50–150 and >150 cm, it appears that, even 20 years after skidding operations, the recovery process will not be complete and there is a significant difference with the control area. The lower number of seedlings in these height classes could be a result of a lack of the time that may be required for seedlings in the understory to grow into seedlings over 150 cm, and, even though there is sometimes an increase in the growth rate of regeneration after logging operations, it is limited to the first few years after operations [15,31,48]. The lack of recovery of soil properties, especially BD and PR, can also have an impact on the regeneration and decrease of growth seedlings, as reported by Picchio et al. [32]. Consistent with the current study, Blouin et al. [49] also stated that the quantity, type and quality of seedlings in skid trails had a significant relationship with changes in soil bulk density. According to Gomez et al. [50], when PR is increased to more than 2 MPa, it may prevent the elongation and rooting growth of tree seedlings. The proportion of seedlings in slopes and traffic differences in height classes of <50 and 50–150 cm showed an approximately uniform distribution. While, in the height class of >150 cm, the proportion of seedlings in low traffic was more than medium and high, but there was no significant difference between two slope classes. Although the establishment of the seedlings in the skid trail was improved compared to undisturbed forests, in the Hyrcanian forests, the establishment of the seedlings is highly influenced by nearby seed trees. However, when seedlings are established, compacted soil can negatively impact their height and diameter growth [31,51,52,53].

5. Conclusions

This study evaluated the recovery of soil physical properties under a combination of treatments, including slope gradient, traffic intensity and skid trail age, in a period of twenty years after the cessation of timber harvesting operations. We found that the slope gradient and traffic intensity were the most effective factors influencing recovery rates. The results revealed that most of the soil disturbances occurred on the slopes above 20%, and at high levels of traffic intensity. Soil bulk density and penetration resistance were still higher than in the control area, with a significant difference as compared to the control area up to twenty years after the cessation of logging operations. Over 20 years after the skidding operation seceded, total porosity was partially recovered, however, some differences were found. The evaluation of seedling density revealed significant differences for the height classes of 50–150 and >150 cm. The number of seedlings per m² was less than in the control area in all skid trails and for all height classes. The proportion of seedlings occupied in low traffic intensities was higher than in medium and high traffic intensities. The recovery of forest soils following ground-based mechanized operations can naturally occur, albeit normally at a slower rate, with the presence of climatic processes and the interaction of roots-soil and activity-soil fauna. Therefore, obtaining an accurate measure of soil disturbance for post-harvesting assessment of skidding operation is necessary. The results of this study can help managers to modify skid trails (which are no longer intended for use) as soon as possible and return them to the forest, in order to rehabilitate the soil’s biological activities and restore regeneration. On the other hand, while the restoration of fixed skid trails in the mid-10-year intervals of non-use of trails is a proper solution, they should be able to facilitate water to penetrate, so as to prevent soil erosion and disrupt the water cycle. Therefore, limiting ground-based skidding operations to fixed trail segments with a gentle slope gradient in medium and high traffic intensities can be an effective way to reduce soil disturbance after timber harvest operations. Also, the evaluation of the performance of each of the common methods of soil compaction correction, including the use of logging residuals, and the creation of diagonal and orthogonal grooves, can be helpful in selecting the most effective method. The importance of studies such as these also lies in the possibility of deploying and updating the guidelines, criteria and indicators for sustainable forest management, as proposed by forest certification schemes and as expected by a conservative management of renewable resources.