Present-day stress field pattern in the Vrancea seismic zone (Romania) deduced from earthquake focal mechanism inversion

Vrancea seismogenic zone in the South-Eastern Carpathians is characterized by localized intermediate-depth seismicity. Due to its complex geodynamics and large strain release, Vrancea represents a key element in the Carpatho-Pannonian system. Data from a recently compiled catalogue of fault plane solutions (REFMC) are inverted to evaluate stress regime in Vrancea on depth. A single predominant downdip extensive regime is obtained in all considered clusters, including the crustal layers located above the Vrancea slab. The prevalent stress regime confirms previous investigations and requires some mantle-crust coupling. The S 3 principal stress is close to vertical, while S 1 and S 2 are horizontal, oriented perpendicularly and respectively tangentially to the Carpathians Arc bend. This configuration is present at any depth level. According to seismicity patterns, there are two main active segments in the Vrancea intermediate-depth domain, at 55 – 105 km and 105 – 180 km, both able to generate major events. The configuration of the tectonic stresses as resulted from inversion is similar in both segments. Also, high fault instability (I > 0.95) is characterizing the segments. The only notable difference is given by the friction and stress ratio parameters which drop down in the bottom segment from μ = 0.95 to μ = 0.55 and from R = 0.51 to R = 0.29. This variation is attributed to possible weakening processes activated below 100 km depth and can explain the intensification of seismicity production as earthquake rate and average energy release in the lower segment versus the upper segment.


Introduction
There are a few areas in the world characterized by localized intermediate-depth seismic activity where earthquakes are generated in subducted lithosphere far from the plate boundaries [Prieto et al., 2012]. Among these we mention Bucaramanga (Colombia), Hindu Kush (Afghanistan) and Vrancea (Romania). In Vrancea, earthquakes The seismicity and stress data are key elements to understand the present-day dynamics of the tectonic system in the Vrancea area.
The anomalous high-velocity body beneath Vrancea is clearly visible in the tomographic images both from teleseismic and local waveforms [Wortel and Spakman, 2000;Martin et al., 2006;Koulakov et al., 2010]. Its highvelocity volume is significantly larger than the seismically active range [Martin et al., 2006]. There is no consensus about the origin of the body -whether it is a piece of the relic continental lithosphere or a relic of oceanic slab and hence about the type of geodynamic processes taking place in the region descending oceanic slab, slab detachment or delamination of the continental lithosphere [for a review of various models of geodynamics and seismicity in Andrei Bala et al. and stress field along the vertical seismically active body beneath the Vrancea region, in correlation with the possible seismicity segmentation. A few investigations were performed on this subject [Oncescu and Trifu, 1987;Oncescu and Bonjer, 1997], but using a reduced amount of data.
The purpose of the present study is to employ a comprehensive set of fault-plane solutions available for the Vrancea intermediate-depth and crustal earthquakes  in order to emphasize the stress characteristics in the subducted lithosphere and overlying crust in connection with the geodynamic processes responsible for generating the seismicity at the South-Eastern Carpathians Arc bend. The dataset comprising of a catalogue of fault-plane for the earthquakes recorded between 1929 and 2012 [Radulian et al., 2018;2019] called REFMC  is much better represented both in number of events and quality (considering the number of stations used but also in some cases the method of determination), comparatively with previous investigations [Radulian et al., 2000;Radulian et al., 2002]. The improved data allow us to estimate the stress within the earthquake-prone descending body and in the overlying crust by inverting focal mechanism data in order to find the horizontal stress and to explore its distribution in depth. in Vrancea zone (VNI) at depths ≥ 50 km and the crustal earthquakes at depths < 50 km located above the subcrustal source, mostly in the southern part of Romania ( Figure 1). An analysis focused on the fault-plane solutions of the events and inversion of focal mechanism for the determination of the horizontal stress in the crust on Romanian territory was carried out by Bala et al. [2020]; this however does not provide many details regarding the special setting of the Vrancea zone.  Bala et al. [2020] together with all other crustal events in Romania. The definition and the abbreviated names of the zones were taken from Bala et al. [2019]. The seismicity in the largest part of the eastern Moesian platform (MO1, MO2) is weak, while to the west of Argeş river, it is almost missing; the only exception is the zone located close to the epicentral area of the Vrancea intermediate-depth earthquakes (MO3 and MO4) which is likely to be in connection with the geodynamics in the mantle. The earthquakes generated here had M w <5.5 and many times they occurred in sequences, especially towards SE of Vrancea, in MO3 [Popescu and Radulian, 2001]. Note also that these sequences follow an alignment NE-SW, which resembles the similar alignment observed for the deep events [Popescu and Radulian, 2001;Popescu et al., 2011;Raileanu et al., 2009]. The seismic activity in NDO is clustered along the major faults crossing the region from Black Sea to the Carpathians (Peceneaga-Camena and Sf. Gheorghe faults). The seismicity is more diffuse in BD and FC. Five shocks of M w between 6.2 and 6.5 were recorded in FC since XVI th century.

Fault plane solutions data and selected groups
In order to invert the fault plane solutions for stress determination we consider two groups of events, first in the slab sinking in the mantle and second in the crust surrounding the Vrancea zone, groups which are investigated separately.
We focus our attention in this work on the seismic activity closely related to the Vrancea zone (study area is simply defined in Figures 2 and 3 as the polygon ABCD). For this we use focal mechanisms of 634 intermediate-depth events with M w between 2.6 and 7.7 and 74 crustal events with M w between 2.1 and 4.6 which occurred in the 1929 -2012 period. The 3-D representation of the hypocenter distribution for the events with reliable focal mechanism (at least 20 well-defined first P-wave polarities) shows interesting and characteristic features ( Figure 2): a clear discrepancy as concerns the energy release and clustering degree between the subducted lithosphere and crustal earthquakes, a relative deficit of earthquakes in the transition layer between crust and subducted lithosphere, a predominant reverse faulting regime in the slab (the same regime can be noticed in the crust as well, but less pronounced). A change is visible in the seismicity regime between the upper segment of the Vrancea seismogenic volume (noted below as VNI_A) and the lower segment (noted below as VNI_B).
Ternary diagrams are plotting the plunge of the intermediate (B) and tension axes (T) of focal mechanisms to illustrate the faulting style of the seismotectonic regimes. The plotting method, as proposed by Álvarez-Gómez [2019] using the FMC script which improves the ternary diagrams of Kaverina et al. [1996], was already used to characterize the faulting style in the crust in the southern Romania by Bala et al. [2019]. Note that the ternary diagrams do not reflect fault orientation, which is an additional criterion for the delineation of the seismotectonic domains [Pondrelli and Morelli, 2008].
In order to apply the inversion of focal mechanism data to evaluate the stress regime, we subdivided the seismogenic volume into different cell zones sufficiently large to include minimum 30 events with computed fault plane solutions, as suggested by Vavryčuk [2014]. The partition for VNI is made following mainly the particular geometry of the seismogenic volume. As such, we defined parallelepiped-like cells of 35x25x25 km dimensions, in NE, SE, and vertical directions, respectively. In this way, we fit the shape of the projection the Vrancea earthquakeprone volume, which is an approximate rectangle of 70-75 km length per 25-30 km width (Figure 3), by 3 cells. Also, the entire depth range (55 -180 km) is fitted within 5 layers, being able to capture the hypocenter transition with depth as reflected by Figure 2. The cell dimensions are also conventionally selected such as to have regular cells fitting the Vrancea earthquake-prone volume and to include in each of them at least 30 earthquakes.
The focal mechanisms for the earthquakes occurring in the crustal seismic zones located directly above VNI (MO4 and MO3) are also examined in order to compare the horizontal stress computed in the subducted lithosphere with that in the overlying crust. For the crustal earthquakes, we used single cells per seismogenic area, as shown in Figure 3. MO4 is practically overlapping the epicentral area of the Vrancea subcrustal earthquakes, while for MO3 a non-regular shape is adopted so that to cover as much as possible the seismicity pattern. With the partition specified in this way, the repartition of the total number of earthquakes with focal mechanism computed (634 events) in VNI is that of Table 1.

5
Stress field pattern in the Vrancea zone Figure 3. Events in the Vrancea intermediate-depth zone (VNI: 55 -180 km depth) projected at surface, in order to reflect the cell clusters in latitude and longitude in which they were separated: VNI NE, VNI middle and VNI SW; these cell clusters can be seen also in 3D in Figure 2 (with A, B, C and D points as reference) and in Figure 6; crustal zones MO3 and MO4 (1 -50 km depth) and corresponding earthquakes can also be seen in purple.

Methodology
Several methods have been developed to determine stress field by inversion of earthquake focal mechanisms. Michael [1984], Gephart [1990], Angelier [2002] contributed to the development of these techniques. These methods usually assume that, for optimum results, a few main conditions should be fulfilled: • tectonic stress is homogeneous in the region; • earthquakes occur on pre-existing faults with varying orientations; • the slip on a fault occurs in the direction of maximum shear stress, the so-called Wallace-Bott hypothesis [see Wallace, 1951;Bott, 1959]; • the earthquakes do not interact with each other and do not disturb the background tectonic stress. Vavryčuk [2014] developed these technique into an algorithm called StressInverse, following the conditions above.
Obviously, these conditions are not entirely satisfied in our case, taking into account that we employ for inversion data from earthquakes with foci spread from near surface down to 180 km depth.
When stress is not uniform over the entire study region, we subdivide it into smaller areas for which the condition of uniform tectonic stress is more easily satisfied. The way is to define subzones as small as possible but containing at the same time a reasonable number of earthquakes (no less than 30 events) (see Table 1 and Figure 3).
The Wallace-Bott assumption of the slip vector parallel to the stress on the fault is valid only in isotropic media.
In anisotropic media, the two vectors are not necessarily parallel and the problem becomes more complicated, particularly when the knowledge of the anisotropy in the focal zone is difficult, which happen to be our case, mostly due to the depth of the hypocenters.
Stress changes due to the occurrence of small and moderate earthquakes are usually negligible, but large earthquakes can significantly affect the background stress field. From this point of view, it is ideally to perform the inversion for earthquakes clustered not only in space but also in time, separately for time intervals between large earthquakes, and after them. For the time being, the dataset available cannot ensure sufficient statistics to obtain a good accuracy of the results on distinct time intervals.
The software StressInverse is applied to estimate the stress tensor. This software is based on Michael's method [Michael, 1984] and incorporates the instability criterion proposed by Lund and Slunga [1999]. It represents a linear, iterative stress inversion based on the four assumptions presented above. Although these assumptions look apparently very restrictive, the analysis of real observations proves that they are well-satisfied in most cases, in particular, for local seismicity consisting of weak to moderate earthquakes (Gorgun et al., 2016;Fojtíková and Vavryčuk, 2018). Details about the method and its accuracy are published in Vavryčuk [2014].
The method proposed by Michael [1987] was developed in several stress inversion codes such as SATSI [Hardebeck and Michael, 2006 The inversion outcome comprises the orientation of the principal stress axes (S 1 , S 2 , S 3 ) and the stress ratio (R): where σ 1 , σ 2 , σ 3 (σ 1 > σ 2 > σ 3 ) are the eigenvalues of the stress tensor and 0 < R < 1.
Both StressInverse and MSATSI algorithms use the same definition of stress ratio R as the one in Equation 1 [see also Vavrickuk, 2014;Martinez Garzon et al., 2014;2015].
The estimation of errors is provided by a repeated stress inversion of focal mechanisms contaminated by artificial noise. In StressInverse the uncertainties are computed as the maximum differences between the results of the inversion for noise-free and noisy data after 1000 noise realizations. This method proves to be more convenient for two reasons: first, it takes into account the case when some nodal planes are more uncertain than the others by specifying differently noise levels for the fault orientations and slip directions. Second, the uncertainties are obtained in a narrow domain. For example, the shape ratio R computed by iterative method has uncertainties between 6 -12% as opposed to those obtained by Michael's method of 20 -44% [Vavryčuk, 2014]. For this study we Andrei Bala et al.
used 100 realizations of random noise in the inversion. The level of noise of 10°-12° corresponds to the estimated accuracy of input focal mechanisms. The inversion process is usually stopped after 5 iterations.
MSATSI is useful as the results can be considered and presented as 1D, 2D or 3D models in some arbitrary cells.
A scheme of damping can be employed between two successive cells in order to smooth the transition of the stress direction from one cell to the next. However, when analyzing real data sets, the uncertainties of the retrieved stress directions, the shape ratio, the fault selection and overall friction on faults are calculated by the bootstrap method as proposed by Michael [1987].
In a recent study about the results obtained with MSATSI, Martinez-Garzón et al. [2016] consider that, if the focal mechanisms sufficiently cover the range of allowed orientations within an ambient stress field with a stress ratio R = 0.5, at least 20 mechanisms are necessary to resolve the stress field orientation with low noise conditions, while at least 40 are necessary for data sets with high noise. In our case we have employed cells defined so that to include at least 30 events with computed mechanism (in eight out of ten cells we have more than 40 events -see Table 1).  Figure 3); when a cell does not contain at least 30 events, we associate it to one or two neighboring cells in the same layer (like C1 in VNI1).

Results
In order to apply inversion procedure, we set the partition of VNI following three rules: using regular cells, The StressInverse algorithm does not require that the events be located in a certain cell. In order to compare the results obtained by inversion using StressInverse and MSATSI algorithms we carried out the computations following the adopted cell partition for both algorithms. The inversion results can be represented in 2D images on different depth layers, or in pseudo-3D images.
Preliminary computation carried out separately on the two crustal seismic areas, MO3 and MO4, shows a similar thrust faulting regime (S 3 vertical) and roughly similar orientation of the principal horizontal stresses (S 1 and S 2 ).
However, taking into account the number of events and the relatively larger errors than for the subcrustal earthquakes, we considered a single area for the two crustal seismic areas, MO3 and MO4 together.
The results of stress inversion after application of StressInverse algorithm are presented in the • in all cases (crust and slab layers) the S 3 principal axis is close to vertical: axis plunge ranges from 71° in the crust and 87° in the bottom layers. This is compatible with a compressive stress regime. Correspondingly, the axes S 1 and S 2 are close to horizontal (plunge angles below 9° in the slab and below 20° in the crust); • stress ratio (R) has values between 0.30 ± 0.1 to 0.84 ± 0.09 (VNI4 and VNI1). The confidence limits are spread because the shape ratio is sensitive to the number of focal mechanisms inverted and to their accuracy; • there are two conjugate focal mechanisms (denoted as "principal focal mechanisms" by Vavrycuk) compatible with the stress configuration. For the subsequent discussions and graphic representations, we chose one of the conjugate focal mechanisms, the one with the parameters (strike/dip/rake) given in Table 2  to the principal mechanism after the values of strike/dip/rake in the Table 2.
Stress field pattern in the Vrancea zone The faulting type represented by a certain DC moment tensor can be determined by the orientation of its eigenvectors, P, B and T, in a sense which may apply to the characterization of the stress regime in an area. Zoback (1992) determined a set of rules for the classification of faulting type depending mainly on the plunge angles of the P and T axes, complemented, in some cases, by the respective plunge of the B axis (Table 3). These also determine the trend of the maximum horizontal stress axis, S Hmax .
There is a region of P and T plunges in which focal mechanisms are not classified in a typical faulting type, but rather as "odd" or "unknown" (U). This is the case of either sub-horizontal faults with horizontal slip or sub-vertical faults with (nearly vertical) dip-slip (both with plunge ≈ 45° or all three axes with 25° < pl < 45°). Such cases are generally rare and may characterize very low-angle normal faults or thrusts where the principal stress field is tilted out of horizontal and vertical planes [Zoback, 1992;Heidbach et al., 2016a]. Our data presented in Table 2 support the conclusion that the results are all included in the line 6 of Table 3 and the tectonic regime is assigned to the category of thrust faulting. The conclusion is valid either for every layer considered in the subcrustal domain from VNI1 to VNI5, or for the principal subcrustal domains VNI_A (55 -100 km depth) and VNI_B (101 -180 km depth).
Somewhat unexpectedly, the crustal part MO3-MO4, which is located above the subcrustal seismically active body, shows the same tectonic regime of thrust faulting.
The seismicity pattern on depth (Figure 2), with a significant deficit at the transition from the crust to the subducted lithosphere, suggests a possible decoupling between the crustal and slab segments. However, the stress regime in the crust is of thrusting type even though it is less constrained than in the subducted lithosphere. For example, the range of the S 3 axis plunge is larger (71° -86°) in the overlying crustal domain ( Stress field pattern in the Vrancea zone

Fault instability
The fault plane instability, I, is an important parameter of evaluating the fault plane proposed by Vavryčuk et al. [2013] and Vavryčuk [2014]. Since differently oriented faults have a different susceptibility to be activated in a given stress field and thus one can define faults more or less unstable for this stress field, Vavryčuk et al. [2013] introduced a quantity to measure this instability defined by the following Equation: ( 2) where τ c and σ c are the shear traction and effective normal traction along the principal fault, and τ and σ are the shear traction and effective normal traction along the analyzed fault; μ is the friction coefficient. Fault instability I ranges from 0 (the most stable faults) to 1 (the most unstable faults). The most unstable fault is the optimally oriented fault for shear faulting called the principal fault.
It can be shown that the fault instability is independent of absolute stress values, and (2) can be expressed as a function of friction μ, shape ratio R (Equation 1) and directional cosines n defining the inclination of the fault plane from the principal stress axes: The directional cosines n 1 , n 2 , n 3 define the inclination of the fault plane from the principal stress axes, n being expressed in the coordinate system of the principal stress directions [Vavryčuk, 2014].
To evaluate the fault instability using Equation (3), a value of friction μ is needed. Friction on faults most often ranges between 0.2 and 0.8, but its value is usually unknown, especially at depth more than 50 km. Numerical tests revealed, however, that the inversion is rather insensitive to μ [Vavryčuk, 2014], so it is often sufficient to assign some mean value to friction during the inversion, for example, μ = 0.6.
Another approach is to run the inversion for several values of friction and adopt the value which produces the highest overall instability of faults for the data inverted. This approach is used in the StressInverse algorithm in order to run synthetic tests as well as in the applications to real data [Vavryčuk, 2014].
Once StressInverse yields a friction coefficient μ, which results in the highest overall instability, we compute fault instability I using (3), (4) and (5). In order to compute the directional cosines (n 1 , n 2 , n 3 ), a useful algorithm is MohrPlotter, which is conveniently transforming all the parameters present in the Mohr circle (MohrPlotter v. 3.0 by Allmendinger, 2020).
An additional parameter that can be derived from the Mohr circle is the angle α:  Table 4. The first two parameters (R and μ) are results of the StressInverse program after inversion; α and I are obtained according to Equations (6) and (3). To test how our results fit the regional stress map, we transformed the orientations of S 1 and S 3 to azimuths of maximum horizontal stress (σ H ) and minimum horizontal stress (σ h ), using the script of Lund and Townend [2007], which is also implemented in the MSATSI code of Martínez-Garzón et al. [2014]. The values obtained are given in Table 4 for crustal segment and the two principal segments in the slab and they are represented in Figure 10. Laboratory measurements indicate for the sedimentary rocks in the upper part of the crust μ values between 0.6 and 1.0 [Byerlee, 1978;Townend and Zoback, 2000]. For deeper regions there are no direct measurements, but some authors consider that static friction might have values close to 1 at these specific depths, deeper than 50 km [Townend and Zoback, 2000]. We should keep in mind also that the values of static friction are significantly higher than the dynamic friction over the earthquake rupture area (friction at the seismic sliding rates of m/s during the earthquake), μ < 0.1, due to various proposed weakening mechanisms [see Marone, 1998].
As shown in Figure 6, the principal stress axes for most domains are nearly horizontal and vertical. The method to distinguish the fault plane from the auxiliary plane is looking for the highest instability following the Mohr-Coulomb criterion.  In all 10 cells from the vertical seismic zones VNI1 -VNI5 we have the same faulting type of thrust faulting determined in the Table 2 and plotted in Figure 4 and Figure 5, with S 3 being very close to vertical and S 1 and S 2 which are having similar orientations, especially at the same level, for example at VNI3 and VNI4. From previous studies  we were aware that most of the earthquakes with mechanism in Vrancea slab are showing the same prevalent faulting type of thrust faulting mechanism (R -yellow or R-SS -orange in Figure 2). In the present study we have obtained thrust faulting and the corresponding tectonic mechanism of radial compression in all the 10 cells that we chose for inverting the event mechanisms in order to obtain the directions of the horizontal stress.

Seismic zone Stress ratio (R) Friction (μ) α = 45° -(tan -1 (μ)/2) Fault instability I S H /S h
The conclusion is that at all subcrustal levels (VNI1 -VNI5) are in the presence of thrust faulting and compressive stress regime in which S H > S h > S v . S Hmax and S hmin are showing similar directions at the same level from VNI1 to VNI5 with only minor changes in direction from one level to the other ( Figure 6 and Figure 8). This is an important evidence of the fact that the entire lithospheric block generating earthquakes beneath Vrancea region is subject of the same kind of forces resulted from a single stress field with homogeneous configuration.

Discussion of the results in correlation with principal boundaries in the crust and subducted slab
The main goal of our study is to evaluate the stress field in the Vrancea and to investigate how this field correlates with the seismicity distribution at crustal and intermediate depths.  Table 1).

Stress field pattern in the Vrancea zone
First, the presence of a transitional layer between the crust and the seismic active body in the upper mantle is important to be analyzed. This layer placed roughly between 40 and 60 km depth seems to be unable to generate significant seismic activity. For this reason, it was considered as a seismic gap area [Fuchs et al., 1979;Hurukawa et al., 2005] and lead to the hypothesis of an oceanic slab detached from the crust and descending in the mantle beneath Vrancea [Fuchs et al., 1979;Oncescu, 1984]. However, Sperner et al. [2001] adjusted this hypothesis considering that the slab is not yet fully detached, because otherwise the high strain rates inside the slab (which indicate strong slab pull forces) are hard to be explained. The transition region was interpreted as a zone of weakened mantle or lower crustal material where slab detachment is presently taking place.  Table 1. Tondi et al. [2009] considered this depth interval as a zone of partial melting which may result from delamination of the European mantle lithosphere and the upwelling of hot asthenospheric material. Also, the refraction seismic data suggest a low-velocity zone at a depth of 47 to 55 km beneath the Vrancea region [Hauser et al., 2007]. Although the depth to Moho boundary is established at 40 -42 km in Vrancea area, with an abrupt increase in Vp from 7.0 to 7.9 km/s, there is another deeper layer, between 45 and 55 km depth, which marks a reversal of Vp velocity from 8.0 km/s to 7.6 km/s. The structure under 55 km depth is characterized by a Vp = 8.5 km/s and is marking clearly the intrusion of a high-velocity slab in the lithosphere under the Moho. Given these data, we fixed the lower limit of the crustal domain to 50 km and the upper limit of the intermediate-depth seismic active body to 55 km depth.
The seismic activity in the bottom side of the seismogenic volume is sharply terminated at around 180 km depth, where probably a critical phenomenon inhibits the instability conditions able to generate brittle-like ruptures.
According to the ROMPLUS earthquake catalogue [Oncescu et al., 1999;INFP, 2020], only four isolated earthquakes have been recorded so far below 180 km, at depths between 186 and 218 km and with magnitudes between 3.2 and 4.1. It is worth mentioning that the seismogenic volume represents a relatively small part of a larger high-velocity body extended in depth to about 350 km, with a rather complicated shaping [Martin et al., 2006].
The last seismicity feature we planned to investigate is the nature of the boundary at 100 -110 km depth which makes transition from an upper active segment (VNI_A) to a lower active segment (VNI_B). According to some authors [Trifu and Radulian, 1994;Enescu and Enescu, 1998;Cărbunar and Radulian, 2011], the two segments are characterized by different seismicity patterns and triggering mechanisms and differences in the number of earthquakes and the energy released in both segments are established [Bala et al., 2001]. Notably, the seismic rate is approximately five times higher in the lower segment than in the upper segment, the average magnitude is systematically higher (by about 0.5 magnitude units) in the lower segment and the major shocks (M w > 6.5) occur alternatively in the two segments. This boundary is also observed by Raykova and Panza [2006] in a study dedicated to the seismic properties of the lithosphere-asthenosphere system for the South-Eastern Carpathians from a large tomography project in the South-Eastern Europe. In the cell 16E of the model of Raykova and Panza [2006], which corresponds roughly with the Vrancea seismic zone, the shear wave velocity in the slab is established between 4.45 This boundary considered by seismological researches at 100 -110 km depth might coincide with the lithosphere/asthenosphere boundary (LAB) as it was proposed by Enescu et al. [1992]. It coincides with the LAB position as considered for the regions surrounding the Vrancea region in more recent studies [Grinč et al., 2014]. The Moesian platform is well delineated by a relatively thick lithosphere, characterized by an E-W trend of lithospheric thickness decrease. In the east, near the Vrancea region, the thickness is about 110 km and in the west it is only 80 km. In the same time the thickness of the lithosphere of the East European Platform is on average more than 120 km, but in the northern part of this area some thicker places can be found Grinč et al. [2014].

Discussion of the results in correlation in terms of field orientation and relative ratio R
The main result of the present study is the homogeneity and consistency of the stress field regime along the entire Vrancea earthquake-prone body (55 -180 km in depth). First, the S 3 axis is close to vertical with a plunge angle of 84.40° ± 2.24° (average of the plunge angles computed by inversion for the five subcrustal layers given in the Table   2). The reduced variability of axis orientation, which is obtained also when considering the results for individual cells instead of layers, is surprising taking into account the inherent errors in determining the fault plane solutions. We can firmly conclude that the stress field in the Vrancea subcrustal seismogenic zone is predominantly downdip extension along the vertically descending body. As a second important result, we note the invariance of the orientation of the principal horizontal axes, with S 1 oriented perpendicularly and S 2 tangentially to the Carpathians Arc bend: 138.20° ± 8.35° (for S 1 ) and 47.60° ± 8.78° (for S 2 ). The average values of σ1 and σ2 are computed using the values on the five layers in the Table 2, converted to the second and first quadrant, respectively. We note that this result remains the same when we use the values computed per 10 cells instead of 5 layers. It is the first time that such a coherent stress field regime is outlined for the Vrancea intermediate-depth source. The result has been anticipated to some extent in other studies as well [for example, Radulian et al., 1999], but either on the basis of Andrei Bala et al.
Stress field pattern in the Vrancea zone empirical inferences than computations, or based on a reduced number of earthquake mechanisms over a large area [Oncescu, 1987]. Oncescu [1987] performed stress inversion in the Vrancea region by applying the method of Gephart and Forsyth [1984] for 27 Vrancea intermediate-depth earthquakes recorded between 1934 and 1986. The results pointed out the vertical orientation of the σ3 axis compatible with a predominant radial compression regime in the Vrancea. This is in agreement with the fault plane solutions showing preferred tension orientation (T axes) consistent with the vertical or close to vertical extension in the slab [Bala et al., 2003]. Whereas the S 3 axis orientation is consistent with the vertical tension stress pattern predicted by our study, the orientation of the horizontal axes, S 1 oriented NE-SW and S 2 oriented NW-SE, is practically rotated by 90 0 as compared with our results and contradicts the most common modelling with principal compressive axis oriented NW-SE [e.g., Radulian et al., 2000;Sperner et al., 2001;Lorinczi and Houseman, 2009]. We can explain this discrepancy in the orientation of the horizontal stresses by the great change in data quality nowadays versus 1986 and also on the reduced number of mechanisms used in the inversion (27) which in the present work is considered as a lower limit.
Interestingly, Oncescu (1987) obtained a low value of stress ratio (R = 0.2) close to the R value from our inversion in the lower Vrancea layer (R = 0.29). We note that from the total number of 27 events considered in Oncescu's analysis, 21 events are located in the lower Vrancea segment (VNI_B). In the same range, Protopopova and Botev [2019] determined for the Vrancea region a compressional stress regime, with R = 0.3.
In another study focused on a larger area centered on the Vrancea region, Müller et al. [2010] considered the maximum horizontal stress (σ H ) given by direct measurements in 94 boreholes with depths down to a maximum of 5.72 km, which were spread around the Vrancea zone. After applying a smoothing algorithm over a grid radius of 40 km on a data set of 84 records, they come up with a mean value of the azimuth of S H = 82.6° ± 55.7°, which is very close to our result of 84°, obtained by inversion applied to the focal mechanisms of 74 crustal events located in MO3-MO4, almost superimposed to the epicentral area of the Vrancea subcrustal earthquakes ( Table   2 and Table 1).
As we know, the value of ratio R expresses whether the value of σ2 is closer to σ1 or to σ3. If we consider the Equation 1 and if we scale the reduced tensor as followed: σ1 = 1; σ2 = 1 − 2R; σ3 = −1. Then we have for VNI_A the values: σ1 = 1; σ2 ~0; σ3 = −1; and for VNI_B: σ1 = 1; σ2 = 0.42; σ3 = −1, which means that the magnitudes of the principal stresses are not close to each other as it was proposed by Muller et al. [2010] for the crustal part.
The difference in R value between VNI_A and VNI_B involves some variation of the tectonic regime among the two segments. In VNI_A we have σ2 which is placed in the middle of the domain established by σ1 and σ3. While in VNI_B the σ2 is closer to σ1 and the tectonic regime is much more homogeneous due to close values of the two horizontal stresses. We conclude that the variation in seismic behavior between the two segments might be due to the difference in tectonic regimes surrounding the slab, not only due to the composition in the slab itself.

Models of the tectonic stress in the slab
The numerical modelling developed by Ismail-Zadeh et al., [2005a;2005b] begin to model the descending slab in the Vrancea zone based on the first model in the area obtained by Martin et al. [2001].
The modelling continued in the next stage [Ismail-Zadeh et al., 2007 starting from the 3D tomographic image obtained by Martin et al. [2006] for the high-velocity slab descending in the mantle. The upper-mantle seismic velocity anomalies are usually associated with significant temperature variations, so that the presence of a highvelocity body down to about 350 km depth (as coming out from seismic tomography experiment) implies lower temperature in the slab down to this depth relative to the neighboring asthenosphere. To evaluate temperature, Ismail-Zadeh et al. [2005b inverted the P-wave anomalies into temperature taking into account the effects of mantle composition, anelasticity and partial melting. The contrast of temperature between the cold high-velocity slab and the warm asthenosphere generates stress acting in the slab region. The numerical computations led to a 3D stress field with the maximum principal stress oriented horizontally and minimum principal stress oriented vertically, and with the highest stress values distributed in the seismically active part of the descending place, localized at depths of about 70-170 km. Therefore, the enhancement in the computed stress field coincides with the seismogenic volume beneath Vrancea, and the stress configuration is consistent with the stress determination based on the focal mechanisms as obtained in the present study. whereas the release of the seismic energy is greater in VNI_B than in VNI_A. To overcome this discrepancy, we assume that the elastic properties of the material in the lower segment VNI_B change significantly as a result of increasing pressure and temperature, some authors suggesting also that faulting due to metamorphic phase transitions [Green and Burnley, 1989] or dehydration-induced embrittlement [Hacker et al., 2003] may also play a role in the regional stress generation and release. This variation in the rheology of the slab is to some extend supported by the variation of stress ratio R and friction coefficient μ parameters as resulted from the stress inversion.

Stress regime modelling at the South-Eastern Carpathians Arc bend
In the Figure 9 we draw a schematic representation of the stress field regime in the study area following the results obtained by the inversion procedure applied in this paper. The first key element of the model is the predominant downdip extension regime with almost vertical elongation (S 3 ) and horizontal compression along the entire segment of the lithospheric slab which is seismically active (55 to 180 km in depth). The second key element is the orientation of the maximum compressive axis (S 1 ) and of the intermediate compressive axis (S 2 ) perpendicular to and respectively parallel to the Carpathians Arc bend. This configuration is present at any depth level with some slight variation in the orientation of the axes.
As regards the stress regime in the crustal area situated above Vrancea subcrustal source, we obtained thrust regime similar with the stress in the slab. A significant variation is observed in the orientation of the horizontal compressive axes (S 1 and S 2 ) which are rotated by about 90° relative to the configuration in the slab.
Andrei Bala et al.
18 This compressive regime is localized close to the Vrancea epicentral area (MO2-MO3 in our analysis) and turns into an extensional regime as we go away from this area.
The study of Tondi et al. [2009] interpreted the isolines of Vp/Vs values, visible between 45 and 180 km depth, in relation with the downgoing slab geometry. Thus, they interpreted where the limit of Vp/Vs = 1.65 which occur at 180 km depth as the possible transition from continental to oceanic crust in the slab itself. The model shows a steep zone delimited by the 1.70 isoline which might be the boundary which we considered at 110 km depth.
Although such a steep boundary could be suggested in our 3D representation of Figure 2, it is hard to be delimited, since it has probably a thickness of a few kilometers only. According to Radulian [2014], the region around 110 km depth might be a region in which mantle material has penetrated and the conditions for occurrence of earthquakes are considerably reduced. The considered region of depth might be a transition zone between two earthquake-prone segments, VNI_A (upper segment) located approximatively between 55 and 105 km depth and VNI_B (lower segment), located approximatively between 106 and 180 km depth, as considered in our paper (Table 2). Although the seismicity rate and magnitude (energy release) differ in the two segments by a factor of ~5 and ~ΔM 0.5 respectively, the orientation of the stress principal axes looks like to be practically the same all along the seismic active volume (within ± 5° which is at the level of errors

Conclusions
In the present study we computed the stress field pattern in the Vrancea area by inverting the earthquake focal mechanisms using an extended and complete set of fault-plane solutions available for the Vrancea intermediatedepth and crustal earthquakes recorded between 1929 and 2012 [Radulian et al., 2018.
The earthquake mechanisms for the Vrancea intermediate-depth earthquakes reveal predominant vertical tension and horizontal compression with variable directions [Heidbach et al., 2007]. However, the fault planes tend to align along NE-SW, following the elongated distribution in the same direction of seismicity, in accordance a NW-SE compression, perpendicular to the Carpathians Arc [Oncescu and Trifu, 1987;Radulian et al., 1999].  Table 2) for the VNI A and VNI B segments.
The data allowed us to perform the investigation on 3D gridding by splitting the seismogenic volume in six layers and 12 cells containing minimum 30 events each. This number threshold is considered adequate for statistical analysis to determine mean horizontal stress directions and stress regimes. For the subcrustal domain (55 -180 km) in all cases (5 layers and 10 cells) the stress pattern shows similar orientation of the stress axes: S 3 -vertical, S 1 -NE-SW oriented, S 2 -NW-SE oriented. Our result assumes a unitary stress behavior of the entire seismic active block descending in the mantle.
The same characteristics are obtained in the overriding crust, but slightly less constrained. Note that this behavior of the tectonic regime is observed in the crust just above and adjacent to the Vrancea slab projection. As shown in Bala et al. [2020], the stress regime changes significantly if we move away from the Vrancea area. Keeping the stress regime in the crustal domain above the sinking slab proves that the stress is transmitted along the slab to the surface and that the seismicity in the overriding crust should be somehow influenced by the seismicity in the slab. Therefore, our results are in favor of a model with certain coupling between the stress regime in the subcrustal domain and that acting in the overlying crust against a model of total break-off.
The regional stress field computed in the last version of the World Stress Map [Heidbach et al., 2016a[Heidbach et al., , 2016b indicates for the SE Romania a strong prevalent sub-horizontal orientation of SH max at crustal level which is close to our results for the MO3-MO4 seismic region: the maximum horizontal stress field oriented E-W to ENE-WSW (see Figure 7).
It is interesting to note that the similarity observed for the fault plane solutions of the major Vrancea earthquakes [e.g., Radulian, 2014] fits the principal mechanisms as they come out from the present inversion approach ( Figure   10). In this figure we selected the fault plane solution associated to the principal fault plane resulted from inversion with the nodal plane dipping toward NW closer to vertical. As shown by previous investigations on the focal mechanism for the Vrancea major events, the real fault plane is the one dipping toward NW (dip of 60°-70°). This is a strong argument to consider that the stress field regime acting upon the entire slab volume generating intermediate-depth earthquakes in Vrancea is compatible with a predominant faulting plane oriented NE-SW and inclined toward NW, plane which defines at the same time the seismicity pattern. In other words, the stress regime inferred by us from inversion of focal mechanisms is controlling the earthquake generation in the Vrancea intermediate depth domain both at small scale (small and moderate magnitude) and at major earthquake scale.
Data and sharing resources. The Romanian Earthquake Focal Mechanism Catalogue (REFMC 1929(REFMC -2012 is the basis for analysis. It can be downloaded and used according to its license using the Mendeley Data repositoryhttps://data.mendeley.com/datasets/mykkx4gygy . For stress inversion we used StressInverse [Vavryčuk, 2014] and MSATSI [Martinez-Garzon et al., 2014] software. For the computation of the directional cosines and fault instability (I) in Table 4 we used MohrPlotter v. 3.0 [Allmendinger, 2020]. Figures 1, 2, 3, 7 and 9 were plotted mostly using the ESRI ArcMap and ArcScene software. For the other figures we used Matlab along with Adobe Photoshop. For focal mechanisms we used the WOLFRAM Demonstration Project "Earthquake focal mechanism" of Scherbaum F., Kuehn N. and Zimmermann B. [2009]. Open content licensed under CC BY-NC-SA. All the computer packages are used only for scientific work which resulted in the present paper.