Earthquake Source Properties of a Lower Crust Sequence and Associated Seismicity Perturbation in the SE Carpathians, Romania, Collisional Setting

Romanian seismicity is mainly confined to the Eastern Carpathians Arc bend (ECAB), where strong subcrustal earthquakes (magnitude up to 7.9) are generated in a narrow lithospheric body descending into the mantle. The seismic activity in the overlying crust is spread over a larger area, located mostly toward the outer side of the ECAB. It is significantly smaller than subcrustal seismicity, raising controversies about possible upper mantle-crust coupling. A significant earthquake sequence took place in the foreland of the ECAB triggered on 22 November 2014 by a mainshock of magnitude 5.7 (the greatest instrumentally recorded earthquake in this region) located in the lower crust. The mainshock triggered a significant increase in the number of small-magnitude events spread over an unusually large area in the ECAB. The paper’s goal is to compute the source parameters of the earthquakes that occurred during the aforementioned sequence, by empirical application of Green’s function and spectral ratio techniques. Fault plane solutions are determined using multiple methods and seismicity evolution at regional scale is investigated. Our results highlight a still active deformation regime at the edge of the EE Craton, while the source parameters reveal a complex fracture of the mainshock and a very high-stress drop.


Introduction
The crustal seismicity developed in front of the SE Carpathians bending zone (SEC) is spread over the entire area between the Neogene Trotuş and Intra-Moesian faults (Figure 1). Crustal earthquakes are generated into a complex structure with significant lateral inhomogeneities involving different kinematics [1]. This is also reflected in the variability of the focal mechanisms that show consistent stress regime with different localization: normal faulting on the eastern and southern flank of the Focşani Basin (FB), thrust faulting in the area adjacent to the Vrancea subcrustal earthquakes and strike-slip solutions preferably observed along the Peceneaga-Camena Fault [2]. The crustal seismicity in the SEC partly overlaps with the intermediate-depth earthquakes generated in a depth range of 60-180 km, beneath the Vrancea Zone (Figures 1 and 2).   The intermediate-depth seismicity is located in a narrow seismogenic volume which descends almost vertically in the mantle, able to generate 2-3 strong events (Mw ≥ 7) per century. The seismicity in the overlying crust is limited to moderate-size events (Mw < 6) and spreads over a larger area following some specific alignments. One of the significant seismically active faults placed in front of the SEC is the Peceneaga-Camena (PCF) fault zone and related deformation, which delineates ( Figure 1) the boundary between the Moesian Platform (MP) and North Dobrogea Orogen (NDO, [5]). Along its path towards NW from the Black Sea Basin to the SEC, the fault shows two segments with distinct seismic activities: (1) from the Black Sea Basin to the Danube River, and (2) from the Danube River to the NE edge of the Vrancea epicentral area. The first segment is less active with a stress regime of thrusting faulting [6,7], while the second segment displays a significant seismic activity with predominant strike-slip faulting along the major fault line and its satellites, oriented SE-NW. Towards the NW edge, where the PCF intersects the Peri-Carpathian Fault, the seismicity configuration becomes more complex (Figure 1).
Several studies [7][8][9][10][11] indicate this area as the place where the descending lithospheric slab started to break off, beneath the SEC, at the time of the continental collision between the Carpathian (as part of Tisza-Dacia) plate overriding the East-European, Scythian and Moesian plates.
In this complex tectonic framework, a significant earthquake sequence occurred in November 2014 close to Mărăşeşti city. According to the Romanian earthquake catalogue ROMPLUS [4] permanently updated by the Romanian Data Centre (RONDC) of the National Institute for Earth Physics (NIEP), the mainshock occurred on 22 November 2014 at 19:14 (GMT), located in the lower crust (H = 41 km), in the range of the Moho discontinuity [12] and had a local magnitude (M L ) of 5.7. The occurrence of this sequence in the Mărăşeşti area is the most significant episode of seismic energy release in this region since instrumental recordings were used to monitor the seismic activity in Romania. A similar event is mentioned in the ROMPLUS catalogue but occurred more than one century ago (01.03.1894, M L = 6.1).
The poly-kinetic character of these events has been pointed out by previous studies [13,14] as well as through the analysis of several significant earthquake sequences [15,16] generated along the Eastern Carpathians [17][18][19] However, the 2014 Mărăs , es , ti seismic sequence differs from past ones, because of the location of a large part of the hypocenters close to the crust-mantle discontinuity and because of the subsequent increase of seismic activity spread over an unusually large area in the adjacent regions (Figures 1 and 2).
The analysis of the source parameters for the earthquakes of the seismic sequence will certainly bring relevant elements to the understanding of the tectonic processes and of the physical mechanism driving seismicity in the region. Moreover, the seismic hazard assessment strongly depends on estimation accuracy. One of the most debated parameters in such analysis is the stress drop. This gives indications about the amount of released energy and also may estimate the destructive potential of an earthquake. A major disadvantage is given by the estimation of accuracy since its uncertainties are rather high, due to the corner frequency determinations that could increase or decrease the stress drop value by at least three times [20].
Large earthquakes with magnitudes of M > 6 occurred in inter-plate or intraplate regions [21] and pointed out stress drops in a range of 3-10 MPa, with an increasing trend for intraplate earthquakes [22,23]. For crustal earthquakes occurring in the Eastern Carpathians' foredeep, previous studies [24,25] highlight stress drops in a range of 2.5-50 MPa. Another notable feature highlighted by previous research [26][27][28] reveals similar properties for small and large earthquakes, changing only in scale, which would imply a constant stress drop for the whole range of magnitudes. However, several studies have shown [20,29] that for small earthquakes (M < 3) there is an increasing trend towards stress drop at the seismic moment. One of the features that could cause this behavior is the strong attenuation effect within the upper crust, which acts as a filter, removing the high frequencies (f > 10 Hz) and limiting the corner frequencies' resolution [30], especially for small earthquakes, leading to an underestimation of the stress drop.
In this study, we determine source parameters of 21 earthquakes that occurred during the 2014 Mărăşeşti sequence (22 November 2014-1 February 2015), as well as for several events that have occurred at higher distances or were delayed in time, but considered as having been triggered by the aforementioned seismic sequence. To achieve this goal, we apply empirical Green's function (EGF) and spectral ratio (SR) techniques, since they allow suitable constraints independent of the structure between the source and receiver, which is difficult to know with sufficient accuracy. In addition, we determine the fault plane solutions for some of these events and investigate the evolution of seismic activity in the region to draw an interpretation in connection with the tectonic structures highlighted by the previous studies.
The previous investigations dealing with the 2014 Mărăşeşti sequence limited themselves to the evaluation of the macro-seismic effects [31] and the computation of the focal mechanisms for the main event and 32 aftershocks using P-wave polarities and amplitude ratios [32].

Sequence Description
The seismic activity in front of the SEC bending zone along the PCF up to the northwestern edge at the contact with the Vrancea epicentral area and Scythian Platform has continuously developed in time with sporadic weak-to-moderate earthquakes or clusters of events (Figure 1). These events are probably related to some extent to the geodynamics in the mantle beneath the Vrancea region. However, this is still a controversial issue [6,7,33] and therefore understanding the interdependence of the tectonic stress in the mantle and coupling and crustal seismicity requires further investigation.
The mainshock of 22 November was followed by 230 aftershocks located by the RONDC using the LOCSAT routine [34] embedded within the Antelope software and using IASP91 as 1-D global velocity model [35]. A noteworthy increase of seismicity in the region was recorded after the mainshock for at least three months which is significantly more than the duration typical for the aftershock activity of moderate-size mainshocks.
The aftershocks were numerous but unusually weak in size relative to the magnitude of the mainshock ( The seismic activity generated after the mainshock occurrence (22 November 2014) manifests through a significant increase in the seismicity rate first in the Mărăşeşti area and second over an extended region to the south and to the north (the Râmnicu Sărat and Adjud areas, see Figure 2). The triggered seismicity developed over a large region and for a long-time interval.

Methods and Data
To determine the source parameters of the earthquakes that occurred during the 2014 seismic sequence close to the Mărăşeşti, and also for the events triggered by this seismic activity, we searched in the ROMPLUS catalogue (November 2014-July 2015) for clusters of events (co-located events with similar waveforms). Each cluster has one or several main events and associated EGFs. The selection of clusters is carried out in two steps. To select the potential EGF events, we started with the mainshock (M1/5.7 in Figure 3) and searched among the events within a magnitude able to solve the source parameters within the frequency range of the available data. First, for each main event (M) we searched in the aforementioned catalogue for potential EGF events if the inter-hypocenter distance is less than 0.3 deg. and the difference in magnitude around one magnitude unit or more. As a second step, we applied the waveform cross-correlation technique to select as EGF events those with large cross-correlation coefficients (CC ≥ 0.7). A high cross-correlation for a pair of events is considered to consistently reflect the waveform similarity among them [28,36]. For the waveform correlation analysis, we used time windows based on the magnitude of the M [21], a length of 0.5-1 s being enough for the small size events. The waveforms are bandpass filtered to an appropriate frequency band, using a two-pole Butterworth filter. We low-pass filter both the M and each potential EGF at a frequency related to the expected corner frequency of the M. This is because large and small earthquakes are not expected to cross-correlate well at high frequencies (e.g., [36]). We applied cross-correlation separately for P and S waves.  Note that some EGF events can be associated with more M. M2, M3 and M5 are selected both as M and EGF (see Table 1).
We found 6 individual master events (M1 to M6 given below) and plot them in Figure  3 as dots filled in different colors and 17 EGF (empty circles, with the outline colors matching one of the master events, see Figure 3), meeting our requirements imposed by the method.   Table 1).
We found 6 individual master events (M1 to M6 given below) and plot them in Figure  To determine the source parameters of the selected events ( Figure 3 and Table 1), only the broadband recordings of the RSN were considered.
The fault plane solution has a significant relevance since it is directly used to determine the stress regime of a region. To obtain a high accuracy in determining the focal mechanisms, we selected only the events whose recordings (broadband or short-period stations) allowed us to read with enough accuracy the P-waves' first motion at a minimum of 12 stations. Two methods provided by the SEISAN algorithm [37] are applied. Both FOCMEC [38] and HASH [39,40] algorithms use the first P-wave polarities together with the amplitudes measured on all three components (vertical, transversal and radial). The FOCMEC method is more rigorous since it removes the erroneous polarities and ratios above an imposed limit (we adopted 0.3 in our case). HASH uses weights for the polarity and amplitude ratio errors and provides averaged values, without eliminating the stations with large errors. For both methods, we filtered the seismograms, removing the frequencies above the corner frequency, and measured the maximum amplitudes at intervals of maximum 2s starting at first P, then respectively S arrivals.
The spectral ratios and empirical Green's function methods have been used in various areas of the world [36,[41][42][43][44][45][46] showing reliable results in retrieving source parameters. These are efficient in eliminating the path, site and instrument effects, which distort to a large extent the ground motion observed at the surface when pairs of earthquakes are available and satisfy the conditions imposed for the selection of the events. Examples of earthquake pairs recorded by the Topalu (TLB) station are given in Figure 4. The SRs are approximated by a nonlinear least-squares procedure with theoretical spectral ratios, depending mostly on the source properties. The fit parameters are the ratios of the low-frequency asymptotes and corner frequencies of the large and small shocks, respectively. The least-squares best fit is estimated by iteration from an initial guess using a simplex algorithm. The relationship which best approximates the Fourier spectral ratio R for a pair of collocated events is: where Ω M 0 , Ω G 0 are the low-frequency spectral amplitudes of main and EGF earthquakes, and f c M , f c G are the corresponding corner frequencies. We used the hypothesis of the Brune's model [47,48] with spectral decay at the high frequencies, γ. The ratio of the low-frequency spectral levels is equal to the ratio of the seismic moments for the collocated earthquakes if they have similar fault plane solutions. The seismic moments for EGF events were determined using the spectral ratio value at low frequency as is given by the following equations: Once the seismic moment and source radius estimates are obtained, we compute the stress drop values using the relation of [47]: To determine the source time function (STF) we applied in parallel the method of empirical Green's functions, deconvolving EGF functions from M events. The deconvolution functions were obtained by spectral division for different time windows (function of the size of the event) for both P and S waves, using a tapering of 10%. To stabilize the spectral division, the amplitude spectrum of EGF was smoothed out, before division, with a five-point running moving window. In addition, to reduce high-frequency instabilities, the deconvolved source time function was band-pass filtered between 0.5 and 10 Hz with The size of the rupture area is computed, taking into consideration the circular fault assumption, directly from the corner frequency [47,48]: where r represents the source radius, K is a constant (K P = 0.32 for P waves and K S = 0.21 for S waves), f c is the corner frequency (for P and S waves) and V is the velocity of P and S waves in the epicentral region. The determinations were performed for both P and S waves using vertical, or respectively horizontal, components of velocity recordings. Because the spectral ratio technique is a relative method, it cannot estimate the absolute values of the seismic moments for the two earthquakes, computing only their ratio. Based on the spectral analysis we determined P and S displacement spectra separately and seismic moments of the main events by using the following equation: where ρ is the density at the source depth, v is the velocity of P or S waves at source depth, Ω 0 is the long period displacement spectral level, R is the hypo-central distance, F (=2) is the free-surface parameter and R θ ϕ is the source radiation pattern (average values of 0.52 for P waves and 0.63 for S waves, according to [49]). The density and velocities adopted for the study area were taken from [12] (V P = 7.0 km/s, V S = 4.04 km/s and ρ = 2.85 g/cm 3 ). The seismic moments determined from P and S wave recordings are rather close, therefore, we estimated M 0 for each station as the mean value and as a final value the median of all selected stations was considered. Based on M 0 we also compute the moment magnitude (M W ) using the equation of [50]: The seismic moments for EGF events were determined using the spectral ratio value at low frequency as is given by the following equations: Once the seismic moment and source radius estimates are obtained, we compute the stress drop values using the relation of [47]: To determine the source time function (STF) we applied in parallel the method of empirical Green's functions, deconvolving EGF functions from M events. The deconvolution functions were obtained by spectral division for different time windows (function of the size of the event) for both P and S waves, using a tapering of 10%. To stabilize the spectral division, the amplitude spectrum of EGF was smoothed out, before division, with a five-point running moving window. In addition, to reduce high-frequency instabilities, the deconvolved source time function was band-pass filtered between 0.5 and 10 Hz with a Butterworth filter [45]. To increase the accuracy of the results, we used EGF pairs recorded by broadband stations by each component and stacked them up. The source rise time and source duration (τ, τ1/2) for the main events were computed from the STFs each time it had a pulse-like shape. The corner frequency (fc) was determined from the source duration using equation (8) from [20] and finally we computed the source radius using formula (2).
In addition, to have a better image of the tectonic processes in the region we computed the Gutenberg-Richter relation (see [51] and Equation (9)) before, during and after the seismic sequence, for an area extending to a distance of 0.75 deg. around the epicenter of the mainshock. The analysis was performed with Zmap software [52] using the maximum likelihood method [53]. For computation, we take into account local magnitude provided by the ROMPLUS catalogue, estimated as a maximum Wood-Anderson amplitude of one of the horizontal components corrected for distance and compiled only for stations with a good signal/noise ratio (SNR ≥ 3).
N(M) is the number of earthquakes with magnitude larger or equal to M, a is a constant and b is the scaling parameter.

Results
The results obtained for fault plane solutions, source parameters and evolution of seismic activity within the study area are in agreement with the findings pointed out by previous studies. By increasing the number of seismic stations in the area, we succeeded in increasing the results accuracy and also in adding complementary information that could improve the seismotectonic characteristics of the region.

Focal Mechanisms
The fault plane solutions determined by applying FOCMEC and HASH techniques are displayed in Figure 5 and listed in Table 2, where we also list the associated errors and the amount of data used. Therefore, for the FOCMEC method the amplitude ratio fit AF, no. of bad polarities NBP and no. of bad amplitude ratios NBAR are provided; while for the HASH routine the fault plane uncertainty FPU, auxiliary fault plane uncertainty AFPU, fit errors F-fit and station distribution ratio STDR are listed. All the fault plane solutions are well constrained by the selected seismic recordings. In addition, the focal mechanisms determined using both methods are comparable, increasing their reliability degree, and are also in agreement with the results pointed out within the study area by past studies [54][55][56]. For the mainshock, the fault plane solution indicates normal faulting with both nodal planes NW-SE oriented, in the direction of the Peceneaga-Camena Fault zone. Our determination shows a good similarity with the solutions computed by several international centers (INGV, USGS, GFZ), as well as with the solution of [32]. The compression axis is plunging almost vertically, while the extension axis is almost horizontal W-E oriented. According to [57], the fault-plane solutions can be grouped into six types of faulting ( Table 2) Table 2. Beach-balls are displayed as a function of magnitude. The circles show earthquake epicenters with outline color as a function of depth.

Spectral Ratios
The spectral ratios method was applied to all earthquake pairs, the main event-empirical Green's function, listed in Table 3. The source parameters in Equation (1) are the asymptote at low frequency (in logarithmic scale), a, the corner frequencies of main (fc M ) and EGF (fc G ) events. The corner frequencies are determined by approximating the computed spectral ratio with the theoretical ratio and averaging the results for all the available earthquake pairs and stations. The analysis was performed for both P and S-waves, taking into account only the recordings where the signal-to-noise ratio was above 2. The results obtained for M and EGF events are given in Table 4 and examples of spectral ratios computed for pairs of earthquakes are displayed in Figure 6. Note that the obtained results show that low-frequency level, a, and the corner frequency of the EGF (fc G ) show a higher variation since they depend on the earthquake pair, while for the main event the corner frequency (fc M ) should not be influenced by the pair of the events if source directivity ef-  Table 2. Beach-balls are displayed as a function of magnitude. The circles show earthquake epicenters with outline color as a function of depth.

Spectral Ratios
The spectral ratios method was applied to all earthquake pairs, the main eventempirical Green's function, listed in Table 3. The source parameters in Equation (1) are the asymptote at low frequency (in logarithmic scale), a, the corner frequencies of main (f c M ) and EGF (f c G ) events. The corner frequencies are determined by approximating the computed spectral ratio with the theoretical ratio and averaging the results for all the available earthquake pairs and stations. The analysis was performed for both P and S-waves, taking into account only the recordings where the signal-to-noise ratio was above 2. The results obtained for M and EGF events are given in Table 4 and examples of spectral ratios computed for pairs of earthquakes are displayed in Figure 6. Note that the obtained results show that low-frequency level, a, and the corner frequency of the EGF (f c G ) show a higher variation since they depend on the earthquake pair, while for the main event the corner frequency (f c M ) should not be influenced by the pair of the events if source directivity effects are negligible. This feature is similar to the results pointed out by previous studies [58].
The corner frequencies ratio of P to S waves (f c P /f c S ) averaged over all stations for the main events varies between 0.94 and 1.06 for the M1 and M2 events, respectively, with a mean value of 1.05 computed for all selected events (see Table 4). These results show similarities with the ones pointed out for various regions of the world [45,46,[59][60][61][62], and follow the circular source model proposed by [63].
Commonly, the P source pulse shows a shorter duration than the S source pulse, which leads to higher corner frequencies for P waves because the propagation time of P waves is faster compared with S waves. Thus, an average value of 1.5 is expected for the f c P /f c S when the waves leave the fault at angles of incidence greater than 30 • from the normal. For cases where the propagation of the rays is close to the normal direction of the fault, the P and S waves leave the source in-phase and interference due to the source limitation is negligible (f c P ≈ f c S ). Using equations from (1) to (7), we determined the seismic moment, source radius and stress drop for the main and EGF events (Tables 5 and 6). The source radius estimated from corner frequency is an average between the estimations using P-wave and S-wave corner frequencies. Since relation (2) assumes a ratio of about 1.5 for f c P /f c S and for our data f c P ≈ f c S , the radius computed from f c P is systematically greater than the radius computed from f c S roughly by a factor of 1.5 (see Tables 5 and 6).  Using equations from (1) to (7), we determined the seismic moment, source radius and stress drop for the main and EGF events (Tables 5 and 6). The source radius estimated from corner frequency is an average between the estimations using P-wave and S-wave corner frequencies. Since relation (2) assumes a ratio of about 1.5 for fc P /fc S and for our data fc P ≈ fc S , the radius computed from fc P is systematically greater than the radius computed from fc S roughly by a factor of 1.5 (see Tables 5 and 6).

Empirical Green's Function Analysis
For the same pairs of events considered in the SR method, we applied in parallel the method of deconvolution with EGF to compute for M the STFs, rise time (τ1/2), source duration (τ), corner frequency (fc) and source radius (r) using Equations (2) and (8). To remove any biases, we computed the relative source time function for a given station by summing the functions resulting from deconvolution from all three components of that station. The obtained values are listed in Table 7 and examples of STFs resulting from the deconvolution method are displayed in Figure 7. The source radius and stress drop values are equivalent to the same values inferred by the SR method. The only exception is the case of the mainshock where the differences are~30% smaller for source radius and~40% larger for stress drop.

Scaling Relationships
The source parameters determined in the present analysis and their scaling relationships are investigated as indicators of geotectonic peculiarities of the Carpathians foredeep region as in previous papers [16,24,25,45,46,64]. Up to now, there has been no systematic investigation of the source parameters for the earthquakes that occurred close to Mărășești city. The scaling of M0 with ML is displayed in Figure 8. The data are well constrained (correlation coefficient of 0.97) by the regression line:

Scaling Relationships
The source parameters determined in the present analysis and their scaling relationships are investigated as indicators of geotectonic peculiarities of the Carpathians foredeep region as in previous papers [16,24,25,45,46,64]. Up to now, there has been no systematic investigation of the source parameters for the earthquakes that occurred close to Mărăs , es , ti city. The scaling of M 0 with M L is displayed in Figure 8 [68,69] pointed out, for these magnitudes, seismic moments with values smaller by up to an order of magnitude.
Scaling of the seismic moment with source radius is shown in Figure 9 for both P and S waves and is approximated by the regression lines: study is similar to those computed by the international seismological centers (GFZ, EMSC, CSM, NEIC). A second possible explanation for this deviation is given by an overestimation of seismic moments of small earthquakes (ML < 4). Several studies achieved in different regions [68,69] pointed out, for these magnitudes, seismic moments with values smaller by up to an order of magnitude. Scaling of the seismic moment with source radius is shown in Figure 9 for both P and S waves and is approximated by the regression lines:    study is similar to those computed by the international seismological centers (GFZ, EMSC, CSM, NEIC). A second possible explanation for this deviation is given by an overestimation of seismic moments of small earthquakes (ML < 4). Several studies achieved in different regions [68,69] pointed out, for these magnitudes, seismic moments with values smaller by up to an order of magnitude. Scaling of the seismic moment with source radius is shown in Figure 9 for both P and S waves and is approximated by the regression lines:    The slope of the regression line is higher than the theoretical value for scaling the seismic source with a homogeneous rupture process (compare dashed with solid lines in Figure 9). Since in all these cases the same relative deconvolution techniques were applied, we assume that the deviation from a theoretical scaling (slope~4 instead of 3) is a consequence of the underestimation of the corner frequency for the smaller earthquakes.  Figure 10 we can notice both for P and S-wave determinations an increasing trend of stress drop with the seismic moment which may suggest a higher slip rather than a higher fault dimension [19]. Acoustics 2021, 3 FOR PEER REVIEW 19 The slope of the regression line is higher than the theoretical value for scaling the seismic source with a homogeneous rupture process (compare dashed with solid lines in Figure 9). Since in all these cases the same relative deconvolution techniques were applied, we assume that the deviation from a theoretical scaling (slope ~ 4 instead of 3) is a consequence of the underestimation of the corner frequency for the smaller earthquakes. The scaling of stress drop with the seismic moment ( Figure 10) indicates an increasing value stress drop over the entire magnitude range. The stress drop values are in the range 1-41 MPa for P waves and 4-101 MPa for S waves. However, the increasing trend of the distribution may be apparent as the decisive element is given by the comparatively higher value of the stress drop in the case of the mainshock of 22 November 2014 (for the rest of the events, if we remove the largest event, a constant stress drop scaling seems to be reasonable). The statistics are too small for larger earthquakes to have an eloquent image of scaling. In Figure 10 we can notice both for P and S-wave determinations an increasing trend of stress drop with the seismic moment which may suggest a higher slip rather than a higher fault dimension [19]. Figure 10. Scaling of stress drop with the seismic moment for P (blue) and S (purple) waves. Figure 10. Scaling of stress drop with the seismic moment for P (blue) and S (purple) waves.

Frequency-Magnitude Distribution
The analysis of the frequency-magnitude distribution at the scale of the entire region (0.75 × 0.75 deg.) outlines a typical distribution with the b slope close to 1 for 'normal' time intervals (see Figure 11a,d).

Frequency-Magnitude Distribution
The analysis of the frequency-magnitude distribution at the scale of the entire region (0.75 × 0.75 deg.) outlines a typical distribution with the b slope close to 1 for 'normal' time intervals (see Figure 11a,d). We performed the analysis for four times intervals, two of them situated outside the perturbed activity before the mainshock of 22 November 201422 November (200522 November -2011, one for a three year interval before the mainshock occurrence (2011-2014), one for 100 days starting with the mainshock occurrence and the last one for a five year interval after the sequence (2015-2020). To compute the b value, we use the maximum likelihood method [53]. The magnitude of completeness (Mc) of 2.1 and the b slope of 1.14 characterize the first interval

Discussion
We discuss in the following the results of the present work taking into account the tectonic and geophysical characteristics of the study region as well as comparison with the results of previous studies.
The space-time and frequency-magnitude analyses on the crustal seismic events that occurred within the study region in 2014-2015, during and after the seismic sequence of Mărășești, show that the seismic activity is distributed in several clusters of earthquakes ( Figure 2). The mainshock (ML = 5.7), located close to the Moho boundary (Figure 13), triggered a significant increase in seismic activity, not only in the source neighborhood but spread over a larger area, as compared with its source dimension. It was also unusual that this intensification consisted of a relatively large number of small-magnitude earthquakes spanning a long-time interval.

Discussion
We discuss in the following the results of the present work taking into account the tectonic and geophysical characteristics of the study region as well as comparison with the results of previous studies.
The space-time and frequency-magnitude analyses on the crustal seismic events that occurred within the study region in 2014-2015, during and after the seismic sequence of Mărăs , es , ti, show that the seismic activity is distributed in several clusters of earthquakes ( Figure 2). The mainshock (M L = 5.7), located close to the Moho boundary (Figure 13), triggered a significant increase in seismic activity, not only in the source neighborhood but spread over a larger area, as compared with its source dimension. It was also unusual that this intensification consisted of a relatively large number of small-magnitude earthquakes spanning a long-time interval. The mainshock occurrence first generated an aftershock activity located in the epicentral area, near Mărășești city. Then an intensification of seismic activity was observed in the Râmnicu Sărat area (~50 km south of Mărășești, see Figures 1 and 2). Apart from the small size seismic activity triggered there, a moderate earthquake occurred in the same The mainshock occurrence first generated an aftershock activity located in the epicentral area, near Mărăs , es , ti city. Then an intensification of seismic activity was observed in the Râmnicu Sărat area (~50 km south of Mărăs , es , ti, see Figures 1 and 2). Apart from the small size seismic activity triggered there, a moderate earthquake occurred in the same area (12 January 2015, 06:08 GMT, M L = 4.2). An increase in seismic activity has also been emphasized towards the north (~30km relative to Mărăs , es , ti), close to Adjud city ( Figure 2). Similar to the previous case, but delayed in time, a moderate size earthquake occurred close to Adjud city (29 June 2015, 22:20 GMT, M L = 4.0). Intensification of the seismic activity was pointed out to the East and South-East of Mărăs , es , ti, but here only small-magnitude events were recorded. Overall, the epicenters distribution within the study region follows a direction parallel to the Eastern Carpathian Bend zone, as was also shown by previous studies [2,7,24,25,73]. However, if the seismicity is investigated at a smaller scale (for specific areas) the shape of the epicenters' distribution may change ( Figure 5). The depth distribution (Figures 2 and 13) emphasizes higher depths for the events that occurred close to the Mărăs , es , ti city with most of them located within the lower crust, while the ones located in the adjacent regions were placed in the middle or upper crust (Figures 2, 3 and 13). The depth estimations of these events are rather confident, even for the low size events, since all the locations were manually performed and the seismic station coverage is good (Figure 1) as a consequence of a continuous improvement of the seismic network in the study area [74] (Figure 12). The depth distribution gives us valuable indications about the depth of the faults. The mainshock of the sequence is located in the Focs , ani Basin, at~5 km, west of the Peceneaga-Camena Fault, a deep crustal fault with a Moho offset of~5 km [75]. It separates the Moesian Platform from the buried part of the North Dobrogea orogen since Neogene [76]. In our analysis, the epicenters distribution seems to follow the characteristics described by [70,77] which showed, based on reflection seismic lines, that in the eastern side of the Focs , ani Basin the deformations occur along numerous normal faults in a wide area with clear topographic expression (see also [1]).
The frequency-magnitude distribution for the earthquakes of the Mărăs , es , ti sequence ( Figure 11c) shows a magnitude of completeness of M L = 1.2 and a slope of the distribution (0.99) similar to the normal value (b~1.0) [78]. This is a good indication that these earthquakes are of tectonic origin. On the other hand, we notice important variations of b value in time, for the first period (2005-2011), the b value is about 1.14, increasing up tõ 1.4 in the second period (2011-2014), and returning close to the normal (b~1.04) for the last period (2015-2020). The changes of b value could be caused by several factors such as geological conditions, degree of heterogeneity of cracked medium, and strain and stress variation in the region [79,80]. The higher b value before the sequence could be explained by a relative blocking of the stress on moderate-size resisting areas (asperities) and intensification instead of small-size earthquakes, taking into account the high heterogeneity of the region [1,81]. The decrease of Mc after 2014 is directly related to the seismic network improvements to monitor seismicity in the region (Figure 12).
Our results for the fault plane solutions reveal three new mechanisms determined with rather good accuracy for events 14, 18 and 19 (see Table 2). The overall mechanisms highlight a wide variety for earthquakes associated or triggered by the 2014 Mărăs , es , ti seismic sequence, although most of these events occurred in the Focs , ani Basin, a thick Neogene foredeep sediment layer as thick as 13 km located near the Carpathian bend zone [70]. Although the region is relatively narrow, this variety of fault plane solutions is not uncommon for the SE Carpathians foredeep area as shown by other studies [7,31,32,54,55] due to the geological complexity of this region, pointed out on one hand by the collision of three major units with different geometries and characteristics, all representing cratonic continental platforms [82,83] and on the other hand by the significant subsidence, although it is still uncertain whether these faults influenced the lateral variations in the thin-skinned thrusting kinematics [83].
We noticed particular characteristics for each group of events. For the moderate-size earthquakes that occurred near Mărăs , es , ti city, normal fault plane solutions are prevalent with the nodal planes oriented in the NW-SE direction, corresponding to the PCF direction. The small size earthquakes that occurred in the same region and in the same depth range have in general the same type of faults, highlighting acceptable deviations of the focal planes, probably caused by a reduced amount of data used in their determination. Conversely, earthquakes that occurred in the upper depth segments (but still at the level of the lower or middle crust) tend to change their fault plane solutions towards a normal left lateral mechanism, which may involve activation of secondary faults of the same system (PCF). Given that the epicenters are distributed parallel with the advancing orogen, a flexural rupture of the lower plate with partial reactivation of inherited normal faults might be an interpretation for the earthquakes located further south of the PCF (Figure 13b). The normal faulting highlighted by our results for this area was also emphasized by the results of previous studies [6,77]. In addition, [1] pointed out a large number of normal faults with offsets in the order of tens of meters at the base Quaternary level which belong to the system of the PCF. The moderate size earthquakes generated near Adjud city (~30 km northern part of Mărăses , ti city) were located in the upper crust (Figures 1, 2 and 13), showing a tendency to change the mechanism to strike-slip. At least one of the nodal planes is oriented in the same direction (NW-SE) with PCF, which may indicate an activation of a secondary faults system located towards its northern side and at the upper crust level dominated by horizontal movements.
The faults pattern observed and described above can be explained by the differences in localization deformation transferred from the advancing Carpathian Orogen to the plate. These differences are generated by the sharp changing rheology and mechanical behavior along inherited lower plate blocks and their obliquity to the advancing orogen [1], such as the case of the major border faults TF and PCF separating the Scythian Platform/NDO from Moesian Platform ( Figure 13). In the case of the TF, which is almost perpendicular to the advancing orogen, the inherited structures are not inverted, but instead form a new left lateral strike-slip system (with associated normal and reverse faults) located above the contact between strong Scythian vs weak Moesian Platform. The geometry of the PCF, orientated roughly at a 45 • angle from the advancing orogen, makes it susceptible for inversion during orogenic loading and flexure, with the normal faults located at the contact between NDO and Moesian Platform.
The earthquakes from the Râmnicu Sărat region (~50 km south relative to Mărăs , es , ti) were also located in the upper crust, with a predominant reverse mechanism, with the nodal planes oriented in the NE-SW direction, parallel to the Eastern Carpathians' bending. This type of mechanism was observed in a large number of earthquakes that occurred in the past in the same area [6,25,84], the existence of a possible fault along the Eastern Carpathians' bending being a debated topic [85].
These NE-SW change in fault type pattern can be explained by the ongoing compressional deformation of the SE Carpathians fold and thrust belt which is responsible for both low angles thrusting on the leading edge orogen fault (Pericarpathian thrust Fault, see Figures 1 and 13) and high angle reverse faults, probably related to the inversion of previous Mesozoic normal faults as suggested by [71]; however, the exact geometry of the faults is not resolvable with the current available refraction seismic data (see also [76]).
It is worth noticing that for the event that occurred in 12012015, 06:08 (GMT), M L = 4.2, our results indicate a reverse focal mechanism, with the nodal planes NE-SW oriented, which contradicts the normal fault plane solution determined by [32]. We assume that the fault mechanism determined by the above study is flawed due to the event location, which in their study is reported as being located~40 km towards the north, and at a greater depth (with~20 km) as compared to the location provided by the ROMPLUS catalogue.
Previous results have shown that source parameters determined using SR and EGF techniques have a better accuracy since the determination uncertainties decrease by averaging the data [23,35,45]. Our data set consists on average of five EGFs for each main event, allowing an uncertainty decrease of up to 20%. The results obtained by using two different techniques (SR and EGF deconvolution) highlight rather close values with an average difference of~30% (for the stress drop) which indicates an adequate degree of accuracy. The largest event in our dataset has an average stress drop of 134 MPa. Similar high-stress drop values were observed in other regions for normal-faulting earthquakes [86] or for the deeper location of the events [87]. At the same time, we noticed an increasing trend of stress drop with seismic moment.
According to our analysis, the STF has a simple pulse-like form for all the studied events. However, the apparent pulse-like form for the source of the largest event as obtained by EGF deconvolution is not so well constrained as for the other smaller main events. This could be attributed to some complexity of the rupture process in the source, visible in the waveform recordings of the mainshock. If we assume a multi-shock release of seismic energy for this event, the resulting stress drop is a dynamic value rather than static, which explains the unusually high-stress drop values obtained for the mainshock (Tables 5 and 7). Source parameter scaling relationships fit well with the results previously obtained for the Carpathian foredeep region [24,25,45,64]. The deviations observed for the seismic moment-source radius scaling (slope~3.9 instead of 3) and seismic moment-stress drop (deviation from a constant stress drop scaling) from the theoretical values could be due to the limitation of the empirical Green's function and spectral ratios techniques to constrain the corner frequency for the weakest earthquakes and the anomalous high-stress drop value of the mainshock. On the other hand, the high attenuation effects within the upper crust might act as a filter, removing the high frequencies content (f > 10 Hz) and limiting the corner frequencies resolution [28] affecting especially small earthquakes and leading to an underestimation of their stress drop. It is also likely that the usage of these methods for the small earthquakes to lead to a saturation of the corner frequency due to the significant increase of the recorded noise at high frequencies. Because of this, we can expect that our analysis is no longer able to discern the higher corner frequencies and, for this reason, the source dimension becomes overestimated and, correspondingly, the stress drops underestimated. We need more high-quality data for different size earthquakes to validate or not the scaling deviations obtained in the present work.

Conclusions
The seismic sequence triggered on 22 November 2014 close to Mărăs , es , ti city offered the possibility of evaluating the source parameters and the focal mechanism solution for the strongest crustal event (M L = 5.7) ever recorded in the area since the instrumental recordings are available.
The focal mechanism solution of the mainshock, obtained from local, regional or global data, indicates normal faulting with the compression axis almost vertical and expansion axis almost horizontal, in agreement with the seismotectonic of the eastern and southern flanks of the Focs , ani Basin. Both nodal planes are oriented NW-SE, in the direction of the PCF. Going toward the south, in the Râmnicu Sărat area, reverse faulting is also obtained, which is characteristic of this area as highlighted in previous investigations.
Applying SR and EGF methods the source parameters for 15 earthquakes of the sequence of November 22, 2014 (2.0 ≤ M L ≤ 5.7) and for five earthquakes produced in the first seven months of 2015 (1.8 ≤ M L ≤4.2) were determined.
The source parameter values determined through SR and EGF methods are close, indicating an increasing trend of stress drop with seismic moment, induced especially by the high value of the mainshock. This seems to explain the not very well constrained pulse-like shape obtained by EGF deconvolution for the mainshock which could indicate a multi-shock release of seismic energy for this event.
Investigation of the frequency-magnitude distribution at regional scale outlines a 'normal' regime with b-slope value around 1, if we are situated sufficiently far from the perturbed activity related to the occurrence of the strong crustal earthquake of 22 November 2014, and a possible anomaly (b = 1.4) before the occurrence of this event. A behavior with pre-shock increase of b from about 1 to 1.5 and sudden decrease during the aftershock activity and further on is typical for acoustic experiment tests and for critical systems in which the population of cracks concentrates after the critical phase along preferential paths.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.