Aftershock analysis and forecasting for the crustal seismicity in Romania

doi:10.21203/rs.3.rs-4799375/v1

Romania is known for its persistent seismicity at intermediate-depths in the Vrancea region, however crustal areas are also a significant source of seismic hazard, although large shallow events are less common. This study is a first attempt to characterize statistically and propose a forecasting model for two recent aftershock sequences occurred at crustal depths in 2014 and 2023, following mainshocks of moderate magnitudes (M_w5.4 for both mainshocks). We apply a robust approach based on a state-of-the-art procedure developed and tested previously for Japan, which is able to determine in quasi real-time the parameters of the Gutemberg-Richter law and Omori-Utsu law for aftershocks and provide probability estimates of larger events, which can be updated in real time. For both the 2014 Vrancea-Marasesti and 2023 Gorj sequences we test several relatively short (hours to day) learning periods and subsequent forecasting periods. Both sequences are characterized by relatively high b-values (~ 1.2), obtained for all tested learning periods, which may point out to the release of stress following the mainshocks or the presence of crustal fluids in the studied regions. The aftershock decay is characterized by a parameter p of around 1.0, commonly observed for crustal aftershock sequences. The c-value, which indicates the onset time of the power-law decay of aftershocks, is on the order of minutes to hours. Although the two studied sequences follow mainshocks having the same magnitude, the probability of larger aftershocks for the 2023 Gorj sequence are larger. The results obtained in this study are encouraging for the development of a real-time monitoring and forecasting system for the Romanian crustal seismicity.

real time aftershock forecasting

seismic sequences

crustal seismicity

Romania

Studying the space-time evolution of earthquake sequences is of importance for understanding their underlying physical processes, like migration of fluids and slow slip (e.g., Shimojo et al., 2021) within the seismogenic zone, as well as for seismic hazard mitigation. The importance of reducing damages and losses after larger mainshocks supports the efforts being made to develop dedicated forecasting tools. Existing methods have proven effective in forecasting seismicity, both at the beginning (Reasenberg and Jones, 1989; Page et al., 2016) and during (Omi et al., 2013, 2014, 2015, 2016, 2019) active aftershock sequences.

In general, earthquake forecasts are given based on the data that is revised and prepared at the time of a study or the data transmitted in real-time into a semi-automatic system. The latter type is more likely to experience forecasting accuracy problems, due to a higher predisposition to errors. These errors may include a lack or reduced number of seismic recording seismic stations in the study area, seismic sensor malfunctioning, overlapping earthquakes on seismograms for certain time periods, especially at short times following a mainshock, miscalculated magnitudes, and errors of epicenter determination (e.g., Omi et al., 2016 and references therein).

Romania has both crustal and subcrustal seismicity. The Vrancea region (Romania), is well-known for its energetic intermediate-depth seismicity (e.g., Enescu et al., 2023), with the largest events having magnitudes above 7.0, like the destructive March 4, 1977 M_w7.5 earthquake (e.g., Fuchs et al., 1979). Although less energetic and less persistent compared to the Vrancea subcrustal earthquakes, the crustal events in Romania are often generated as moderate seismic swarms and foreshock-mainshock-aftershock sequences (Popescu and Radulian, 2001, Popescu et al., 2011, Popescu et al., 2012), which, due to their relatively shallow hypocentral locations pose a considerable seismic hazard.

In the last decade, the National Seismic Network of Romania has benefited from an increasing number of seismic stations, which led to a significant increase in the number of earthquakes detected and located on the territory of Romania and neighboring regions (Marmureanu et al., 2021), as recorded in the ROMPLUS seismic catalog (e.g., Oncescu et al., 1999; Popa et al., 2013), prepared and provided by National Institute for Earth Physics (NIEP).

The purpose of this study is twofold. First, we focus on the statistical analysis of two crustal seismic sequences associated with the largest recent events that occurred in the recent years on the Romanian territory, in the Vrancea-Marasesti area and Gorj-Oltenia region, respectively (Fig. 1). The mainshocks in both regions had a moment magnitude of Mw5.4 and occurred at shallow depths, as described in detail in the next section. A related purpose of this study is the forecasting of the aftershock activity in the two regions and the development of an aftershock forecasting system that would be the first of its kind for the territory of Romania. The presented system for automated aftershock forecasting in Romania targets generating results from the first few hours to the first few days after a mainshock, using previously proposed statistical techniques that have been applied successfully to the seismicity of Japan (Omi et al., 2013, 2014, 2015).

Earthquake data for the Marasesti and Gorj earthquake sequences

One of the largest crustal events, M_W5.4 (M_L5.7) instrumentally recorded in Romania occurred in the Vrancea-Marasesti (or simply Marasesti) area (Fig. 2a), at a depth Z = 40.9 km, on November 22, 2014, at 21:14 (Ghita et al., 2021) and was followed by an aftershock sequence consisting of 222 recorded earthquakes, with M_L ≥ 0.1, that spanned a period of about 70 days from the mainshock. The sequence occurred in the Focsani basin, part of the Moesian Platform (a major structural unit of the Carpathian and Balkans foreland, e.g., Seghedi et al., 2005), and is related to the normal fault system associated to the major Peceneaga-Camena fault, which separates the Moesian Platform from the North Dobrogea promontory (Craiu et al., 2019).

Another intense earthquake sequence occurred on February 2023 in the Gorj area, northwest of Targu Jiu city (Fig. 2b), situated in the northwestern part of Central-South Charpatian zone (Oros et al., 2023), at the contact between the Getic Depression, to the South, and the Southern Carpathians orogen, to the North. The sequence started with a magnitude M_W4.8 (M_L5.2) at a depth Z = 18.3 km foreshock on February 13, 2023, at 16:58:09, followed by the M_W5.4 (M_L5.7) mainshock on February 14, 2023, at 15:16:52, the largest earthquake ever recorded or reported in the area, at a depth Z = 14 kmApart from this seismic event, the second largest magnitude (M_W5.2) earthquake occurred in this area is the one from northwest of Targu Jiu, at a depth of 9.9 km, on June 20, 1943. It is worth noting that another sequence was recorded in the study area in 2012 (Radulian et al., 2014), with smaller magnitude earthquakes (maximum magnitude of M_W3.8), which occurred on reverse faults located eastward from the current sequence (Fig. 1). Starting February 22, 2023, due to the installation of three mobile seismic stations in the Gorj epicentral zone: BVPR, CPSR, DBRR (Fig. 2b), a better location of the aftershocks was achieved, together with the decrease of the magnitude of completeness, Mc, in the area.

For the analysis, we used the real-time data-portal earthquake catalog (https://dataportal.infp.ro) and the ROMPLUS earthquake catalog (Oncescu et al., 1999) for the Gorj and Marasesti areas, respectively. Our analysis and forecasts are based on the early stages of aftershocks sequences. The real-time system component of the data-portal is based on the Antelope data acquisition and processing software for handling large amounts of real-time recorded data (Popa et al., 2015). Because the data-portal real-time catalogue has the magnitudes given in M_L, from now on, in our study we will use only the M_L magnitude. Local time is used throughout the study.The forecasting parameters for the Marasesti area were estimated using seismic events from the ROMPLUS catalog (Oncescu et al., 1999; Popa et al., 2013), between November 22, 2014, at 21:14, and February 1, 2015, at 23:31. ROMPLUS has been used since the data-portal was not fully developed at the time of the sequence. The Gorj area dataset contains the location parameters and magnitudes of the foreshock, mainshock and other 2175 earthquakes with M_L ≥ 1.2 recorded from February 13, 2023, at 14:58:09 to March 24, 2023, at 4:02, as retrieved from the data portal.

For the Vrancea-Marasesti sequence (Fig. 2a), we refer to the crustal events with depth 18.6 ≤ Z ≤ 50.0 km, and magnitudes 0.1 ≤ M_L ≤ 5.7, recorded from 22 November 2014 to February 1, 2015, in the area bordered by 45.67^o N – 46.17^o N latitude and 26.83^o E – 27.5^o E longitude (Fig. 2a) (Ghita et al 2019).

For the Gorj sequence (Fig. 2b), we refer to the crustal events with depth 1.4 ≤ Z ≤ 26.3 Km, and magnitudes 0.1 ≤ ML ≤ 5.7, recorded from February 14, 2023, to March 20, 2023, in the area bordered by 44.93o N – 45.27o N latitude and 22.80o E – 23.40o E longitude (Fig. 2b, with a large potential of significant regional earthquake activity, as documented in previous studies (Placinta et al., 2016; Ghita et al., 2019). In the last years, studies of Bala et al. (2020) and Oros et al. (2023) have referred to the earthquakes in this area as belonging to a crustal seismogenic source named Central-South Carpathians (CSCSZ)(Fig. 2b).

In the present study we estimate the Omori-Utsu law and Gutenberg-Richter law (see below for their definition) parameters for the two earthquake sequences introduced in the previous sub-section, using various time periods as learning periods, and forecast the aftershock activity for periods that follow the learning periods. The approach used in this paper follows closely that proposed by Omi et al. (2013).

The Omori-Utsu law (Utsu, 1961; Utsu et al., 1995) is an empirical relation that describes the decay of the aftershock rate following a mainshock, according to a non-stationary Poisson process, as expressed by the formula:

\(\:\text{N}\left(\text{t}\right)=\:\frac{\text{K}}{{(\text{t}+\text{c})}^{\text{p}}}\) (1),

where N(t) represents the aftershock rate function of the time t (in days) after the mainshock, while K, c and p are parameters dependent on the earthquake sequence. The parameter K is a constant that sets the scale of the aftershock sequence (i.e., describes the productivity of the sequence), c is the time delay before the onset of the power-law decay of aftershock activity with a slope given by the parameter p, representing how fast the aftershock activity is decaying with time, t, from the mainshock.

The Gutenberg-Richter law (Gutenberg and Richter, 1944), which expresses the frequency-magnitude distribution (FMD) of earthquakes, can be defined as:

\(\:\text{m}\left(\text{M}\right)=\text{A}\times\:{10}^{-bM}\propto\:\text{e}\text{x}\text{p}(-\beta\:M)\) (2),

where m(M) is the frequency of earthquakes (normalized number of events) of magnitude M and A and b (or β = b*ln10) are constants. In a statistical framework, m(M) is the intensity of events with magnitude M.

The joint occurrence rate of aftershocks, λ, as a function of time from the mainshock, t, and magnitude M is represented by the product of N(t) and m(M). So, by combining N(t) and m(M) we get:

\(\:{\lambda\:}\left(\text{t},\:\text{M}\right)=\frac{{K}^{{\prime\:}}}{{(t+c)}^{p}}\beta\:{e}^{-\beta\:(M-{M}_{0})}\) (3),

where K’ = β⁻¹KA exp (−βM₀), and M₀ is the magnitude of the mainshock used to adjust the scale of K’ (Utsu 1970). Hereafter, to simply notations, we use K instead of K’ (e.g., Morikawa et al., 2021).

Note that the occurrence rate of aftershocks for earthquakes with magnitudes equal or larger than a threshold magnitude Mt can be obtained by integrating (3) as follows:

\(\:{\lambda\:}\left(\text{t}\right)={\int\:}_{{M}_{t}}^{\infty\:}\lambda\:\left(t,M\right)dM\) (4).

Many small aftershocks are incompletely recorded in seismic catalogs at early times after the mainshock (e.g., Utsu et al., 1995; Kagan, 2004; Peng et al., 2007; Enescu et al., 2009). In order to address the incompleteness issue, the detection rate as a function of magnitude (e.g., Ogata and Katsura, 1993, Omi et al., 2013) is used:

\(\:{\upvarphi\:}\left(\text{M}\right|{\mu\:},\:{\sigma\:})=\:\frac{1}{\sqrt{2{\pi\:}{{\sigma\:}}^{2}}}\:{\int\:}_{-\:{\infty\:}}^{\text{M}}\text{d}\text{x}\:\text{e}\text{x}\text{p}\:[-\frac{{(\text{x}-\:{\mu\:}(\text{t}\left)\right)}^{2}}{2{{\sigma\:}}^{2}}]\) (5),

where µ represents the magnitude with a 50% detection rate and σ represents a the scale for a range of magnitudes around m where earthquakes are more or less frequently detected. In order to obtain the „corrected”, real frequency-magnitude distribution (i.e., the distribution that would be observed if there were no detection issues), the observed frequency-magnitude distribution m(M) (Eq. 2) is multiplied with the detection rate function φ(M|µ,σ) (Eq. 4). The parameters β, µ and σ can be simultaniously determined (e.g., Ogata and Katsura, 2006) by fitting the magnitude data of all detected earthquakes to a probability density that is the normalized version of m(M)*φ(M|µ,σ), through a maximum likelihood procedure. Note however that the detection rate strongly depends on the elapsed time since the mainhock (e.g., Enescu et al., 2009; Peng and Zhao, 2009).

Omi et al. (2013) have developed a state-space model for an objective Bayesian inference of the time-dependent parameter µ(t) and employed a Gaussian prior probability distribution of the b-value for robust estimation using short data from an early period. Given the estimate of the b-value and detection rate function, one can next estimate the parameters K, c and p of the Omori–Utsu formula for the underlying aftershocks by maximizing a log-likelihood function (for details we refer to Omi et al., 2013). Finally, we can calculate the expected number of underlying aftershocks, Nt, with magnitudes equal or above a threshold magnitude Mt, as well as the predictive interval (Omi et al., 2013) for a forecast window spanning from t1 to t2, by integrating (4):

\(\:{N}_{t}={\int\:}_{{t}_{1}}^{{t}_{2}}\lambda\:\left(t\right)dt\) (6).

In our computations, we have chosen Mt as the magnitude of completeness, Mc, for the data in the learning period. Note that while the magnitude of completeness, Mc, is not explicitly given by (5), it can be estimated using the formula:

\(\:{M}_{C}\left(n\right)=\:\mu\:+n\sigma\:\) (7),

where n indicates the confidence level: n = 0 means that 50% of the events are detected above M_c, while n = (1, 2, 3) means that 68%, 95% and 99% of the events are detected (e.g., Mignan and Woessner, 2012). In this study we chose n = 2 in order to assure a feasible estimation.

Parameter estimations and afterhocks’ number forecasting for the 2014 Vrancea-Marasesti earthquake sequence

We present first the Omori-Utsu and Gutemberg-Richter laws parameter estimations for the Vrancea-Marasesti earthquake sequence, as well as the aftershock forecasts for this sequence. To perform calculations that are meaningful for real-time forecasting, we chose relatively short learning periods of 3 hours (Fig. 3), and 6, 9, 12 and 24 hours (Additional File 1 – Figures S1-S8) after the mainshock. We forecasted the cumulative number of aftershocks for the same number of respective hours following the learning periods.

By using a learning period [0 3h], we estimated the following Omori-Utsu law (Eq. 1) parameters: K-value = 8.17*10 − 4, c-value = 0.0276 and p-value = 1.03. The β-value of the Gutenberg-Richter law (Eq. 2) is 2.70.

The b-value calculated from the β-value estimated for the first 3 hours of aftershock activity in the area is approximately equal to 1.17, which is slightly higher than the average value of ~ 1.0, observed for worldwide seismicity, indicating a higher proportion of smaller aftershocks that may reflect the release of stress in the region. Indeed, an increased b-value after the mainshock is reported for some well observed aftershock sequences and interpreted in terms of decreased differential stresses after the mainshock (e.g., Tormann et al., 2015, Gulia et al., 2018). Alternatively, a higher b-value could indicat the presence of fluids in the area (Tormann et al., 2015). The c-value of 0.0276 days indicates that the power-law onset of the aftershock decay starts after approximately 40 minutes from the mainshock, while the p-value of 1.03 indicates a commonly observed decay rate of aftershocks (e.g., Utsu et al., 1995).

Figure 3a shows the cumulative number of aftershocks for the learning and forecasted windows, as well as expected number of afterhocks during the forecasted period. It can be noticed that the observed aftershock number is within the 95% confidence interval and close to the forecasted aftershock number. The frequency-magnitude distribution during the forecasted period is also in good agreement with the estimations based on the learning period.

Figure 4a presents the scatter plot matrix for the estimated parameter sets, in the case of the Vrancea-Marasesti sequence. Figure 4b shows a slow decrease of the detection magnitude µ(t) (i.e., the magnitude for which 50% of the earthquakes are detected) during the 3 hours learning period, from M1.6 to M1.3.

Several other learning periods have been used for the parameter estimation and forecasting and the results are summarized in Table 1. In general the parameter estimations prove to be stable, without large variations from case to case. The forecasted number of aftershocks is in general within the 95% confidence limits. Table 3 summarizes the detection parameters obtained for the same learning periods.

Table 1

Summary of parameter estimations for the Vrancea-Marasesti sequence, including their errors, according to the corresponding learning periods.
Learning period and related figs.	K	p	c [days]	β	b-value
0–3 h (Figs. 3 and 4)	\(\:8.17\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:1.11\:\times\:\:{10}^{-4})\)	\(\:1.03\) (\(\:\pm\:0.12\))	\(\:2.76\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:1.66\:\times\:\:{10}^{-2})\)	\(\:2.70\)	\(\:1.17\)
0–6 h (Figs. S1, S2)	\(\:6.83\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:1.14\:\times\:\:{10}^{-4})\)	\(\:1.03\) (\(\:\pm\:0.11\))	\(\:2.75\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:1.77\:\times\:\:{10}^{-2})\)	\(\:2.75\)	\(\:1.19\)
0–9 h (Figs. S3, S4)	\(\:6.39\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:0.69\:\times\:\:{10}^{-4})\)	\(\:1.08\) (\(\:\pm\:0.11\))	\(\:2.56\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:1.35\:\times\:\:{10}^{-2})\)	\(\:2.71\)	\(\:1.17\)
0–12 h (Figs. S5, S6)	\(\:6.29\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:0.64\:\times\:\:{10}^{-4})\)	\(\:1.09\) (\(\:\pm\:0.10\))	\(\:2.54\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:1.19\:\times\:\:{10}^{-2})\)	\(\:2.71\)	\(\:1.17\)
0–24 h (Figs. S7, S8)	\(\:4.93\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:0.80\:\times\:\:{10}^{-4})\)	\(\:1.12\) (\(\:\pm\:0.10\))	\(\:2.57\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:1.15\:\times\:\:{10}^{-2})\)	\(\:2.75\)	\(\:1.19\)

Table 2

Summary of parameter estimations and their uncertainties for the Gorj-Oltenia sequence according to the corresponding learning periods.
Learning period and related figs.	K	p	c [day]	β	b-value
0–3 h (Figs. 5, 6)	\(\:3.69\:\times\:\:{10}^{-3}\:\) \(\:(\pm\:0.37\:\times\:\:{10}^{-3}\:)\)	\(\:1.01\) \(\:(\pm\:0.14)\:\)	\(\:4.67\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:4.22\:\times\:\:{10}^{-2}\:)\)	\(\:2.77\)	\(\:1.2\)
0–6 h (Figs. S9, S10)	\(\:3.53\:\times\:\:{10}^{-3}\:\) \(\:(\pm\:0.44\:\times\:\:{10}^{-3}\:)\)	\(\:1\) \(\:(\pm\:0.12)\)	\(\:8.32\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:7.08\:\times\:\:{10}^{-2}\:)\)	\(\:2.84\)	\(\:1.23\)
0–9 h (Figs. S11, S12)	\(\:2.58\:\times\:\:{10}^{-3}\:\) \(\:(\pm\:0.28\:\times\:\:{10}^{-3}\:)\)	1.01 \(\:(\pm\:0.11)\)	\(\:7.93\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:3.67\:\times\:\:{10}^{-2}\:)\)	2.91	1.26
0–12 h (Figs. S13, S14)	\(\:2.18\:\times\:\:{10}^{-3}\:\) \(\:(\pm\:0.23\:\times\:\:{10}^{-3}\:)\)	\(\:1.0\) \(\:(\pm\:0.12)\)	\(\:8.44\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:3.96\:\times\:\:{10}^{-2}\:)\)	\(\:2.98\)	\(\:1.29\)
0–24 h (Figs. S15, S16)	\(\:6.37\:\times\:\:{10}^{-4}\:\) \(\:(\pm\:4.95\:\times\:\:{10}^{-4}\:)\)	\(\:1.04\) \(\:(\pm\:0.10)\)	\(\:2.76\:\times\:\:{10}^{-2}\:\) \(\:(\pm\:16.68\:\times\:\:{10}^{-2}\:)\)	\(\:3.30\)	\(\:1.43\)

Table 3

Mean detection magnitudes (µ), standard deviations (\(\:\sigma\:\)) and threshold magnitude (M_th) for the Vrancea-Marasesti and Gorj-Oltenia sequences.
Seismic sequence	Learning period	µ	\(\:\sigma\:\)	M_th
Vrancea -Marasesti	0–3 h	1.56	0.08	1.8
	0–6 h	1.57	0.06	1.7
	0–9 h	1.57	0.06	1.7
	0–12 h	1.57	0.06	1.7
	0–24 h	1.56	0.06	1.7
Gorj -Oltenia	0–3 h	2.03	0.03	2.1
	0–6 h	1.96	0.06	2.1
	0–9 h	1.98	0.05	2.1
	0–12 h	1.79	0.05	1.9
	0–24 h	1.96	0.06	2.1

For the learning period of 6 hours ([0 6h]), the β-value has a small increase to 2.75, leading to a b-value of approximately 1.19, once again suggesting that small aftershocks occur more frequently in the Vrancea-Marasesti area. Like for the previously used learning period [0 3h], aftershock activity decreases following a power-law with an onset of approximately 40 minutes from the mainshock (c = 0.0275 days), at a „regular” decay rate, given by p = 1.03 (Figure S1b). In addition, for the forecasted period [6h, 12h], the difference between expected and observed cumulative aftershock numbers is significantly larger than for the [3h, 6h] forecasted period, the lower limit of the 95% confidence interval being close to the observed aftershock counts (Figure S1a). Figure S2b shows a slow decrease of the detection magnitude µ(t) during the 6 hours of the learning period, from an initial estimate M1.6 to M1.3, then an increase to M1.4. Such variations may reflect detection and location issues dependent on the number of stations that recorded the ongoing sequence.

The learning period of 9 hours from the mainshock consists of 117 events and is characterized by a β-value of 2.71, equivalent to b-value = 1.17, reconfirming the higher frequency of smaller aftershocks for this sequence. The aftershock decay parameters are c-value = 0.0256 and p-value = 1.08 (Figure S3b). The lower limit of the 95% confidence interval for the forecast period [9h, 18h] is also close to the observed cumulated aftershocks, while the difference between forecasted and observed aftershock counts (Figure S3a) is smaller than for the [6h, 12h] forecasted period.

For all tests on the Vrancea-Marasesti sequence, the estimated c-value is relatively small, consistent with well-observed aftershock sequences (e.g., Peng et al., 2007; Enescu et al., 2009) and may reflect early aftershock incompleteness. The p-value remains close to the average of ~ 1.0, which indicates typical stress relaxation. The b-value is relatively large (~ 1.2), which might hint to the release of stress (e.g., Gulia et al., 2018) and/or the presence of crustal fluids in the area, as they are often associated with an increased b-value (e.g., Tormann et al., 2015). Note that fluids might be present at crustal depths in both study areas (e.g., Oros et al. 2023).

The K-value of the Omori-Utsu law has the same order through all estimations (Table 1), although it slightly decreases with time. Longer estimation periods should lead to more robust estimates.

Parameter estimations and afterhocks’ number forecasting for the 2023 Gorj earthquake sequence

The earthquake catalog for the Gorj sequence is retrieved from the data-portal and consists of 2068 events. Among these, 66 aftershocks correspond to the first 3 hours of activity after the mainshock and produce a β estimation of 2.77, resulting in a b-value approximately equal to 1.2, similar to the β and b-value estimations made for the Vrancea-Marasesti sequence. The afteshock activity decreases at a commonly observed decay rate (p ~ 1.0), but from an onset time of about 67 minutes from the mainshock (c = 0.0467 days) (Fig. 5b). A threshold magnitude of 2.1 (determined as the magnitude of completeness during the learning period) produces an expected cumulative aftershock number that is relatively close to the observed one, above the lower limit of 95% confidence interval (Fig. 5a).

Figure 6a presents the scatter plot matrix for the estimated parameter sets, in the case of the Gorj sequence. Figure 6b shows a slow decrease of the detection magnitude µ(t) during the 3 hours of the learning period, from M2.3 to M1.9.

In the same way as for the Vrancea-Marasesti sequence, we have made parameter estimations for longer learning periods (Table 2) and and forecasted the number of aftershocks for subsequent forecasting periods (Additional File 1 – Figures S9-S16).

The learning period [0 6h] has 108 events and is characterized by a β-value of 2.84 (b-value of 1.23). The aftershock decay is characterized by p = 1.01, while the c-value gets slighlty larger (c = 0.0832 days (119.8 minutes)). A threshold magnitude of 2.1 produces an expected cumulative aftershock number that is close to the observed one, being above the lower limit of the 95% confidence interval (Figure S9a). Figure S10b shows a slow decrease of the detection magnitude µ(t) during the 6 hours of the learning period, from M 2.0 to M 1.7.

When analysing the case of the [0, 9h] learning period, the number of aftershocks is 141, chaeacterized by a β-value of 2.91 (b-value of 1.26). The aftershock decay is characterized by p = 1.01 and c = 0.0793 (114 minutes) (Figure S11b). A threshold magnitude of 2.0 produces expected cumulative aftershock numbers that are similar to the observed ones (i.e., above the lower limit of 95% confidence interval, Figure S11a). Figure A12b shows the detection magnitude µ(t) remains constant and equal to threshold magnitude M = 2.0 for approximately 0.05 days (72 minutes) after the mainshock, then decreases to 1.7 at the end of the 9 hours interval.

Due to catalog incompleteness and the lack of mobile seismic stations until February 22, 2023, the productivity, K-value, for the Gorj-Oltenia sequence presents a relatively large variation when comparing the values for the 12-hour and 24-hour learning periods (Figures S13 and S16, respectively). In addition, the b-value is larger than for the Vrancea-Marasesti area, hinting at an overall higher frequency of smaller earthquakes (Table 2) and probably the involvement of fluids. The p-values estimated for all learning periods are around 1.0 and c-values are on the similar order as for the Vrancea-Marasesti sequence. The results are relatively robust relative to the choice of the learning period and the forecasted number of earthquakes agree well with the observed numbers for the respective forecasting periods. Values utilized for threshold magnitude calculation, along with the errors, are highlighted in Table 3.

Aftershock occurrence probability estimations for the two aftershock sequences

In this section we estimate the probability of aftershocks for both regions, function of magnitude. For our analysis, we chose learning periods of activity after the respective mainshocks of 3 h in Figs. 7, and of 6, 9, 12 and 24 hours in the Additional File 1 – Figures S17 – S20. The forecasting was performed for the next 3, 6, 9, 12 and 24 hours after the learning periods, respectively.

We found that overall aftershock probabilities were lower for the Vrancea-Marasesti sequence for all forecasting periods. Aftershocks of M_L5.7 or larger in the Vrancea-Marasesti region were forecasted with probabilities close to 0.1%, while the same magnitudes were determined to have probabilities ranging between 1% and 10% for the Gorj sequence, after the first 3 hours from the mainshock.

Aftershock probabilities remain lower for the Vrancea-Marasesti region for all forecasting periods, with magnitudes equal or above the mainshock magnitude M_L5.7 having occurrence probabilites close to zero, while the same range of magnitudes has a probability of at least 0.1% for the Gorj region.

Resulting probabilities for the 9-hour learning period are similar to the ones obtained for the 6-hour learning period for both areas, with the Gorj area remaining more susceptible to stronger aftershocks.

The results obtained in this study highlight the importnnce of using real-time forecasting for the Romanian aftershock sequences. Despite having the same mainshock magnitudes, the two sequences evolve rathet differently, thus pointing out that a generic model would not be appropriate to characterize and forecast the two sequences.

Limitations of the approach proposed in this study and how these could be tackled in real-time applications

Figure 8 shows schematically the concept of our forecasting system for the crustal sequences in Romania. In the present configuration of the system, the Omori–Utsu law is used as a temporal model to issue short-term forecasts immediately after a mainshock. While this simple configuration proves rather robust, it does not capture the secondary or higher-order triggering of aftershocks. To solve this in the future, we will examine the performances of the ETAS model for real-time data and consider replacing in the future the Omori–Utsu modeling with the ETAS model for longer forecasting periods. We may also consider extending our forecasting model to a space–time model to provide spatial forecasts that would be useful for aftershock sequences, as in the case of the 2023 Gorj sequence.

We have proposed for the first time a relatively simple, but robust, analysis and forecasting workflow for the Romanian seismicity, in particular aftershock sequences that occur in crustal regions on the territory of Romania, based on an approach that is currently used for quasi real-time monitoring and forecasting of aftershock sequences in Japan. This approach has been applied to two recent aftershock sequences occurred in 2014 and 2023, in Vrancea-Marasesti and Gorj regions, respectively. Both sequences show a relatively high b-value (around 1.2) of the Gutenberg-Richter law, which may indicate the presence of crustal fluids or that the stress has been released in the studied regions. In both cases the estimated and observed number of aftershocks, for the chosen forecasting windows, are in good agreement. We also present probabilities for the ocurrence of aftershocks equal or larger than a certain magnitude. Such estimations may help the risk assessment efforts and contribute to better informing the population on future aftershock hazards. Our results are encouraging for the real-time monitoring and forecasting of shallow aftershock sequences in Romania.

NIEP stands for the “National Institute for Earth Physics”; ROMPLUS stands for the Romanian earthquake catalogue created and maintained by NIEP; CSCSZ stands for the “Central-South Carpathians”.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable.

Availability of data and materials

The dataportal belongs to the National Seismic Network (RSN) and can be accessed at https://dataportal.infp.ro. The automatic catalog is provided by the National Institute for Earth Physics. The revised ROMPLUS catalog is available at the NIEP website (https://www.infp.ro). They were last accessed on March 31, 2024. We used the following GitHub codes for aftershock forecasting: https://github.com/omitakahiro/AftFore. We also used of the Zmap software (Wiemer, 2001) for some of the computations.

Competing interests

The authors declare that there are no competing interests regarding this manuscript.

Funding

This research is supported by the JSPS KAKENHI Grant Number JP21H05206 in Transformative Research Areas (A) “Science of Slow to Fast Earthquakes”. It also benefited from funding by the AFROS Project PN-III-P4-ID-PCE-2020-1361, 119 PCE/2021, supported by UEFISCDI, Romania, the NUCLEU Program SOL4RISC, Project no PN23360201, supported by MCID, and PNRR-DTEClimate Project nr. 760008/31.12.2023, Component Project Reactive, supported by Romania - National Recovery and Resilience Plan. The funding bodies provided had no specific role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript should be declared.

Authors' contributions

Cristian Ghita: Data curation, Formal Analysis, Writing - Original Draft, Visualization, Data acquisition; Bogdan Enescu: Conceptualization, Formal Analysis, Writing - Review & Editing, Supervision, Funding acquisition; Alexandru Marinus: Formal Analysis; Iren-Adelina Moldovan: Writing - Review & Editing, Supervision; Constantin Ionescu: Supervision; Eduard Gabriel Constantinescu: Help with the preparation of map figures. All the authors discussed and approved the final version of this manuscript.

Acknowledgements. We thank our NIEP colleagues for providing the catalog data used in this work and fruitful discussions. We are grateful to Dr. Takahiro Omi for providing his aftershock analysis and forecasting packages nd guidance in using his routines.

Bala A, Radulian M, Toma-Danila D (2020) Crustal stress partitioning in the complex seismic active areas of Romania. Acta Geodaetica et Geophysica 55:389–403. doi:10.1007/s40328-020-00299-0.
Bora DK, Borah K, Mahanta R, Madhab J (2018) Seismic b-values and its correlation with seismic moment and Bouguer gravity anomaly over Indo-Burma ranges of northeast India: Tectonic implications. Tectonophysics 728–729:130-141. doi:10.1016/j.tecto.2018.01.001.
Craiu A, Ghita C, Craiu M, Diaconescu M, Mihai M, Ardeleanu L (2019) The source mechanism of the seismic events during the sequence of the moderate-size crustal earthquake of November 22, 2014 of Vrancea region (Romania). Ann Geophys 61. doi:10.4401/ag-7617.
Enescu B, Ito K (2001) Some premonitory phenomena of the 1995 Hyogo-ken Nanbu earthquake: seismicity, b-value and fractal dimension. Tectonophysics 338:297-314. doi:10.1016/S0040-1951(01)00085-3.
Enescu B, Ghita C, Moldovan IA, Radulian M (2023) Revisiting Vrancea (Romania) Intermediate-Depth Seismicity: Some Statistical Characteristics and Seismic Quiescence Testing. Geosciences 13(7):219. doi:10.3390/geosciences13070219.
Fuchs K, Bonjer KP, Bock G, Cornea I, Radu C, Enescu D, Jianu D, Nourescu A, Merkler G, Moldoveanu T, Tudorache G (1979) The Romanian earthquake of March 4, 1977 ii. Aftershocks and migration of seismic activity. Tectonophysics 53(3–4):225-247. doi:10.1016/0040-1951(79)90068-4.
Ghita C, Diaconescu M, Moldovan IA, Oros E, Constantinescu EG (2020) Spatial and temporal variation of seismic b-value beneath Danubian and Hateg-Strei seismogenic areas. Rom Rep Phys 72:704.
Ghita C, Diaconescu M, Raicu R, Moldovan IA, Rosu G (2021) The analysis of the seismic sequence started on November 22, 2014 based on ETAS model. Rom Rep Phys 73:708.
Gulia L, Rinaldi AP, Tormann T, Vannucci G, Enescu B, Wiemer S (2018) The effect of a mainshock on the size distribution of the aftershocks. Geophys Res Lett 45. doi: 10.1029/2018GL080619.
Gutenberg B, Richter CF (1944) Frequency of Earthquakes in California. Bull Seismol Soc Am 34:185-188. doi: 10.1785/BSSA0340040185.
Hajra S, Hazarika D, Shukla V, Kundu A, Pant CC (2022). Stress dissipation and seismic potential in the central seismic gap of the north-west Himalaya. J Asian Earth Sci 239:18-51. doi: 10.1016/j.jseaes.2022.105432.
Hussain H, Shuangxi Z, Usman M, Abid M (2020) Spatial Variation of b-Values and Their Relationship with the Fault Blocks in the Western Part of the Tibetan Plateau and Its Surrounding Areas. Entropy (Basel). doi: 10.3390/e22091016.
Iwata T (2008). Low detection capability of global earthquakes after the occurrence of large earthquakes: Investigation of the Harvard CMT catalogue. Geophys J Int 174:849–856. doi: 10.1111/j.1365-246X.2008.03864.x.
Main IG, Meredith PG, Jones C (1989) A Reinterpretation of the Precursory Seismic b Value Anomaly from Fracture Mechanics. Geophys J Int 96:131-138. doi:10.1111/j.1365-246X.1989.tb05255.x.
Marmureanu A, Ionescu C, Grecu B, Toma-Danila D, Tiganescu A, Neagoe C, Toader V, Craifaleanu IG, Dragomir CS, Meita V, Liashchuk A, Dimitrova L, Ilies I (2021) From National to Transnational Seismic Monitoring Products and Services in the Republic of Bulgaria, Republic of Moldova, Romania, and Ukraine. Seismol Res Lett 92. doi: 10.1785/0220200393.
Mignan A, Woessner J (2012) Estimating the magnitude of completeness for earthquake catalogs. Community Online Resource for Statistical Seismicity Analysis. doi: 10.5078/corssa-00180805.
Nanjo KZ, Hirata N, Obara K, Kasahara K (2012) Decade-scale decrease in b value prior to the M9-class 2011 Tohoku and 2004 Sumatra quakes. Geophys Res Lett 39. doi: 10.1029/2012GL052997
Ogata Y, Katsura K (1993) Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues. Geophys J Int 113:727–738. doi: 10.1111/j.1365-246X.1993.tb04663.x.
Olsson R (1999) An Estimation of the Maximum b-Value in the Gutenberg-Richter relation. Geodynamics 27:547-552. doi: 10.1016/S0264-3707(98)00022-2.
Omi T, Ogata Y, Hirata Y, Aihara K (2014) Estimating the ETAS model from an early aftershock sequence. Geophys Res Lett 41:850–857. doi: 10.1002/2013GL058958.
Omi T, Ogata Y, Hirata Y, Aihara K (2013) Forecasting large aftershocks within one day after the main shock. Sci Rep 3:2218. doi: 10.1038/srep02218.
Omi T, Ogata Y, Hirata Y, Aihara K (2015) Intermediate-term forecasting of aftershocks from an early aftershock sequence: Bayesian and ensemble forecasting approaches. J Geophys Res 120:2561–2578. doi: 10.1002/2014JB011456.
Omi T, Ogata Y, Shiomi K, Enescu B, Sawazaki K, Aihara K (2016) Automatic aftershock forecasting: A test using real-time seismicity data in Japan. Bull Seismol Soc Am 106:2450–2458. doi: 10.1785/0120160100.
Omi T, Ogata Y, Shiomi K, Enescu B, Sawazaki K, Aihara K (2019) Implementation of a Real-Time System for Automatic Aftershock Forecasting in Japan. Seismol Res Lett 90. doi: 10.1785/0220180213.
Oncescu MC, Marza VI, Rizescu M, Popa M (1999) The Romanian earthquake catalogue between 984–1997. Vrancea Earthquakes: Tectonics, Hazard and Risk Mitigation. doi: 10.1007/978-94-011-4748-4_4.
Oros E, Popa M, Paulescu D, Placinta AO, Ghita C (2023) A refined catalogue of focal mechanisms for the Intra-Carpathian region of Romania: implications for the stress field and seismogenic features assessment. Ann Geophys 65. doi: 10.4401/ag-8387.
Page MT, van der Elst N, Hardebeck J, Felzer K, Michael AJ (2016) Three ingredients for improved global aftershock forecasts: Tectonic region, time-dependent catalog incompleteness, and intersequence variability. Bull Seismol Soc Am 106:2290–2301. doi: 10.1785/0120160073.
Placinta AO, Popescu E, Borleanu F, Radulian M, Popa M (2016) Analysis of source properties for the earthquake sequences in the South-Western Carpathians (Romania). Rom Rep Phys 68: 1240–1258.
Popa M, Chircea A, Dinescu R, Neagoe C, Grecu B, Borleanu F (2022) Romanian earthquake catalogue (ROMPLUS). Mendeley Data 2022 V2. doi: 10.17632/tdfb4fgghy.2.
Popa M, Radulian M, Ghica D, Neagoe C, Nastase E (2015) Romanian Seismic Network Since 1980 to the Present. Nonlinear Mathematical Physics and Natural Hazards. doi: 10.1007/978-3-319-14328-6_9.
Popescu E, Borleanu F, Rogozea M, Radulian M (2012) Source analysis for earthquake sequence occurred in Vrancea (Romania) region on 6 to 30 September 2008. Rom Rep Phys 64
Popescu E, Neagoe C, Rogozea M, Moldovan IA, Borleanu F, Radulian M (2011) Source Parameters For The Earthquake Sequence Occurred In The Ramnicu Sarat Area (Romania) In November December 2007. Rom J Phys 56
Popescu E, Radulian M (2001) Source characteristics of the seismic sequences in the Eastern Carpathians foredeep region (Romania). Tectonophysics 338:325-337. doi: 10.1016/S0040-1951(01)00087-7.
Reasenberg PA, Jones LM (1989) Earthquake hazard after a mainshock in California. Science 243:1173-1176.
Ringdal F (1975) On the estimation of seismic detection thresholds. Bull Seismol Soc Am 65:1631–1642. doi: 10.1785/BSSA0650061631.
Scholz CH (2015) On the stress dependence of the earthquake b value. Geophys. Res. Lett. 42: 1399 – 1402. doi: 10.1002/2014GL062863.
Seghedi A, Vaida M, Iordan M, Verniers J (2005) Paleozoic evolution of the Romanian part of the Moesian Platform: an overview. Geologica Belgica, 8/4, 99–120.
Scholz CH (1968) The Frequency-Magnitude Relation of Micro-Fracturing in Rock and Its Relation to Earthquakes. Bull Seismol Soc Am 58:399-415. doi: 10.1785/BSSA0580010399.
Shimojo K, Enescu B, Yagi Y, Takeda T (2021) Nucleation process of the 2011 northern Nagano earthquake from nearby seismic observations. Sci Rep 11:8143. doi: 10.1038/s41598-021-86837-4.
Tormann T, Enescu B, Woessner J, Wiemer S (2015) Randomness of megathrust earthquakes implied by rapid stress recovery after the Japan earthquake. Nat Geosci 8:152–158. doi:10.1038/ngeo2343.
Utsu T (1965) A method for determining the value of “b” in a formula log n = a-bM showing the magnitude-frequency relation for earthquakes. Geophys Bull Hokkaido Univ 13:99–103.
Utsu T (1961) A statistical study on the occurrence of aftershocks. Geophys Mag 30:521-605. doi: 10.4401/ag-5094.
Utsu T, Ogata Y, Matsu'ura RS (1995) The centenary of the Omori formula for a decay law of aftershock activity. J Phys Earth 43:1–33. doi: 10.4294/jpe1952.43.1.
Wiemer S (2001) A Software Package to Analyze Seismicity: ZMAP. Seismol Res Lett 72: 373-382. doi: 10.1785/gssrl.72.3.373.
Wiemer S, Wyss M (1997) Mapping the Frequency-Magnitude Distribution in Asperities: An Improved Technique to Calculate Recurrence Times? J Geophys Res 102:15115-15128. doi: 10.1029/97JB00726.
Wiemer S, Wyss M (2000) Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the Western United States, and Japan. Bull Seismol Soc Am 90:859–869. doi: 10.1785/0119990114
Woessner J, Wiemer S (2005) Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. Bull Seismol Soc Am 95:684–698. doi: 10.1785/0120040007.

Aftershock analysis and forecasting for the crustal seismicity in Romania

Status:

Version 1

Abstract

Figures

Introduction

Earthquake data for the Marasesti and Gorj earthquake sequences

Methods

Results and Discussion

Parameter estimations and afterhocks’ number forecasting for the 2014 Vrancea-Marasesti earthquake sequence

Parameter estimations and afterhocks’ number forecasting for the 2023 Gorj earthquake sequence

Aftershock occurrence probability estimations for the two aftershock sequences

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1