Multiscale Organization of Landscape Structure in the Middle Taiga of European Russia

Dominant landscape-ecological models either focus on the hierarchical organization of a single phenomenon or describe relations at a single hierarchical level. We proposed the tool MALS (Multiscale Analysis of Landscape Structure) to reveal multiple independent hierarchies based on the interactions between properties of relief, soils and vegetation and tested it on the example of the middle-taiga landscape in European Russia. Morphological properties of soils and abundance of plant species were measured in operational territorial units. Multidimensional scaling was used to reveal ecological drivers. Combinations of landforms from DEM were used to describe spatial heterogeneity in the higher-order geosystems. Response surface regression was applied to relate soils and vegetation to each other and to relief of several hypothetic higher-order geosystems. Spatial extent of a higher-order geosystem was determined from the series of equations. Then we compared contributions of external (inter-level) and internal (intra-level) interactions to spatial variability of soils and vegetation. Herbs, low shrubs, and morphologic soil properties turned out to be controlled mainly by the geosystems with the linear size 1200 m, while trees, shrubs, and sediments – by the geosystems with size 2000 m. From 2 to 5 levels of the higher-order geosystems should be considered in order to obtain the proper explanation of spatial heterogeneity.


Introduction
The Interest in the multiscale organization of landscapes is encouraged by the necessity to translate information among hierarchy levels (O'Neill 1988;Wu & David 2002;Burnett & Blaschke 2003;Zhang et al. 2013).The key idea of the concept is that each property of a landscape reflects the superposition of effects generated at various hierarchical levels.Hierarchy theory suggests that multiple scales of pattern will exist in landscapes because of the multiple scales at which processes are acting (Turner & Gardner 2015).Since landscape-ecological studies traditionally focus on the linkages between both natural geocomponents (e.g.bedrock, soil, water, air, vegetation, animals) and spatial units it is crucial to reveal interactions between spatial and nonspatial effects.Limitations upon the interactions between geocomponents are imposed by the differences both in characteristic space scale and time scale inherent to phenomena and processes.Characteristic time was defined as the span of time during which particular processes occur or, in the case of self-regulating systems, the time required for a system to return to a state of equilibrium (Armand & Targul'yan 1976).Theoretically, interactions can occur only between those natural phenomena that have comparable time and space scales (Delcourt et al. 1983;Puzachenko 1986;Shugart 1999).To avoid uncritical application of the famous metaphor "everything is linked to everything" the concept of partial geocomplexes, or partial geosystems, was developed in physical geography (Neef 1967, Sochava 1978, Solon 1999) as a tool to consider strong linkages between groups of certain properties with similar space and time scales.
The definition of landscape as a geosystem is adopted mainly by physical geographers and landscape ecologists (Bastian et al. 2015).Solnetsev (1948) defined landscape as a "genetically uniform territory, with regular and typical repetition of some interrelated combinations of geological structures, landforms, surface and groundwater, microclimates, soil types, phytocoenoses and zoocoenoses".Geosystems form hierarchies but the specific spatial and temporal scales at which interrelationships operate must be recognized.A hierarchy is defined as a system of interconnections wherein the higher levels constrain the lower levels to various degrees, depending on time constraints of the behavior (Turner & Gardner 2015).Since scale multiplicity is inherent in spatial heterogeneity multiscale analysis is imperative for understanding the structure, function and dynamics of landscapes (Wu et al. 2000).
The studies of multiscale organization of landscapes in recent decades has relied upon the concepts of general system theory, theory of complex adaptive systems, hierarchy theory (Carlile et al. 1989;Vasconcelos et al. 1993;Haila 2002;Perry & Enright 2002;Hall et al. 2004;King et al. 2004;Yao et al. 2006;Cushman 2016).Evidences for the multiplicity of scales affecting ecological processes were found on the examples of dead wood patterns and dynamics (Kennedy et al., 2008), spatial patterns of vegetation in wetlands (King et al. 2004), variability of soil properties and processes (Kachanoski 1988), soil humidity as a factor of plant cover pattern (Lookingbill & Urban 2004), influence of topographybased hydrologic features on patterns of woody plant encroachment in savanna landscapes (Wu & Archer 2005), shrub-encroachment due to ecogeomorphic feedbacks (Turnbull et al. 2010), spatial shifts of alpine treelines (Zeng & Malanson 2006), the relationship between the forest vegetation and relief (Chang et al. 2006), the linkage between runoff and vegetation in semiarid landscapes (Wilcox et al. 2003), chemical soil properties as related to soil biomass in mountainous tundra (Oline & Grant 2002), soil attributes across an alpine topographic/snow gradient (Litaor et al. 2002), species distribution in connection to position in a landscape (Fletcher et al. 2016).
The most critical issue in multiscale landscape studies is the quantitative evaluation of the contributions from each scale level to the spatial variability of landscape attributes (Borcard & Legendre 2002;Cushman & McGarigal 2002;Yao et al. 2006;Khoroshev et al. 2007).To compare the contributions of ecological non-spatial and spatial factors a number of techniques have been applied involving autocorrelation, canonical ordination, Fourier-analysis, and general additive mixed modeling (Lookingbill & Urban 2004;Legendre 1993;Couteron et al. 2006;Musio et al. 2007;Cross & Perakis 2011).
The most commonly used concept that relates pattern to process in landscape ecology, the patchcorridor-matrix model, perceives the landscape as a planimetric surface (Hoechstetter et al. 2008).The need to include vertical dimension to landscape pattern models requires consideration for the role of topographic and geological structures (Bolstad et al. 1998;Dorner et al. 2002;Sebastiá 2004;Khoroshev & Aleshchenko 2008;Dragut et al. 2010;Bastian et al. 2015).The focus on the critical significance of abiotic environment in shaping landscape structure has been instrumental for the Central and East European schools in landscape ecology and physical geography, even since the earliest stages (Berg 1915;Solnetsev 1948;Neef 1967;Angelstam et al. 2013;Bastian et al. 2015).The landscape's relief can be interpreted in various ways: (i) as a legacy of former-time processes (e.g.sedimentation in Pleistocene), (ii) as an indicator of geological structure affecting nutrient supply and (iii) as a binding factor for the present-day matter flows (e.g.erosion, water migration, seeds dispersion etc.).The patterns in a landscape surface that are of interest to landscape ecologists may also be interpreted as emergent properties of particular combinations of surface heights and slopes across the study area (McGarigal & Cushman 2005).Weaver & Perera (2004) criticized the models simulating the fate of each pixel independently (termed pseudospatially explicit by Malanson 1996) and argue that accounting for the spatial dependence create more reliable output for analyzing spatial patterns and relating those patterns to ecological processes.
The purpose of our research was to evaluate the contributions of the emergent effects generated within the higher-order geosystem and the effects of internal interactions within the focus-level geosystem.Correlating properties of soil and vegetation were hypothesized to vary in space within the constraints imposed by combination of landforms in some neighborhood.As McGarigal et al. (2016) wrote, determining the right neighborhood size is a major focus of current multiscale habitat selection modeling.We argue that the statistically significant relations within a set of correlating soil-vegetation properties as well as the spatial patterns of the higher-order geosystems may indicate the present-day or former-time processes that govern spatial heterogeneity.According to Malanson et al. (2017), for the analysis of changes in scale, two approaches are common: multilevel (defined as including multiple hierarchical levels of observation in a single model) and multiscale (multiple extents in separate analyses).In this paper, the multiscale analysis is performed to explain the soil-vegetation relationships at the focus level by the processes operating at the higher level of landscape organization.The dominant landscape-ecological models focus either on hierarchical organization of one natural geocomponent (in most cases plant cover, land use or relief, which can be easily detected from DEM and satellite imagery) or describe relations at a single hierarchical level.In contrast, we propose the tool to reveal multiple independent hierarchies based on the interactions between properties of natural geocomponents.If the rule that relates geocomponents to each other is uniform across the entire space of a heterogeneous geosystem, this indicates the emergent property.To prove the existence of such a rule one needs to express the interdependency between geocomponents by any statistically significant quantitative model for a spatial series of landscape units (Khoroshev & Aleshchenko 2008).By this we test the hypothesis as follows: if the combination of spatial units in some neighboring area changes, the properties of the focus unit will change as well (Khoroshev et al. 2007).The size of a neighboring area that affects processes in a focus unit is a matter of analysis.Hence, we face the need to compare the quality of statistical models designed for several hypothetic higher-order geosystems.Thus, the paper focuses on statistical and cartographic methods that could provide answers to the following questions: How to evaluate the contributions to spatial variability created by radial interactions between natural geocomponents and those imposed by higher-order geosystems?How to range the spatial and nonspatial factors affecting low-order landscape units?

Study area
The research was performed in the middle taiga of East-European plain (the southern Arkhangelsk region of Russia) between 60.83N -60.94N and 43.00E -43.38E (Fig. 1).The study area (the Zayachya river basin, 154 km 2 ) is located within the Ustyanskoye plateau and composed of Permian sedimentary rocks.Elevations range from 100 to 175 m a.s.l.Physical environment of the landscape was shaped by morainic and limnoglacial accumulation in the Riss period of the Pleistocene.Flow directions of the first-and second-order streams are affected strongly by the system of lineaments stretching northeastward and northwestward (Khoroshev 2003).Development of the Zayachya terraces with alluvial deposits overlaying Permian marlstones and morainic loams dates back to the late Würm (Avessalomova et al. 2016).
Humus accumulation in soils indicates high content of base cations due to close to surface carbonate morainic loams or Permian marlstone, in most cases -on the valley slopes.Podzolization develops on the poorer substrates where loams are covered by 30-70 cm thick sandy layer deposited in the Pleistocene glacial lakes.Peat accumulation occurs mainly in the oligotrophic mires at the central sections of the flat interfluve areas.Pinetum eriophoro-sphagnosum communities on Histic Gleysols or Histosols on the flat interfluves are gradually replaced by Piceetum myrtillosum forests on Haplic Podzols and further towards the valley slopes -by Piceetum oxalidosum forests on Umbrisols or Rhendzic Leptosols.In the secondary forests Populus tremula dominates after clear-cutting on the most nutrient-rich soils, Betula pendula -on the soils with medium nutrients supply, and Pinus sylvestris -on the poorest soils with thick sandy layer.During recovery succession all of them are replaced by Picea abies, except for the units on sandy terraces where pine preserves its domination.

Data
Field data were collected at 184 forested sample plots 20x20 m distributed relatively evenly across the Zayachya river basin in accordance with proportions of various landforms (flat interfluves, slopes, terraces, floodplains), and variety of succession stages, water and nutrients supply level.At each plot we described landforms (genesis, shape, slope angle, aspect), soil and vegetation cover.Phytocoenosis was described by 5 groups of attributes: abundance of species in the layers of trees (10 variables), shrubs (10 variables), low shrubs (9 variables), herbs (50 variables) and mosses (4 variables).Soils were characterized by three groups of attributes: thickness of genetic horizons (11 variables), Munsell color -hue, value, and chroma -measured over the interval 5 cm (30 variables), and texture measured over the interval 5 cm as a ranks of clay content (10 variables).
We used topographic map (1: 50,000) to compose digital elevation model (DEM) using triangulation technique.We deliberately focused on low-resolution (400 m) DEM to omit from examination fine details of relief (small rills, oxbows, debris slopes etc.) and to concentrate on mesoscale landforms (e.g.valley slopes, ravines, terraces, floodplains) and the related contrasts in soil and vegetation cover.Each sample plot was georeferenced to a DEM pixel which was used as a square operational territorial unit (OTU).ArcView 3.2a and Fracdim (Khoroshev & Aleshchenko 2008) software were applied to calculate 4 morphometric features of relief characterizing intensity of lateral matter transfer: standard deviation of elevations (vertical dissection, VD below), total length of valleys (horizontal dissection, HD below), vertical curvature (VC), and horizontal curvature (HC).These variables were calculated in a moving window with linear dimension ranging from 1200 to 6000 m with step size 800 m, the results being centered on the focal OTU provided by field data.

Multiscale analysis of landscape structure
We elaborated the special procedure called "Multiscale Analysis of Landscape Structure" (MALS) which includes seven steps.
At step 1 reduction of dimensionality of field data was performed.Strong deviation of raw quantitative data from normal distribution required applying non-parametric techniques.Multidimensional scaling was chosen as a method with no restrictions for normality and non-linearity (Cox & Cox 2001), more precisely its nonmetric version (NMDS below) (Legendre & Legendre 1988).For each group of attributes, we calculated non-parametric Gamma correlations and converted correlation (r) matrix to distance (d) matrix by equation d=1-r.We applied multidimensional scaling to distance matrix to calculate the measure of sensitivity of each variable (a 1i ...a 4i ) to the axes (dimensions it terms of NMDS).
To derive appropriate number of axes (NA) we plotted the stress value against numbers of dimensions and analyzed the scree plot obtained.Then we calculated the coordinates of the sample plots (i.e.OTU provided by field descriptions) on the axes of ecological factors from the system of equations: where y i j =a 1i x 1 j +a 2i x 2 j +a 3i x 3 j +a 4i x 4 j (1) -known value of the variable i measured at the sample plot j, a 1i ...a 4i -sensitivity coefficients for the variable i in relation to the axes 1, 2, 3 or 4; x 1 j ...x 4 j -coordinate value for the sample plot j.
The axes were rationalized as the geocomponents properties, that is sensitivity of soils or phytocoenoses to ecological factors according to the requirements of horizons or species, respectively.The coordinates of the sample plots on NMDS axes x 1 j ...x 4 j calculated from (1) are below referred to as D rg , where r -rank, g -geocomponent, e.g.herbs, trees, soil horizons etc.).
At step 2 the data base was composed and processed in order to distinguish the spatial and non-spatial effects in interactions between the geocomponents.The data base included the coordinates of each sample plot on the axes D rg of ecological factors (i.e.properties of soils or phytocoenoses) and morphometric relief features for the hypothetic higher-order geosystems with the abovementioned linear dimensions (1200…6000 m).Normal distribution of the plots coordinates on the NMDS axes values allowed performing principal components analysis (PCA) (Davis 2002) to reduce dimensionality and distinguish the groups of properties with various contributions of the internal (focus-level) and external (higher-level) factors of variability.The appropriate number of principal components was determined from the plot of Eigenvalues vs. Number of principal component considering the sharpest decrease of Eigenvalue.The purpose was to derive two groups of the orthogonal "super-factors".To distinguish them, we analyzed factor loadings for each principal component.The first group of "superfactors" was expected to control the properties D rg that vary in concordance with the relief properties of one or several higher-order geosystems (i.e.D rg and relief properties have factor loading far from zero and close to either maximum or minimum values).This would mean that the properties are sensitive to the external spatial influences generated at the higher level (i.e.inter-level interactions).The second group of "super-factors" was expected to control the properties D rg that vary in concordance with each other (soil-vegetation, sediments-soils relationships etc.) but independently of the higherorder geosystems which was indicated by close to zero factor loadings for relief properties.This would testify the result of the internal non-spatial relations and self-development at the focus scale level (i.e.intra-level interactions).Of course, each property to some extent can be sensitive to both groups of super-factors.
By comparing coefficient of determination (r 2 ) and statistical significance (p-value) of equations (2.1, 2.2) we separated contributions of external and internal factors and concluded whether it is necessary to consider relief of surrounding landscape in order to explain spatial variability of a property under consideration.
Step 4 is aimed at the determination of relevant scale level of a higher-order geosystem to explain variability of the topography-driven attributes.RSR was used to relate the NMDS coordinates to the morphometric features of relief in a square neighborhood of those pixels in which field description was performed.RSR equations were composed separately for each hypothetic size of a higher-order geosystem (namely, 3, 5, …, 15 times as large as the linear dimension of OTU): where D rg -value of a property (coordinates of a sample plot on the NMDS axis), X n -morphometric features of relief in a higher-order geosystem (VD, HD, VC, and HC) (Fig. 2).
Comparison of equations composed for various sizes of a higher-order geosystem provides the opportunities as follows: • to choose the equation with the highest r2 and, hence, to determine "resonance scale level" of a higher-order geosystem that affects the OTU state1 ; • to clarify whether one or several scale levels of the higher-order geosystems are critical for the OTU state; • to identify the set of topographic variables that serve as statistically significant predictors for the OTU state; • to determine whether the topographic predictor correlates with the dependent property positively or negatively.
At step 5 we tested the hypothesis that the combined effect of several higher-order geosystems (step 3) contributes to spatial variability of OTU attributes more effectively than the individual effects of each higher-order geosystem (step 4).It was assumed that we inevitably lost some amount of information about relief at step 2. In contrast, at step 4 we did not lose any information while composing equations separately for each hypothetical "resonance scale level".Suppose that equation for the combined effect of several higher-order geosystems (step 3) explains more variance of a property as compared to any of the equations explaining the individual influence of some higher-order geosystem.This would testify that, despite the loss of information, the combined effect of several higher-order geosystems is more significant than the individual contribution of any higher-order geosystem.If this is not true, one can conclude that no emergent effect of a higher-order geosystems exists.
At step 6 the purpose was to create a series of cartographic models that would express both the interrelationships of geocomponents and the possible multiplicity of higher-order geosystems.We performed classification of pixels (OTUs) by a set of four abovementioned morphometric relief features in a "resonance" square environment determined at step 4. Multistructural organization of a landscape forces us to compose a series of maps, not a single one.Each of these maps follows the deterministic logic: "under the given combination of landforms in surroundings the unit has the soil-vegetation class X".A map in this case shows partial geosystems, or geocomplexes (e.g.water-sensitive, nutrientsensitive etc.), shaped by specific system-forming process with varying degree of manifestation across a landscape.
After that, the OTUs described in field were classified by those attributes of soils and vegetation that appeared to be sensitive to the same ecological factor.Separate classifications were performed for the water-sensitive and nutrient-sensitive attributes.Then, class membership was used in discriminant analysis as a grouping variable.Morphometric features of relief in several "resonance" square environments were taken as response variables.Statistically significant morphometric features were selected using Forward stepwise method (F=1) in Statistica 7.0 software.By this, we predicted the soil-vegetation class with the highest membership probability for each OTU.At this step we followed the probabilistic logic: "under the given influence of the higher-order geosystems the soil-vegetation class X has the highest probability of occurrence".
Step 7 included quantitative evaluation of uncertainty in predicting soil-vegetation classes, since we can easily suppose landforms combinations that permit occurrence of several soil-vegetation classes with similar probabilities.To calculate uncertainty of class prediction we applied the Shannon formula: Landscape Online -supported by the International Association for Landscape Ecology and its community

Khoroshev
Landscape Online 66 (2019) -Page 8 where H -uncertainty of class membership, p i -probability that class i will occur under current combination of landforms in significant neighborhoods.

Driving ecological factors
For most groups of variables NMDS showed smooth decrease of stress values which level off to the right of number of dimensions "4".For example, for the tree layer (10 variables) under 3, 4, 5 and 6 dimensions D-star raw stress accounted for 0.714242 0.2323432, 0.0575684 and 0.0000051, respectively.Therefore, for each group of variables, we derived four NMDS axes (NA=4) to ensure the same number of degrees of freedom in RSR equations.Similarly, four morphometric features of relief were calculated for each size of square neighborhood.RSR modeling indicated that 60-80% of variance of most raw field data was explained by the coordinates on axes nos.1…4.The 1 st and 2 d axes for most geocomponents were interpreted as the sensitivity to either nutrients or water supply.For example, if the thickness of humus horizon was scoring low on the axis and that of podzolic horizon scoring high, the axis was interpreted as a degree of nutrients supply depending strongly on soil-forming sediments.
We found evidence that the nutrients supply has the highest significance for trees, shrubs, and herbs, while water supply -for low shrubs, mosses, thickness of soil horizons, and soil colors.Spearman correlations indicated significant relationships between the D rg values for various geocomponents.For example, the values of the 1 st axis for trees (Pinus sylvestris vs. Picea abies) correlate negatively with the 2 d axis for soil horizons (humus vs. podzolic horizon).This corresponds to response of trees to a gradient from nutrient-poor sandy soils to nutrientrich loamy ones.

Number and size of higher-order geosystems
VD, HD, VC, and HC as the predictors in Eq. ( 3) were repeatedly computed at a range of spatial extents around each location using a moving window with linear dimensions 1200, 2000, 2800, 3600, 4400, 5200, and 6000 m.Coefficients of determinations r 2 in RSR models relating the NMDS axes to relief see in Fig. 3. Evidently, most properties of vegetation and soils are scale-sensitive and exhibit the best response to the relief of geosystems with linear dimension either 1200 m or 2000 m.Higher-order geosystems with linear dimension 1200 m provide constraints for the properties of herbs, shrubs, low shrubs, thickness and color of soil horizons (Fig. 4,  A).Tree and shrub layers and, to some extent, soil texture (i.e., more precisely, genesis of soil-forming sediments) demonstrated the most obvious response to the properties of the larger "2000 m" geosystems which show explicitly step-like organization of relief (Fig. 4, B).

Contributions of internal and external effects: relief-sensitive vs. relief-insensitive properties
Eight mutually independent "super-factors" explained 54% of variance for the whole set of the NMDS axes.Among them the 1 st , 2 d , 6 th and 7 th superfactors (referred to as "external" below) turned out to be highly sensitive to the relief of the higher-order geosystems.For example, super-factor 1 reflects variability of VD and HD in distant environment (2800-6000 m) and indicates the high-order neotectonic blocks, while super-factor 2 controls curvatures in the closest environment (1200-2800 m).The same super-factors describe variability of some properties of plant and soil cover, e.g. the 1 st and the 2 d axes for the herb layer being indicative of nutrient and water supply, respectively.
In contrast, four other super-factors (nos.3, 4, 5, 8) were insensitive to the relief of surroundings.They have high factor loadings for the 1 st axis for trees, the 2 d -for low shrubs, the 1 st axis for mosses, and the 3 d axis for herbs.From this fact we conclude that the internal interactions between these attributes of vegetation layers in situ have higher significance than the broad-scale processes.These superfactors are referred to as "internal" below.Several properties demonstrated sensitivity both to external and internal super-factors.
For each axis we compared proportion of variance explained by four relief-sensitive super-factors vs. four relief-insensitive ones (Table 1).For example, percentage of variance of the 2 d axis for trees (sensitivity to water supply) explained by the reliefsensitive external super-factors accounted for 31%, while relief-insensitive super-factors independently explained 40%.For most axes super-factors explain much less than 100% of variance, 86% for the 1 st axis for herbs being the maximum.
Remember that at previous steps we determined individual contribution of each level of the higherorder geosystems.Percentage of variance of the 1 st axis for herbs (sensitivity to nutrient supply) explained by properties of geosystems with linear dimensions 1200, 2000, 2800, 3600, and 6000 m accounted for 40, 34, 28, 25, and 26% respectively.However, the combined emergent effect of the external superfactors accounted for 59% (Table 1).This greater figure integrates not only the influence of broad-scale The internal super-factors contributed to the spatial variability of the properties linked with feedbacks.For example, the 3 d axis for herbs (Majanthemum bifolium, Oxalis acetosella, Linnaea borealis vs. Heracleum sibiricum, Carex canescens, Equisetum sylvaticum) is slightly sensitive to relief of the "2000 m" geosystem but experience no any emergent effect of higher-order geosystems.At the same time emergent effect of the four internal super-factors is large.The property has poor correlation with abiotic environment but strong -with the forest stand age, canopy cover, canopy height, and tree diameter.
Obviously, the abundance of this group of herbs is related to self-regulation in phytocoenosis during the recovery succession.

Probabilistic cartographic models of partial geosystems
To compose the series of maps of partial watersensitive (WS) and nutrient-sensitive (NS) geosystems we applied the probabilistic approach.We performed classification of units based on watersensitive properties that were linked by strong Spearman rank correlation.8 soil-vegetation WSclasses were identified.Then we used discriminant analysis to test the hypothesis that each WS-class occurs under the specific combination of relief attributes imposed by the higher-order geosystems as well as by slope gradient in the OTU.Forward stepwise method (F=1) was used to identify the significant relief morphometric features that were able to distinguish 8 soil-vegetation WS-classes (Table 2).Only 32% of the units provided with field descriptions were correctly classified (Wilks' Lambda = 0.41441).This is not a surprise since most water-sensitive properties depend much more on the internal factors than on the external ones.The map in Fig. 5 shows the location of 8 soil-vegetation WS-classes that were predicted based on the highest probability of occurrence under the given relief conditions.However, it shows location more or less perfectly only for well-discriminated WS-classes 2, 6, 7, 8 which can occur under the specific relief conditions that are not inherent for any other class.
For example, WS-class 2 (Fig. 5) is typical for the areas with the highest VD due to the outcroppings of Permian marlstones close to water streams.In the well-drained habitats secondary aspen forests dominate, soils have thick humus and very thin podzolic horizons.Absence of Vaccinium myrtillus (commonly, typical for the taiga) indicates low level of groundwater.Hence, water regime resembles the southern taiga rather that the middle taiga.processes (e.g.relief-dependent water migration) but also the non-linear relationships between the relief-dependent properties of soil and plant cover.
In contrast, the single-level dependence holds true for several axes responsible for variability of texture and soil colors (the 1 st and 4 th axes).
Landscape Online -supported by the International Association for Landscape Ecology and its community   Analogously, we identified 8 soil-vegetation NSclasses based on the set of the nutrient-sensitive attributes and performed discriminant analysis in relation to the relief attributes.Percentage of correct classification accounted for 43% (Wilks' Lambda=0.25127).Since most attributes are governed by external factors, the quality of discrimination was higher than that for the watersensitive attributes.The most probable NS-class for each pixel is shown in Fig. 6.NS-classes 1, 3, 6, 7 exhibited the highest percentage of correct discrimination.
For example, NS-class 1 (Fig. 6) to a large extent overlaps with the WS-class 2 (Fig. 5) for the watersensitive properties.However, its areal is less perforated by the other NS-classes.The maps of uncertainty (Fig. 7) show in which locations a set of external and internal factors allows occurrence of several possible soil-vegetation WSand NS-classes.Less uncertainty was detected at the map of nutrient-induced patterns (Fig. 7, B), the lowest uncertainty occurring in the most dissected and the least dissected areas.These locations correspond, respectively to the maximum and minimum possible influence of marlstone.

Discussion
Advance in the formalized techniques of landforms delineation provides the opportunities to predict habitats, forest types, soil patterns etc. based on the models of their various relationships with abiotic template.Our results at step 2 of MALS showed the efficiency of distinguishing two groups of landscape attributes.The first group is governed by the internal interactions between the geocomponents within a geosystem; the second one -by broad-scale processes in the higher-order geosystems.In order to compose the topography-based landscape maps (step 6) correctly, we identified the topographysensitive properties and determined the relevant scale levels of their manifestation (step 4).In our research hierarchical levels were revealed by evaluating linkages between the properties of the focus unit and spatial emergent properties of embracing higher-order geosystem.Each attribute of soil and vegetation can receive significant signal from one or several rank-orders of geosystems simultaneously.Our findings showed evidence that similar relief conditions and corresponding water-and nutrientinduced patterns allow multiplicity of combinations of soil and vegetation properties.Percentage of unexplained variance in the range of 20-40% for most axes (step 3), most likely, indicates that many processes are operating at the other scales, e.g.matter redistribution among micro-landforms.Our previous investigation in the same study area for the water-sensitive properties of plant cover showed that at least 30% of variance unexplained by coarse-scale model (DEM resolution 400 m) can be explained by regression model based on calculations from the more detailed DEM (30 m) (Khoroshev et al. 2013).
The nutrient-induced patterns are more heterogeneous in the interior interfluve areas (Fig. 6) than the water-induced ones (Fig. 5).The latter form a series of belts from the water divides towards valleys.In the deeply dissected sectors, vice versa, nutrient-induced patterns are less diverse.It follows, that neither of two factors can be assigned the higher importance for landscape mapping.Multilevel cartographic models are less uncertain and more effective for the nutrient-sensitive properties than one-level ones (step 7).We identified three principal ways of subordination for the properties of soil and vegetation cover (step 5).First, a property can undergo priority influence of the other geocomponents independently on external factors.Second, a property can be influenced by ecological processes of a single higher-order geosystem that impose strict constraints on the possible range of values.Third, a property can be influenced by an emergent effect of several higher-order geosystems, i.e. by a set of broad-scale processes.
Our research confirmed the hypothesis that the combined effect of several higher-order geosystems provides emergent effect for the low-order landscape unit.The proportion of oligotrophic and megatrophic herb species is governed, to some extent, on the largest higher-order geosystems (3600-6000 m) which were shaped by different rates of Quaternary cover removal depending on rate of neotectonic uplift and resulting in soil enrichment by base cations from underlying marlstones.The lower-level geosystems (1200-2000 m) contributed by means of the present-day migration of nutrients with surface and subsurface waters along the slopes and thalwegs.This kind of results follows the idea that cross-scale interactions can generate emergent behavior that cannot be predicted based on observations at single or multiple independent scales and the interactions may produce nonlinear dynamics with thresholds (Peters et al. 2004) which is one of the most important attributes of the complex adaptive systems (Messier & Puettmann, 2011).
The information translated from the higher-order geosystems is manifested mainly in the drainage conditions that depend on landforms pattern in certain neighborhood and control soil and vegetation properties.Most soil and vegetation properties are governed by not a single scale level but by a simultaneous action of a variety of processes at different spatial and temporal scales, number of which ranges from 2 to 5. Our multiscale cartographic landscape model with the chosen OTU size was more relevant to describe nutrient redistribution.To detect water redistribution in details the chosen OTU size is too coarse.Most likely, it means that present-day redistribution follows to lesser extent the geological patterns but to a much greater extent -fine-scale flows among micro-landforms.Probabilistic landscape mapping (step 7) showed the areas with perfect adaptation of soils and vegetation to abiotic environment as well as areas with high possibility of several stable states.This way of mapping provides the promising opportunity to flatten the contradiction between the concepts of discrete and continual organization of nature.

Conclusions
The MALS procedure showed that in the investigated middle-taiga landscape herbs, low shrubs, and morphologic soil properties are controlled mainly by the geosystems with approximate linear size 1200 m.Trees, shrubs, and to some extent sediments are subject to influence of more broad-scale phenomena in larger geosystems with linear size 2000 m.Though for most properties one of scale levels contributes more than the others, from 2 to 5 levels of higherorder geosystems should be considered in order to obtain proper explanation of spatial heterogeneity.This information provides rationales to objectively delimit landscape units following the criterion of both geomorphologic and soil-vegetation homogeneity.
From the concept of scale multiplicity, it follows that uniform rules of establishing boundaries for the whole set of properties hardly exist.We do not reveal relief hierarchical levels based on DEM beforehand as it is commonly performed in most studies.Instead, we tested a series of hypotheses about possible hierarchical levels of landscape organization using the criterion of statistically significant linkages between soil/vegetation property and relief of a higherorder geosystem.If a researcher finds evidence that the property is controlled by relief-dependent processes it is reasonable to turn to the advanced steps or research focusing on modeling gravityinduced matter and energy flows.If it is not the case the researcher should focus either on modeling dependence on succession stage or on geology.
If one misses the stage of separation of factors contribution he/she faces the risk to mix up effects of independent factors.The proposed methodology is believed to be scale-invariant.It could be applied to any size of OTU and any neighborhood size.

Figure 1 :
Figure 1: Location (left) and digital elevation model (right) of the study area.

Figure 2 :
Figure2: Procedure for determination of relevant size of the higher-order geosystems that control spatial variability of a geocomponent property.The value of property is measured in operational territorial unit (OTU) and expressed as coordinate on the axis of multidimensional scaling (NMDS)

Figure 3 :
Figure 3: External effects of linkages between soil-vegetation properties D rg and morphometric features of the higherorder geosystems relief.Percentage of variance explained by the Response regression surface equation (3).Red outline: internal non-spatial factors are stronger than external spatial ones.Yellow outline: external spatial factors are more significant.

Figure 4 :
Figure 4: Maps of partial geosystems based on the deterministic model: "Properties of geosystems (D rg ) are strictly determined by landforms combination in the square neighborhood".Colors correspond to the types of units embedded into the higher-order geosystems with linear dimensions: A -1200 m (properties of herb, low shrub, shrub layers, soil horizons), B -2000 m (properties of tree and shrub layers).Bold black curve -boundary of the Zayachya river basin Relief properties: 1 -rugged terrains with dense erosion network; 2 -flat terraces and marginal parts of narrow flat interfluves; 3 (B) -slightly concave catchment areas with flat slopes; 4 (A) and 5 (B) green -gentle slightly dissected slopes; 6 (A) -Slightly convex narrow interfluves with gentle slopes; 7 (B) -flat interfluves surrounded by the upper reaches of the valleys; 8 -wide flat poorly drained interfluves; 9 (A) -flat drained interfluves; 10 -narrow flat poorly drained interfluves

Figure 5 :
Figure 5: Multiscale cartographic model of the most probable soil-vegetation WS-classes based on water-sensitive properties.Bold black curve -boundary of the Zayachya river basin Main features of soil-vegetation WS-classes: 1 -dominance of spruce, Vaccinium myrtillus, podzolic horizon in soils with long-term gleyization.2 -dominance of aspen and alder, absence of Vaccinium myrtillus, well-drained soils with humus horizon.3 -high abundance of pine and Pleurozium schreberii on drained sandy soils with short-term gleyization in spring with legacy of ploughing.4 -high abundance of pine, Vaccinium myrtillus, Vaccinium vitis-idaea, on drained sandy soils with temporary gleyization and thick podzolic horizon.5 -high abundance of spruce, Pleurozium schreberii on well-drained soils without gleyization.6 -high abundance of alder and herbs, absence of mosses on well-drained soils with legacy of ploughing.7 -high abundance of spruce, pine, Sphagnum mosses on poorly-drained soils with peat and gley horizons.8 -dominance of pine, Vaccinium myrtillus, Vaccinium vitis-idaea, and Sphagnum mosses on poorlydrained soils with peat, podzolic and gley horizons.

Figure 6 :
Figure 6: Multiscale cartographic model of the most probable soil-vegetation NS-classes based on nutrient-sensitive properties.Classes 1…8 -see in text.Bold black curve -boundary of the Zayachya river basin Main features of soil-vegetation NS-classes: 1 -dominance of alder, aspen, and megatrophic herb species or cultivated areas on humus-rich weakly alkaline soils. 2 -high abundance of mesotrophic and megatrophic herb and shrub species on nutrient-rich weakly acid soils.3 -dominance of pine, Vaccinium vitis-idaea, and boreal herbs on nutrient-poor strongly acid soils.4 -elevated abundance of spruce, pine, boreal mesophytes, and hygrophytes on nutrient-poor weakly acid soils.5 -high abundance of pine, boreal mesophytes and hygrophytes on weakly acid soils with medium nutrients supply.6 -dominance of spruce, pine, and typical boreal species on nutrient-poor acid soils.7 -elevated abundance of aspen, Vaccinium myrtillus, and boreal herbs on weakly acid soils with medium nutrients supply.8 -high abundance of pine and boreal species or hygrophytes on strongly acid nutrient-poor soils.

Figure 7 :
Figure 7: Uncertainty of class membership for the multiscale cartographic models of soil-vegetation classes based on water-sensitive (A) and nutrient-sensitive (B) properties.Bold black curve -boundary of the Zayachya river basin

Table 2 :
Forward stepwise selection of the relief morphometric features significant to distinguish 8 classes of units identified by either water-sensitive or nutrient-sensitive properties Landscape Online -supported by the International Association for Landscape Ecology and its community