Advanced search

Geophysical analysis of landscape polystructures

Full Text:


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.

For citation:

Sysuev V.V. Geophysical analysis of landscape polystructures. GEOGRAPHY, ENVIRONMENT, SUSTAINABILITY. 2020;13(1):200-213.


The identification of multi-scale polystructural geosys­tems and the boundaries between them is among the prin­cipal 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 gradi­ents; 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 geosys­tems are determined by the magnitude and sign of flow divergence. For example, if we consider the behavior of el­ementary 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 co­efficients 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 hier­archy 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, 2006) considers three mechanisms of landscape struc­turing - 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 manage­ment. The development of territorial planning assumes the need to use different methods for selection of spatial units, depending on the objectives of environmental man­agement. For example, the maps showing the structure of biocentricity are necessary to embed nature reserves; those representing positional-dynamic and morphologi­cal structures are essential for industrial facilities (Pozach- enyuk 2006); for agroforestry purposes the catenary dif­ferentiation 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-morpholog­ical 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 in­dicators 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-territori­al 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 (ceterisparibus) could be improved by increasing the num­ber of taxonomic units. Moreover, the landscape typolog­ical approach is the most appropriate for the large-scale work (Tkachev and Bulatov 2002).

Modern landscape ecology is based on the patch mo­saic 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 in­tuition. In addition, the patch mosaic model is consistent with well-developed and widely understood quantitative statistic techniques designed for discrete data (e.g., anal­ysis 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 land­scape structure based on continuous spatial heterogene­ity. 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 etc.), 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 mo­saics, 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 oth­er hand) of radial and lateral processes. This interpretation turned out to be especially productive for regions highly transformed by anthropogenic activity, where the zon­al landscape was preserved in the form of a few "islands” From the point of view of modeling the landscape struc­ture based on structure-forming processes these meth­odological 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-mathe­matical models, and the empirical and semi-empirical pa­rameters need to be introduced into rigorous descriptions of the processes of transfer of matter and energy to over­come 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 parameter­ization, approximation, and similar fitting methods are re­quired.

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 relative­ly 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 im­plemented.
  • Development of regional norms or recommendations for planning spatial relationships between the main land­scape 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 land­scape structure, which is interpreted in the classical defini­tion 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 nu­merical modeling to describe the landscape polystructure by the parameters of major continuous geophysical fields.


The elaboration of any physical-mathematical mod­el 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 physi­cal laws. It is essential that physical laws and their param­eters be applied relevant to the description of the struc­ture-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 val­ues 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 digi­tal 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 pro­cedures. 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 construct­ed, the pixel size could be 10x10 m. However, the pixel size depends on the resolution of aerial photo or satellite im­age 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 pa­rameters of the state of elementary material objects (pix­els) to distinguish uniform areas. Theoretical description of the geostructure, i.e. stationary (for a certain time inter­val) state of a dynamic geosystem, begins with identifying morphometric parameters (MP) which describe the force fields that determine the main structure-forming process­es. Morphometric formalization of the Earth's surface in the gravitational field was systematized in (Shary 1995). Logi­cally meaningful association of morphometric parameters includes three groups characterizing: 1) the distribution of solar energy - the dose of direct solar radiation (daily, an­nual), aspect and illumination of slopes; 2) the distribution and accumulation of water under the influence of gravity - the specific catchment area and the specific dispersive area, depth of β-depressions and the height of β-hills, the slope gradient; 3) the mechanisms of matter redistribution un­der 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-com­posite) 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 geo­physical fields and structure-forming processes.

The physical meaning of the morphometric parame­ters is quite clear. For example, the dose of radiation char­acterizes the potential input of direct insolation. Slope as­pect 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 con­vexity/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 conserva­tion 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, hu­midity 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 gravitation­al, electromagnetic, and other parameters, are very prom­ising 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 pa­rameters of the gravity and insolation fields acquires a fun­damental geophysical meaning. In this case, the concept of polystructural landscape organization is absolutely logi­cal: 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 dif­ferentiation 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 classifica­tion of landform elements according to a known landscape structure. On the other hand, the change in the set of pa­rameters and their numerical values makes it possible to model landscape structure changes under the influence of climate change, neotectonic events, etc. A rigorous land­scape approach is needed for such modeling that allows identification of the main factors of differentiation and ex­clusion 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 be­longs to the end moraine belt of the last Valdai (WGrm) Ice Age in the northwestern East-European Plain. The loamy moraine deposits with residual carbonates reach a thick­ness 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 san­dy 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 bio­logical 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 x 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 geomorpho- metric 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 corre­spond 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 multi­tude of DEM pixels. Geometrically, the closer two such vec­tors in the parameter space, the less different the param­eter 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 geo­metrically 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 iden­tified. 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 au­tomatic classification of objects is reduced to the problem of diagonalizing such a communication matrix. The au­tomatic classification can be understood as a geometric task of distinguishing "dense” concentration of points in a certain space. Such geometric approach allows elaborat­ing the methods of solving the task of automated classi­fication 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 decod­ed 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 de­scriptions, 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.)


Fig. 1. Distribution of the stand volume (m3/ha) along the transects crossing typical landscapes of the Valdai National Park: I - landscape of the ridge-hilly moraine-kame plain on carbonate moraine loams often covered with silt- sandy loam deposits of low thickness; II - landscape of ridge-hilly kames-eskers plain on sandy-loam sediments; III - landscape of outwash glacio-fluvial plain with sand ridges. The alphabetic indices characterize the localities ("mestnost"). In the lower part - elevations a.s.l. (m) and a schematic lithological section: 1 - moraine deposits, locally overlaid by kames, 2 - glacio-fluvial sands, 3 - peat deposits


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 per­formed using the discriminatory methods.


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




Table 1. Identification of decoded classes basing on the results of the continuous forest inventory along the transect

Class (Fig. 2, B)

Number of sites




Open water



Water plants



Meadow swamp



Closed spruce forests



Spruce-pine and pine forests



Mixed and spruce-pine forests



Boggy pine forests



Small-leaved and boggy pine forests



Bogs with stunted pine trees

Correlation of the compiled stand map and the calcu­lated type of site conditions with scale 1:10000 forest com­partment 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 di- chotomous grouping of landform elements (DEM pixels) based on the parameters of geophysical fields and the state of landscape cover. Independent morphometric pa­rameters were chosen based on preliminary analysis (dig­ital modelling): the annual dose of direct solar radiation, elevation, slope gradient, horizontal and vertical curvature, specific catchment area, as well as numerical data of Land- sat 7 spectral channels and the NDVI index.


Fig. 3. The structure of the NTC based on the relief classification using the geophysical field gradients and the space image data Landsat-7.

1 - moraine ridges and kame hills with loamy Umbric Albeluvisols (sod-podzolic soils) under Piceetum oxalidosum forests; 2 - summits of kame hills and ridges with sandy Umbric Podzols under Pinetum cladinosum, Pinetum cladinoso-hylocomiosum and Pinetum herbosum forests; 3 - foot slopes of the hills and flat concave hollows with Umbri-Endogleyic Albeluvisols and Umbric Albeluvisols under mixed forests; 4 - river and lake terraces with Endogleic Umbrisols and Distri-Fibric Histosols under spruce forests and mixed forests; 5 - dune ridges and sandy hills with Umbric Albeluvisols under pine forests; 6 - flat and convex upland bogs with deep Histosols under sparse pine forests; 7 - river floodplains with Endogleic Umbrisols under flooded meadows; 8 - steep slopes of hills with Umbrisols under coniferous forests; 9 -anthropogenically modified and anthropogenic landscapes (roads, power lines, quarries, farmland, nurseries and residential areas)


Classification results depend significantly on the num­ber 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­tions and a substantial number of preliminary numerical experiments (Sysuev 2003, 2014). For example, if the heat (energy) supply of the territory is modeled with the same weight coefficients for all insolation parameters, the ob­tained groups do not satisfy the landscape structure re­vealed in the field studies. The problem is that at first and subsequent levels of the dichotomic classification the class­es of relief surface are distinguished, first, against the expo­sure 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 veri­fication of the classification of relief elements with the se­lected 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. Similar­ly, 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 identifi­cation 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 signifi­cance 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 ex­perimental 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 creat­ed a landscape map basing on the classical method. The comparison showed that the landscape map resulting from numerical experiments is sufficiently accurate in re­producing the boundaries of the NTC independently ob­tained 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 integrat­ed 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 de­termined 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 parame­ters, namely slope gradient, specific catchment area, hori­zontal 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-men­tioned parameter to select cells with values exceeding a predetermined threshold that are considered to be poten­tial source points. At the second step, channels from the given sources are drawn, and the sources which have tran­sit flow from higher elevations are removed. At the third step, channels smaller than a certain minimum length are cut off (Fig. 4).


Fig. 4. Automated procedure for cutting off the drainage channels shorter than the critical length (catchment basin of the Loninka River, the Valdai National Park) in MapWindow GIS using the TauDem module (Tarboton et al. 1991)


The process can be easily adjusted by changing the lim­it values of the catchment area and the minimum length of the drainage channels. The resulting array of morphomet- ric characteristics with geographic reference to watersheds of various orders is a characteristic of landscape-hydrologi­cal systems (LGS), or catchment geosystems (Sysuev 2014).

Application of typological approach to obtain land­scape 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 fluviogla- cial deposits. Modern drainage streams are mostly tempo­rary, with poorly developed alluvial relief and poorly pro­nounced valley forms. They occupy the drained inter-ridge depressions, which are ancient valleys of the glacial melt- water runoff. The inter-ridges depressions lacking the ac­tive 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 al­gorithm for constructing a drainage network in SAGA num­bers 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 meth­od 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 forma­tion noticeably reduced.


Fig. 5. The structure of watercourses on a digital elevation model (A) and catchment geosystems (B) of the Upper Mezha River (Central Forest Reserve, Tver Region). The higher is the territory, the lighter is the DEM tone; a specific catchment area (SCA) was selected as a critical value for calculating the order of watercourses (A). Watersheds as delineated by zero SCA values (B)


The geosystem order is the same that the order of a watercourse (Fig. 5B). A special category - zero order geo­systems - was introduced for lower rank complexes that do not have a pronounced drainage watercourse. About 400 such geosystems have been identified within the ter­ritory under study. The maximum fourth order geosystem is the catchment of the Mezha River in the vicinity of the Fedorovskoye village (Table 2).


Table 2. The main statistical parameters of catchment geosystems

Order of catchment

Number of catchments

Average area, km2

Minimum area, km2

Maximum area, km2


Average discharge, June, l/s




































The dependence of the average values of a catch­ment area Y on its order X is described by the equation Y=b0*(X+1)b1, where the parameter values are b0=0.419665, b1=2.526742, and the model reliability is R = 0.99977. Areas of zero order geosystems have the approximately lognor­mal distribution.

The qualitative and quantitative characteristics of runoff for any watercourse depend on the physical and geograph­ical features of its catchment. To identify these dependen­cies, a map of the landscape structure of catchment geosys­tems within the river basin was compiled (Fig. 6).


Fig. 6. Map of the landscape structure of catchment geosystems. The table legend is below the map. The color indicates the type of vegetation cover, the intensity of color shows the moisture gradient distribution. Dots indicate the sites of hydrological and hydrochemical testing


The methodology for compiling the map of the land­scape 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 struc­ture 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 imag­ery. Eleven classes were also identified, and two of them were then combined to reflect anthropogenically modi­fied 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 cal­culated 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 sim­plify the map compilation, classes with the area exceeding 10% of a geosystem total area were left for further consid­eration. 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 water­sheds boundaries.

The technique simplifies the application of landscape characteristics and differs from methods used for identify­ing 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 calcu­lating 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).


Fig. 7. The discharge dynamics for the measuring sites (red points) in the Upper Mezha River basin (the Tver region) during 2001 summer, and the decrease in active runoff area as a result of the gradual depletion of spring flood water. Schematic hydrographs of runoff are shown below in small graphs; the black area is the period of discharge averaging for Table 3


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 de­pressions are usually not well pronounced in the relief and their depth could be only dozens of centimeters. Howev­er 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 geo­systems 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 wa­ter period associated with the low drainage from soil and ground, the flow continued only in the largest streams and the Mezha River itself.


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)

Order catchment

Number of measuring sites

Discharge, l/s


Electrical conductivity, μC/cm.





























No runoff





No runoff





No runoff
















No runoff





No runoff





No runoff














Analysis of the catchments parameters and hydrological measurements showed a close relationship between the structure and functioning of geosystems. This provides op­portunity to calculate the water flow discharge basing exclu­sively on the geosystems structure and precipitation data. Hydrological functioning and water protection zoning of geosystems

The calculation of surface runoff from a priori topo­graphic DEM data can be performed in different GIS sup­porting hydrological procedures. For example, in order to calculate the water flow in SAGA (Olaya 2004), a sufficiently large number of individual catchment parameters, such as the Shezi-Manning coefficient of surface roughness ("Man­ning's n” - MN) and the coefficient of soil influence on the intensity of surface runoff ("Curve number” - CN), are re­quired in addition to general parameters (elevation, slope gradient, specific catchment area, etc.).

The values of parameters must be assigned to each pix­el of the model, which is objectively possible only with re­liance 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 pa­rameters of the hydrological model aimed at calculating the runoff rates within the Loninka River basin.


Fig. 8. Distribution of hydrophysical parameters in the basin of the Loninka river for modeling water flow on the basis of the landscape structure (SAGA GIS). А - Curve number (CN); B - Manning's n (MN)


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 sur­face, 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 ob­served 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 chan­nels with water flow.


Fig. 9. Calculation of surface runoff rates in the Loninka River basin by means of the SAGA GIS. The precipitation intensity is 10 mm/h, CSS = 1, MFT = 180, CDT = 360


The results of runoff simulation using various parame­ters 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 embank­ment (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 precipita­tion 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 val­ues 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 pre­cipitation intensity. Thus, the different location of the ob­servation 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 sug­gest 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.


Fig. 10. Comparison of the measured and calculated flow rates in the Loninka River at the gauging stations.

1- field measurements; 2 - model calculation, precipitation intensity 0.0 mm/h; 3 - model calculation, precipitation intensity 0.66 mm/h; 4 - model calculation, precipitation intensity 10.0 mm/h; 5 - model calculation, intensity of precipitation 100.0 mm/h. The abscissa shows the numbers of gauging stations downstream the river sources


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. Never­theless, 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 iso- chrones of flow delay time were calculated using the SAGA GIS (Fig. 11, A). However, this method of calculating is diffi­cult 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 devel­oped 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 sec­ond-order catchments, and so forth (Fig. 11, B).


Fig. 11. Running time of surface water runoff, hours. A - to the control section of the Loninka River, according to the SAGA GIS algorithm; B - to the river for each particular point of the channel, according to the cascade algorithm



The formation of landscape structures is described in traditional empirical concepts using the geomorphomet- ric parameters of geophysical fields, i.e. gravity and inso­lation. The concept of landscape polystructure becomes physically defined: by choosing the main structure-form­ing 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 typo­logical approach makes it possible to obtain a hierarchy of natural-territorial complexes (facies - urochishche - mest- nost - 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 cur­vature, specific collection area and specific dispersive area; B-depression depth) and insolation (direct solar radiation dose; aspect and illumination of slopes) fields are consid­ered 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 gra­dients; horizontal/ planar curvature is the divergence of streamlines; vertical curvature is a derivative of the steep­ness factor along the streamline; the dose of direct solar ra­diation 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 defini­tions and preliminary numerical experiments. The need for professionally correct and justified selection of the physical state parameters, i.e. principal structure-forming process­es, 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 redis­tribution of water over the surface in the gravitational field (slopes, specific catchment area, horizontal and vertical curvature). Such classification makes it possible to identi­fy the contours of various-order catchment geosystems in accordance with the Horton-Strahler-Tokunaga scheme.

The structural parameters obtained from the typolog­ical description of landscapes allow simulating the hydro­logical functioning of catchment geosystems with satis­factory accuracy for particular types of water protection zoning. 


1. Akbari Kh.Kh., Bondar Yu.N., and Sysuev V.V. (2006). Indicative properties of the stand in the landscapes of the edge zone of the Valdai glaciation. Proceedings of Moscow University, series 5 Geography (6), 59-66. (in Russian with English summary)

2. Antipov A.N. and Fedorov V.N. (2000). Landscape and hydrological organization of the territory. Novosibirsk: Nauka. (in Russian with English summary)

3. Armand A.D. (1988). Self-organization and self-regulation of geographical systems. Moscow: Nauka. (in Russian)

4. Beven K. J. (2012). Rainfall-runoff modelling: the primer. Chichester: Wiley-Blackwell.

5. Chow V.T., Maidment D.R., and Mays L.W. (1988). Applied Hydrology. New York: McGraw-Hill.

6. Chow V.T. (1959). Open-channel hydraulics. New York: McGraw-Hill. Dodds P.S. and Rothman D.H. (1999). Unified view of scaling laws for river networks. Phys. Rev. E, 59(5), 4865–4877.

7. Dyakonov K.N. and PuzachenkoYu.G. (2004). Theoretical positions and directions of modern landscape studies. In: K.N. Dyakonov, ed. Geography, society and environment. V. II. Functioning and present-day state of landscapes. Moscow: Gorodets, 21–36. (in Russian).

8. Dyakonov K.N. (2008). Basic concepts and intention of landscape studies. In: N.S. Kasimov, ed. Geographical scientific schools of Moscow University. Moscow: Gorodets, 348-381. (in Russian) Forman R.T.T. (2006). Land mosaics: the ecology of landscapes and regions. Cambridge University Press, Cambridge.

9. Kleidon A. (2010). Life, hierarchy, and the thermodynamic machinery of planet Earth. Physics of Life Reviews, 7(4), 424–460.

10. Khoroshev A.V. (2012). Geographical concept of landscape planning. Proceedings of the RAS. Geographical series (4), 103-112. (in Russian with English summary)

11. Khoroshev A.V. (2016). The large-scale organization of the geographical landscape. Moscow: Fellowship of scientific publications KMK.

12. Landscape Planning: Principles, Methods, European and Russian Experience. (2002). In: A.N. Antipov and A.V. Drozdov, ed. Irkutsk: Publishing House of the Institute of Geography SB RAS. (in Russian with English summary)

13. Likens G.E. and Bormann F.H. (1995). Biogeochemistry of a forested ecosystem. New York: Springer-Verlag.

14. McGarigal K., and Cushman S.A, (2005). The gradient concept of landscape structure. In: J. Wiens and M. Moss, ed. Issues and perspectives in landscape ecology. Cambridge: Cambridge University Press, 112–119

15. McGarigal K., Tagil S., Cushman S.A. (2009). Surface metrics: an alternative to patch metrics for the quantification of landscape structure. Landscape Ecol (24), 433–450. DOI 10.1007/s10980-009-9327-y

16. Olaya V. (2004). A gentle introduction to SAGA GIS. Edition 1.2.

17. Orlova I.V. (2014). Landscape-agroecological planning of the territory of the municipal district. Novosibirsk: Publishing House of the SB RAS. (in Russian)

18. Pozachenyuk E.A. (2006). Territorial planning. Simferopol: Publ. House of the TNU. (in Russian)

19. Rulev A.S. (2008). Landscape planning of agroforestry land improvement. Environmental planning and management, 1(6), 25-28. (in Russian with English summary)

20. PuzachenkoYu. G. (2014).Organization of a landscape. In: K.N. Dyakonov, V.M. Kotlyakov, T.I. Kharitonova, ed. Issues in Geography. Vol. 138. Horizons of landscape studies. Moscow: Kodeks, 35-64. (in Russian)

21. Reteyum A. Yu. (1975). Physical-geographical zoning and the allocation of geosystems. Issues in Geography, 98, 5-27.(in Russian)

22. Shary P.A. (1995). Land surface in gravity points classification by a complete system of curvatures. Mathematical Geology, 27(3), 373–390.

23. Solntsev N.A. (1948).The natural geographic landscape and some of its general rules. Proceedings of the Second All-Union Geographical Congress. Moscow: State Publishing House for Geographic Literature, 258-269. (in Russian). See also in: J.A. Wiens, M.R. Moss, M.G. Turner, D.J. Mladenoff, ed., (2006). Foundation papers in landscape ecology. New York: Columbia University Press, 19-27.

24. Solntsev V.N. (1997). Structural landscape studies. Moscow: Faculty of Geography MSU. (in Russian)

25. Sysuev V.V. (2003). Morphometric analysis of geophysical differentiation of landscapes. News of the Academy of Sciences. Series geographical, (4), 36-50 (in Russian with English summary)

26. Sysuev V.V. (2014). Basic concepts of the physical and mathematical theory of geosystems. In: K.N. Dyakonov and V.M. Kotlyakov, ed. Issues in Geography. Vol. 138. Horizons of landscape studies. Moscow: Kodeks, 65-100. (in Russian with English summary)

27. Sysuev V.V., and Solnetsev V.N. (2006). Landscapes of the edge zone of the Valdai glaciation: classical and morphometric analysis. In: K.N. Dyakonov, ed. Landscape science: theory, methods, regional studies, practice. Proceedings of the XI International landscape conference. Moscow: MSU Publishing House, 249-252. (in Russian)

28. Sysuev V.V., Sadkov S.A., Erofeyev A.A. (2011). Basin principle of functional zoning: modeling of the structure and drainage of catchment geosystems based on a priori data. In: K.N. Dyakonov, ed. Landscape science: theory, methods, regional studies, practice. Proceedings of the XI International landscape conference. Moscow: MSU Publishing House, 101-105. (in Russian)

29. Tarboton D.G., Bras, R.L., Rodriguez-Iturbe, I. (1991). On the extraction of channel networks from digital elevation data. Hydrological Processes, 5(1), 81-100.

30. Tkachev B.P., and Bulatov V.I. (2002). Small rivers: state-of-the act and ecological problems: Analit. review. Novosibirsk: Publishing House of the SB RAS.

31. Troendle C.A. (1985). A variable source area model for stormflow prediction of forested watersheds. In: M.G. Anderson and T.P. Burt, ed. Hydrological forecasting. Chichester: John Wiley & Sons, 431-495.

32. Turner MG (2005). Landscape ecology: what is the state of the science? Annu. Rev. Ecol. Evol. Syst. (36), 319–344

33. Turner M., Gardner R.H. (2015). Landscape Ecology in Theory and Practice. Pattern and Process. New York: Springer.

34. Viktorov A.S. (1986). Landscape pattern. Moscow: Mysl. (in Russian with English summary).

About the Author

Vladislav V. Sysuev
Lomonosov Moscow State University
Russian Federation

Faculty of Geography, 


For citation:

Sysuev V.V. Geophysical analysis of landscape polystructures. GEOGRAPHY, ENVIRONMENT, SUSTAINABILITY. 2020;13(1):200-213.

Views: 262

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 License.

ISSN 2071-9388 (Print)
ISSN 2542-1565 (Online)