Biosystems Diversity

Recreation is an important cultural ecosystem service and is able to significantly affect soil heterogeneity and vegetation functioning. This study investigated the role of the relief and tree stand density in the apparent soil electrical conductivity variation within an urban park. The most suitable variogram models were assessed to evaluate the autocorrelation of the regression models. The map of the spatial variability of apparent soil electrical conductivity was built on the basis of the most suitable variogram. The experimental polygon was located in the Botanical Garden of Oles Honchar Dnipro National University (Dnipro City, Ukraine). The experimental polygon was formed by a quasi-regular grid of measurement locations located about 30 m apart. The measurements of the apparent electrical conductivity of the soil in situ were made in May 2018 at 163 points. On average, the value of soil apparent electric conductivity within the investigated polygon was 0.55 dSm/m and varied within 0.17–1.10 dSm/m. Such environment predictors as tree stand density, relief altitude, topographic wetness index, and potential of relief to erosion were able to explain 48% of the observed variability of soil electrical conductivity. The relief altitude had the greatest influence on the variation of soil electrical conductivity, which was indicated with the highest values of beta regression coefficients. The most important trend of the electric conductivity variation was due to the influence of relief altitude and this dependence was nonlinear. The smallest values of the soil electrical conductivity were recorded in the highest and in lowest relief positions, and the largest values were detected in the relief slope. Recreational load can also be explained by the geomorphology predictors and tree stand density data. These predictors can explain 32% of the variation of recreational load. The variogram was built both for the soil apparent electrical conductivity dataset and for the residuals of the regression model. As a result of the procedure of the models’ selection on the basis of the AIC we obtained the best estimation of the variogram models parameters for the electrical conductivity and for the regression residuals of the electrical conductivity. The level of recreation was correlated statistically significantly with the apparent soil electrical conductivity. The residuals of regression models in which geomorphological indicators and tree stand density were used as predictors have a higher correlation level than the original variables. The soil electrical conductivity may be a sensitive indicator of the recreation load.


Introduction
Urban ecosystems provide a combination of ecosystem services such as provisioning, regulating, and habitat cultural services (Brouwer et al., 2013). Soils of urban parks influence carbon and nitrogen pools and fluxes (Raciti et al., 2011). The topsoil layer, forest litter, and vegetation cover have the key function of preventing soil erosion (Zuazo & Pleguezuelo, 2008;Kunah et al., 2019). Urban soils and vegetation are very different from natural ones due to the recreation impact (Levin et al., 2017). Recreation is an important cultural ecosystem service (Chiesura, 2004;Balzan & Debono, 2018). Urban soils are subjected to high anthropogenic influence (Scalenghe & Marsan, 2009;Yorkina et al., 2019). Recreation is able to significantly affect soil heterogeneity and vegetation functioning (Kutiel et al., 1999). The increasing of recreation loading was shown to induce the decrease of vegetation cover and biological diversity (Sujetovienė & Baranauskienė, 2016;Yorkina et al., 2018). Urbanisation transforms the vegetation and soil cover on a city territory, which in turn leads to the change of parameters of the carbon cycle (Svirejeva-Hopkins et al., 2004). Soils and the aboveground structure of urban systems are highly spatially heterogeneous and interaction of various kinds of heterogeneity in urban systems is an open question (Pickett et al., 2008). Urban soils are very variable in space and time (Vasenev et al., 2014).
It has been proposed to categorize the soils of anthropized areas according to the ecosystem services they provide in urban areas (Morel et al., 2014). Urban green spaces support conservation of the biodiversity in urban areas (McKinney, 2006;Lepczyk et al., 2017). The heterogeneous microenvironment structure of parks promotes the preservation of natural vegetation (Sarah et al., 2015;Gritsan et al., 2019). It has been shown that the measurable biological indices may be applied for assessment of ecological, environmental-regulating, and productive functions of urban soils (Vasenev et al., 2012). Soil biota has a considerable utility for estimation of the ecological potential of urban soils (Maltsev et al., 2017). The diversity of soil biota is important to many environmental functions such as water depollution, biochemical cycles, fertility and carbon storage (Guilland et al., 2018). The variability of soil properties of urban parks affects the growth and development of plants (Pregitzer et al., 2016). Significant variation of the soil properties was found in a distance gradient of measurements taken around selected individual trees affecting the quality and quantity of understorey vegetation in park forest stands (Sikorski et al., 2013). The variability of soil properties promotes the maintenance of biodiversity in urban areas (McKinney, 2006). The apparent soil electrical conductivity (ECa) is a useful and express measure of soil variability (Corwin et al., 2003;Yorkina et al., 2018). This index is related to soil properties affecting ecosystem primary production (Corwin & Lesch, 2005). The characterization of soil spatial variability using ECa may be used for soil quality assessment (Corwin, 2005). The spatial variability of the apparent soil electrical conductivity was modeled on the basis of regression dependencies from remote sensing predictors (Zhukov et al., 2016).
The principal issue of the ecological modeling is the precise assessment of the spatial variability of soil properties (Shit et al., 2016). An inverse distance weighting or ordinary kriging are the effective approaches for interpolation of the spatial patterns of soil properties (Uygur et al., 2010;Göl et al., 2017;Tang et al., 2017). The efficiency of predicting spatial variability of soil properties was proposed to be improved by a combination of regression and spatial interpolation (Hengl et al., 2004). This approach was called regression-kriging (Kumar et al., 2012;Peng et al., 2013;Mondal et al., 2017). Regression-kriging is one of the most popular, practical and robust hybrid spatial interpolation techniques for modeling of the soil distribution patterns at multiple scales in space and time (Keskin & Grunwald, 2018).
The objectives of this study were (a) to investigate the role of the relief and tree stand density in the apparent soil electrical conductivity variation within an urban park, (b) to assess the most suitable variogram models to evaluate the autocorrelation of the regression models, and (c) to map spatial variability of the apparent soil electrical conductivity.

Methods
The experimental polygon was located in the Botanical Garden of Oles Honchar Dnipro National University (Dnipro City, Ukraine) ( Fig. 1). The climate at the experimental polygon is temperate-continental. According to statistics from 1998 to 2018, the average yearly precipitation was approximately 565 mm. The average temperature was the highest in August at 25.7 °C, while the lowest was in January at -10 °C. There were two soil types within the experimental polygon: calcic chernozem (siltic, tonguic, upland and hillside positions) and technosol (formed on construction wastes in the lowland positions). To measure the electrical conductivity of the soil in situ, the sensor HI 76305 was used (Hanna Instruments, Woonsocket, R. I.). This sensor works in conjunction with the portable device HI 993310 (at the depths of 0-5 cm in a three-times repetition). The experimental polygon was formed by a quasi-regular grid of measurement locations with distance between them about 30 m. The measurements were made in May 2018 at 163 points (Fig. 1). The investigated area of the experimental polygon was 115,296 m 2 . The sampling point coordinates were measured using a GPS-navigator. The number of tree trunks was determined within a 5 m radius of the soil penetration resis-tance measurement point. In the forest plantation the presence of 19 tree species was revealed, among which the most common were Robinia pseudoacacia L., Ulmus glabra Huds., Populus nigra L., and Acer campestre L.
Digital elevation model (DEM) is a presentation of the Earth's surface in numerical format (Dowman, 1999). The Advanced Land Observation Satellite -ALOS (www.eorc.jaxa.jp/ALOS/en/index.htm) data were used to generate a digital elevation model. Spatial resolution for the study area is 10 meters, nominal vertical accuracy and nominal horizontal accuracy is 5 meters. By means of kriging procedure, DEM was resampled to a resolution of 10 m (Susetyo, 2016). The kriging procedure also made DEM suitable for calculation of derived layers -topographic wetness index and erosion factor (Hojati & Mokarram, 2016).
The concept of the topographic wetness index (TWI) was first proposed by Beven & Kirkby (1979) and may be calculated by the formula: TWI = ln(a/tanβ), where a is the upslope contributing area per unit contour length, β is the local slope. TWI is unitless. High values of TWI indicate an area high with increased accumulated runoff potential (Kunah & Papka, 2016a, b).
Potential of relief to erosion (LS) is one of the components of the Universal Soil Loss Equation (USLE). LS is the product of L-and Sfactors. The L-factor defines the slope length, and the S-factor is slope steepness. The LS-factor is dimensionless, having values equal to and greater than 0 (Panagos et al., 2015). In this study TWI and LS-factor were calculated with the aid of the SAGA (Olaya & Conrad, 2008).
Kriging is common technique in geostatistics (Minasny & McBratney, 2005). The variogram is a central concept in geostatistics and knowledge of the precise mathematical form of the variogram enables us to determine a spatial variation (McBratney & Pringle, 1999). The intercept of the variogram model curve is specified as the nugget (τ 2 ), the difference between the asymptote and the nugget as the sill (σ 2 ), and the distance at which the theoretical variogram curve reaches its maximum as the range. For models with an infinite range, the value at which the variogram reaches 95% of the asymptote is called the practical range. Commonly used variogram models (spherical, exponential and Gaussian) are characterized by lack of flexibility (Stein, 1999). As an alternative, one can consider the Matern variogram class of models (Matern, 1986). Matern models have considerable flexibility for modeling the spatial covariance and are able to describe a wide variety of local spatial processes. Based on this, the Matern model is proposed to be used as a general approach for the simulation of soil properties (Minasny & Mc Bratney, 2005). Matern isotropic covariance function has the form (Handcock & Stein, 1993;Stein, 1999): , where h is the separation distance; K ν is modified Bessel function of the second kind of order κ (Abramowitz & Stegun, 1972), Г is the gamma function, φ is the range parameter (φ > 0), which measures how fast correlation decays with distance; κ is the smoothness parameter. The Matern model has a high flexibility compared with traditional geostatistical models in view of the smoothing parameter κ. When the parameter κ = 0.5, the Matern model fully corresponds to an exponential model. When κ → ∞, the Matern model corresponds to a Gaussian model. If κ = 1, it corresponds to a Whittle's function (Whittle, 1954;Webster & Oliver, 2001;Minasny & McBratney, 2005). If the range parameter r is large (r → ∞), then the spatial process is approximated by the power function when κ > 0, and a log function or de Wijs function if κ → 0 (de Wijs, 1951(de Wijs, , 1953. The spatial dependence level (SDL) or nugget to sill ratio is an indicator of the strength of the spatial autocorrelation (Cambardella et al., 1994;Sun et al., 2003;Zhukov et al., 2019a, b). Calculations were made using geoR library (Paulo et al., 2016).
Regression kriging is a spatial interpolation technique that combines a regression of dependent variables on predictors with kriging of the prediction residuals (Hengl et al., 2004): where is the fitted deterministic part, is the interpolated residual. Thus, the first part of the right-hand side of the equation represents the regression and the second part represents the kriging of the residual.
To measure the accuracy of differential entropy maps we use crossvalidation procedure and consequently we compute normalized root mean squared error (NRMSE), mean error (ME) and mean squared deviation ratio (MSDR) (Vašát et al., 2013). Mean squared error (RMSE) was calculated as follows: . Normalized root mean squared error (NRMSE) was calculated as follows: . Mean squared deviation ratio (MSDR) was calculated as follows: . where x 1 is a prediction of the variable X; x 2 is a measure of that variable; n is the number of records; var is a kriging variance. The smaller the NRMSE values, the more accurate the map. The MSDR indicates whether the variance of measurement data is well reproduced with the kriging interpolation and ideally it equals to 1 (Vašát et al., 2013). The R-squared of the regression between observed and predicted after cross validation values was used as they are very intuitive. Cross-validation procedure was performed using function xvalid from package geoR library (Paulo et al., 2016). Spatial variation of predictors and regression models of the soil mechanical impedance residuals was displayed using the Surfer ® 12 from Golden Software, LLC (www.goldensoftware.com).
As a measure of the recreation loading the Strava data were used (www.strava.com). Strava (San Francisco, CA, USA) is a widespread social network for cyclists and runners. Strava consists of a mobile app and a website. The Strava app records distance, time, average speed and route (GPS trajectory) of each activity. Strava's database comprises nearly a trillion GPS points globally and is growing by over 8 million activities every week (Sun & Mobasheri, 2017).

Results
On average the value of soil apparent electric conductivity within the investigated polygon was 0.55 ± 0.01 dSm/m and varied within 0.17-1.10 dSm/m ( Table 1). The tree density reached 8 trees in a circle with radius of 5 meters. The study area was influenced by a significant recreational load (median of 43.2%). The main load was concentrated along the walkways, but the areas with natural vegetation were also impacted by the significant recreational exposure (Fig. 2). Relief altitude ranged 133.8-153.2 m above sea level. The topographical wetness index varied within 5.81-12.76, and the potential of relief to erosion (LSfactor) was 0.03-1.32. The predictors were able to explain 48% of the observed variability of soil electrical conductivity ( Table 2). The relief altitude had the greatest influence on the variation of soil electrical conductivity, which was indicated with the highest values of beta regression coefficients. The topographical wetness index and LS-factor were found to be the statistically significant predictors of the soil electrical conductivity variation. The tree stand density was not a statistically significant predictor of the soil electrical conductivity on the spatial level investigated.
Recreational load can also be explained by the geomorphology predictors and tree stand density data. These predictors can explain 32% of the variation of recreational load. As in the case with soil electrical conductivity, the relief altitude was shown to be a leading predictor. The tree stand affected the recreation. This effect was non-linear. The most favourable park areas for recreation occupied the places with moderate tree density stands. Areas of the park with too dense tree stands or with too sparse tree stands were less attractive for recreation. The negative correlation between recreational load and topographic wetness index indicated that the areas of the park with a tendency to waterlogging were less attractive for recreation. Fig. 2. The spatial variation of recreation loading within the studied polygon according to Strava Global data (www.strava.com): data were rescaled to a range: 0 -no loading; 100% -the maximum loading The variogram was built both for the soil apparent electrical conductivity dataset and for the residuals of regression model (Fig. 3). The nugget effect was searched in first stage of the assessment of the best values of variogram model parameters for fixed values of Kappa = 2 and the starting value of Phi = 7 m (Fig. 4). The variation of soil electrical conductivity was characterized by the much smaller spatial dependence than the variation of regression model residuals in terms of the AIC criterion (Fig. 4a). An increase in the fixed values of the nugget effect was accompanied by a monotonic growth of the optimal value of variogram range parameter when Kappa parameter was fixed (Fig. 4b). In the next stages the parameters that were found to be optimal in the previous stage were selected as starting parameters. An experimental increase of Kappa parameter for electrical conductivity demonstrates the existence of optimal value, which is equal to 1.5 (Fig. 5). The AIC change is monotonous with the increase of the Kappa parameter for variogram of the residuals of regression, indicating its best value, which approaches to ∞.
As a result of the conducted procedure we obtained the best estimation of the variogram models parameters for the electrical conductivity and for the regression residuals of the electrical conductivity (Table 3). The parameter Phi and the practical range were found to be much greater for OK than these parameters for RK. In turn, the spatial component of the variation characterized by an index of SDL is much higher for RK than for OK. The optimal value of the Kappa parameter (Kappa → ∞) for OK of the Matern model turns it to the Gaussian model. a -x-axis is the fixed value of the nugget-effect represented for comparability by the SDL value, %; y-axis is the AIC values for ordinary kriging (left axis) and for regression-kriging (right axis); for ordinary kriging the AIC has a minimum in x (nugget) = 0.0103 (corresponding SDL = 27.2%), for regression-kriging the AIC has a minimum in x (nugget) = 0.0028 (corresponding SDL = 9.6%)

Fig. 5.
The dependence of the AIC on the model parameter Kappa (Nugget = 0.0103 for OK and Nugget = 0.0028 for RK; Phi = 17.3 m for OK and Phi = 5.3 for RK): x-axis is the Kappa parameter, y-axis is the AIC values for ordinary kriging (left axis) and for regressionkriging (right axis); the AIC has a minimum value at kappa is 1.5 for the ordinary kriging and the AIC demonstrates an increase with the growth of parameter Kappa (Kappa → ∞) for regression-kriging The best model for the RK is slightly smoother than the Whittle model (Whittle, 1954;Webster & Oliver, 2001;Minasny & McBratney, 2005), for which Kappa = 1. The lower NRMSE value indicates a more accurate map of electrical conductivity in the case of regression-kriging (Fig. 6). The more precise reproduction of the electric conductivity variations is in the case of RK, as indicated by the higher value of the MSDR statistic. Notes: Phi -variogram range (the distance at which the theoretical variogram curve reaches its maximum as the range); Practical Range -the value at which the variogram reaches 95% of the asymptote; Sill -the difference between the asymptote and the nugget; Nugget -the intercept of the variogram model curve; SDLnugget to sill ratio as an indicator of the strength of the spatial autocorrelation; Kappa -Matern model smoothing parameter; Regression Radj2-adjusted R 2 of the regression model with terrain and tree stand variables as predictors; NRMSEnormalized root mean squared error; MSDR -mean squared deviation ratio. The level of recreation was correlated statistically significantly with an apparent soil electrical conductivity (Fig. 7). However, as was shown earlier, both of these parameters were dependent on the geomorphology predictors. The residuals of regression models in which geomorphological indicators and tree stand density were used as predictors have a higher correlation level than the original variables.

Discussion
The soil electrical conductivity is an express parameter which can be easily measured in large quantities for the spatial analysis. It can be used as a direct indicator of the soil condition including the influence of recreational load Sarah et al., 2015). Also, information on soil electrical conductivity can be used to design an optimal placement of test polygons for the spatial modeling of the other soil and ecological properties whose number of samples is limited due to the complexity of carrying out field studies (Siqueira et al., 2016). The results obtained indicated an important role of relief predictors for explanation of the spatial variability of soil electrical conductivity. The most important trend of the electric conductivity variation was due to the influence of relief altitude and this dependence was nonlinear. The smallest values of the soil electrical conductivity were recorded in the highest and in lowest relief positions, and the largest values were detected in the relief slope. It should be noted that the highest or lowest relief positions were the most favourable for recreation. The soil properties and herbaceous vegetation characteristics were revealed to be affected by human activities. In turn, the above characteristics were affected by natural factors mainly in the microenvironments which were subjected to low levels of recreation loading (Sarah et al., 2015). The litter layer and soil organic horizon are most significantly affected by recreation (Amrein et al., 2005;Brygadyrenko, 2015;Faly et al., 2017;Faly & Brygadyrenko, 2018). Recreation leads to soil compac-tion (Özcan et al., 2013) and reduces the grass cover and litter layer, resulting in a deteriorating water regime of soils (Oral et al., 2013), which may be manifested as a reduction of the soil electrical conductivity. This interpretation also explains the fact that we have not found the statistically significant influence of the tree stand density on the soil electrical conductivity. This is due to the fact that a significant recreational load is characteristic for areas both with dense tree stand and for areas without forest cover. This result confirms the assumption that the soil electrical conductivity may be a sensitive indicator of the recreation load. a b Fig. 7. Correlation between recreation load and apparent soil electrical conductivity (a) and correlation between residuals of regression models of recreation and apparent soil electrical conductivity with geomorphological characteristics and tree stand density as predictors: a -absciss axis is a recreation load (%), ordinate axis is an observed electrical conductivity (dSm/m); b -absciss axis is the residuals of regression model of recreation load (%), ordinate axis is the residuals of regression model of apparent soil electrical conductivity (dSm/m) A positive relationship between the soil electrical conductivity and the topographic wetness index and a negative relationship with LS-factor are logical. These predictors are factors in the variation of the soil electrical conductivity, which form the natural background of this indicator. It is obvious that the estimation of the recreational component in the variation of the soil electrical conductivity is possible after extraction of the underlying variability of this indicator induced by the relief factors.
The variation of soil electrical conductivity is characterized by the presence of a significant spatial component. Hengl et al. (2004) introduced the process of using the regression-kriging (RK) method for spatial prediction of soil variables. The ordinal kriging and regression-kriging comparison indicates that the relief predictors contribute to the formation of a large-scale component of the spatial variation of soil electrical conductivity. An extraction of the influence of relief predictors after regression procedure allows one to obtain a fine-scale component variation of the soil electrical conductivity with a much larger spatial autocorrelation. That is why regression-kriging allows one to produce a more detailed map of the spatial variation of soil electrical conductivity. The knowledge of the precise mathematical form of the variogram enables one to predict the soil properties on a local or regional level (Minasny & McBratney, 2005). The procedure for searching the best variogram model parameters based on the AIC has shown that the Gaussian model is the best for the regression residuals. There is not any suitable model from the set of the traditional models for OK, so the Matern model with the parameter Kappa = 1.5 is the most appropriate model.

Conclusion
The recreation load and apparent soil electrical conductivity are influenced by geomorphological properties of relief and tree stand density. The electrical conductivity can be used as an indicator of the soil cover transformation under the influence of recreation. However, in order to correctly estimate the level of recreation through electrical conductivity of the soil, it is necessary to make a preliminary extraction of the variation component due to other environmental factors. This procedure also has an impact on the geostatistical characteristics of the spatial pattern of the apparent soil electrical conductivity. The traditional variogram models are not suitable for spatial modeling of the apparent soil electrical conductivity. The Matern model is the most flexible and allows one to obtain a more accurate model of the spatial process. The variogram of the residuals of regression model of apparent soil electrical conductivity with geomorphological properties and density of tree stand as predictors is characterized by a smaller practical range, which also indicates a possible recreational component of the formation of spatial patterns of this indicator.