GEOPHYSICAL ANALYSIS OF LANDSCAPE POLYSTRUCTURES

. The objective identification of landscape cover units is very important for sustainable environmental management planning. The article proposes a method-algorithm for describing the formation of landscape structures, which is based on the classic landscape analysis and applies the parameters of geophysical fields. The main driving forces of all structure-forming processes are the gradients of gravitational and insolation fields, parameters of which were calculated using the digital elevation models and the GIS-technologies. A minimum number of principal parameters are selected for typological and functional classification of landscapes. The number and importance of parameters were identified basing on the results of numerical experiments. Landscape classifications elaborated on the basis of standard numerical methods take a fundamental geophysical value. In this case, a concept of polystructural landscape organization is logical: by selecting different structure-forming processes and physical parameters, different classifications of landscapes could be elaborated. The models of geosystem functioning are closely related to their structure through boundary conditions and relations between parameters. All models of processes and structures are verified by field experimental data obtained under diverse environmental conditions.


INTRODUCTION
The identification of multi-scale polystructural geosystems and the boundaries between them is among the principal problems of landscape research. The fundamentals of non-equilibrium thermodynamics show the principles on which the classifications of natural-territorial complexes (NTC) should be based. In accordance with the Onsager bilinear equation, classifications should account for: 1) the system-forming flows; 2) the force fields and their gradients; 3) the phenomenological coefficients of generalized conductivity. The fields of gravity and insolation are the most common for any geosystem. Selection (classification, integration) of geographical objects by the parameters describing the geophysical fields and their gradients leads to the identification of geosystems according to the flows of matter and energy. This is a functional approach to the identification and investigation of geosystems, which is being developed in the works by Armand (1988), Reteum (1975), and other. In such a case the boundaries of geosystems are determined by the magnitude and sign of flow divergence. For example, if we consider the behavior of elementary water volumes in a geopotential field, we obtain a hierarchy of catchment geosystems (river basins) which corresponds to the formalized Horton-Strahler-Tokunaga schemes. The drastic change in the phenomenological coefficients is the basis for the classification of NTC according to the principle of homogeneity (typological approach in accordance with N. A. Solntsev's theory (Solntsev 1948)). Considering the spatial distribution of plants and animals in both geopotential and other physical fields -insolation, chemical, thermodynamic etc. allows obtaining the hierarchy of ecosystems (biogeocenotic systems) and their spatial distribution. Such approaches to classification of geosystems are mutually complementary, and should not be contrasted. For example, V.N. Solntsev (Solntsev 1997(Solntsev , 2006 considers three mechanisms of landscape structuring -geostationary, geocirculatory and biocirculatory, which could operate individually or simultaneously. The objective identification of landscape cover units is essential for the planning of sustainable nature management. The development of territorial planning assumes the need to use different methods for selection of spatial units, depending on the objectives of environmental management. For example, the maps showing the structure of biocentricity are necessary to embed nature reserves; those representing positional-dynamic and morphological structures are essential for industrial facilities (Pozachenyuk 2006); for agroforestry purposes the catenary differentiation should be considered, and different types of units, such as urochishche, mestnost etc., and catchments should be used (Rulev 2008). In landscape-agroecological planning the preference is given to genetics-morphological approach, the principal units being the urochishche and a group of urochishches (Orlova 2014). An example of the use of landscape planning is water protection zoning (Landscape Planning… 2002). Complicated environmental situations and a variety of conflicts between land and water users are characteristic of the water protection zones. On the other hand, the most complex landscape-hydrological systems (LHS) are presented in the areas adjacent to water bodies. The combination of landscape and hydrological indicators at the basin level is used to refine the calculation of hydrological parameters or to assess the distribution of LHS (Antipov and Fedorov 2000). The concept of LHS as applied to landscape hydrology is a set of natural-territorial complexes (NTC) similar in runoff formation conditions. The NTC only indirectly characterizes the catchment area with similar capacitive features; therefore the experimental observations in each small river basin are necessary for an accurate estimate of the LHS area. However, in most cases LHS are formed from a set of NTC. The calculation accuracy (ceteris paribus) could be improved by increasing the number of taxonomic units. Moreover, the landscape typological approach is the most appropriate for the large-scale work (Tkachev and Bulatov 2002).
Modern landscape ecology is based on the patch mosaic paradigm, in which landscapes are conceptualized and analyzed as mosaics of discrete patches (Forman 2006;Turner et al. 2005).The strength of the patch mosaic model lies in its conceptual simplicity and appeal to human intuition. In addition, the patch mosaic model is consistent with well-developed and widely understood quantitative statistic techniques designed for discrete data (e.g., analysis of variance). Developing this rather limited approach to environmental considerations McGarigal and Cushman (McGarigal and Cushman 2005) introduced the «landscape gradient» model, as a general conceptual model of landscape structure based on continuous spatial heterogeneity. Based on the continuous characteristics of the earth's surface, obtained from DEM (slope, topographic wetness index, topographic position index, normalized difference vegetation index NDVI etс.), the «landscape gradient» model is constructed using the statistical characteristics of patch mosaic (patch density, largest patch index, edge density, mean patch area, area-weighted mean patch area, coefficient of variation in patch area, mean patch shape index, etc.). Statistical characteristics of patches are called landscape metrics (McGarigal et al. 2009). However, these are metrics of just the sizes and shapes of patches in mosaics, and not of complex geosystems, such as landscapes of any dimension and hierarchy. The structure (pattern) is understood, first of all, as a combination of interacting spatial elements with their area, configuration, orientation, neighborhood, connectedness or fragmentation (Turner and Gardner 2015), i.e. close to the concept of "landscape pattern" (Viktorov 1986). The structure is interpreted as an indicator (on the one hand) and a condition (on the other hand) of radial and lateral processes. This interpretation turned out to be especially productive for regions highly transformed by anthropogenic activity, where the zonal landscape was preserved in the form of a few "islands". From the point of view of modeling the landscape structure based on structure-forming processes these methodological approaches are not entirely correct. Amount of empirical data, composite indices with a fuzzy (intuitive) physical meaning cannot be used directly as parameters of the equations of mathematical physics. As a result, a gap arises between landscape-ecological and physical-mathematical models, and the empirical and semi-empirical parameters need to be introduced into rigorous descriptions of the processes of transfer of matter and energy to overcome it. For example, GIS SAGA software (Olaya 2004) is supplemented by the TOPMODEL hydrological unit (Beven 2012), which describes the migration of moisture based on the Darcy equation, with a significant number of empirical parameters that are difficult to determine, so parameterization, approximation, and similar fitting methods are required.
Planning decisions based on the landscape-ecological analysis could become more reasonable if the following problems are solved (Landscape Planning… 2002): • Identification of quantitative indicators describing both the structure and the functioning and development of a landscape. We need objective indicators which are relatively easy to calculate. • Development of classifications of landscapes according to their sustainability, vulnerability, suitability and capacity for particular types of environmental management. • Identification of quantitative characteristics of the land-scape self-organization, or at least qualitative description of the regional developments of this process. • Search for relationships between the spatial structures of natural and socio-economic systems. • Determination of minimum natural ranges within which ecological stabilization of cultural landscapes could be implemented.
• Development of regional norms or recommendations for planning spatial relationships between the main landscape elements (by area and configuration).
These tasks are currently relevant. In our work, possible ways of implementing some of these tasks are given.
The aim of the work is to justify the choice of the least number of objective parameters characterizing the landscape structure, which is interpreted in the classical definition of the Moscow State University School of Landscape Sciences. In fact, this is a synergistic task of determining the main parameters of structure-forming processes. The method proposed in the article allows application of numerical modeling to describe the landscape polystructure by the parameters of major continuous geophysical fields.

COMMON PRINCIPLES OF THE LANDSCAPE STRUCTURE MODELS
The elaboration of any physical-mathematical model begins with basic axioms and postulates. The principal point is the identification of elementary material objects (particles, points) forming the system and the assignment of independent variables and functions of the system's states. Further, it is necessary to adopt a number of binding postulates so that it is possible to apply particular physical laws. It is essential that physical laws and their parameters be applied relevant to the description of the structure-forming landscape processes (Sysuev 2014).
To describe geosystems, it is first of all necessary to substantiate the potentials of the main geophysical fields that determine the structure-forming processes, and then to formalize the description of elementary geosystems and hierarchical invariants of geosystems. The quantitative values of spatially distributed physical parameters of the state of landscapes could be obtained: 1) from digital elevation models (DEM) -morphometric parameters describing the gradients of the gravity and insolation fields; 2) from digital remote sensing data -parameters of the Earth's surface cover; 3) from field and laboratory measurements, and 4) during special experiments.
The space of geographical coordinates is provided by the construction of a digital elevation model (DEM). Pixels of the 3D DEM are elementary material points (similar to the material points in theoretical mechanics), from which the NTC structure is synthesized using the formalized procedures. DEMs are constructed to achieve the maximum resolution of a particular hierarchical level of geosystems. For example, if a regular-grid DEM based on the contours from a detailed topographic map M 1:10 000 is constructed, the pixel size could be 10x10 m. However, the pixel size depends on the resolution of aerial photo or satellite image as well. So, the resolution of Landsat images (30x30 m) allows us to identify the NTC of just urochishche level.
Differentiation of geographical space could be realized using various mathematical methods (cluster analysis, neural networks, etc.). However, we need numerical parameters of the state of elementary material objects (pixels) to distinguish uniform areas. Theoretical description of the geostructure, i.e. stationary (for a certain time interval) state of a dynamic geosystem, begins with identifying morphometric parameters (MP) which describe the force fields that determine the main structure-forming processes. Morphometric formalization of the Earth's surface in the gravitational field was systematized in (Shary 1995). Logically meaningful association of morphometric parameters includes three groups characterizing: 1) the distribution of solar energy -the dose of direct solar radiation (daily, annual), aspect and illumination of slopes; 2) the distribution and accumulation of water under the influence of gravitythe specific catchment area and the specific dispersive area, depth of B-depressions and the height of B-hills, the slope gradient; 3) the mechanisms of matter redistribution under the influence of gravity -horizontal, vertical and mean curvature, slope gradient, height (Sysuev 2014). It should be noted that the minimum number of simple (non-composite) state parameters is selected that independently describe the gradients of geophysical fields generating the main structure-forming processes. Thus, the selection (classification) of geosystems is not carried out according to relief elements, or elementary locations, or patches, or other spatial units, but directly by the parameters of geophysical fields and structure-forming processes.
The physical meaning of the morphometric parameters is quite clear. For example, the dose of radiation characterizes the potential input of direct insolation. Slope aspect and gradient are the components of the geopotential vector gradient. Horizontal curvature is responsible for the divergence of flow lines. Vertical curvature is the derivative of the steepness factor and characterizes the slope convexity/concavity. These parameters are directly included in the equations of mathematical physics. Specific catchment area shows the area from which suspended and dissolved substances could be transported to a surface element. It is a substance balance parameter (included in conservation laws), and a component of a number of related indices (water flow capacity index, erosion index, etc.).
The state of the Earth's surface (vegetation cover, snow cover, soil cover, etc.) is detected from the digital data of space spectral image and from the related indices (e.g., the normalized difference vegetation index, snow index, humidity index -NDVI, NDSI, NDWI, etc.). The most important sources of data are field studies, which also allow verifying the interpretation of the state of the covers. In addition to the traditional complex methods of landscape science, it is necessary to use automated complexes for recording geophysical parameters of the lowest atmospheric layers, natural waters and soils. Methods of applied geophysics, based on measuring the spatial distribution of gravitational, electromagnetic, and other parameters, are very promising as well.
The parameters for describing the structure are chosen in accordance with the classical approaches of landscape studies. All formal algorithms for selecting the smallest and the higher order units of relief surface based on the parameters of the gravity and insolation fields acquires a fundamental geophysical meaning. In this case, the concept of polystructural landscape organization is absolutely logical: by selecting different structure-forming processes and physical parameters, different classifications of landscapes could be elaborated. Let us demonstrate the approach through particular cases below.

Typological model of the landscape structure
Typological approach allows obtaining a hierarchy of classical NTC (facies -urochishche -mestnost -landscape) in accordance with N.A. Solntsev's theory (Solntsev 1948). The parameters are chosen in accordance with the wide-ly-known definitions. For example, "An elementary NTC -facies ... is confined to one element of mesorelief; this territory is homogeneous in terms of its three principal characteristics: the lithological composition of the rocks, the slope aspect and gradient. In this case, the total solar radiation and atmospheric precipitation coming to the surface are the same within any part of it. Therefore, one microclimate and one water regime are formed, ... one biogeocenosis, one soil unit and a uniform complex of soil mesofauna" (Dyakonov and Puzachenko 2004). As it follows from the definition, elementary NTCs could be selected by the parameters of solar radiation and water distribution over the surface; more precisely, by the distribution of the gradients of insolation and gravity fields. Thus, the classical definition already requires the description of the NTCs differentiation using the theory of field and the morphometry of the Earth's surface. The classification results essentially depend on weight values and number of parameters. By changing the latter, it is possible to optimize the classification of landform elements according to a known landscape structure. On the other hand, the change in the set of parameters and their numerical values makes it possible to model landscape structure changes under the influence of climate change, neotectonic events, etc. A rigorous landscape approach is needed for such modeling that allows identification of the main factors of differentiation and exclusion of derivative or dependent variables. Automatically obtained classes of landscape cover require identification and verification of their physical content.
The investigated territory of the Valdai National Park is located in the central part of the Valdai Upland, which belongs to the end moraine belt of the last Valdai (Würm) Ice Age in the northwestern East-European Plain. The loamy moraine deposits with residual carbonates reach a thickness of 25 m in the ridges of the Crestets end moraine belt and overlie the glacio-fluvial sands of preceding stages. Locally, the moraine is covered by kame silty-sandy loam sediments. The fluvioglacial plains are covered with peat bog sediments, the massifs of which are separated by sandy eskers. Peatland systems are connected by streams and the Loninka and the Chernushka rivers. Such a variety of landforms and sediments causes the high degree of biological and landscape diversity within the study area. The digital elevation model (DEM) was constructed on horizontally a detailed topographic map with a scale of 1:10 000 using the regular grid method with 28 × 28 m pixel size attached to the Landsat-7 DTM+. The pixel size makes it possible to reliably distinguish NTCs of the locality (mestnost) rank in the study area with dimensions of about 10x10 km. Based on the values of the height and size of pixels (vertical and horizontal steps), the main geomorphometric parameters were calculated, i.e. slope, dose of direct solar radiation, aspect and illumination of slopes, specific watershed area and dispersive area, B-depression depth, B-hill height, horizontal and vertical curvature. GIS ECO (P. Shary), GIS FractDim (G. Aleschenko, Yu. Puzachenko), GIS DiGem (O. Conrad), and GIS SAGA (V. Olaya) were used for calculations. The best result of building a smoothed DEM and morphometric parameters of the studied territory was shown by GIS ECO.
Next, a matrix (database) was built: the rows correspond to the relief surface elements (DEM pixels), and the columns to the parameters (MP) describing the state of an element (height, geomorphometric parameters, as well as digital values of the brightness of the Landsat-7 DTM+ channels and NDVI). The parameters describing the same surface element have different physical meanings and are not comparable in dimensions and size. They are therefore normalized and reduced to a standard form. The resulting matrix is ready for any algebraic operations.
The row vectors of the data matrix characterize a multitude of DEM pixels. Geometrically, the closer two such vectors in the parameter space, the less different the parameter values for both objects. This suggests that the closer two vectors in the parameter space, the more "similar" and less distinguishable the corresponding objects in many of their other properties, not only for those included in the data matrix. Therefore, if it is possible to distinguish a geometrically sufficiently isolated "group" of vectors close to each other within the set of all object vectors, then a class of objects with similar internal properties could be identified. The Euclidean metric of the distance between the corresponding object vectors was used as a measure of the proximity of two objects. If the proximity measure function is selected and calculated, the matrix of communication between objects is thus constructed, and the task of automatic classification of objects is reduced to the problem of diagonalizing such a communication matrix. The automatic classification can be understood as a geometric task of distinguishing "dense" concentration of points in a certain space. Such geometric approach allows elaborating the methods of solving the task of automated classification of spatially distributed objects, such as elements of relief surface in the DEM, remote sensing data, etc. The number of objects in the geographic data matrix is very large and could reach dozens and hundreds of thousands. Efficient recurrence algorithms could be applied in these cases which provide that the calculations are performed with only one successive object (or with one row of the corresponding communication matrix) at each step. In our work, the FractDim software was used to classify the relief according to the MP matrix.
The classes of vegetation cover automatically decoded from space spectral image were verified on the basis of field data (Akbari et al. 2006). Some results are shown in Fig. 1. The investigated transect (integrating a series of analogous transects) was about 5 km long; leveling was performed at 5 m interval; integrated descriptions were made at sample plots 20 m apart. Along with complex descriptions, a complete forest inventory was carried out at 239 plots of 20x20 m area. The forest inventory included strip enumeration of trees with measurement of standard parameters of each tree and the sample plot (composition, stratum, height, diameter, age, crown attachment height, stocking, canopy density, underwood, advance growth, grass-shrub cover, type of soil, type of station, etc.) The map of automatically decoded classes (Fig. 2, A) compiled from a priori geophysical data (annual dose of direct solar radiation, slope gradient, numerical data of Landsat 7 spectral channels and the NDVI) was refined on the basis of field data geo-referenced to space image. The map of vegetation cover (Fig.2, B) shows classes according to the field-based verification. Interpolation of continuous forest inventory data on the studied territory was performed using the discriminatory methods.
Correlation of the compiled stand map and the calculated type of site conditions with scale 1:10000 forest compartment map for this territory showed that the simulation results are significantly more detailed compared to the standard forest inventory data.
The scheme of the geosystems structure (NTC at the urochishe level) is shown in Fig. 3. We used successive dichotomous grouping of landform elements (DEM pixels) based on the parameters of geophysical fields and the state of landscape cover. Independent morphometric parameters were chosen based on preliminary analysis (digital modelling): the annual dose of direct solar radiation, elevation, slope gradient, horizontal and vertical curvature, specific catchment area, as well as numerical data of Landsat 7 spectral channels and the NDVI index.
Classification results depend significantly on the number of parameters and their weight values. By changing them, it is possible to optimize the classification of relief elements according to the known (assumed) landscape structure. The selection of parameters was carried out in accordance with the classical landscape science defini-GEOGRAPHY, ENVIRONMENT, SUSTAINABILITY 2020/01

Fig. 2. Identification of the physical content of landscape cover classes in the Valdai National Park according to a priori information (A) and the vegetation cover identified by continuous forest inventory along the transect (B). Legend see in Table 1. The dots show the transect location
Class (Fig. 2 (Sysuev 2003(Sysuev , 2014. For example, if the heat (energy) supply of the territory is modeled with the same weight coefficients for all insolation parameters, the obtained groups do not satisfy the landscape structure revealed in the field studies. The problem is that at first and subsequent levels of the dichotomic classification the classes of relief surface are distinguished, first, against the exposure of slopes, which is not true for a significant part of the territory, which is a swampy landscape of the fluvioglacial outwash plain. On the other hand, the cooling properties of swampy landscape are not taken into account, and even vast massifs of upland bogs with lakes are not identified. To obtain more correct classification, the weight coefficient of the slope parameter was increased, which improved the classification of the relief surface. At each stage, the verification of the classification of relief elements with the selected values of weight coefficients was checked by the method of discriminative analysis. At all levels, according to the values of F-criterion, the distinguished classes differ statistically reliable. Thus, the leading role of waterlogging of the territory, obvious to landscape scholars, could be numerically expressed by the value of weight coefficients of slope, i.e. the parameter of gravity field gradient. Similarly, the significance of other geomorphometric parameters was substantiated. Because of numerical experiments, the level of required numerical classification was also objectively revealed. The 5-level dichotomy adopted at the beginning with identification of 32 classes turned out to be excessive. The size of urochishche (simple and complex) corresponding to the main mesostructures of the relief of the studied territory is the most adequately displayed at the 4th level of classifica-tion. Moreover, as can be seen in Fig. 3, a good number of identified classes are automatically combined into larger groups according to close-valued colors of a palette.
The need for an objective justification of the significance of the state parameters and the level of numerical classification dichotomy is emphasized by the crucial role of the landscape approach, which makes it possible to single out the role of individual factors (structure-forming processes) of NTC differentiation in specific geographical conditions.
Finally the distinguished classes were characterized using the parameters obtained during field study of experimental landscape transects. The characteristics of the grass-shrub and soil cover, and the lithological structure were extrapolated by discriminative methods.
In addition, an independent team of researchers created a landscape map basing on the classical method. The comparison showed that the landscape map resulting from numerical experiments is sufficiently accurate in reproducing the boundaries of the NTC independently obtained by the classical field methods (Sysuev and Solntsev 2006).
Revealing spatial structures in such a way is a process of synthesis, since the material points (pixels) are integrated into elementary natural-territorial complexes (at a given hierarchical level) according to the selected parameters of main geophysical fields and the state of covers.

Functional model of geosystem structure
The functioning of low-order geosystems is largely determined by water flows. Hence, the classification is aimed at the construction of a hierarchy of catchment geosystems according to the morphometric parameters describing the redistribution of water in the gravity field. These parameters, namely slope gradient, specific catchment area, horizontal and vertical curvatures, determine the boundaries of the divergence zones and the convergence of streamlines. The hierarchy of catchment geosystems is determined in accordance with the Horton-Strahler-Tokunaga scheme (Tarboton et al. 1991;Dodds and Rothman 1999). The automated algorithm for identifying drainage channels on the basis of the GIS raster layers involves three main steps. First, we use the digital map of an above-mentioned parameter to select cells with values exceeding a predetermined threshold that are considered to be potential source points. At the second step, channels from the given sources are drawn, and the sources which have transit flow from higher elevations are removed. At the third step, channels smaller than a certain minimum length are cut off (Fig. 4).
The process can be easily adjusted by changing the limit values of the catchment area and the minimum length of the drainage channels. The resulting array of morphometric characteristics with geographic reference to watersheds of various orders is a characteristic of landscape-hydrological systems (LGS), or catchment geosystems (Sysuev 2014).
Application of typological approach to obtain landscape characteristics of catchment geosystems could be demonstrated by the example of the Upper Mezha River basin (Central Forest Reserve, Tver Region). The southern taiga spruce forests with broad-leaved species, shrubs and nemoral herbs dominate this southwestern part of the Valdai Upland (the East European Plain). The territory is a combination of flat ridges and inter-ridge depressions, a weakly dissected plain composed of glacial and fluvioglacial deposits. Modern drainage streams are mostly temporary, with poorly developed alluvial relief and poorly pronounced valley forms. They occupy the drained inter-ridge depressions, which are ancient valleys of the glacial meltwater runoff. The inter-ridges depressions lacking the active runoff are occupied by upland and transitional bogs. The upper reaches of the Mezha River are located in the study area. The river channel is winding within the valley with a wide floodplain (50-100 m) winding channel. The methods for the DEM construction and calculation of parameters are similar to those described above. The algorithm for constructing a drainage network in SAGA numbers the selected segments of watercourses in the order of R. Shreve. To obtain data on the orders of catchments according to the Straler classification, the channels of the first, second, and subsequent orders were sequentially cut off, until maps of watercourses with a breakdown in order were received (Fig. 5 A). An analysis of the use of the method of the specific catchment area critical values in different landscape showed that as the age of the relief increases (in the series: secondary glacial plain → periglacial plain of the marginal zone of the Valdai IceAge → end moraine zone of the Valdai Ice Age) the average area of watercourse formation noticeably reduced.
The geosystem order is the same that the order of a watercourse (Fig. 5B). A special category -zero order geosystems -was introduced for lower rank complexes that do not have a pronounced drainage watercourse. About 400 such geosystems have been identified within the territory under study. The maximum fourth order geosystem is the catchment of the Mezha River in the vicinity of the Fedorovskoye village ( Table 2).
The dependence of the average values of a catchment area Y on its order X is described by the equation Y=b0*(X+1) b1 , where the parameter values are b 0 =0.419665, b 1 =2.526742, and the model reliability is R = 0.99977. Areas of zero order geosystems have the approximately lognormal distribution.
The qualitative and quantitative characteristics of runoff for any watercourse depend on the physical and geographical features of its catchment. To identify these dependencies, a map of the landscape structure of catchment geosystems within the river basin was compiled (Fig. 6).
The methodology for compiling the map of the landscape structure of geosystems is as follows. Maps of relief structure and vegetation cover were created for the study area by means of classification analysis of the digital relief model and satellite imagery of Landsat 7. The relief structure map was created using "with training" classification according to the digital relief model. Using the K-means method, eleven classes were distinguished, which reflect the differentiation of the territory into flat sections and slopes with convex (spurs) and concave (hollows) sections. A map of the structure of vegetation cover was created through the interpretation of the Landsat 7 satellite imagery. Eleven classes were also identified, and two of them were then combined to reflect anthropogenically modified territories (village, roads, fields and hayfields, deposits etc.). At the next stage the number of pixels corresponding to a particular class of vegetation and topography was calculated for each geosystem. As a result, the areas occupied by the main classes of relief and vegetation within each geosystem were estimated. These data were summarized in a single table. Then, the percentage of the area occupied by each class was calculated for each geosystem. To simplify the map compilation, classes with the area exceeding 10% of a geosystem total area were left for further consideration. A generalized matrix was then obtained and the map of the landscape structure of catchment geosystems of various orders was compiled (Fig. 6). Thus, the landscape structure is described within physically determined watersheds boundaries.
The technique simplifies the application of landscape characteristics and differs from methods used for identifying the landscape-hydrological systems (LGS). According to Antipov and Fedorov (2000), the area of LGS varies from year to year, from season to season, and from day to day. Therefore, the selection of NTCs that have similar state in relation to runoff in particular time period (for example, floods) is not simple.
In mathematical modeling of surface runoff, there is a clearer concept of an "active runoff area" that changes during the process (Troendle 1985). We suggested calculating the active runoff area by the value of the territory elevation excess above the mouths of the particular order watercourses for a certain period of time (Fig. 7).
The relationship of functioning and the structure of geosystems could be analyzed through the gradual decline in flow hydrographs (Fig. 7, Table 3). During spring floods, and after heavy prolonged rains, practically all first order catchments function, i.e. a characteristic surface runoff along the hollow-like depressions is observed. These depressions are usually not well pronounced in the relief and their depth could be only dozens of centimeters. However they are well marked by moist grasses and humus-gley soils. So for the beginning of June (the final stage of 2001 spring flood), the active runoff area for the first-order geosystems with an excess above the watercourses mouths of <1.0 m was 0.2 to 1.7 ha. As a rule, such areas accounted for about 60-70% of the total area of a geosystem (Fig. 7). For some geosystems, particularly those with flat surface, the active areas could not estimated. During the low water period associated with the low drainage from soil and ground, the flow continued only in the largest streams and the Mezha River itself.
Analysis of the catchments parameters and hydrological measurements showed a close relationship between the structure and functioning of geosystems. This provides opportunity to calculate the water flow discharge basing exclusively on the geosystems structure and precipitation data. Hydrological functioning and water protection zoning of geosystems The calculation of surface runoff from a priori topographic DEM data can be performed in different GIS supporting hydrological procedures. For example, in order to calculate the water flow in SAGA (Olaya 2004), a sufficiently

Table 3. Average water discharges and mean values of some hydrochemical indicators for the sections of different order catchments in the Mezha River basin (the Tver region)
large number of individual catchment parameters, such as the Shezi-Manning coefficient of surface roughness ("Manning's n" -MN) and the coefficient of soil influence on the intensity of surface runoff ("Curve number" -CN), are required in addition to general parameters (elevation, slope gradient, specific catchment area, etc.). The values of parameters must be assigned to each pixel of the model, which is objectively possible only with reliance upon the information on the landscape structure. In Fig.8, the numerical values of MN, taken from the standard Chow tables (Chow 1959), are depicted in accordance with the typological structure of landscape (see Fig. 3). These data are the basis for setting the spatially distributed parameters of the hydrological model aimed at calculating the runoff rates within the Loninka River basin.
The average precipitation intensity was assumed 0.0, 0.66, 10.0, or 100.0 mm/hour in numerical experiments for the calculation of runoff. In addition, Channel Site Slope (CSS) parameters, runoff characteristics (slope surface, Mixed Flow Threshold -MFT, and Channel Definition Threshold -CDT) and some other parameters were subject to changes in calculations. Numerical modeling has shown that even tabular values of MN, CN, CSS, MFT and CDT, not adapted to taiga wetlands, reveal significant features in the distribution of surface water runoff in various geosystems (Fig. 9). Extremely low runoff values (<0.01 m/s), were observed for most of the basin. Higher rates are characteristic only for the channels of streams and rivers (0.025-0.2 m/s) and the runoff increases up to 2 m/s within certain sections of the Loninka River. The pattern of runoff rates distribution is quite realistic, since the catchment of the Loninka River is a flat swamped hummocky sandy plain, cut by rare channels with water flow.
The results of runoff simulation using various parameters of Average Rain Intensity (ARI) and Channel Site Slope (CSS) revealed some regularities. In all cases, the increase in ARI caused higher runoff, for example, at the source point located near the drainage pipe under the railway embankment (the source of the Loninka River). This site is highly modified by human activities, and, consequently, it has low MN values (0.025) contributing to the surface runoff, and high CN (98) impeding water infiltration into the soil. Thus, the changing intensity of precipitation successively leads to the change in surface runoff rates. That is, in this section of the river the sensitivity of runoff to the precipitation intensity is great, although the flow rates are not very high. Most other observation sites are located in natural forest and mire landscapes. These sites are characterized by high MN values (0.5-0.9) and low CN (30-40). High values of MN prevent surface runoff, while low CN favor active infiltration. As such, the runoff rates decrease substantially at these sites and weakly respond to changes in the precipitation intensity. Thus, the different location of the observation sites in terms of landscape structure results in significant differences in runoff characteristics. Low regular runoff rates at zero precipitation intensity confirm the high capacity of flat over-moisturized catchment to accumulate water and regulate runoff in geosystems.
Verification of model calculations was carried out in the field. Experimental measurements of flow rates and discharge of the Loninka River at the gauging stations suggest the following. In all cases, the predicted rates differ from the measured ones, but the calculated values were not so far from the real ones as it was expected (Fig.10). The closest results were obtained for the precipitation intensity of 10 mm/h, although lower-intense precipitation is more probable.
More accurate simulation results could be obtained by adjusting the values of model coefficients (tabular values for non-waterlogged rivers were used for calculations). A more detailed DEM could also be useful. It turned out that the channel width of 1.0-1.5 m, and the 30 m pixel size does not allow delimiting valleys, as well as the microrelief which is very important for runoff from the flat plains. On the other hand, errors in the measurement of flow rates are quite possible for flat, boggy meandering channels, often blocked by forest debris and beaver dams. Nevertheless, under the lack of information, the values obtained during GIS modeling could become a basis for predicting runoff values in areas where the direct measurements are labor-consuming or otherwise impossible.
Let us demonstrate the possibility of water protection zoning of geosystems based on modeling the structure of catchment basins and the runoff from their areas using a priori data. An important environmental characteristic of the processes in the catchment basin is the delay time of water flowing to the river or control stations. The isochrones of flow delay time were calculated using the SAGA GIS (Fig. 11, A). However, this method of calculating is difficult to use to predict the time of pollutants arrival from side streams, which is important for taking measures to localize pollution before it gets into the river channel. We developed a modified cascade algorithm (Sysuev et al. 2011) to calculate the running time to first-order channel for each first-order catchment, combine them into second-order catchments, and then calculate the running time for each second-order catchment before merging together all second-order catchments, and so forth (Fig. 11, B).

CONCLUSION
The formation of landscape structures is described in traditional empirical concepts using the geomorphometric parameters of geophysical fields, i.e. gravity and insolation. The concept of landscape polystructure becomes physically defined: by choosing the main structure-forming processes and their principal parameters different classifications of landscapes could be elaborated. Formal mathematical algorithms of selecting surface relief units acquire fundamental geophysical meaning if combined with the state parameters. Implementation of the typological approach makes it possible to obtain a hierarchy of natural-territorial complexes (facies -urochishche -mestnost -landscape); implementation of the approach of the hydrological functioning of the landscape results in the hierarchy of catchment geosystems; implementation of the classification approach for parameters and normalized coefficients of remote sensing data produces the structure of vegetation cover.
Geomorphometric values describing the gradients of gravity (height, slope, horizontal, vertical and average curvature, specific collection area and specific dispersive area; B-depression depth) and insolation (direct solar radiation dose; aspect and illumination of slopes) fields are considered to be the parameters of physical state of individual units of relief surface, i.e. the DEM pixels, which form the geosystems. The state parameters were preferred due to their simple form and direct description of physical fields For example, slope is the modulus of the geopotential gradients; horizontal/ planar curvature is the divergence of streamlines; vertical curvature is a derivative of the steepness factor along the streamline; the dose of direct solar radiation is the relative amount of incoming energy, etc. The state parameters are also independently included into the description of structure-forming processes. Digital remote sensing data are also physical parameters of the state of individual units of relief surface and geosystems.
Parameters of the typological model of the landscape structure are selected in accordance with classical definitions and preliminary numerical experiments. The need for professionally correct and justified selection of the physical state parameters, i.e. principal structure-forming processes, as well as their weights suggests the crucial role of the landscape approach.
The functional model of the landscape structure is based on morphometric parameters describing the redistribution of water over the surface in the gravitational field (slopes, specific catchment area, horizontal and vertical curvature). Such classification makes it possible to identify the contours of various-order catchment geosystems in accordance with the Horton-Strahler-Tokunaga scheme.
The structural parameters obtained from the typological description of landscapes allow simulating the hydrological functioning of catchment geosystems with satisfactory accuracy for particular types of water protection zoning.