RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 12, ES4001, doi:10.2205/2012ES000516, 2012

*P. Haderka ^{1}, A. N. Galybin^{2}*

^{1}Escad Design GmbH, Donauwörth, Germany

^{2}Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia

This paper presents an investigation of the applicability of the stress trajectories concept and the stress trajectories – slip lines alternations method to geomechanical problems. We extend our approach introduced for the stress analysis of two-dimensional plastic bodies to the problem of the stress reconstruction in plastic regions of the lithosphere. The method is developed for the Cauchy boundary value problem and utilizes the data on principal directions as one of the boundary conditions. For this purpose the first order stress indicators of the World stress map (WSM) project database (release 2008) are utilized in computations. The set of considered boundary conditions is supplemented by the normal derivatives of the stress orientations. Complete formulation of the problem involves a yield condition. Although the general approach is not limited to a specific yield criterion, present calculations are performed for the Mohr-Coulomb criterion. Applications of the method include the stress reconstructions in three regions of the Earth's crust (Swiss Alps, Tibetan plateau and Eastern Anatolia). The continuous boundary conditions are derived by an averaging method applied to the discrete data in immediate vicinity of the starting boundary. Thereafter, for the chosen strength parameters of the Mohr-Coulomb theory (friction angle and cohesion), the unique grids of stress trajectories and slip lines are determined. These fields are further compared against the WSM data available inside the regions. The computations are made for different strength parameters in order to provide the best fit to the data. The results of the analysis are presented as two plane fields: the map of normalized mean stresses and the grid of corresponding trajectories of principal directions. The normalization parameter is unknown (it represents an initial value of the mean stress in a single node of the boundary), which is a consequence of non-uniqueness of the stress reconstruction problem based on the data on stress orientations alone. The reconstructed stress orientations are compared with the observations from the WSM database.

Identification of stress fields inside tectonic plates is an essential part of any global relative-plate-motion model [

Methods widely accepted for stress modeling are based on finite element analyses in view of gravitational potential energy (GPE) predictions [

However, model predictions resulting from various estimates of driving forces lead to inconsistent and contradictory results [

In contrast to FEM based GPE prediction methods, we suggest using the WSM stress indicators as boundary conditions in the direct formulations of boundary value problems (BVP) and hence overcoming the necessity for estimation of boundary conditions crucial in inverse problems. We propose to employ a plastic model with friction to describe the state of failure equilibrium of intraplate continental crust. In this case, the governing equations consist of the differential equations of equilibrium (DEE) and a criterion that imposes certain algebraic relationships on the stress components. The assumption of plastic rheology has been accepted in [

Fundamentals of the non-elastic stress modeling are developed in works of

In contrast to the SL approach we propose to build a solution to the problem by using the trajectories of principal stress. The stress trajectories, according to

We further use both approaches in turn and introduce the ST-SL alternations method (for detail see [

Applications of the developed approaches are presented for the problems of stress field identification in plastic regions of the lithosphere. Following

The DEE for a plane problem can be written in terms of the mean stress $P$ and the maximum shear stress $\tau$ in the following form:

\begin{eqnarray*} \partial_x P + \cos(2\theta)\partial_x \tau + \sin (2\theta)\partial_y \tau - 2 \tau \sin (2\theta)\partial_x \theta + 2 \tau \cos (2\theta)\partial_y \theta = 0, \\ \end{eqnarray*} \begin{equation} \tag*{(1)}\label{1} \partial_x P + \sin(2\theta)\partial_x \tau - \cos (2\theta)\partial_y \tau + 2 \tau \cos (2\theta)\partial_x \theta + 2 \tau \sin (2\theta)\partial_y \theta = 0. \end{equation}Here, symbols $\partial_x$, $\partial_y$ stand for differential operators. The stress functions $P, \tau$ are functions of the Cartesian coordinates $x,\: y$. They represent the sum and the difference of the principal stresses $\sigma_1$, $\sigma_2$, respectively:

\begin{equation} \tag*{(2)}\label{2} P= \frac{\sigma_1 + \sigma_2}{2}, \;\; \tau = \frac{\sigma_1 - \sigma_2}{2}, \;\; (\tau \geq 0, \; \sigma_2 \leq \sigma_1) \end{equation}Angle $\theta$ denote the orientation of the major principal stress $\sigma_1$ measured counter-clockwise with respect to the positive direction of the $x$-axis (it is accepted that compression is negative).

In problems of plane plasticity a yield criterion has to be specified additionally to equation (1). In the works by

where the maximum shear stress is expressed as a linear function of the mean stress via the angle of internal friction, $\mu$, and cohesion, $C$. It should be noted that yielding criterion (3) is valid if the maximum shear stress in (3) does not vanish and the normal stresses on the slip lines are compressive (negative).

Integration of (1) in view of (3) and the boundary conditions leads to four characteristic equations; two of them for identification of the slip lines:

\begin{equation} \tag*{(4)}\label{4} dy = dx \tan (\theta \pm \varepsilon) \end{equation}and two for the relationships valid along the slip lines:

\begin{equation} \tag*{(5)}\label{5} dP^* \pm 2 P^* \tan\mu d\theta =0 \end{equation}Equations (5) define the relationships between the mean stress and the principal directions valid along the slip lines ($P^*$ represents the mean equivalent stress, $P^* = P - C \cot \mu$). Equations (4) define two families of isogonal lines inclined to the direction of the major principal stress $\sigma_1$ at the angles:

\begin{equation} \tag*{(6)}\label{6} \pm \varepsilon = \frac{\pi}{4} + \frac{\mu}{2} \end{equation}From the latter relationship it is evident that this inclination depends on the angle of internal friction.

Figure 1 |

For completeness it has to be mentioned that if the angle $\theta$ in (4) is:

\begin{equation} \tag*{(7)}\label{7} \theta \pm \varepsilon = \frac{m}{2} \pi, \;\; m= 0,\:1,\:2, \ldots \end{equation}then the SL approach fails as the nodal points can not be identified. Since (7) is valid everywhere along the slip lines, none of these lines can be used as a stand-alone starting boundary.

The present approach is based on utilizing the stress trajectories instead of the slip lines in solving the BVP introduced in the previous section. The relationship between the orthogonal families of stress trajectories (specified by the tangents $t_s$, $t_q$) and the slip lines (specified by the tangents $t_{\alpha}$, $t_{\beta}$), are given as:

\begin{equation} \tag*{(8)}\label{8} \hat{\alpha} = \theta-\varepsilon, \;\; \hat{\beta} = \theta+\varepsilon \end{equation} where $\hat{\alpha}, \hat{\beta}$ represent the tangent angles between the slip lines $\alpha$ and $\beta$ and the $x$-axis, respectively (see also Figure 1). The system in (1) can be rewritten along the stress trajectories, hence representing the Lamé-Maxwell form of the DEE (e.g. in [

The arc lengths $s, \:q$ are measured along the stress trajectories of the first and second family, respectively. The coefficients and are functions of spatial coordinates; for the Mohr-Coulomb yield criterion they have the form:

\begin{equation} \tag*{(10)}\label{10} a = \frac{1}{2} \frac{1- \sin \mu}{C \cos \mu - P \sin \mu}, \;\; b = \frac{1}{2} \frac{1+ \sin \mu}{C \cos \mu - P \sin \mu} \end{equation}Equations (9) in view of (10) form a closed system and, together with the boundary conditions, constitute the BVP for plastic media.

A method for solution of (9) along the stress trajectories for a general, explicit yield criterion has been developed by

Similarly to SL, the solution is performed numerically by the finite difference approach and it is decoupled into two steps for the determination of the stress trajectories and the relationships along them. Based on the definition of the stress trajectories, the nodal points are sought as the points of intersections of two different, orthogonal families which emanate from neighboring nodes on the boundary. The curvilinear families of trajectories are approximated by polygons, hence the used piecewise linear approximation introduces some errors in the identification of the intersection points of trajectories. However, these errors are not significant and they are of the same order as in the SL approach where the same approximation is accepted (see (4)). The mean stresses and the orientations are calculated at the new nodes by continuation of their boundary values by using expansions into Taylor's series along the stress trajectories (derivatives of the order greater than one are neglected). The introduced directional derivatives along the trajectories ($\partial_s$ and $\partial_q$) are then transformed into the tangential and normal derivatives with respect to the contour ($\partial_t$ and $\partial_n$), respectively). Once the derivatives of the orientations are known ($\partial_t \theta$ can be retrieved from the information along the starting boundary and $\partial_n \theta$ is specified in boundary conditions) the derivatives of the stresses are determined from (9). Hence, the determined derivational terms satisfy the DEE and are consequently used to recalculate the sought mean stresses and the orientations at the new nodal points. However, it has to be noted that the mean stress $P$ can be determined only in a normalized form, thus its initial value remains unknown. With a sufficiently dense stress trajectories network the set of values ($P_{j,k}$ and $\theta_{j,k}$) at every point $z_{j,k}$ presents an approximate solution to the BVP. An overview of the equations obtained by the described numerical integration of the BVP and the inherent limitations are presented in Appendix A.

In order to illustrate that both approaches produce close results we conduct numerical tests for the Mohr-Coulomb criterion. Boundary conditions are given in terms of orientations and their normal derivatives. It is assumed that they both have linear distributions along the starting boundary, i.e. $\theta = ax +b$ and $\partial_n \theta = cx +d$ on a part of the $x$-axis. Following the classical approaches it is also assumed that all the variables having the dimension of stresses are normalized with respect to cohesion, $C$.

Figure 2 |

Table 1 |

The results obtained by the two approaches are in good agreement, however, the deviation of the results may increase (depending on the boundary conditions) when the nodes in the ST approach cross the boundaries of the domain of dependence (specified by the slip lines). These observations suggest that, while the results of the ST approach are unique, the identified stress field might comprise regions (lying outside of the domain of dependence of the boundary segment) where the real stress field is no longer dependent on the boundary conditions along the starting boundary. Therefore, the solution in these regions, obtained by the ST approach, can be considered as an extrapolation of the real solution (see [

The method of alternations is developed for the starting boundaries along which relationship (7) are valid in the case of the slip lines or (A.8) for the case of the stress trajectories. It is shown that the mean stresses and the stress orientations identified along these lines can be utilized for further extension of the stress field (beyond the boundaries of the reconstructed triangles). For this purpose it is assumed that the whole domain is in state of limiting equilibrium and that the boundary conditions allow for such extension. In the following examples boundary conditions along the starting boundary are specified according to Figure 2.

Figure 3 |

Maps of the mean stress in the extended region and inside the characteristic triangle built from the starting boundary are presented in Figure 3 (left) and (right), respectively. The contours of triangles obtained by the different methods are emphasized with full thick lines. The results obtained for different starting boundaries are tested for accuracy as described in the following sub- section.

In the case when the starting boundary is an existing slip line, both numerical procedures, developed in previous sections for the ST and SL approaches, fail to give reliable results. However, the ST approach can be modified following the idea mentioned in

Figure 4 |

Slip lines between the layers $k$ and $k-1$, are given by the nodes $z_{j, k-1}$, $z_{j+1, k-1}$ and $z_{j, k}$, $z_{j+1, k}$ available from the previous calculations (they belong to the characteristic triangle). The nodes $z_{j, k+1}$, $z_{j+1, k+1}$ and the slip line $\alpha_{k+1}$ which passes through them are unknown. Because of the expression in (8) the position of the stress trajectories $s_{j,k}$, $q_{j,k}$ is known at these nodes. Therefore, the nodal points $z_{j, k+1}$, $z_{j+1, k+1}$ can be found as the intersections of $s_{j,k},\: \beta_{j+1}$ and $q_{j+1,k},\: \beta_{j}$, respectively. After identification of $z_{j, k+1}$, $z_{j+1, k+1}$, the mean stresses and the orientations are sought at these points by employing a finite difference form of DEE (9). Details of the numerical analysis together with the corresponding relationships are presented in Appendix B.

Figure 5 |

Table 2 |

The results are in good agreement in the areas where the identified stress trajectories and the slip lines grids overlap.

In this section the developed approaches are applied to the problem of stress reconstruction in some areas of the Earth's crust. Three particular regions are considered and investigated in view of the WSM data. As there is no information available regarding the stress magnitudes, the boundary conditions are posed in terms of $\theta$ (WSM database release 2008 is considered, [

Computational areas chosen for the analysis are the Alps region, a region in Tibet (both parts of the Eurasian plate) and a region in Eastern Turkey (a part of the Anatolian and Arabian plate). Lithospheric body forces, from which the stresses arise, are a consequence of high topography of these regions, as each of these areas belongs to a mountain massif. Therefore, the methods in plasticity can be considered as more adequate. Moreover, following the research of

Application of the method of alternation is performed for straight starting boundaries. This simplification is not vital but rather used because of the nature of the data (which are scattered in the whole region). An averaging method is introduced to transform the discrete data into continuous boundary conditions. For this purpose the moving window method is used. Thus, an elongated rectangle is introduced at the vicinity of the starting line and all WSM data falling within this rectangle are selected. Consequently, the mean of the orientations inside the area is associated with the boundary condition $\theta$ and the standard deviation is associated with the normal derivative of the principal orientations, $\partial_n \theta$.

The numerical calculation begins with the identification of the stresses along the starting boundary from (A.3) and (A.4). Here, the strength parameter $C$ is chosen to be equal to unity during the whole analysis whereas $\mu$ varies between $5\mbox{°}$ and $40\mbox{°}$ in order to achieve the best correlation of the results. The free parameter in (A.3); namely the mean stress $P_{j, k}$, can be chosen arbitrarily (in the calculations the value is set to unity).

Solution proceeds in accordance with relationships (A.2)–(A.5) till the stress trajectory grid is built. Using the idea of alternations, the identified legs of the stress trajectories triangle are taken as starting boundaries for the SL approach. Next, using the slip lines reconstructed by the SL approach, the stress field is extended by the ST approach as described in sub-section 3.2 (similarly to Figure 5). The results of the calculations for different regions are presented in the following sub-sections.

Figure 6 |

Figure 7 |

The stress trajectories grid, presented by dots in Figure 7 (top, left), is built from the starting boundary. Nodes of the slip lines grid are presented by circles. It can be noted that the starting boundaries of the SL approach coincide with the stress trajectories. Finally, information along the slip lines is used to extend the solution area once more. These grids are presented by the nodes in the form of small triangles. The map of the mean stresses is shown in Figure 7 (bottom, left). The pictures on the right (Figure 7) present two inserts that show the predicted and observed stress orientation data used for the analysis of the accuracy in the stress field reconstruction.

Figure 8 |

Figure 9 |

The Eastern Anatolia region represents an example of a continental collision zone [

Figure 10 |

Figure 11 |

As evident from the figures in this example, in some regions the results correlate to a good degree of accuracy whereas in others the accuracy of the predicted orientations is fair. However, it follows directly from the nature of the observed orientations that a good correlation cannot be achieved in this area because of the randomness present in the data. One way to improve the certainty of the results is to use the data of the highest quality (WSM rank A) in this region. However, it has to be noted, that the inaccuracies can also stem from the choice of the constitutive law (rheology), adoption of the 2D model or other, including geological, effects.

Applicability of the ST-SL alternation method, developed for the limiting equilibrium analysis of two-dimensional plastic bodies, has been studied to address the geomechanical problem of stress reconstruction in the Earth's crust. Namely, identification of stress fields in mountain regions of the crust has been of the main concern. Discrete data on stress orientations taken from the WSM database have been utilized for numerical analysis, the stress magnitudes remained unknown. The data have been treated by a specially developed procedure to obtain continuous functions of the orientations of the principal stress and their normal derivatives along the starting boundary. As a result, the proposed numerical procedure allowed for the unique identification of the stress trajectories grid (or the slip lines grid). The mean stresses are found in the normalized form such that the normalization parameter remains unknown, which emphasizes the non-uniqueness of the stress reconstruction problem based on stress orientations alone.

The numerical scheme has been developed for the Mohr-Coulomb criterion which is widely accepted for modeling of stress fields in plastic regions of the lithosphere. The criterion contains two parameters controlling the stress fields; the angle of internal friction and cohesion. The latter has been used for normalization of stresses, while the former can be estimated from the best fit of the predicted stress orientations and the external data, i.e. for data that have not been involved in the identification of the boundary conditions by the moving window method. It should be noted that the Mohr-Coulomb yield condition does not limit the ST approach and the ST-SL alternations method. It has been proved in previous studies that some other criteria are also admissible.

The method has been applied to reconstruct plastic stress fields in three different regions of the Earth's crust. In particular, the best-fit model of the Swiss Alps region predicts the NNW stress orientation with a NW trending near Lake Maggiore. For the vast part of the Tibetan plateau, the maximum compressive stresses are oriented in North direction with the exception of the south-eastern part where a NNW trend has been observed. Stress orientations in the eastern Anatolia and Arabian plate show NNE to NE directions. In general, in all the regions analysed, the numerical modeling satisfactorily predicts the data on the WSM stress indicators (calculated relative errors are within the expected errors inherent in the data). However, the predictions are less accurate for the north-eastern Anatolian region which can be explained by the high level of WSM data scattering.

The analysis performed on the basis of plastic rheology has shown satisfactory agreement between the predicted and observed stress orientations in some mountain regions of the lithosphere. This justifies the existing opinion [

Figure 12 |

The fact that both relationships in (A.1) should specify the same nodal point allows for the determination of the unknown distances $\delta_{j,k}^{(s)}$, $\delta_{j+1,k}^{(q)}$ as follows:

\begin{eqnarray*} \left( \begin{array}{c} \delta_{j,k}^{(s)} \\ \\ \delta_{j+1,k}^{(q)} \end{array} \right) = \frac{1}{\cos(\theta_{j,k} - \theta_{j+1,k})} \times \end{eqnarray*} \begin{equation} \tag*{(A.2)} \left( \begin{array}{cc} \cos \theta_{j+1,k} &\sin \theta_{j+1,k} \\ \\ \sin \theta_{j,k} & - \cos\theta_{j,k} \end{array} \right) \left( \begin{array}{c} {\rm Re} (z_{j+1, k} - z_{j, k}) \\ \\ {\rm Im} (z_{j+1, k} - z_{j, k}) \end{array} \right) \end{equation}

The accuracy of the utilized piecewise linear approximation is illustrated in Figure 12 (right). The mean stresses and the principal directions at the new nodal points $z_{j, k+1}$ are calculated by means of the Taylor's expansion as:

\begin{equation} \tag*{(A.3)} \left( \begin{array}{c} P_{j,k+1} \\ \\ \theta_{j,k+1} \end{array} \right) = \left( \begin{array}{c} P_{j,k} \\ \\ \theta_{j,k} \end{array} \right) + \delta_{j,k}^{(s)} \left( \begin{array}{c} \partial_s P_{j,k} \\ \\ \partial_s \theta_{j,k} \end{array} \right) \end{equation}

Derivatives of the order higher than one are neglected. Derivational terms on the right-hand side of (A.3) are unknown and are yet to be determined. Based on the boundary conditions, posed in terms of orientations $\theta$ and their normal derivatives $\partial_n\theta$, the numerical procedure utilizes the transformation formulas between the directional derivatives along the stress trajectories $(s,q)$ and the tangential and normal directions to a contour ($t,n$). DEE (9) in view of this transformation (DEE in terms of normal and tangential derivatives) can be written in the form:

\begin{eqnarray*} \left( \begin{array}{c} \partial_t P_{j,k} \\ \\ \partial_n P_{j,k} \end{array} \right) = \left( \begin{array}{cc} a_{j,k}\cos \xi_{j,k} & a_{j,k} \sin \xi_{j,k} \\ \\ -b_{j,k} \sin \xi_{j,k} & b_{j,k} \cos\xi_{j,k} \end{array} \right)^{-1}\times \end{eqnarray*} \begin{equation} \tag*{(A.4)} \left( \begin{array}{cc} + \sin \xi_{j,k} & - \cos\xi_{j,k} \\ \\ - \cos\xi_{j,k} & - \sin \xi_{j,k} \end{array} \right) \left( \begin{array}{c} \partial_t \theta_{j,k} \\ \\ \partial_n \theta_{j,k} \end{array} \right) \end{equation}

or equivalently:

\begin{eqnarray*} \left( \begin{array}{c} \partial_n P_{j,k} \\ \\ \partial_n \theta_{j,k} \end{array} \right) = \left( \begin{array}{cc} \sin \xi_{j,k} & - \cos \xi_{j,k} \\ \\ -b_{j,k} \cos \xi_{j,k} & a_{j,k} \sin \xi_{j,k} \end{array} \right)^{-1}\times \end{eqnarray*} \begin{equation} \tag*{(A.5)} \left( \begin{array}{cc} - a_{j,k} \cos \xi_{j,k} & \sin \xi_{j,k} \\ \\ b_{j,k} \sin \xi_{j,k}& - \cos\xi_{j,k} \end{array} \right) \left( \begin{array}{c} \partial_t P_{j,k} \\ \\ \partial_t P_{j,k} \end{array} \right) \end{equation}It should be noted that the tangential derivatives in the equations above can be expressed via the values of $\theta$ and $P$ on the $k$-th layer by the forward finite differences:

\begin{equation} \tag*{(A.6)} \left( \begin{array}{c} \partial_t P_{j,k} \\ \\ \partial_t \theta_{j,k} \end{array} \right) = \frac{1}{|z_{j+1,k} - z_{j,k}|} \left( \begin{array}{c} P_{j+1,k} - P_{j,k} \\ \\ \theta_{j+1,k} - \theta_{j,k} \end{array} \right) \end{equation}

Notations $a_{j,k}, \; b_{j,k}$ refer to the values at the node $z_{j,k}$ of the corresponding coefficients introduced in (10). It can be readily observed that by substitution of the second equation in (A.6) into (A.4) the tangential and normal derivatives of $P$ can be found directly. These derivational terms (along the tangent and the normal) specify the directional derivatives in the Taylor's expansion (A.3). As the result, the unknown stress field $P$ is found in the normalized form such that the initial value of $P_{j,k}$ remains unknown.

Since (A.4) and the boundary conditions now allow for recalculation of the stresses from (A.3) along the starting boundary the solution further continues with equations (A.5) and (A.6) (where, the unknown normal derivatives are sought as functions of the known values along the boundary).

One obvious limitation of the method is the existence of the inverse matrix that becomes singular if its determinant vanishes, which yields:

\begin{equation} \tag*{(A.7)} a_{j,k} b_{j,k}=0 \Leftrightarrow\cos \mu =0 \end{equation}in (A.4) or:

\begin{equation} \tag*{(A.8)} a_{j,k} - b_{j,k} - (a_{j,k} + b_{j,k}) \cos 2 \xi_{j,k} =0 \end{equation}in (A.5). Note, that the restriction in (A.7) implies non-physical values for the friction angle that usually varies in the range $(0, \pi/4)$.

The finite difference form of equation (9) written with respect to the point $z_{j,k}$ and along trajectories $s_{j,k}$, $q_{j+1,k-1}$, yields:

\begin{eqnarray*} a_{j,k} \frac{P_{j+1,k+1} - P_{j,k}}{\Delta s_{j,k}} + \frac{\theta_{j,k} - \theta_{j+1,k-1}}{\Delta q_{j+1,k-1}} =0, \end{eqnarray*} \begin{equation} \tag*{(B.1)} b_{j,k} \frac{P_{j,k} - P_{j+1,k-1}}{\Delta q_{j+1,k-1}} + \frac{\theta_{j+1,k+1} - \theta_{j,k}}{\Delta s_{j,k}} =0. \end{equation}Here, $\Delta s_{j,k} = | z_{j+1, k+1} -z_{j,k}|$ and $\Delta q_{j+1,k-1} = | z_{j, k} -z_{j+1,k-1}|$ determine the distances between the given nodes (see Figure 4). Similarly this scheme can be applied for the trajectories $ s_{j,k-1}, \; q_{j+1, k}$ at the point $z_{j+1, k}$ which gives:

\begin{eqnarray*} a_{j+1,k} \frac{P_{j+1,k} - P_{j,k-1}}{\Delta s_{j,k - 1}} + \frac{\theta_{j,k+1} - \theta_{j+1,k}}{\Delta q_{j+1,k}} =0 \end{eqnarray*} \begin{equation} \tag*{(B.2)} b_{j+1,k} \frac{P_{j,k+1} - P_{j+1,k}}{\Delta q_{j+1,k}} + \frac{\theta_{j+1,k} - \theta_{j,k-1}}{\Delta s_{j,k-1}} =0 \end{equation}Because of the piecewise linear approximation adopted for the ST approach, some errors are introduced. The averages of the solutions in (B.1) and (B.2) at the node $z_{j+1, k+1}$ have been used to improve accuracy.

*Tectonophysics, 29,* 223–49, doi:10.1016/0040-1951(75)90148-1.

*Tectonophysics, 305,* 153–164, doi:10.1016/S0040-1951(99)00004-9.

*Geophys. J. Int., 160*, 634–650, doi:10.1111/j.1365-246X.2004.02528.x.

*J. Geophys. Res., 94*(B4), 3967–3990.

*Computers & Geosciences, 25, * 383–394, doi:10.1029/JB094iB04p03967.

*J. of Nepal Geological Soc., 36*.

*Tectonophysic, 485,* 164–177, doi:10.1016/j.tecto.2009.12.012.

* Geophys. Res. Lett., 12,* 77–80, doi:10.1029/GL012i002p00077.

*Geology, 22*, 831–834, doi:10.1130/0091-7613(1994)022<0831:TSITAP>2.3.CO;2

*Earth Planet. Sci. Lett., 133,* 299–309, doi:10.1016/0012-821X(95)00084-P.

*Proc. R. Soc. A: London, 254*, 1–45.

*Geophys. J. Int., 101*, 1990 425–478, doi:10.1111/j.1365-246X.1990.tb06579.x.

*Geophys. J. Int., 161*, 445–468, doi:10.1111/j.1365-246X.2005.02598.x.

*Bull. Fac. Sci., Univ. Ryukyus, 89*, 5–25.

*Science, 278*, 647–649, doi:10.1126/science.278.5338.647.

*J. of Geophys. Res., 110*, B12401, doi:10.1029/2004JB003541.

*J. Geophys. Res., 106*(B8), 16,435-16,460, doi:10.1029/2001JB000208.

*Geophys. J. Int., 169*, 866– 896, doi:10.1111/j.1365-246X.2007.03274.x.

*Structural Geology,* Cambridge University Press, Cambridge, doi:10.1017/CBO9780511777806.

*Photoelasticity, Vol. 1,* Wiley, New York.

*Int. J. of Sol. and Struc., 41*, 5125–5142, doi:10.1016/j.ijsolstr.2004.04.007.

*Geophys. J. Int.*, 1–26, doi:10.1111/j.1365-246X.2009.04326.x.

*Nature, 404,* 139–140, doi:10.1038/35004703.

*Tectonophysics, 266*, 11–24, doi:10.1016/S0040-1951(96)00180-1.

*J. Geophys. Res., B97-11*, 805–820.

*Int. J. of Sol. and Struc., 48,* 450–462, doi:10.1016/j.ijsolstr.2010.10.016.

*The World Stress Map database release 2008,* (available online at www.world-stressmap.org), doi:10.1594/GFZ.WSM.Rel2008.

*J. of Geophy. Res., 113*, B07404, doi:10.1029/2007JB004953.

*The Mathematical Theory of Plasticity,* Clarendon Press, Oxford.

*Rev. Geophys., 45*, RG3001, doi:10.1029/2005RG000181.

*The Geological Society of America, Special Paper, 430*, 693–722, doi:10.1130/2007.2430(32).

*Geophys. J. Int., 154,* 5–34, doi:10.1046/j.1365-246X.2003.01917.x.

*Photoelastic stress analysis,* Wiley, London and New York.

*Geophys. J. Int., 170*, 1328–1335, doi:10.1111/j.1365-246X.2007.03468.x.

*Proc. 11th Int. Conf. of IACMAG: Turin, 4,* 441–450.

*Geophys. J. Int., 155, * 126–138, doi:10.1046/j.1365-246X.2003.02023.x.

*J. of Geology, 5*, 212–216, doi:10.1130/0091-7613(1977)5<212:ROTTOE>2.0.CO;2.

*Physics of the Solid Earth, 41*, 902–915.

*J. of Geophys. Res., 97-B8,* pp. 783–803, doi:10.1029/91JB01096.

*Tectonophysics, 275*, 199–219, doi:10.1016/S0040-1951(97)00021-8.

*Tectonophysics, 484*, 27–35, doi:10.1016/j.tecto.2009.08.029.

*J. of Geodyn., 27*, 1–21, doi:10.1016/S0264-3707(97)00024-0.

*Statics of Granular Media,* Pergamon Press, Oxford.

*Ann. Rev. Earth Planet. Sci., 27*, 417–462, doi:10.1146/annurev.earth.27.1.417.

*Science, 258*, 1328–1332, doi:10.1126/science.258.5086.1328.

*Earth Planet. Sci. Lett., 186*, 75–91, doi:10.1016/S0012-821X(01)00229-1.

*Antarctica: A keystone in a changing world, Cooper A. K. et al., eds., Proceedings of the 10th International Symposium on Antarctic Earth Sciences*, The National Academies Press, Washington, DC, 115–124.

*J. Geophys. Res., 112*, B01401, doi:10.1029/2005JB004244.

*J. Geophys. Res., 112,* B11403, doi:10.1029/2006JB004807.

*Earth Planets Space, 58*, 1463–1473.

*Neotectonics of North America, Slemmons, D. B., Engdahl, E. R., Zoback, M. D., Blackwell, D. D. (eds.),* 339–366.

*International Geology Review, 44*, 383–401, doi:10.2747/0020-6814.44.5.383.

*Earthquake Seismology – Treatise on Geophysics, 6*, 253–274.

Received 15 May 2012; accepted 17 May 2012; published 5 June 2012.

**Citation:** Haderka P., A. N. Galybin (2012), Plastic stress field reconstruction based on stress orientations data, *Russ. J. Earth Sci., 12*, ES4001, doi:10.2205/2012ES000516.

Copyright 2012 by the Geophysical Center RAS.