This article presents a review of the long-term efforts and results of studies of linear fold kinematics, initiated by Pr Vladimir Beloussov (who was born in 1907, and died in 1990). These methods are based on the kinematic mechanisms of structure formation at different scales, from individual grains up to the scale of the entire fold belt system. The set of these methods is called multi-rank strain analysis (Yakovlev, 2008a). The models link the geometry of natural objects to variation laws of an object's geometry, including its formation mechanisms and value of strain. Obviously, having the possibility to distinguish among various mechanisms and measuring geometric parameters, we can define the value of spatial shortening. It is important that a solution is searched from small structures to large ones, at any scale and without omissions. The final result is the reconstructed geometry of large structures, which is independent from any geodynamic concept. Within the hierarchic system of objects, several levels are recognized (Rebetsky et al., 2004), which differ by their number of layers.

Level 1: intralayer strain of grains and inclusions; 2: separate folds (layer and pairs of competent and incompetent layers); 3: domains (package of layers, see example on Fig. 1C); 4: structural cells (term by Goncharov (1988)) or tectonic modules; 5: tectonic zone; 6: folded system (such as the Greater Caucasus); 7: the whole folded thrust belt. In this hierarchic system (Rebetsky et al., 2004), the sedimentary cover is usually at level 4, crust at level 5, and the whole lithosphere at the highest level 7.

Boundaries between objects in the hierarchic system are given in such way that it remains coherent within the limits of acting mechanisms. The relationships between the different hierarchic levels will be given in the article for different scales, using a wide range of structures from the Greater Caucasus example (Fig. 1).

Single folds are recognized at level 2 of the hierarchic system. They occur both at small and large scales. Two mechanical models are used to determine the value of shortening within a single fold:

• • from a mechanics point of view, folding of single viscous layers incorporated in a less viscous environment has a clear theory of formation, which is used for computation within the finite-element model (Hudleston and Stephansson, 1973). This model predicts the whole geometry of the fold and makes it possible to construct a diagram, linking geometric parameters to the amount of shortening during the process (Yakovlev, 1978, 2008b). Measuring the same parameters in a natural structure (Fig. 2) allows one to define the value of shortening and to derive the viscosity contrast from the diagrams. Application of the method to 78 structures of the Chiauri zone in the Greater Caucasus gave 56% of average shortening (−ɛ = (L₁−L₀) × 100/L₀) with data scattering from 25 to 82%;
• • for multilayer folds a kinematic iterative model of shape transformation was built. This model combines buckling and flattening for the competent layers and shearing with rotation and flattening for adjacent incompetent layers. Based on this model, inclination of the layers in the flanks and ratio of competent layer thickness between the flanks and hinge of a given fold allow one to estimate the shortening in the direction perpendicular to the axial plane (Yakovlev, 2002). Application of the iterative model to 36 folds of the Chiauri zone gave results (average shortening 57%, scatter from 27 to 83%) similar to those obtained from the modeling of single layers observed in the same area. For 8 local structures the two types of folds were found, which makes it possible to compare results of the two methods. The comparison demonstrates that the kinematic model of multilayers folds is, in general, close to the model of formation of single viscous layer folds.

Multilayers folds were found in the Sochi region in the North-West Caucasus in the footwall of the Vorontsov nappe (Yakovlev et al., 2007, 2008) near the foreland/hinterland boundary (location on Fig. 1B). Two possible deformation models of the structure were tested, corresponding to different geodynamic situations: pure shear (Ramsay and Huber, 1987) with an horizontal shortening linked to the Greater Caucasus development and simple shearing along the sole of the nappe. Two parameters, i.e., the inclination of folds axial plane and the value of shortening perpendicular to the axial plane, have shown significantly different trends of development within the two models (Fig. 3).

These parameters were measured in 39 folds of the Vorontsov nappe and are close to the simple shear model prediction. A difference in 20 degrees between regression lines could be explained by the nappe body dipping 20 degrees to the north at the measurement spot. The initial general dip should have been to the south, which would agree with a gravitational origin for the nappe emplacement. Additional detailed analysis of sedimentary sequences in adjacent tectonic zones and the analysis of the geological map of the Caucasus have revealed the following sequences for the major tectonics events in this region: (1) main folding inside the nappe itself occurred 35–25 Ma ago; (2) southward gravitational sliding of the nappe toward the basin over gently dipping submarine relief took place 22–15 Ma ago; (3) main uplift of the mountain range has started at 14 Ma (Yakovlev et al., 2008). This kinematic modeling at the hierarchic level 2 has proved that gravitational sliding is the most likely cause for the Vorontsov nappe formation, which is therefore not directly related to crustal shortening and mountain building processes. For this reason the 15 km long southward displacement of this 1.5-km thick plate should not be included in the overall shortening value of the Greater Caucasus, which must be deduced from higher hierarchic levels.

The hierarchic level 3 (Rebetsky et al., 2004) integrates several folds which have similar parameters in a unique domain: (1) axial plane dip; (2) value of shortening; and (3) local dip of the envelope of the small scale folds. It is important that characteristics of folds from a single layer in a limited volume can be extrapolated to a pack of layers (0.3–1.5 km thick), which make it possible to describe the strain of the whole volume using the concept of the strain ellipsoid (Fig. 4). The boundaries of individual tectonic domains are defined parallel to the axial plane of folds. These boundaries are located at the place where fold parameters change or where sharp changes of lithology occur. Fault plains are using as boundaries also, if significant displacement of blocks takes place.

This hierarchic level (Rebetsky et al., 2004) is a major one to incorporate the local scale natural observations (Fig. 1C) in the regional cross-sections crossing the tectonic zones. This work step is implicitly or explicitly used during any cross-section balancing procedure (Boyer and Elliott, 1982). The value of shortening of folds in individual tectonic modules should be taken from results of the methods described above. However, mostly this parameter is defined from the fold interlimb angle, an assumption also developed in the works (Suppe, 1983) on kink-bends folds and Ghassemi et al. (2010) for “similar folding”. When doing this it is assumed that layer length at flanks is constant, which is in general proved in studies of natural structures in the Caucasus at least. From the three above retained parameters, kinematical models were issued allowing the realization of diagnostic diagrams (Yakovlev, 2001).

For studying natural linear folding it is very important to have detailed structural cross-sections. The North-West Caucasus was studied in details by T. Giorgobiani and Y. Rogozhin (Giorgobiani and Zakaraya, 1989; Sholpo et al., 1993). Eleven cross-sections compiled by them served as the basis for several types of studies (Fig. 1B, C). In these cross-sections nearly 250 domains were defined. Three afore-mentioned geometrical parameters of the structures were measured (Fig. 4B), i.e., the dip angle of axial planes (AX), the dip of the envelope plane of the small scale folds (EN), and the value of fold shortening (SH). Using diagnostic diagrams that take into account the three measurements (Yakovlev, 2001), it was found that near-fault fold measurements make a specific cluster (Fig. 5B). Such a structure could be linked to a superposition of simple shearing along inclined zone and horizontal pure shear (Fig. 5C). Two kinematical models were compiled with different dips of the shear zone (45° and 20°) and different model increments concerning the amount of simple shear and pure shear. It was found that the best fit for the morphology of those local structures is obtained with a shear zone that initially dips at 45° and a succession of 6° shearing along the shear zone and 1% of horizontal pure shear (horizontal shortening) (Yakovlev, 2005). Distribution on a map of such structures with northern or southern vergence demonstrates regular features and gives a key for tectonic zoning of the region. This study implies that a study at hierarchic level 3 can provide a key to understanding the geodynamics at a higher hierarchic level. Small-scale folding is probably related to pairs of shear zones with standard shear angle, occurring under horizontal compression, rather than slip along tilted bedding (20° as initial tilt) or detachment operating at the top of the basement. It should be noted that the structures of high hierarchic level cannot be easily identified during geological mapping or as a result of structural section compilation.

The description of finite deformation in individual tectonic domains by a strain ellipsoid makes it possible to reconstruct the initial position of a single domain and of the whole cross-section (Yakovlev, 2002, 2009). To achieve this, we add to three measured geometrical parameters the next two parameters (Fig. 4B), i.e., (4) the length of the domain section line; and (5) the dip of section line. By so doing, the final state of the domain can be described as a strain ellipsoid (aliquot ellipse for 2D linear folding), the orientation of which is constrained by the initial geometry of the layers. Its volume is limited horizontally by a specific segment of the profile, which is characterized by a given length and a given dip.

The main problem to reconstruct the pre-folded state of a given domain relates to the selection of the backward kinematic operations, which will bring the fold architecture back to horizontal layers and the aliquot ellipse back to a circle. We have used 3 distinct operations: (a) a rotation back up to the horizontal position of the fold envelope plane; (b) a horizontal shearing to restore the fold axial planes back to a vertical position; and (c) a horizontal extension up to fold straightening, when the ellipse will turn back into a circle. During these operations both the length and the orientation of the profile segment are changed. For instance, in the pre-folded state, this line will have a new length and a new inclination relative to the horizontal layering (Fig. 6). It is known (Ramsay and Huber, 1987) that the order of the three operations will influence the result. However, the same result of reconstruction can be obtained with another order of operations. For instance, we can operate first the rotation back up to obtain a vertical position of the axial plane, followed by a vertical simple shearing up to horizontal position of envelope plane and finally a stretching; but, of course, other parametric values should then be used for these operations.

Another very important parameter is the “initial depth”, that is, the location of the tectonic domain in the stratigraphic column. Each domain extends between two initial depths along the stratigraphic column. Fault surfaces are considered as planes belonging to one of the domains and initial pre-folded orientation of this plane is calculated using the same sequence of the three basic kinematic operations. Then, the “initial depth” difference between two adjacent blocks will define the amplitude of the vertical offset whereas the inclination angle of the plane in its pre-folded state will define the amplitude of the horizontal displacement. A full restoration of the pre-folded state of the whole cross-section is then possible by merging the un-folded states of individual domains and displacing them along the faults.

The basic concept here is similar to the goals of the balanced cross-section procedure (Dahlstrom, 1969), but we propose a systematic procedure to incorporate the data from one hierarchic level to the other. Furthermore, we would like to outline on the following modelling example that restoration of the geometry of sedimentary cover structures is not restricted to line length balancing methods. It requires also involving the deformation obtained at a smaller hierarchic level and to assess correctly the tectonic (external) component of the shortening.

At the scale of a structural cell (hierarchical level 4) Goncharov (1988) developed a pure kinematic model for calculations of field displacement within a convective structural cell based on hydrodynamic equations. Such a model makes it possible to calculate strains (strain ellipses) for small volumes (domain, hierarchical level 3) from the intensity of convection and of the external shortening of the whole structural cell. This study makes it possible also to model kinematically the combination of folding of the sedimentary cover together with internal deformation within individual domains, and to solve a problem which does not have analytical solution in mechanics (Fig. 7). In particular, we can get a model for quasi-buckling at the scale of the sedimentary cover by calculating the combination of parameters which keep constant the whole length of the central layer of the structural cell (la-lb-lc on Fig. 7). If we trace the changes of the length of segments formed by several domains (such segments not being located this time in the middle of the structural cell), their shortening values will strongly differ from each other. Therefore, if we want to evaluate tectonic shortening within complex folded structure, we have to study the structural location of the structural cells in order to integrate the shortening of these objects in the whole structure.

Information on thickness of the cover together with data on pre-folded length of all domains makes it possible to define the initial geometry of these structural cells. Both the pre-folded state coordinates of structural cells boundaries and the initial thickness of sedimentary cover are required to build series of 2D restored profiles.

For the example of the North-West Caucasus region, data on the thickness of all stratigraphic horizons from the Lower Jurassic to the Upper Eocene were collected for the 244 tectonic domains of the 11 cross-sections (Fig. 1B). This allows one to reconstitute the whole sedimentary cover. Interpolation between the entire set of cross-sections provided a 3D description of the architecture of the sedimentary cover of the North-West Caucasus.

Another important problem relates to the onset and duration of the folding process and its relation to the mountain building process for the Greater Caucasus, the precise timing of which is still being debated. The general reviews of the information on the whole structure of the Greater Caucasus can be found in Leonov (2007), and Saintot et al. (2006). It is often considered that the main folding episode occurred before or during the Oligocene and that mountain topography has started to be formed during the Sarmatian. According to the age of angular unconformities on the edges of the Greater Caucasus and by the age of the first conglomerates, these two processes may be indeed diachronous. The lack of refolding of erosional surfaces confirms this idea (Nesmeyanov, 1992). However, there is more information, which demonstrates the long duration of folding formation and its simultaneity with the relief formation at least in a few local structures.

To restore the former architecture of the Greater Caucasus, it is necessary to have some ideas about its tectonic evolution. As a first hypothesis, we have considered here three successive evolutionary stages: Stage 1 is characterized by the end of the deposition of Pre-Oligocene series; during stage 2, still before the Sarmatian, the whole folded structure is formed without any erosion, whereas these folded structures became subsequently eroded during mountain growth processes. Accordingly, stage 3 relates to the modern post-orogenic structure. At present it is not possible to calculate the amount of erosion, which could have operated during the Oligocene episode of folding or later, nor to offer more precise model based on it. In this context, the second stage is rather artificial and is mainly used for the purpose of numeric simulations.

The main parameters for the restoration are deduced from the analysis of the cross-sections. They comprise the pre-folded length of sections, ranging from 40–50 to 70–80 km, as well as the initial thickness of the sedimentary cover, ranging from 7 to 17 km, with a mean value of 13 km. Only 3 to 5 structural cells are sufficient to describe the Greater Caucasus during its first stage of evolution. Depths of the basement top for the second evolutionary stage (folding stage) of the model were calculated from value of shortening and structural (post-folded) thickness of the sedimentary columns for 42 structural cells (definitions are analogue to those of Hossak (1979)). This 3D post-folded cover model is derived from the present boundaries of structural cells. The mean value of shortening of structural cells is 35%, with only a few values as low (and even elongation 10%) and a few others up to 67% of shortening. The mean structural thickness of the sedimentary cover is 22 km with extreme values comprised between 9–15 and 45–49 km, respectively. By comparing the post-folded cover model and the present-day geometry, which constitutes the post-orogenic stage 3 (Fig. 8), we can assess the former thickness of the presently eroded sedimentary column. The mean depth of the basement top in the 42 modules is 13 km. Actually, the depth to the basement has changed through time, from 2–5 to 25–30 km, thus outlining a regular pattern of uplift and subsiding blocks. Part of the present-day 3D structure is shown in Fig. 8. Within the frame of restored structures, it is almost impossible to trace a continuous detachment plane near top of the basement that would be expected for a mode of folding associated to A-subduction. We should also note once more that the most reliable data relate to the first evolutionary stage, i.e., the pre-folded geometry of the former sedimentary basin, and to the third evolutionary stage, i.e., the present-day post-orogenic structure of the Greater Caucasus.

Estimates of the initial thickness of the presently eroded sedimentary cover of the North-West Caucasus constitute another important result. The intensity of the erosion is related to the uplift amplitude and to change in the altitude of the topography.

The distribution of the amplitudes of erosion along strike shows an increase of erosion from the west toward the central part of the Caucasus (mean values are 2.6, 3.9, 6.6, 10.0 and 15.2 km for the first 5 sections, respectively) and some decrease further to the east (12.2, 10.7, 6.7 km for sections 7, 8 and 10, respectively). In almost all the dip sections there is a smooth increase of erosion toward the axis of the folded system (Yakovlev, 2008c, 2008d). Such regularities are expected for such a natural process. From our point of view this observation confirms the accuracy of the results and the reliability of the method and of the input data. The mean value of erosion for the whole structure is about 10 km. Its maximal values reach 20–22 km. Such a quantitative approach has been compiled here for the first time for the entire sedimentary pile of a hinterland domain (Yakovlev, 2007, 2008d).

Certainly, this data set is somewhat unexpected. If we compare these amplitudes of erosion with values of neotectonic uplift obtained from geomorphologic criteria (Nikolaev, 1978), the latter is 5–10 times smaller. Possibly this discrepancy relates to different physical boundary conditions controlling these parameters, but some corrections of the estimates of the geomorphologic amplitudes of uplift must probably be made.

The reconstruction of the initial architecture of the Chiauri flysch zone in the central sector of the Greater Caucasus was also performed by using the balanced cross-section method (Fig. 1A). This reconstruction makes it possible to test a model of the structure of the transition zone between the hinterland and the foreland extending between the Greater Caucasus and the Trans-Caucasus massif, where the same values of shortening occur in the basement and in the sedimentary cover.

The detailed structural profile comprises, from south to the north, the Middle Jurassic deposits of the foreland, the regional-scale Racha-Lechkhum Fault (RLF), the Upper Jurassic to Upper Cretaceous terrigenous-carbonate flysch deposits of the Chiauri zone, and the Middle Jurassic deposits which are located in the hangingwall of a large thrust farther north in another tectonic zone. This section is representative of the entire structure of the Chiauri zone without gaps. Twenty-six tectonic domains were defined in the cross-section (including two domains in the adjacent zones), in which all requested parameters were measured. After unfolding the present-day section which is only 29 km wide, we could measure the initial length of the section, which amounts to 65 km. The values of shortening measured from south to the north for four structural cells are 44, 58, 60 and 59%, respectively (Yakovlev, 2008e). In the pre-folded geometry of the section, the sedimentary deposits of 8-km thick were observed. According to the literature on Upper Jurassic sediments there should be 5–6 km of Middle and Lower Jurassic slates, whereas the upper part of the cross-section should also comprise up to 1–2 km of Paleocene and Eocene flysch sequence, the total thickness of the sedimentary cover being thus comprised between 12 and 15 km.

We should note that in the Trans-Caucasus massif to the south 4–6 km of Lower-Middle Jurassic deposits were initially accumulated, before being gently folded (before the Upper Jurassic). These Middle Jurassic deposits were overlain by no more than 1 km of almost undeformed Upper Jurassic and Cretaceous platform carbonates. The Racha-Lechkhum Fault must have been a normal fault, the top of the metamorphic basement in its northern hangingwall, in order to allow the accumulation of such sedimentary thickness in the Chiauri flysch zone. We first assumed that the shortening occurred with the same value in the basement and in the sedimentary cover. Then, the new post-folded sedimentary cover thickness was calculated for four structural cells based on initial thickness of 13.5 km and on values of shortening. New sedimentary column was placed in vertical direction so that certain age sediments have present-day hypsometric level in section. Depths to the base of the sedimentary column in the four structural cells are 19, 22, 24, and 17 km, respectively (Fig. 9). We can make three main conclusions from this structural scheme: (1) the deep structure of the Racha-Lechkhum Fault must accommodate 10–15 km of normal offset for the basement top; (2) the same value of shortening is geometrically possible in the Greater Caucasus for both the basement and the sedimentary cover; (3) in this interpretation the normal throw of the Racha-Lechkhum Fault increased during the Greater Caucasus structure development with thrusting in the upper level and normal motion in the lower level. This scenario accounts for a structural inversion of the former basin (Gillcrist et al., 1987), but it is sounds quite unusual because normal faulting and shortening would have been simultaneous.

These results have to be commented. Thrusting of the hinterland over the foredeeps of the Greater Caucasus is an undeniable fact. However, in the proposed structural sketch, the overall thickness of one tectonic block has increased. It means that the displacement amplitudes have been different for various horizons. This is the reason why normal offset of the basement top does not exclude at all overthrusting in the upper part of the structure. Therefore this overthrusting should be regarded as a local structure, which is not representative of the main tectonic style of this orogen. It is clear that local overthrusting can occur in a much easier way within the down-faulted blocks of the foreland (Yakovlev, 2008e). The main diagnostic thin-skinned tectonic features of accretionary prisms (Dotduev, 1997) do not apply to the real basement-involved structure of the Greater Caucasus (Yakovlev, 2011). However, this specific case does not deny the possibility of the implementation of thin-skinned models to other regions.

In structural geology, the existing morphological classifications of folds and folded structures use the principle that one class of fold must relate to one mechanism of folding. As a result, the origin of a given structural feature, i.e. a fold, can be directly derived from its classification division. For example, folds of 1B type according to the J.G. Ramsay classification (Ramsay and Huber (1987), p. 349) occur under buckling mechanism. However, to explain the wide range of 1C type folds we need to combine buckling and pure shear deformation (Srivastava and Shah, 2008). The respective order of these two components is unknown. Moreover, the two processes could develop in several steps. The mechanical explanation of any sequence is difficult, and there is no independent method to choose the correct model of fold formation within the frames of such a speculative, a priori and non tectonophysical approach. Only tectonophysic modelling and strain-stress analysis could make it possible to compare quantitatively natural structures with key features deduced from rock mechanics. For instance, it sounds possible to distinguish large folds formed by lateral buckling from transversal bending ones, using the paleostress fields (Gzovsky, 1963).

We do not use the morphological classifications and a priori related mechanisms within the proposed approach. Instead, we propose to address a correct description of the deformations in a hierarchic frame of kinematic models with a sufficient data set to constrain specific kinematic mechanisms. We face two methodological possibilities. The first possibility occurs when a correct kinematic model is used as an integrative system, in which we have to settle the studied structure (folds of single viscous layer, Fig. 2). In this case, we get parameters characterizing the deformation mechanisms and strain value. The second possibility occurs when we study several pilot mechanisms in a natural case study using the same key parameters describing deformation. This makes it possible to choose the most realistic mechanism for the studied structure (for example, finding the mechanism of gravitational sliding for the Vorontsov nappe, as in Section 3 of this article).

The possibility of getting a modern balanced geometry of large structures of linear folding (3D model of the North-West Caucasus) within the whole sedimentary cover is an essential step forward. First, such reliable structure makes it possible to sort out the different types of geodynamic model. Second, they serve as an initial material for elaboration of new geodynamic scenarii.

Linear folding in thin-layered deposits of slates and flysch of the hinterland of the Greater Caucasus is a complex hierarchic structure. Mechanisms, which form different rank structures, cannot be restricted to the simplistic contradiction between lateral buckling or transversal bending processes.

Within the frame of the proposed methods, the boundaries of different rank modules to which we apply the kinematic models, are chosen in such way that specific fold forming mechanisms can be applied to these modules. It means that the geometric restrictions of folded objects and of kinematic mechanisms of its formation should be the same. The selection of module boundaries makes the basis of hierarchic systems of objects of linear folding.

The complex kinematic models used to describe the formation process linking the geometry of individual structures and the amplitude of acting model components (elementary mechanisms) are the basis of the methods that we propose for the study of different scales of structural features. This makes it possible to give a quantitative description of strain for natural objects, including the value of shortening for various structures at different scales. This set of methods that we call multi-scale strain analysis of linear folding comprises its main hierarchic levels.

The tectonophysical approach used here makes it possible to find differences in the deformation pattern of theoretical structures with different origins. The geometry of natural objects and kinematic evolutionary models allow one also to detect exactly which formation mechanism operated in each case. The present-day and reconstructed geometries of the whole sedimentary cover enable us also to test various geodynamic models.