Soil Natural Recovery Process and Fagus orientalis Lipsky Seedling Growth after Timber Extraction by Wheeled Skidder

The growth and quality of natural seedlings are important goals of forest management in uneven-aged high stands. In this study, the recovery process of the physical and chemical properties of soil, as well as growth variables of beech seedlings on the skid trails after skidding operations in three time periods (10, 20, and 30 years) were investigated in the Hyrcanian forests of Iran. Results revealed that the soil chemical properties (organic C, total N, and pH) were fully recovered after 20 years, while soil physical properties (bulk density, porosity, and penetration resistance) were not yet fully recovered even after 30 years. The values of growth characteristics (stem and root system) and the quality index of beech seedlings were not statistically different from those of the control area ones after 20 years. According to our findings, the process of recovering soil characteristics after forest operations is long and takes more than 30 years. Considering the effect that soil compaction can have on seedling establishment and growth, proper planning of the forest operation is needed to limit the negative effects of forest operations, which meets the requirements of sustainable forest management. We suggest designing skid trails with a minimum distance of 140 m from each other and with a longitudinal slope of less than 15%, in compliance with the principles of best management


Introduction
In recent decades, the use of powerful and heavy machinery in forest management has increased exponentially [1]. In addition, anthropogenic activities and climate changes, could have caused an increase in the risk of shallow landslides and erosion over the last decade [2][3][4]. Mechanized forest operations can have detrimental effects on soil physiochemical properties and regeneration of forest stands [5,6]. Many studies have in fact made evident that the main impact of logging operations on forests is ranges of soil compaction [7][8][9]. On the other hand, appropriate technical options can be taken to reduce this type of damage to soil and reduced regeneration capacity [10,11].
Previous studies have reported that the recovery process of the machine-induced compacted soil to return to its initial state may take from a few years to several decades [5,6,30,31]. However, soil-root interactions, soil biota activities, and freezing-thawing of water in soil are the main drivers that accelerate the recovery of soil properties under natural conditions [5]. For example, Meyer et al. [5] reported that soil properties were recovered by planting black alder (Alnus glutinosa L.) trees in the skid trails. In the Hyrcanian forest, soil physical properties recovered on skid trails, but 25 years was not long enough for the full recovery [5]. Similarly, Ebeling et al. (2016) found that the recovery of compacted forest soils physical properties was partially completed after 20 years over the sites with high biological activity and high clay content [31]. On the other hand, Von Wilpert and Schäffer (2006) found in a study regarding the ecological effects of compaction found that at least 10 years are needed to recover pre-harvesting soil porosity [32]. Although soil compaction has become a common problem in current forestry, little is known about the potential to promote the structural recovery of compacted forest soils [5].
Severe soil compaction caused by the forestry machinery traffic drastically affected seedling growth in both root and shoot variables [28,33,34]. Previous studies have revealed that root elongation and branching patterns can be negatively influenced by increasing soil bulk density and subsequent decreasing pore spaces [27,[35][36][37]. For example, Jourgholami et al. (2020) demonstrated that the greatest and the lowest of seedling root mass ratio (RMR) and root to shoot ratio (R/S) were in the undisturbed (control) area and skid trails with high traffic intensity [34]. Likewise, Cambi et al. [35] concluded that shoot/root ratio and main root length of Quercus robur L. seedlings decreased as soil bulk density increased. Furthermore, Jourgholami et al. [36] reported that soil compaction had an adverse impact on the size and biomass of Quercus castaneifolia C.A. Mey.seedling.
However, in previous studies, the effect of soil compaction on beech seedling growth has been cross-sectional, but changes in growth over a long period have not been studied on skid trails.
The Hyrcanian forests of Iran have a substantial ecological, environmental, and socioeconomic importance. They were ascribed to the World Heritage List in 2019. Ground based skidding is the most common logging system in these forests and soil compaction is the most critical consequence of this method. As both environmental aspects and engineering aspects are of interest, the research focused on ecological engineering and thus the aims of the present study were: (1) to evaluate natural recovery of soil physico-chemical properties after 10, 20, and 30 years from the timber extraction by wheeled skidders, (2) to determine the effects of soil compaction on beech seedling morphology, growth, and architecture, and (3) to identify the growth variables that are the most responsive to the soil compaction level on the skid trails.

Study Area
The study was conducted in the Hyrcanian forests of Iran ( Figure 1A). The study area was located in the Nav forests (latitude: 37 • 38 34" to 37 • 42 21" N; longitude: 48 • 48 44" to 48 • 52 30" E) in the Guilan province, north of Iran. The elevation in the study area ranged from 1150 to 1550 m a.s.l. The site was oriented towards the north. The mean annual precipitation is approximately 950 mm and the mean annual temperature is 9.1 • C. The soils are Alfisols (brown forest soil) according to the United States Department of Agriculture (USDA) soil taxonomy, and the soil texture varies from clay to clay loam. The bedrock type is silt-stone and limestone, which belong to the upper Jurassic and lower Cretaceous periods. The average soil depth to the bedrock ranged from 60 to 90 cm, which is welldrained. The forest type ( Figure 2) is uneven-aged mixed hardwood high stand, dominated by beech (Fagus orientalis Lipsky) and hornbeam (Carpinus betulus L.). Other accompanying tree species are caucasian alder (Alnus subcordata C.A. Mey.), maple (Acer velutinum Boiss.), Land 2021, 10, 113 3 of 17 cappadocian maple (Acer cappadocicum Gled.), Norway maple (Acer platanoides L.), velvet maple (Acer insigne Boiss.), wych elm (Ulmus glabra Huds.), lime tree (Tilia begonifolia Stev.), caucasian elm (Zelkova carpinifolia Diopp), chestnut-leaved oak (Quercus castaneifolia Gled.), and common ash (Fraxinus excelsior L.).
The canopy cover ranges from 80% to 100%. The tree density ranges from 280 to 340 stem ha −1 , while the standing volume from 218 to 296 m 3 ha −1 , and the basal area of trees is in the range between 16.8 and 22.5 m 2 ha −1 . In these forests, the main silvicultural method is selection cutting. The usually applied harvesting system is cut to length (CTL) by ground-based skidding, in fact the wheeled winch skidder is the most common machinery used for timber extraction. Skid trail density in the study area ranges from 23.4 to 29.5 m ha −1 . Marked trees were felled using chainsaws and logs were extracted using wheeled skidders. The weight of the used skidders ranges from 8.8 to 10.1 Mg (55% on the front and 45% on the rear axle), and their width and length range from 3.5 to 3.8 m and 6.1 to 6.5 m, respectively, with engine power ranging from 125 to 132 kW. They were equipped with a blade for light pushing of obstacles and stacking of logs. The skidders have generally been fitted with size 24.5 to 32 tires inflated to 220 kPa on both front and rear axles. The width of the skid trails was 4 m, average their distance from each other was 140 m, and their density ranged from 24.5 to 29.5 m ha −1 . The canopy cover ranges from 80% to 100%. The tree density ranges from 280 to 340 stem ha −1 , while the standing volume from 218 to 296 m 3 ha −1 , and the basal area of trees is in the range between 16.8 and 22.5 m 2 ha −1 . In these forests, the main silvicultural method is selection cutting. The usually applied harvesting system is cut to length (CTL) by ground-based skidding, in fact the wheeled winch skidder is the most common machinery used for timber extraction. Skid trail density in the study area ranges from 23.4 to 29.5 m ha −1 . Marked trees were felled using chainsaws and logs were extracted using wheeled skidders. The weight of the used skidders ranges from 8.8 to 10.1 Mg (55% on the front and 45% on the rear axle), and their width and length range from 3.5 to 3.8 m and 6.1 to 6.5 m, respectively, with engine power ranging from 125 to 132 kW. They were equipped with a blade for light pushing of obstacles and stacking of logs. The skidders have generally been fitted with size 24.5 to 32 tires inflated to 220 kPa on both front and rear axles. The width of the skid trails was 4 m, average their distance from each other was 140 m, and their density ranged from 24.5 to 29.5 m ha −1 .

Experimental Design and Data Collection
To evaluate the soil recovery process and the growth response of beech seedlings to soil compaction, three age classes of skid trails (10, 20, and 30 years from harvesting) and three replicates for each was identified and selected in the study area ( Figure 1B). In total, tests were performed on nine abandoned skid trails (Table 1). In all the skid trails, the skidding direction was down-slope. The length of skid trails ranged from 1053 to 1270 m and had a wide range of longitudinal slopes (0 to 36%). In each skid trail, two longitudinal slopes of less than 15% (low slope class) and equal to or more than 15% (steep slope class), which had the same traffic intensity (15 passes), were

Experimental Design and Data Collection
To evaluate the soil recovery process and the growth response of beech seedlings to soil compaction, three age classes of skid trails (10, 20, and 30 years from harvesting) and three replicates for each was identified and selected in the study area ( Figure 1B). In total, tests were performed on nine abandoned skid trails (Table 1). In all the skid trails, the skidding direction was down-slope. The length of skid trails ranged from 1053 to 1270 m and had a wide range of longitudinal slopes (0 to 36%). In each skid trail, two longitudinal slopes of less than 15% (low slope class) and equal to or more than 15% (steep slope class), which had the same traffic intensity (15 passes), were identified to analyze the effect of slope on soil recovery and seedling growth. Five sample plots (10 m × 4 m) with random starting points and a regular distance of 5 m intervals were taken on the skid trails for each slope class (10 sample plots for each skid trail). Sample plots were placed to cover the width of the skid trail (4 m) and 10 m along the skid trail ( Figure 1C). For each skid-trail sample plot, one sample (2 m × 2 m) was taken at a 50 m distance from the skid trail in the undisturbed interior stand as control plots (10 control sample plots for each skid trail). In each sample plot, five sampling lines were designed at a distance of 2 m from each other and perpendicular to the skid trail, of these, three lines were randomly selected for sampling. On each skid-trail sample, nine soil core samples (three samples on the left wheel tracks (LWT), three samples on the right wheel tracks (RWT), and three core samples on the between wheel tracks (BWT)), were taken to measure soil bulk density. In addition, one core sample from the center of the control plot was taken to measure soil bulk density. In total, 90 soil core samples were taken from each skid trail (45 samples on slopes < 15%, and 45 samples on slopes ≥ 15%); and 10 soil core samples from the undisturbed (control). Soil parameters were derived using ASTM D854-00 2000 soil laboratory measurement standards. Surface litter and duff were removed before sampling. The soil samples of 10 cm in depth were collected with a soil hammer and rings (diameter 5 cm, length 10 cm), put in polyethylene bags, and immediately labeled. Samples were weighed on the day of collection. The soil samples were dried in an oven at less than 105 • C for 24 h to obtain dry bulk density [38]. The dry bulk density was calculated from the Equation (1): where BD is the dry bulk density (g cm −3 ), WD is the weight of the dry soil (g), and VC is the volume of cylinder (cm 3 ). The total soil porosity was calculated as Equation (2): where TP is the apparent total porosity (%), BD is the bulk density (g cm −3 ), and 2.65 g cm −3 is the soil particle density measured by a pycnometer on the same soil samples used to determine the bulk density [39]. Macroporosity and microporosity were calculated using Equations (3) and (4), respectively.
where MIP is the microporosity (%), BD is the dry bulk density (g cm −3 ), and θm is the water content on a mass basis (%) [15,40].
where MP is the macroporosity (%), TP is the total porosity (%), and MIP is the microporosity (%). On each skid-trail sample, the soil penetration resistance (PR) was measured by 18 points (six points on left wheel track, six points on right wheel track, and six points on log skidded routes, between the two-wheel tracks) using a hand-held penetrometer (Eijkelkamp, Zevenaar, Netherlands) at depths of 0, 5, and 10 cm. In addition, PR was measured by eight points for each control plot. Soil pH was determined using an Orion Analyzer Model 901 pH meter in a 1:2.5, soil/water solution, soil organic C (OC) was determined using the Walkley-Black technique [41], and total N (TN) using a semi-micro-Kjeldahl technique [42]. The hydrometric method was used to determine soil texture. Seedling response to different soil properties was measured by several variables. The traditional measure is height growth. Forms of below-ground growth and biomass are also significant measure growth responses. On each sample plot, two normal seedlings of beech tree (nearest to the soil bulk density samples) with stem lengths between 30 and 60 cm and seven years of age of were selected (20 seedlings in each skid trail). Seedling age was determined by phyllotaxy, and the following parameters were measured for each seedling: Morphological parameters: stem length (SL), stem diameter (SDM), stem height (SH), main root length (MRL), main root diameter (MRD), lateral root length (LRL), and root penetration depth (RPD); growth parameters: total dry biomass (TDB), stem dry biomass (SDB), and root dry biomass (RDB); architectural parameters: ratio of lateral to main root Land 2021, 10, 113 6 of 17 length (RLM), root mass ratio (RMR; ratio of RDB to TDB), stem mass ratio (SMR, ratio of SDB to TDB), ratio of main root length to stem length (RRS), and the ratio of root penetration to main root length (RPL). Stem and root diameter were measured by a vernier caliper (Insize Model 1205, Austria), 5 cm above and below the soil surface, respectively. The vertical distance from the tip of the root to the soil surface was measured by a metal ruler as the root penetration depth. The dry weights (biomass) of seedlings were obtained after drying at 70 • C until a constant weight was reached. Seedling quality index (SQI) was calculated by Equation (5) [43]:

Data Analysis
Data analyses were performed using SPSS version 19 (Chicago, IL, USA) software. The normality of the distribution of the data was checked using the Kolmogorov-Smirnov test (p = 0.05). The homogeneity of variance among treatments was tested and approved with Levene's test (p = 0.05). The two-way analysis of variance (ANOVA) was used to compare recovery of soil properties and seedling characteristics among treatments (SKT10, SKT20, SKT30, and UNCO). Duncan's test was applied to find differences among treatment means at p ≤ 0.05. The independent samples t-test was applied to the comparison means of soil properties and seedling characteristics between slope classes. The relationship between soil properties and seedling characteristics were determined using Pearson correlation. Principal component analysis (PCA) was applied to investigate any possible linear correlations between the recovery levels of soil physical and chemical properties, as well as on morphology and growth of seedlings among treatments (SKT10, SKT20, SKT30 and UNCO, all replicated for two slope conditions LS and HS). The main criteria used were eigenvalue (>1), scree plot (retaining all components within the sharp descent), loading score for each factor (± 0.10), and meaningfulness of each dimension. To minimize the scaling effect due to different measurement units, the data corresponding to each independent variable were standardized using the Box-Cox lambda. All statistical tests were performed by the SPSS software package (release 17.0; SPSS, Chicago, IL, USA).

Soil Recovery
All the soil physico-chemical properties were significantly different among the varying ages of the skid trails ( Table 2). The slope of the skid trail also had a significant effect on the bulk density, total porosity, macro porosity, penetration resistance, and organic C, while it did not have any significant effect on the micro porosity, moisture content, total N, and pH. Interaction of age and slope of skid trail had significant effects only on BD, TP, and MPI. The values of bulk density, total porosity, and micro porosity were not fully recovered after 30 years from skidding and were significantly different from the undisturbed area (Table 3). The values of macro porosity and penetration resistance were fully recovered on low slopes, while they were not fully recovered on steep slopes. The values of moisture content, organic C, total N, and pH were fully recovered. The values of bulk density in the steep slopes were significantly higher than those in the low slopes in the SKT10 and SKT20, while they were not significantly different in the SKT30 and in the undisturbed area. The values of total porosity in the low slopes were significantly higher than those in the steep slopes in the SKT10 and SKT20, while they were not significantly different in the SKT30 and in the undisturbed area. The values of macro porosity in the low slopes were significantly higher than those in the low slopes in the SKT10, while they were not significantly different in the SKT20, SKT30 and in the undisturbed area. The values of micro porosity, moisture content, total N, and pH they were not significantly different in low and steep slope classes. The value of penetration resistance in the steep slope was significantly higher than in the low slope in SKT10. The values of organic C in the low slopes were significantly higher than those in the steep slopes in the SKT10, while they were not significantly different in the SKT20, SKT30 and in the undisturbed area.

Morphology and Growth of Beech Seedlings after Logging
All the properties of seedling morphology and growth were significantly different among the varying ages and slopes of skid trails, but their interaction did not have a significant effect on seedling morphology and growth (Table 4).
All properties of morphology and growth of seedling grown in the SKT10 and SKT20 were significantly different from the undisturbed area, while they were not significantly different between the SKT30 and undisturbed area ( Table 5). All properties of morphology and growth of seedling grown in the SKT30 were significantly different from SKT20, and all properties of morphology and growth of seedling grown in the SKT20 were significantly different from SKT10 (except stem height and stem length in low slopes).  All properties of morphology and growth of seedling grown in the low slopes were significantly higher than those values in the steep slopes in the SKT10. In the SKT20, low slopes had significantly higher values of main root length, lateral root length, root penetration depth, stem dry biomass, and total dry biomass than those values in steep slopes. In the SKT30 and undisturbed area, they were no significant differences between low and steep slopes.

Architecture of Beech Seedlings after Logging
All the architectural characteristics of seedling were significantly different among the ages of the skid trails. The slope of the skid trails had a significant effect only on RLM (ratio of lateral root length to main root length), and the interaction of age and slope of the skid trail had no significant effect on seedling architectural characteristics (Table 6). All the architectural characteristics of seedlings grown on the SKT10 and SKT20 were significantly different from the undisturbed area (except root mass ratio (RMR), in the low slope class) ( Figure 3). All of the architectural characteristics of seedlings grown on the SKT30 were not significantly different from the undisturbed area. The value of RLM (ratio of lateral root length to main root length) decreased by the increasing age of the skid trail, and the undisturbed area had the lowest value of RLM. The value of RLM in SKT10 and SKT20 in steep slopes was significantly higher than low slopes, while there were no significant differences between low and steep slopes in SKT30 and the undisturbed area. Values of RMR were not significantly different between age classes or slope classes of the skid trails. The value of stem mass ratio (SMR) increased by increasing age of the skid trail, and the undisturbed area had the highest value of SMR. The value of main root/stem length (RRS) increased by increasing age of the skid trail, and the undisturbed area had the highest value of RRS. The value of ratio of root penetration depth to main root length (RPL) decreased by the increasing age of skid trail, and the undisturbed area had the lowest value of ratio of main root length to stem length (RRS).

Seedling Quality Index (SQI)
Results of ANOVA revealed that the age class of the skid trail had a significant effect on the mean values of SQI in the both the low slope (F = 1021.119, p < 0.001) and the steep In SKT10 and SKT20, the values of ratio of lateral root length to main root length (RLM) in the steep slope were significantly higher than those values in the low slope. In SKT20, the value of RRS in the low slope was significantly higher than in the steep slope.

Seedling Quality Index (SQI)
Results of ANOVA revealed that the age class of the skid trail had a significant effect on the mean values of SQI in the both the low slope (F = 1021.119, p < 0.001) and the steep slope (F = 2293.410, p < 0.001) classes (Figure 4). The mean of SQI increased by the increasing age of the skid trail in both slope classes of skid trails. The highest mean value of SQI was in the UNCO (75.22 in LS, and 73.66 in HS), and the lowest was in the SKT10 (36.50 in LS, and 24.59 in HS). The results of Duncan's test indicated that there were significant differences in the mean values of SQI between SkT10 and UNCO, between SKT10 and SKT20 in both low and steep slope classes, while differences between SKT30 and UNCO was significant only in the steep slope class.

Relationship between Characteristics of Beech Seedlings and Soil Properties
Results (Table 7) showed that SH (stem height) was not correlated with soil physicochemical properties, SL (stem length) significantly correlated with MC (moisture content) in steep slope, SDM (stem diameter) significantly correlated with bulk density, main root length, lateral root length. Root dry biomass significantly correlated with all the soil physico-chemical properties, main root diameter significantly correlated with bulk density and moisture content, root penetration depth significantly correlated with all the soil physico-chemical properties (except pH), stem dry biomass significantly correlated with bulk density, and total dry biomass significantly correlated with bulk density, penetration resistance, and moisture content.  Results of the t-test showed that the mean values of SQI in steep slope was significantly lower than the low slope classes in all the age classes of skid trails (t = 15.786, p < 0.001 in SKT10; t = 14.774, p < 0.001 in SKT20; and t = 8.7384, p < 0.001 in SKT30), while there was no significant difference in UNCO (t = 1.927, p = 0.061 in SKT10).

Relationship between Characteristics of Beech Seedlings and Soil Properties
Results (Table 7) showed that SH (stem height) was not correlated with soil physicochemical properties, SL (stem length) significantly correlated with MC (moisture content) in steep slope, SDM (stem diameter) significantly correlated with bulk density, main root length, lateral root length. Root dry biomass significantly correlated with all the soil physico-chemical properties, main root diameter significantly correlated with bulk density and moisture content, root penetration depth significantly correlated with all the soil physico-chemical properties (except pH), stem dry biomass significantly correlated with bulk density, and total dry biomass significantly correlated with bulk density, penetration resistance, and moisture content. Table 7. Pearson correlation between characteristics of beech seedlings and soil properties in two slope classes of skid trails after logging operation. Two principal component analyses (PCAs) were carried out for the soil and seedlings conditions, in which the principal components, PC1 and PC2, explained 95.6% and 2.2% of the total variance, respectively (Figures 5 and 6). The PC1 and PC2 scores for the six treatments (SKT10, SKT20, and SKT30, replicated for two slope conditions LS and HS) and the two controls (UNCO-LS and UNCO-HS) are shown in Figure 6. In general, a positive trend is shown from SKT10 to SKT30 and in detail the slope influence is low (LS and HS). In particular, SKT30 (LS and HS) showed a similar situation and close to UNCO (LS and HS), while SKT10 (LS and HS) and SKT20 (LS and HS) were clearly distinct from UNCO (LS and HS), based on the score plot depicted in Figure 6. treatments (SKT10, SKT20, and SKT30, replicated for two slope conditions LS and HS) and the two controls (UNCO-LS and UNCO-HS) are shown in Figure 6. In general, a positive trend is shown from SKT10 to SKT30 and in detail the slope influence is low (LS and HS). In particular, SKT30 (LS and HS) showed a similar situation and close to UNCO (LS and HS), while SKT10 (LS and HS) and SKT20 (LS and HS) were clearly distinct from UNCO (LS and HS), based on the score plot depicted in Figure 6.

Soil Recovery
Our results revealed that natural recovery of all the soil physico-chemical properties (BD, TP, MP, MIP, PR, MC, OC, TN, and pH) were influenced by the age according to the time since the forest utilization, while only natural recovery of BD, TP, MP, PR, and OC were dependent on skid trail slope. In line with the current results, Venanzi et al. (2016) demonstrated that physical, chemical and biological soil features can be strongly impacted by harvesting operations and, consequently, certain soil processes may be influenced at

Soil Recovery
Our results revealed that natural recovery of all the soil physico-chemical properties (BD, TP, MP, MIP, PR, MC, OC, TN, and pH) were influenced by the age according to the time since the forest utilization, while only natural recovery of BD, TP, MP, PR, and OC were dependent on skid trail slope. In line with the current results, Venanzi et al. (2016) demonstrated that physical, chemical and biological soil features can be strongly impacted by harvesting operations and, consequently, certain soil processes may be influenced at the imposed compaction levels [44]. Our study indicates that none of the soil properties were fully recovered after 20 years. Even after 30 years, the soil BD, TP, and MIP were not fully recovered, while MP and PR in low slope, OC, TN, and pH were fully recovered. Similar to our findings, Ezzati et al. (2012) [45] and Jaafari et al. (2014) [46] showed that soil properties were not fully recovered during a long-term period (20 years) after skidding operations in the Hyrcanian forests of Iran. Moreover, Webb et al. (2012) [47] and Rab (2004) [48] demonstrated that more than 20 years are required for BD to recover in the upper 10 cm of soil. Accordingly, Sohrabi et al. (2020) found that the soil bulk density, the penetration resistance, and the microporosity of a 25-year-old skid trail were 8.4% to 27.4% and 50.4% greater. Total porosity, macroporosity, and soil moisture were 1.9% to 17.1% and 4.6% lower than the undisturbed area [49]. Such results were confirmed also by , who indicated that 20 years after skidding operations, the micro-morphological properties of compacted soil on the skid trails had not yet recovered [33].

Seedling Response
Our results revealed that all the seedling growth characteristics (SH, SL, SDM, MRL, MRD, LRL, RPD, SDB, RDB, and TDB) were dependent on both age and slope classes of skid trails. Furthermore, all the architectural properties of seedlings (RLM, RMR, SMR, RRS, and RPL) were dependent on the age class of the skid trail, while only RLM was dependent on the slope class of skid trail. Consistently with the current study, Jourgholami et al. (2016) studied the effects of soil compaction on seedling morphology, growth, and architecture of chestnut-leaved oak in the Caspian forests of Iran [36]. Their results are in line with our results, showing that seedling characteristics, including size and biomass, were negatively affected by soil compaction. Furthermore, Jourgholami (2018) determined that all morphological responses of Cappadocian maple except lateral root length showed a negative quadratic curve as soil penetration resistance increased [27]. All the average values of seedling growth characteristics increased by the increasing age of the skid trail. Moreover, growth and architectural characteristics (except RMR in low slope) of seedlings grown on SKT10 and SKT20 were significantly different from seedlings grown in UNCO, while seedlings grown on SKT30 were not significantly different from UNCO. In fact, the effect of the age of the skid trail on the reduction of all underground (root system) and above ground (stem) seedling growth characteristics was the same and continued for up to 20 years. The average values of SH, SL, SDM, MRD, and RDB in steep slope classes were significantly lower than those in the low slope classes only in the SKT10; the average values of MRL, LRL, RPD, SDB, and TDB in steep slope classes were significantly lower than those in the low slope classes in SKT10 and SKT20; there were no significant differences between the growth characteristics in high and low slope classes in the SKT30 and UNCO. In fact, the effect of the slope of the skid trail on the reduction of the underground growth profile of the seedlings (root system) lasted longer than its effect on the above ground growth profile of seedlings (stem). Again consistent with the current study, Cambi et al. (2017) showed a higher soil compaction (+27%) reduced the Pedunculated Oak (Quercus robur L.) seedling performance in terms of height (−26%) and root growth (−24%) growth [35]. Similarly, Picchio et al. (2019) monitored and evaluated physical soil properties and their effects on maple and beech seedlings on 10 year old skid trails in the Hyrcanian forests [28]. Their results revealed that, for both seedling species, the morphological parameters, root length (RL) and root penetration depth (RPD) showed the highest decreasing response to soil compaction (decreased RL: 32% in maple, and 27% in beech; decreased RPD: 42% in maple, and 49% in beech); reduction in root growth (28% in maple and 33% in beech) was the main growth response to soil compaction. The negative effects of soil compaction on the morphological attributes can be explained by the physiological status of the seedlings as reported by Cambi et al. (2016) [35]. Our results showed that after 30 years of skidding operations, none of the architectural characteristics of the seedlings were significantly different from the UNCO. In fact, after 30 years of skidding operations, all the architectural features of the seedlings had been restored. The slope of the skid trail had a significant effect only on RLM in SKT10 and SKT20, and on RRS in SKT20. Plant growth is usually negatively affected by soil compaction [50] due to the resistance of the soil to be penetrated by the radicle [21,51,52]. Several studies have shown that soil compaction impacts tree metabolism such as by decreasing photosynthesis, reducing foliage area, and tree growth [52,53]. Furthermore, our results indicated that skidding operations significantly reduced seedling quality index (SQI) in comparison to UNCO 20 years after skidding operations. The values of SQI in the skid trails reached those in the UNCO after 30 years from skidding operations. In compacted soil, seedlings and saplings are more susceptible to drought, as the root system is prevented from quickly reaching deep and moist horizons [52,54,55].

Relationship between the Soil Properties and Seedling Parameters
In this study, three physical characteristics of soil (bulk density, porosity, and penetration resistance) were measured, which are usually considered to assess soil compaction [35], and are good indicators of the response of the plant root system to soil compaction [1,21,28]. Our results indicated that growth characteristics of seedlings were significantly correlated with the soil BD and MC. Root growth characteristics (MRL, MRD, LRL, RPD, and RDB) were significantly correlated with soil PR, OC, TN, and pH. In line with our results, Ponder et al. (2012) showed that compaction reduced the total stand biomass of Aspen trees, but had no effect on other species [56]. From the biplot analysis, in Figure 5, it can be affirmed that the main soil physical parameters were associated with the PC1 scores and soil chemical parameters as well as seedling characteristics, which were closely associated with PC2. Ten years after logging operations in the skid trails, the recovery rates of soil and seedlings characteristics are low ( Figure 6). From the 20th year after logging, the recovery is more evident and reaches a complete recovery in the 30th year after logging operations ( Figure 6). This finding was also confirmed by , who showed that the root length of maple seedlings was inversely related to compaction [33]. Root length decreased consistently with increasing traffic frequency on all trail gradients and with increasing trail gradient at all traffic frequencies, with the greatest reductions at the highest traffic frequency and steeper gradients. In accordance with the current results, Jourgholami et al. (2020) reported that stem and main root length, total, stem and root dry masses, RMR, and R/S were positively correlated with total porosity and macroporosity, and negatively correlated with soil bulk density and penetration resistance [34]. Our results are also in line with Picchio et al. (2019), who showed that different soil physical properties can lead to differences in beech and maple seedling growth between the skid trails and non-skid trails [28]. Furthermore, our results revealed that higher soil bulk density was associated with thinner roots, shorter lateral roots, and lower penetration depth. Similar results were found by Mosena and Dillenburg (2004) [57], Jourgholami (2018) [27], and Picchio et al. (2019) [28]. Forest management should thus primarily aim to reduce the degree of compaction as well as the soil area influenced by machine traffic. Therefore, it remains crucial that reduced impact logging techniques are used, such as permanent skid trails, reduced soil contact pressure of machines, and a brash mat. [53].
Notably, this paper is limited to physical compaction on plant (beech seedling) response, but future research could assess the impacts of soil compaction on factors including food webs, understory species, and mycorrhizal fungi. Mycorrhizae, impacted soil compaction, benefit plants (herb layer species) mostly by enhancing their nutrient access and stress tolerance [58][59][60].

Conclusions
This research investigated the natural recovery of soil physico-chemical properties and seedling growth parameters under the skid trail surface after a different number of years from skidding operations conducted in the Hyrcanian forest. The results showed that after 20 years from the time of skidding, none of the soil characteristics were fully recovered, and the seedlings grown in these skid trails were significantly different (reduced growth parameters) from the seedlings grown in the undisturbed area. After 30 years, soil chemical properties (organic C, total N, and pH) were fully recovered, while soil physical properties (bulk density and porosity) were not fully recovered. Growth parameters and quality of seedlings grown in these skid trails showed no differences from seedlings grown in the undisturbed area. Moreover, the results of this study showed that the soil recovery rate on low slope skid trails was faster than steep slope skid trails. Sustainable timber production needs the continuous establishment of natural tree regeneration and adequate growth of seedlings in natural high forests. From the results of this study, it can be concluded that the process of recovering the physical characteristics of the soil in skid trails is long and takes more than 30 years. Therefore, an accurate execution of skidding operations, including management before, during and after skidding operations according to the guidelines of Best Management Practice (BMP) is necessary to reduce damage to the soil. The design of skid trails with a minimum distance of 140 m from each other and with a longitudinal slope of less than 15% should be considered as important principles to reduce the area and the severity of forest soil disturbance.