Co-variation relations of physical soil properties and site characteristics of Finnish upland forests

Physical soil properties have a marked influence on the quality of forest sites and on the preconditions for forest growth and management. In this study, water retention characteristics (WRC) and related physical soil properties in addition to vegetation coverage and tree stand data were studied at upland forest sites in Finland. Fixed and mixed models between soil and site characteristics were formed to estimate physical and hydrologic soil characteristics and the site quality with indirect co-varying variables. In the present data, the site quality index (H100) shows a high coefficient of determination in respect to the temperature sum. It is also related to soil fine fraction content, topsoil pH and water retention at field capacity. The thickness of the humus layer is predictable from the pH and cover of xeric and mesic plant species. The soil fine fraction content (clay + silt) is closely related to water retention at field capacity, the soil layer and site type, and without WRC to the temperature sum and site index and type, as well as the slope angle. The soil bulk density is related to organic matter, depth (layer) or alternatively to organic matter, slope and field estimated textural class (fine, medium, coarse). Water retention characteristics were found to be best determinable by the fine fraction content, depth and bulk density. Water content and air-filled porosity at field capacity are closely related to the fine fraction. This study provides novel models for further investigations that aim at improved prediction models for forest growth, hydrology and trafficability.


Introduction
Accurate spatial information on the quality (productivity) and properties of forest sites is essential for forestry management.Site quality indicating the growth potential of the forest soil is expressed in Finland typically using forest site types based on vegetation (Cajander 1949) and the site index (H100) based on the dominant height of the trees with largest diameter at 1.3 m (DBH) at the projected age of 100 years (Gustavsen 1980;Vuokila and Väliaho 1980).Forest stands and the ground vegetation tend to form according to the climatic and edaphic conditions (Cajander 1949;Childs and Flint 1990;Sims et al. 1996;Wang and Klinka 1996;Salemaa et al. 2008).Thus, the nutritional, thermal and hydrological regimes of the soil markedly affect the site quality.The hydrological conditions of a site depend not only on the physical properties of the soil, but also on the topographical location and the ambient weather conditions (Childs and Flint 1990;Nyberg 1996;Mäkitalo 2009;Campbell et al. 2013).
High-resolution data on weather conditions (precipitation, temperature, evapotranspiration) are of major importance in estimating the annual growth changes in forests (Mäkinen et al. 2002;Henttonen et al. 2014).In the Nordic countries, the temperature and precipitation regimes are predicted to change in the future due to the climate change whereupon the soil water content during the growth and dormant seasons, as well as the amount and timing of snow cover, soil frost and the soil bearing capacity may alter markedly.This means the preconditions for forest growth and forest management operations will undergo changes as well.Recent studies have indicated that, in addition to high-resolution weather data, accurate soil hydrological information are also needed to describe site quality and forest growth prerequisites more accurately (Henttonen et al. 2014;Peltoniemi et al. 2015).
The hydrological soil properties and conditions affect the terrain trafficability, which is one of the key requirements for the use of machinery in forest harvesting, silvicultural operations, soil preparation and planting for reforestation (Uusitalo et al. 2015(Uusitalo et al. , 2018;;Niemi et al. 2017).Knowledge of the soil strength and bearing capacity is required, e.g., for wood procurement planning and for the prevention of subsequent soil damage (Niemi et al. 2017).The soil strength and bearing capacity is commonly determined as the penetration resistance, which depends mainly on the soil bulk density, in addition to the clay and water content (Ayers and Perumpral 1982;Vaz et al. 2011;Pirnazaro et al. 2012;Uusitalo et al. 2015Uusitalo et al. , 2018)).Thus, physical and hydrological soil data are needed also for accurate determination of terrain trafficability.
In Finland, forest soil types and mean particle size classes are assessed by visual and tactile evaluation in the Finnish National Forest Inventory (NFI), which is based on a sampling method carried out on field sample plots (Tomppo et al. 2011).The soil type is classified in the following categories: organic soil, bedrock, stony soil (or boulder field), glacial till, and sorted soil.For glacial till and sorted soil, three main categories are determined according to the median particle size (Mälkönen 2003).In fine-grained soils, the median grain size is smaller than 0.06 mm, reflecting a high silt and clay content.Medium-grained soils are sandy (0.06-2 mm), and the median particle size of coarse-gravelly soils is 2-20 mm.
The proportion of forest land on mineral soils is estimated to be 72% in southern Finland and 81% in northern Finland (Tomppo et al. 2011).The most common mineral forest soil in Finland is a glacial till with a medium texture, which cover 63% and 75% of the forest land area on mineral soils in southern and northern Finland, respectively.Sorted soil (such as alluvial or aeolian soils with narrow particle size distribution) is the next most common mineral soil class covering 23% of the area of mineral soil forests in southern Finland and 17% in northern Finland.The most frequent soil types (according to World reference base soil classification system, WRB) are Podzols (50%), Histosols (25%), Arenosols (11%) and Leptosols (9%) (Tamminen and Tomppo 2008).Finer-textured soils, Cambisols (1.9%), Gleysols (1.4%) and Regosols (1.2%), have only a small proportion.
Despite the fact that the Finnish NFI covers all the forest area of Finland, it does not contain information of many edaphic determinants, which are important for tree growth and practical forestry management and are difficult to detect indirectly without laboratory analyses, such as soil physical and chemical or climatic and topographic attributes.For example, no measurements of soil water content or water retention characteristics are made in the NFI.Thus, there is a lack of feasible information on the relations between soil and site properties for forestry management, such as silviculture, trafficability and forest growth estimation.
The aim of this study is firstly to examine the co-variation relations of physical soil and site characteristics.Secondly, these characteristics are examined if they can provide novel relevant information and input parameters from indirect explanatory variables for operational and ecological modelling and practice in forestry management in Finland.The focus was therefore to determine predictive models of co-variation relations rather than to discuss in detail the causalities and ecological meanings of all factors involved.Thus, a chosen available set of soil water retention characteristics, related physical soil properties and site characteristics for upland forest sites were studied on a subset of permanent sample plots of the 8th Finnish National Forest Inventory (NFI).Transfer functions (fixed and mixed models) between the soil and site characteristics were formed to achieve relevant models from indirect explanatory variables.These functions and models can be applied on all sites measured in NFI and on other data containing the required information.Possibilities using our models as improved prediction tools of the changes in forest growth, hydrological and trafficability conditions are discussed.

Field data collection
Under the EU BioSoil project, upland forest soils in Finland were inventoried in the 2006-2007 period (Tamminen and Ilvesniemi 2013).The plots measured in the BioSoil project were a subset of the permanent plots of the 8th NFI (done in 1985-86) sites.The sampling network was relatively dense in southern Finland (16 km × 16 km), but was sparser in Lapland (24 km × 32 km).In the BioSoil project, 636 plots (each 400 m 2 in area) of the total of about 3000 NFI plots were measured, and over 6000 plots in all in Europe.
The rod penetration method (Viro 1952) was applied to estimate the volumetric percentage of stones in the upper 30 cm of the soil in the study plots.The particle size was sampled from two soil layers (10-20 cm and 20-40 cm) for each sample plot.Besides the BioSoil data, a set of undisturbed mineral soil samples have been previously collected for determining the water retention characteristics and related physical soil properties.On a representative subset of plots with low to moderate stoniness, cylindrical soil cores (d = 58, h = 60 mm) were taken from the depths of 0-6 cm and 30-35 cm (sampling depths differ from the NFI soil sample depths; partly due to the physical dimensions of the sampling cores).Altogether 82 plots with complete soil and vegetation data were available for the study (Fig. 1; Supplementary file S1, available at https://doi.org/10.14214/sf.9948).
Dominant heights (m) for forest stands were estimated at the age of 100 years (i.e.H100), which represents the site quality and fertility (Gustavsen 1980;Vuokila and Väliaho 1980).The stand age and height of the 100 largest trees at breast height per ha were determined regardless of the tree species and then estimated for the age of 100 year using equations for height development of the dominant trees.The site type was determined separately using the site-type classification (Cajander 1949;Miina et al. 2009) including herb-rich forests (Oxalis-Maianthemum and Oxalis-Myrtillus type), mesic heath forests (Myrtillus type), sub-xeric heath forests (Vaccinium type), xeric heath forests (Calluna type) and barren heath forests (Cladina type) (Table 1).The coverage of each plant species was determined for each sample plot.
At each plot, the terrain topography was defined as a slope, hilltop, flat or depression.The mean topsoil particle size class (fine-, medium-or coarse-grained) was assessed in the field by a visual and tactile evaluation.The precise GPS co-ordinate point of each sample plot was determined to also enable the characterization of the location of the site within the larger hydrological continuum.

Laboratory analyses
The particle size was analysed in soil layers (10-20 and 20-40 cm) using the laser-diffraction method and for 60 samples from the layer at a depth of 10-20 cm (containing also the samples with the highest clay contents) using the sedimentation method (Table 1).The sedimentation method was applied after organic matter had been removed from the samples using H 2 O 2 (Elonen 1971).
The values from the laser-diffraction method were calibrated with unpublished test data to the average level of the sedimentation method using 259 soil samples analysed using both methods.
The approximate organic soil matter content was determined based on loss on ignition at 550 o C for two hours (=100% minus percentage ash content).The soil pH in water (1:5) was estimated by a regression equation based on 33 soil samples [pH(CaCl 2 ) (in the range 2.9-5.8):pH(H 2 O) = 0.723 + 0.972 × pH(CaCl 2 )].
The volumetric water-retention characteristics of undisturbed soil samples were measured by means of a pressure plate apparatus (Soilmoisture Equipment Corp., Santa Barbara, CA, USA) at decreasing matric potentials (i.e.under desorption) (Klute 1986;Heiskanen and Mäkitalo 2002).Soil-filled metal cylinders (d = 58, h = 60 mm) were first saturated overnight, then allowed to drain freely (to -0.3 kPa matric potential) and then exposed to successive matric potentials down to -1500 kPa.The bulk density was determined as the ratio of the dry mass (dried at 105 °C) to volume at -0.3 kPa.The particle density (Dp) was estimated using an average density of 2.65 g cm -3 applied to mineral and 1.5 g cm -3 to organic components (Heiskanen 1992).The total porosity (TP) was estimated to be (Dp-Db) Dp -1 , where Db is the bulk density.The air-filled porosity (AFP) was estimated as TP-WC, where WC refers to the volumetric water content.The water content at field capacity was estimated as the water retention at a matric potential of -10 kPa (Al Majou et al. 2008).

Statistical analysis and models
We modelled the soil properties using linear mixed effects models and a linear model without a random factor.Sample plots were considered as a random variable in the mixed effects models.The number of sample plots (level j in the models) was 82 all over the country (Fig. 1).The number of observations between the models varied a little in the two hierarchy levels due to some missing observations in the data.Soil properties measured from two layers, near the soil surface (0-6 cm) and from the depth of 30-36 cm (level i in the models) were considered to be the correlated observations within the sample plots.The explanatory variables represented both levels in the mixed effects models.Log-and square-root transformations were tested in the response variables (y ij , y j ) to normalize the residuals.A square-root transformation was selected for all of the models where the transformation was needed to achieve adequate residuals.The mixed effects models can be written as: where: , = fixed coefficients that have been measured in sample plot (j) and at soil layer level (i) , = fixed variables that have been measured in sample plot (j) and at soil layer level (i) μ 0j = random sample plot effect (random intercept) ε ij = error variance k 1 = the number of fixed variables representing the soil layer effects k 2 = the number of fixed variables representing the sample plot effect.
The coefficients of determination (R 2 ) for the mixed effects models were computed using the R software package piecewiseSEM (Nakagawa and Schielzeth 2013;Lefcheck 2015).The mixed models were constructed using lme4 (Bates et al. 2015) R packages.Because of the large number of achieved models, we present only the effects of the model predictors in the effect plots of the most relevant variables.
The model for the thickness of the humus layer was a general linear model without the random sample plot effect, because the response variables were only measured at the independent sample plot level.The soil characteristics of the soil layers were averaged for use as the explanatory variables in the models.Alternatively, the soil properties were tested in the models separately, but the averaging produced almost similar results.The general formula for the linear model could be written as: where: β 0 = the intercept β j = the coefficients of the model, measured at the sample plot level (j) x j = variables of the model, measured at the sample plot level (j) ε j = error variance.
There were about 300 plant species, cover percentage of which had been measured on the sample plots.The species were considered as potential explanatory variables for the soil properties, not in a causal sense, but merely as indicators of the soil properties.However, testing the effect of each species would have been too laborious a task.Thus, we used two approaches to study the effects of the species.First, groups of closely related species were formed using Item Cluster Analysis (ICA), which is an advanced alternative e.g. to principal component or factor analysis (Revelle 2017).However, the groups had to be ecologically relevant.Sum variables (the cover of the species) were formed based on the ICA.A Spearman's rank order correlation matrix of the species with a cover of at least 0.5% was used as an input to the ICA (the threshold of 0.5 was used because the correlation matrixes could not be computed for all the measured species).The predictions for the explanatory variables in the models, based on the fixed effects were computed and plotted using the R software package effects (Fox 2003).
The second approach was to analyse which single species were strongly related to the gradients of the soil properties to select these species for the explanatory variables for the models.A Constrained Correspondence Analysis (CCA) with gradient vectors combined with a CCA-plot was used in this approach.The R software package, vegan (Oksanen et al. 2017), was used in the CCA analysis.The species that were selected based on the CCA (situated near the ends of the gradient vectors) were tested in two ways: based on their cover and based on their existence (whether they exist or not in the sample plot).Low frequencies and low cover with the relatively small number of sample plots may have affected non-significant effects of these variables.The defined vegetation groups are described in Table 1.

Results
In the present data covering whole Finland, the atmospheric temperature sum has the strongest correlation with the site index H100.The H100 index is also closely related to soil attributes and is positively related to the topsoil pH and soil water retention at field capacity (WC10) (Table 2).Almost equally good coefficients of determination can also be achieved with the fine fraction content and site type as explanatory variables (Tables 2, 3).The tree species composition showed no significant effect on H100.The coverage of plant species (Table 1) had significance in models only as vegetation groups which indicate site quality (see below).
The thickness of the humus layer decreases with increasing topsoil pH (Fig. 2).It is well predicted from the pH and cover of xeric and mesic plant species (Tables 2, 3).The fine fraction content (clay + silt) of the soil decreases with temperature sum and site type quality (Fig. 3) and is also closely related to water retention at field capacity (WC10), as well as the soil layer (Tables 2, 4).Without water retention characteristics, it is related to the temperature sum, site index and type as well as the slope angle (Tables 2, 3; Fig. 3).The soil bulk density is negatively related to organic matter and slope angle (Fig. 4; Table 3).Alternatively the bulk density relates to organic matter, the slope and field estimated soil textural class (fine, medium or coarse) (Table 2).
The water retention characteristics are determinable by the fine fraction content, the soil depth (layer) and the bulk density (for details see Tables 2, 4).The water content (WC10) and air-filled porosity (APF10) at field capacity are especially closely related to the fine fraction content (Fig. 5).3).The shaded area denotes a 95% confidence interval.3).The shaded area and vertical bars denote a 95% confidence interval.
Table 4. Estimate statistics of mixed effects models for the soil's fine fraction content and bulk density as well as for water retention (WC10) and air-filled porosity (AFP10) at -10 kPa matric potential (SE for fixed and SD for random effects).Estimates for categorical variables refer to the zero reference level (i.e.site types vs. reference category of herb-rich site type and deeper soil layer vs. upper layer and field estimated soil particle size vs. reference category of fine sites).3).The shaded area denotes a 95% confidence interval.4).The shaded area denotes a 95% confidence interval.

Discussion
In the present data, the forest site index (H100) was found to be most closely related to the atmospheric temperature sum, which suggests the effects of latitude and altitude on temperature and preconditions for growth.The topsoil pH and soil water retention at field capacity (WC10) were significant predictive soil attributes for the H100 index; the effects of which presumably vary along the transect from south to north due to climate (temperature, precipitation, transpiration).The soil productivity is affected by the nutritional, thermal and hydrological regimes of the site and varies from fertile to poor sites and from south to north in Finland (Cajander 1949;Tamminen 1993;Mäkinen et al. 2002;Henttonen et al. 2014).There is a dominant relationship between understorey vegetation and the site type and productivity in the boreal forest ecosystem (Cajander 1949).
In this study, the thickness of the humus layer was found to be closely related to the topsoil pH and cover of xeric and mesic plant species.The thickness of the humus layer typically varies between 1-15 cm in Finland (Tamminen 1993(Tamminen , 2000;;Westman 1990).The thickness of the humus layer and the age of the tree stand reduce the pH in the humus layer while the soil's fine fraction content and the dominance of deciduous tree species increase it.In the mineral soil layers, the soil pH is higher on average on stoneless, fine-grained sites with a thin humus layer and no Norway spruce (Tamminen 2000).Organic matter accumulation usually results in a reduction in soil pH, which is caused by partial microbial degradation and production of organic acids (Skyllberg 1990;Westman 1990;Ilvesniemi 1991;Tamminen 2000;Abakumov et al. 2012).
In general, humus-layer development is influenced by the vegetation type, forest-management practice, climate (temperature), geological substrate, soil faunal and microbial activity, and the stage of forest succession.Spatial variability of the humus layer thickness can have important impacts on soil water dynamics, nutrient storage and availability, as well as plant growth (Hokkanen et al. 1995;Ponge et al. 2002;Bens et al. 2006;Laamrani et al. 2014).The structure and thickness of the humus layer can also affect the trafficability of the site (Haarlaa 1973;Wästerlund 1989).
In the present site data, the soil's fine fraction content (clay + silt) was found to be closely related to water retention at field capacity (WC10), soil layer and site type, and without water retention characteristics, it was found to be related to the temperature sum, site index and type as well as slope angle.Water retention characteristics are determinable by the soil's fine fraction content, soil depth (layer) and bulk density.The water retention (WC10) and air-filled porosity (APF10) at field capacity are especially closely related to the soil's fine fraction.The soil's bulk density was found to be related to organic matter, soil depth (layer) or alternatively to organic matter, slope, and the field estimated soil textural class (fine, medium, coarse).
The hydrological regime of a site depends not only on the physical properties of the soil, but also on the topographical location and the ambient weather conditions (Childs and Flint 1990;Nyberg 1996;Heiskanen and Mäkitalo 2000).In a study in northern Finland, for example, the soil's fine fraction content was found to be dependent on the topographical position (18% on summits, 61% on toe-slopes) (Mäkitalo 2009).Forest site types and the dominance of tree species can differ markedly in the soil texture, hydrological regime and soil water retention characteristics of the site (Sepponen et al. 1979;Heiskanen 1988;Heiskanen and Mäkitalo 2002).
The soil's bulk density is correlated to the organic matter content and particle size (Westman 1990;Tamminen 1993;Tamminen and Starr 1994).In Finland, the soil bulk density rapidly increases with depth at the surface but remains uniform at depths over 20 cm (Tamminen and Starr Silva Fennica vol. 52 1994).The soil's organic matter content decreases the bulk density significantly.The fine fraction (<0.06 mm) content tends weakly to decrease the bulk density and the gravel (2-20 mm) content to increase it.
On average, the soil water retention characteristics (WRC) are mainly dependent on the soil particle size, bulk density and organic matter content (Gupta and Larson 1979;Heiskanen 1988;Childs and Flint 1990;Heiskanen and Mäkitalo 2002;Mecke et al. 2002;Wall and Heiskanen 2003;Jauhiainen 2004;Heiskanen et al. 2016).As in this study, water retention especially near field capacity has also previously been found to correlate well with the fine fraction content in forest soils (Heiskanen 1988;Sepponen et al. 1979;Hewelke et al. 2015).If the soil's fine fraction content is too high this may result in excessive soil water and poor soil aeration (AFP10 below 10%) (Sepponen et al. 1979;Heiskanen and Mäkitalo 2002;Wall and Heiskanen 2003;Mäkitalo 2009; cf.Fig. 5).
The mean actual soil water content (SWC) is about the same as the water retention at field capacity or at -10 kPa matric potential (i.e.WCfc ≈ WC10) (Heiskanen 1988;Al Majou et al. 2008).The actual SWC is similarly dependent on the soil properties as the WRC but is also dependent on the depth of the ground water as well as transient water storage changes in the weather (precipitation/ transpiration) (Heiskanen and Mäkitalo 2002;Mäkitalo 2009;Henttonen et al. 2014).The SWC can also be related to topographic slope data (TWI, topographic wetness index) or the distance to ground water (DTW, depth-to-water index, m) (Mäkitalo 2009;Campbell et al. 2013;Niemi et al. 2017).In northern Finland, with an increasing TWI, the soil's fine fraction content, organic matter content and WC10 increases, while the AFP10 decreases significantly (Mäkitalo 2009).
Soil type, along with traffic intensity, affects significantly the soil bearing capacity and rut formation in logging operations (Ayers and Perumpal 1982;Vaz et al. 2011;Mohtashami et al. 2017;Uusitalo et al. 2018).The soil bearing capacity (expressed as penetration resistance, PR) is mainly dependent on the soil bulk density, clay content and SWC.In fully saturated silty and clayey soils, the PR in the uppermost 15 cm can be as low as 600-800 kPa, whereas in dry soils the PR varies between 1500-2500 kPa (Uusitalo et al. 2018).The PR of the humus soil is typically as low as 500-600 kPa, which indicates that the whole humus layer is compressed or moved aside forming deep ruts under heavy logging operations.Campbell et al. (2013) have shown in Canada that an explanatory model for the PR (cone index) using elevation, slope and distance to ground water level (DTW index) (R 2 = 0.50) can be improved by including the field-measured SWC to the model (R 2 = 0.66).DTW alone may not predict rut formation in forest logging (Mohtashami et al. 2017).

Conclusions
The present data and models for the site index, water retention characteristics and related physical soil properties form a novel basis and input parameters for further investigations, which aim to more accurately predict the relevant soil properties for practical forestry management using topographical and forest inventory data.The soil bearing capacity is related to the soil penetration resistance, which is commonly determined by the soil bulk density, clay and soil water contents.Thus, the estimated soil bulk density, soil water content at field capacity (WC10) and fine fraction content of this study can, for example, be used as explanatory soil data in a geographic information system for estimating the average soil bearing capacity locally and extensively for the Finnish forest upland sites.However, the present data cannot fully reveal the effects of physical soil properties on the site characteristics, hydrological regime or trafficability at all locations because of its relatively limited sample and incomplete cover of the forest land area.For example, western and southern Finnish coastal areas, where forest soils contain more silt and clay, are scarcely represented.Therefore, further verification studies are needed in the next phase with more comprehensive data sampling together with more detailed location-based data analyses.This would achieve a location-based data analysis resulting in high-resolution data which would contribute to achieving more accurate prediction models for forest growth, hydrology and trafficability in various sites.

Fig. 2 .
Fig. 2. Predicted relation of the thickness of the humus layer to topsoil pH (0-6 cm) (see the fixed model in Table3).The shaded area denotes a 95% confidence interval.

Fig. 3 .
Fig. 3. Predicted relations of the soil's fine fraction content to the temperature sum and site type (see the fixed model in Table3).The shaded area and vertical bars denote a 95% confidence interval.

Fig. 4 .
Fig. 4. Predicted relations of the soil bulk density to the terrain slope angle and soil organic matter by soil layer (see the mixed effects model in Table3).The shaded area denotes a 95% confidence interval.

Fig. 5 .
Fig. 5. Predicted relation of the water content (WC10) and air-filled porosity (AFP10) at -10 kPa matric potential to the soil's fine fraction content (clay + silt) (see the mixed models in Table4).The shaded area denotes a 95% confidence interval.

Table 2 .
Summary statistics for the selected major response variables with the three most significant (Anova F-test) explanatory variables.Site type and soil layer are categorical variables.
Db = bulk density, Org.= organic matter (%), Fines = soil fine fraction, * the cover percentages of plant species indicating site quality (see Table1), ** the fine fraction as a mean of layers.

Table 3 .
Estimate statistics of fixed effects models for the H100 site index, thickness of humus layer (cm) and soil fine fraction content (%) (without hydraulic explanatory variables).Estimates for categorical variables refer to the zero reference level (i.e.site types vs. reference category of herb-rich site type).