Relocation of Seismic Events and Validation of Moment Tensor Inversion for SENTINELS Local Seismic Network

Abstract

Seismic activity of Podhale (Poland) and Spiš (Slovakia) regions has been recognized for years. The first information about tremors from this area comes from the XVIII century. Four earthquakes, with intensity over VI in the MSK scale, were reported in 1643, 1724, 1840 and 1901. Since 1960’s an instrumental measurements have been conducted. However, until 2011 no stronger tremors have been recorded in the area of artificial water reservoir Czorsztyn lake located in extremely complex geotechnical condition, between border zone of Inner and Outer Carpathians separated by Pieniny Klippen Belt. Before Czorsztyn 2D seismic survey, knowledge of tectonic boundaries and velocities was very limited. The seismic survey confirmed assumptions of a flower type faults system and provided the first 3D velocity model for this area. After a series of earthquakes, in 2013 a SENTINELS network started its operation and since that time it recorded almost 200 events. The focal mechanisms were calculated for around 20 of them. Both location and moment tensor are crucial in the investigation of the origins of seismogenic process related to industrial operations. Therefore the relocation of the events and validation of the moment tensor solutions for the SENTINELS network were conducted with the use of a 3D velocity model. The validation of mechanisms was conducted with the use of the synthetic tests based on the 1D velocity model derived from the 3D velocity model, taking into account its lateral velocity anisotropy. It was based on the synthetic amplitudes generated with the assumed normal and strike-slip faulting type, similar to the obtained solutions. The validation proved that focal mechanisms are reliable even in a sparse focal coverage and noise not exceeding 40% of the P-wave amplitude. Most of the events are normal or strike-slip with nodal planes striking NW–SE or NNE-SSW. The latter ones are in agreement with the orientation of the main discontinuities in this area.

1 Introduction

Poland is located in a stable, continental part of the Eurasian plate. The seismicity in Poland exemplifies intraplate earthquake activity. The largest seismic activity in Poland is observed in the Polish Carpathians: Orawa–Nowy Targ basin, Podhale area and Pieniny Klippen Belt (PKB). The first information about earthquakes in this region dates back to the XVIII century. Other earthquakes were observed in 1840, 1901, 1935, and 1966. More recent seismic events occurred in the area of Krynica in 1992–1993 and in the Podhale region in 1995 and 2004 (Guterch 2009).

Czorsztyn Lake is situated at the eastern border of Orawa–Nowy Targ Basin in a border zone of two major Carpathian tectonic units: inner and outer Carpathians. PKB has complicated geological–tectonic structure; it is cut by a series of young faults mostly of N–S strikes. It is mostly built of carbonate rocks (limestone and marlstone), shale, and claystone. Its narrow width and the significant length suggest correlation with dislocation or a system of dislocations extending along the tectonic boundary of Western Carpathians (Birkenmajer 1960). Compared to negligible seismicity of Poland, from 1998 to 2010, around 50 earthquakes were recorded in a direct area of artificial Czorsztyn Lake (filled in 1997), which makes, on average, about one event in 3 months. Unexpectedly, a significant increase in seismic activity was observed in 2011. More than 60 seismic events were recorded close to Czorsztyn Lake in November 2011. Most of them were weak. The strongest earthquake, which occurred on 5 November 2011, had a magnitude M_W 2.6. After November 2011, the activity returned to its previous low level. However, since January 2013 a gradual increase in activity was noted. Within the period of time covered by this work (from August 2013 to December 2018) around 200 events with M_L from − 0.7 to 2.9 were recorded by the local seismic network. The catalog completeness magnitude was M_C = 0, calculated according to Goodness of fit test at 95% confidence bounds (Mignan and Woessner 2012). There were 155 events above completeness level in this period. Seismicity of the area is considered as a delayed seismic response type of triggered seismicity—the significant increase of seismic activity started 14 years after the lake’s impoundment. Czorsztyn Lake is shallow and the effect of the water load at depth is negligible. Temporal water level changes were not correlated with seismic activity (Fig. 1). The highest water level change was about 20 m during winter 2015/2016. Following refilling of the lake to its regular level the seismic activity slightly increased to the level similar to what it was before, but it did not reach the activity from the 2013/2014 highest rate. Additionally, both the groundwater level and precipitation cannot be linked with seismic activity changes. The seismicity is considered as a response to fault weakening due to diffusion effects. The earthquake triggering in this region would indicate the existence of a fault or fault system being close to failure (Białoń et al. 2015). Location and moment tensor are crucial in the investigation of the origins of seismogenic process related to industrial operations such as artificial reservoir exploitation. In Czorsztyn lake area the proper location of seismic events and focal mechanisms may indicate the active fault orientations.

Fig. 1

Location of seismic events (circles) and stations (black triangles SENTINELS network, grey triangle NIE broadband station. Lower panel: seismic activity and Czorsztyn lake water level from August 2013 and February 2018

The location of the hypocentres was determined with the use of LocSAT application based on ray tracing approach and simple 1D velocity model (Bratt and Bache 1988). The question is how accurate these locations are, do they fit local tectonic boundaries? The software using probabilistic Bayesian approach to the relocation of earthquakes from local seismic network in the area of Czorsztyn lake in Podhale region has been used to address these issues. The relocation was performed using a 3D velocity model of the studied area, which included boundaries between Pieniny Klippen Belt, Magura Flysh, and Inner Carpathians. This velocity model was based on the Czorsztyn 2D seismic survey. The main part of the work will focus on earthquake relocation with the use of TRMloc software (Debski and Klejment 2016). The selection and analysis of all necessary parameters for TRMloc will be presented. Relocation results may influence the MT inversion, therefore the reliability of MT solutions was also tested taking into account both relocation results and the velocity model influence on take-off angles as well as the noise in the data. SENTINELS (SEismic Network for TrIggered and Natural Earthquake Location and Source determination) network allowed for estimation of the moment tensor (MT) from P-wave amplitudes of the recorded events (Lizurek 2017). Focal mechanisms of 23 events were calculated. Most of the events are strike-slip or normal faults with the strike-slip component. Since the moment tensor is an important feature in discrimination between triggered and induced anthropogenic seismicity (Dahm et al. 2013) it is important to validate the reliability of the full MT in terms of the non-DC components as well as the nodal planes and principal axes orientation (Lizurek 2017). The validation of the moment tensor reliability was performed using synthetic tests. It was based on the synthetic amplitudes generated with assumed normal and strike-slip faulting similar to the obtained solutions. Take off angles and incidence angles were calculated with HypoDD ray tracer using 1D velocity model derived from the 3D velocity model (Kwiatek and Martinez-Garzon 2016). The validation proved that the focal mechanisms are stable with RMS error of amplitude fitting usually below 0.3 and dominance of DC MT component even in a sparse focal coverage and noise not exceeding 40% of the P-wave amplitude.

2 The Outline of Pieniny Klippen Belt Geological Structure

The Pieniny Klippen Belt in the northwestern Carpathians is a very complex geological structure. It is approximately 600 km long, and at the same time very narrow, from several 100 m to 20 km, stretching from Vienna in Austria to northern Romania. The PKB marks the suture zone between the Inner Carpathians (part of ALCAPA plate) and the Outer Carpathians, PKB separates Outer Carpathians composed of Upper Jurassic to the Lower Miocene deposits deformed during Miocene to the north and the Inner Carpathians composed of deformed in Cretaceous, Paleozoic and Mesozoic rocks covered by weakly deformed post-nappe deposits of the Central Carpathians Paleogene Basin to the south (e.g. Birkenmajer 1960; Plašienka et al. 1997). The PKB involves predominantly Jurassic, Cretaceous and Paleogene sediments with variable lithology and intricate internal structure. It has been subdivided into numerous lithostratigraphic and tectonic units (e.g. Plašienka 2012). The boundaries of PKB display tectonic character, represented by sub-vertical faults and shear zones. The PKB originated 20–14 million years ago as a flower structure limited by two parallel faults merging into a single vertical fault, which currently reaches the North European Platform. The Inner Carpathians are thrust over the North European Plate (Golonka et al. 2018).

3 SENTINELS Network Observations

Instrumental seismicity measurements in Orawa-Nowy Targ Basin began in the 1960s when Institute of Geophysics, Polish Academy of Sciences (IG PAS), had installed Su-69 seismometer in Niedzica castle basement. When Niedzica dam constructions works had started in the late 80 s, ZEW Niedzica S.A. (the state-owned company managing the facility), built a new bunker for the seismometer, which was located on the hill, 2 km to the west from Niedzica castle, to provide a proper signal to noise ratio for seismic measurements. Digital measurement started parallel to analog one in the 90 s, and both are carried on until now. The station in Niedzica (NIE) is part of Polish Seismological Network system PLSN. In 2008–2012, the project Monitoring of Seismic Hazard of Territory of Poland (MSHTP) was launched by the IG PAS. The Western Carpathians region was monitored by five stations. Podhale seismic monitoring continued after the end of the MSHTP project. Due to increased seismicity rate near Czorsztyn Lake observed since 2011, IG PAS in cooperation with Niedzica Hydro Power Plants Company (ZEW Niedzica SA) launched local seismic network SENTINELS (SEismic Network for TrIggered and Natural Earthquake Location and Source determination) in August 2013. The first configuration of network contained seven short period Lennartz 3D Lite seismometers with NDL recorders. At the turn of 2013, under the IS-EPOS project (IS-EPOS: Digital Research Space of Induced Seismicity for EPOS Purposes), the network has been reorganized. Lennartz equipment was replaced by GeoSig short period (5 s) VE-53 seismometers with GMSplus3 recorders. One station, due to noise, was relocated, and three new stations were founded. Since that time network consists of 10 short period stations, located around Czorsztyn Lake in facilities of ZEW Niedzica SA, Krościenko administration of forestland and Nowy Targ municipality, and in private land. The short period stations are supported by broadband station NIE from PLSN. SENTINELS network was part of IS-EPOS project (see “Acknowledgments”). SENTINELS seismic network has recorded almost 200 seismic events since August 2013. Figure 1 shows the locations of stations and seismic events determined by LocSat software. The main cluster of events was located below Czorsztyn lake. For all recorded events a local magnitude M_L was determined, for stronger events also moment magnitude M_W was calculated. The strongest—Mw 3.1, earthquake recorded by SENTINELS network so far took place on 1st April 2014 and was located in the main cluster. Magnitude completeness level for the catalog is M_L = 0. Since SENTINELS began operating, two series of increased seismicity periods could be distinguished. The first occurred between August to September 2013, with 35 earthquakes (Fig. 1). That temporal cluster cumulated on 5 September 2013, when 11 events occurred with the strongest one of M_L 1.5. The second period of increased activity occurred between February and April 2014, and contained 43 seismic events. In the most active day (23.03.2014) 10 events occurred, with the strongest event of M_L 2.1. The strongest event of the whole activity was M_L 2.9 1st April 2014 (Fig. 2). For the rest of the time, seismic activity rate recorded by SENTINELS network was low and stable, with just few earthquakes per month. The data set used for relocation contains 195 seismic events from 25.08.2013 to 15.11.2017.

Fig. 2

Magnitude release in time and activity rate presented as a number of events per month since SENTINELS network operation began

4 Geophysical Data and Velocity Model Derived from the Czorsztyn 2D Seismic Survey

Podhale, Pieniny and PKB, due to complex tectonics, have been an area of many seismological research. The geodynamic research, after building a dam in Czorsztyn, was investigated by Pachuta et al. (2010) and Walo et al. (2016). They focused on leveling, gravimetric and, geodetic observations. The study yields 7 mm vertical displacement and changes in gravity (0.1 mGal over 30 years) after filling the reservoir. The horizontal observations show less than 1 mm/year linear trend in the northeast direction. A similar study was made by Łój et al. (2009), the research covered geodynamic process monitoring, but the two observation lines were skipping Czorsztyn lake area. Data mentioned above, the small volume of the Czorsztyn Lake, as well as small water level changes during the normal exploitation of the reservoir, are all supporting the thesis that the effect of the water load at depth is negligible. There are no deep boreholes in the study area of Czorsztyn lake. The nearest one, in Bańska Niżna, is placed almost 20 km from the area and was drilled for the geothermal purpose. The deep seismic sounding, made during CELEBRATION 2000 project under the leadership of Janik et al. (2009), shows only the general seismic velocities and tectonic structures. There were also no commercial 2D or 3D seismic profiles in this area. The first high-resolution deep seismic reflection survey Czorsztyn 2D was carried out in the framework of the IS-EPOS project. The aim of Czorsztyn 2D seismic survey was to define geology, tectonics and rocks properties in the area of Czorsztyn lake, especially area of Pieniny Klippen Belt, and Inner and Outer Carpathian. Four seismic lines (01-01-15 K, 02-01-15 K, 03-01-15 K, and 04-01-15 K) of the total length of 50 km, were designed to provide the seismic imaging as deep as possible. The recording was provided in all four lines simultaneously, which gave 2004 active channels. The sampling rate was 2 ms, and record length was 8 s. The distance between receiving point was 25 m, at each point were 12 geophones in the circle patch set. The shot interval was 50 m. Seismic sweep signal was generated with the length of 16 s, and bandwidth 6–80 Hz. For most of the time, four vibroseis machines were running. Low Velocity Zones (LVZ) were measured, in 15 places located in nearby of seismic line with the use of 200 m refraction survey. The field work took place in February 2015, lasted about 2 weeks and was performed by Geofizyka Kraków SA. Four archival seismic lines from 1987 (24-5-87 K, 24A-5-87 K, 26-5-87 K, and 28-5-87 K) were reprocessed to provide a good tie to the nearest borehole. In the area of seismic survey, data from following boreholes were available: Bańska IG-1, Bańska PG-1, Bukowina Tatrzańska PGP-1, Biały Dunajec PAN-1, Maruszyna IG-1, Nowy Targ IG-1. Unfortunately, most of them have only lithostratigraphic profiles, and only two of them: Maruszyna IG-1 and Nowy Targ IG-1 has seismo-acoustic profiles. Czorsztyn 2D seismic survey provided the first high-resolution velocity model of the Pieniny Klippen Belt. The velocity model was developed with the use of seismic inversion technique. Seismic inversion is a reverse modeling, where the properties of rocks and tectonic structures are recreated, based on seismic processed data from the surface survey (after migration) and wellbore data. Czorsztyn 2D was processed with the use of Model-Based Inversion. The model was the base of the synthetic seismic section, then the iteratively synthetic section was fitted to the measured data, an adjustment factor was higher than 0.98 and the fitting error was less than 0.05. For 171 super bins (their central positions called further CMP are represented as the violet dots at Fig. 3), in eight processed profiles, interval velocities were estimated, which was the base for average velocity estimation and spatial velocity estimation, binned into Cube 3D. Figure 3 shows an average velocity cross section in a map view at the depth of 3000 m below sea level, based on Cube3D (Dec 2015).

Fig. 3

A cross-section of velocity at the depth of 3000 m below sea level. The 171 violet points represent super bins positions (CMP), where the mean velocity field was estimated (Dec 2015)

Czorsztyn 2D seismic survey allowed to create 1D and 3D velocities models for SENTINELS purpose. The points from reprocessed profiles have information about velocities to the maximum depth of about 10 km below sea level. The points from newly made profiles, due to different designing, have information of depth to about 22 km, from sea level. Unfortunately, Czorsztyn 2D did not cover all area of SENTINELS network, due to the necessity of tying it to the nearest wellbore located 15 km from the research area, and limited funding for that purpose in the IS-EPOS project. Nevertheless, Czorsztyn 2D covered the western and central part of Czorsztyn lake. The line 04-01-15 k runs almost directly through the area of seismic events cluster, located with the use of LocSat. The velocity curves start from about 700 m above sea level and end at the depth of about 21.5 km below sea level. Velocity field could be separated into two main parts. The first one, from the surface to the depth of about 8 km, has a bigger velocity range, but the increase of velocity is less rapid than in the second part. Figure 4 shows all CMP values plotted in one graph, the black bold line represents a gradient average model used in the first relocation attempt with LocSat. A velocity of S wave was assumed as ${1 \mathord{\left/ {\vphantom {1 {\sqrt 3 }}} \right. \kern-0pt} {\sqrt 3 }}$ of P wave velocity.

Fig. 4

A cloud of CMP curves are marked as thin lines, the black bold line shows a gradient 1D velocity model for SENTINELS purpose

Based on the Czorsztyn 2D a 1D velocity model was made as a gradient model, best fitting to average velocity value of CMP points. The 3D velocity model was constructed with all available CMP data, plus 15 points of Low Velocity Zone (LVZ) refraction profiles. As mentioned previously, velocity data do not cover all research area. The span of SENTINELS network, limited by the stations and events locations, is higher than 50% of available data. In Fig. 5, a cross-section of calculated velocity model with the location of CMP, LVZ points, earthquakes, and stations are shown.

Fig. 5

A Cross section of 3D velocity model prepared for use in TRMloc software, at the depth of 3 km below sea level, covering seismic network area. The CMP points are marked as red dots, violet points represent LVZ, the SENTINELS seismic stations are marked as green triangles, green circles are representing seismic events calculated with the use of LocSat software. The contour mark the Czorsztyn lake, solid black line denotes the profile location used in the interpretation of relocation results

TRMloc application requires using 3D velocity model. The range of the model is approximately 30 × 20 km in the X–Y plane and 20 km in Z plane. It is expanded by default 3% in each plane to prevent a computational error near the edges. The algorithm to calculate the model uses extrapolation to fulfill missing data, it also smooths model to prevent fast velocity changes in small areas. The grid density can be chosen by the user. The grid is automatically adjusted to the shape of the square. TRMloc uses a Cartesian coordinate system, in this case, all coordinates were converted to XY2000 coordinate system.

5 Methods

Location method with probabilistic (Bayesian) inversion technique was applied for Czorsztyn seismic data. According to TRMloc software’s authors (Debski and Klejment 2016), used in this work, the software approach relies on exploring the space of all possible source locations and assigning to each point in this space (e.g. each possible location) a probability of being the true hypocentre location. This a posteriori probability can then be used for any analysis of location errors using standard statistical methods. To obtain uncertainty the two meta-characteristics are used: the evidence (Z) and entropy (H). The likelihood function is actually the convolution of statistics of observational and modeling errors (Debski 2015a). If we assume that these joint errors can be modeled by the Gaussian statistics so the choice of the quadratic misfit function is justified and we only need to set the covariance matrix defining misfit function. If we set Cp (variance of the distribution) significantly larger than the variance of the real errors we end up with the likelihood function which is apparently consistent with the statistic of real errors (somehow encapsulates it). However, if Cp is taken too small, the likelihood function will fail to describe the real errors becoming ‘too narrow’ compared to statistics of the real errors. Thus, a dependence of some properties of the a posteriori solution on Cp could provide us information when our assumption about the Gaussian nature of the joint errors breaks down. Cp parameters influence mainly shape of the a posteriori probability distribution rather than its particular feature—the position of its maximum, then H, Z defined by Eq. (5) in Dębski (2015b) are the simplest meta-characteristics of the a posteriori PDF.

The probabilistic approach to hypocentre location provides the most complete information about the solution, however, the method is computationally demanding as it requires an exhaustive exploration of the 3 or 4-dimensional model space (3 or 2 spatial dimensions and time for origin time calculation, respectively). TRMloc uses the time reversal and reciprocity principles of the wave equation together with finite-difference type forward modeling solvers. The Fast Sweeping Method (FSM) eikonal solver, applied in TRMloc enables very efficient calculation of the wavefront positions in the entire 3D domain for a general velocity model. Nevertheless, the algorithm has some limitations, the eikonal solver provides solutions only for the first arriving seismic phases (or the first arriving P or S waves in case of elastic waves). Another limitation is related to the fact that the spatial resolution achieved by the algorithm is limited by the grid size used by the forward modeling algorithm. Achieving higher resolution requires a finer spatial grid but this increases the computation time linearly. The next issue is connected to the spatial grid in eikonal solver used by TRMloc, the first-order differential solver requiring quite fine grid for high numerical accuracy. Analysis with TRMloc should be run when velocity models are relatively smooth, which means that the lateral and vertical velocity ratio between the grid cells should be around 0.2. In case of very high-velocity contrast, the eikonal solver may become numerically unstable. TRMloc returns two types of final solutions: Rml and Rav. The Rml is a maximum likelihood solution when Rav is an average solution. The importance of the average model Rav comes from the fact that it provides not only information on the best fitting model but also includes information about other plausible models from the neighborhood of the “best” model Rml. Comparison between them provides a qualitative evaluation of the reliability of the inversion procedure: the more average solution differs from maximum likelihood one, the more complex and non-bell-shaped is the form of the a posteriori pdf. In such a case more care must be taken when interpreting the inversion results, especially inversion uncertainties.

Focal mechanisms were estimated with the use of moment tensor inversion (MT). Moment tensor was calculated using inversion of the P wave amplitudes in the time domain (Wiejacz 1992; Kwiatek and Martinez-Garzon 2016). The synthetic tests of the noise influence and take off angle uncertainties on the results of full MT solutions were carried out to investigate the reliability of the MT results. Synthetic amplitudes and polarities were generated for a priori known fault orientations. The mechanism for the synthetic data was assumed as pure shearing (99.99% of the DC full MT component). The synthetic input data were amplitude and polarities of P-wave displacement following the proper focal sphere quadrants to the configuration with the amplitudes following the P and T axes orientations for the assumed fault. Bootstrap amplitude resampling tests of noise and take off angle uncertainty influence were conducted (Kwiatek and Martinez-Garzon 2016). In amplitude resampling procedure, the random noise was introduced to input amplitude data. The amount of noise added to the original amplitude was specified as a random number drawn from Gaussian distribution with mean 0 and standard deviation of 0.1, 0.2 and 0.3. The noise to input amplitude data reached a factor of 20, 40, and 60% respectively. In a similar manner, the take-off variability was set up to 15° due to lateral anisotropy of the 3D velocity model and limitations of 1D velocity model used in routine MT inversion. A bootstrap procedure was similar in the take-off angle influence test. Take-off angle was drawn randomly from the range of 15° around the initial take-off angle value for every station and then the MT was calculated including original amplitude input data, such procedure was repeated 100 times.

6 Relocation of the Events

Location is a basic procedure in the seismic data processing. There are few ways to determine the event location. The fastest and most common is the ray tracing method used for example in LocSat (Bratt and Bache 1988). This software is designed to use the 1D velocity model in global seismic observation, which may be inappropriate for local scale networks like SENTINELS or mine networks. Also, it is hard to estimate errors made during the location. The second type of software is based on the probabilistic approach. This type of software is designed to use a 3D velocity model, in local and global applications, returning full calculations errors, for further analysis. The disadvantages of this software include the calculation time, strongly depending on the 3D calculations grid. Examples of this type of software are TRMloc (Debski and Klejment 2016) and NonLinLoc written by Lomax et al. (2000). The third type of software is based on the double differences method, but it was not considered in this paper.

Determining the location and focal mechanisms of the events is part of a routine operations. Before Czorsztyn 2D seismic survey there was no reliable velocity model for the Czorsztyn lake area. Location of the event was calculated using LocSAT program of SeisComP3 with IASPEI-91 velocity model, which in fact is the isotropic model for the shallow subsurface. The data set used for relocation contains 195 seismic events from 25.08.2013 to 15.11.2017, however, after taking into account the completeness level of the catalog, 155 events left. Locations determined in this way have some errors because of a general regional velocity model applied for the complex local geological situation. After Czorsztyn 2D seismic survey the new 1D average model from 3D (Fig. 4) was used in LocSat program to determine events location and in FOCI to calculate moment tensors (Kwiatek and Martinez-Garzon 2016). All events locations were recalculated. To compare the obtained solutions the residuum normalized by a number of stations (normalized RMS) was calculated following Aster et al. (2013). There was not much difference in epicentre location of the main cluster, recalculated events were slightly moved to the NE (Fig. 6), with average location change around 2 km. Less than 20% of the events above completeness level had the location change higher than 5 km with small fraction of events with location change between 5 and 8 km, and with several events outlying extreme changes over 10 km (Fig. 6). The biggest differences were in depth the average depth change was around 5 km. Moreover, normalized RMS did not differ significantly (Fig. 6).

Fig. 6

Relocation results on a map view (green dots denote relocated events, grey dots denote primary location). Histograms of normalized RMS errors for primary locations obtained with IASPEI velocity model (upper plot) and histogram of normalized RMS errors for relocation (lower plot) with use of average 1D velocity model derived from a 3D velocity model

Next, the 3D velocity model was used in TRMloc, which is a software developed to a spatial location with probabilistic, also called Bayesian, inversion technique. The TRMloc application requires few input files: a priori solution, wave origin time, velocity model and control parameters. In this paper, a location from seismic catalog determined by LocSat program with 1D velocity model was used as an a priori solution. To recalculate events from SENTINELS network P waves have been chosen, because these waves are better recognizable in the seismograms. Information about P wave arrival times were taken from the local bulletin of SENTINELS network. The 3D model was used in calculations, it covered all SENTINELS network area with a 3% margin. The range of 3D model includes 23 × 21 × 19 km and has a size of 229 × 207 × 171 points. The X–Y, and also Z step was equal to 100 m, so it fulfills the necessary condition of the proportional model grid. These grid parameters were chosen after testing various grid sizes with different steps, the chosen parameters were the best trade-off between accuracy and time of the estimation. The control parameters file contains information about all necessary input and output parameters, file names, data formats, etc. Altogether there are 30 parameters that can be defined. Most important are: ‘Cp’, ‘Cme’ and ‘Cmz’. The detailed information about all parameters, input/output data can be found in TRMloc Reference manual by Dębski (2015a). The Cp parameter declares likelihood variance. The diagonal elements of the a posteriori covariance matrix Cp are convenient estimators of the inversion uncertainties. The default value of Cp parameter is equal to 0.1 s. In the calculation, the range of Cp parameter was changing in the range of 0.001 to 1 s. As it was shown by Dębski (2015b) the wide range of Cp is necessary to find the optimum Cp value for the located event. If too narrow Cp range is used, it might be impossible to find out optimum Cp and then it has to be set arbitrary. The parameter Cme is a postulated uncertainty for the horizontal coordinates X and Y of the a priori model, by analogy Cmz is a postulated uncertainty for the depth coordinate Z. The value of these parameters is given in meters. A large value of Cme and Cmz parameters are causing a minimum influence of the a priori location to the inversion results. Analogously a small value of Cme and Cmz strongly influence inversion results limiting the space of possible locations.

The next step was to determine the residual for each station. TRMloc gives back residual value for each station as Rml (a maximum likelihood solution), and Rav (average solution). To achieve that goal the first iteration of TRMloc was made. At that moment only residual values for the stations were important. Velocity model was set as a 0D model with the P wave velocity equal to 5300 m/s, the parameters Cme, and Cmz were set to 20000 m and 5000 m respectively. There were 19 seismometers available, located at 12 points, including POTK station from MSHTP network and NIE station from PLSN network (Table 1). If the residual value was bigger than 0.2, the station was rejected from further processing. That situation may occur for few reasons, for example, the recorded signal was very weak and the noise was picked as P wave onset, or there were problems with velocity model, or timing for the station was incorrect. In general, only in few cases, a station was rejected due to the residual value. An example of residual values for a station located in Kacwin is shown in Fig. 7. There were two stations operating in this place. The station KACW recorded 124 events and was picked up during network upgrading, the station KAC2 recorded 19 events. Figure 7 shows the residuals of Rml, Rav, and both plotted on one picture for station KACW.

Table 1 SENTINELS stations and other available stations used in the relocation, stations with similar codes were used when located in the same place but the Lenartz equipment was replaced with GeoSig

Fig. 7

Comparison between Rml and Rav residual values for a station located in Kacwin

The Rml residuals of KACW station were in the range of − 0.01 to 0.02, only few Rml values were close to 0.1. The Rav residuals were in the range between − 0.1 to 0.2 with a maximum at 0.05. Two values were rejected from the future calculation due to an exceedance of 0.2 value. The plot of both residuals shows that the Rml is gathering over the 0 points, and is much more narrow than the average residual distribution. The Rav values are shifted to 0.1 and the spectrum of values is wider. Other stations show similar characteristics as described above. The Rml residuals are mostly concentrated around 0 when Rav is shifted from the 0 about ± 0.05 to 0.1. In the case of stations in Harklowa, Maniowy, Ochotnica, Niedzica, and Potok both residuals are distributed very similarly. The more data was recorded by station the more residual distribution tends to follow a normal distribution. The next step in the relocation procedure was to apply the 3D model. In this case, the 2D pdf distribution of events was taken into account in the validation of the results, checking if it has a bimodal distribution. The Cp parameter was set to default 0.1 value, the Cme and Cmz parameter was set to 20000 m and 5000 m respectively, the same as in the previous step. There were two main cases: first with events that occurred inside the network, and second with events that occurred near the edge or outside of the network. The minimum number of the station to gain location is four, but as much as possible is preferred. When the event occurred inside the network and was recorded by numerous stations, a point solution may be expected. The probability distribution takes Gaussian bell-shaped, and epicenter location is very accurate and the depth location uncertainty is small (Fig. 8). This case covers only a part of events with higher magnitudes. If the event was weaker and recorded by few stations, the solution is less accurate, especially in case of depth. An example probability distribution of magnitude M_L 1.4 event recorded by 9 stations is shown in Fig. 8, which could be treated as an example of a well-located event. A case of the weaker magnitude M_L = 0.1 event recorded inside the network by five stations is plotted in Fig. 9. It can be clearly seen that the weaker event has a much wider probability distribution of the epicenter location. The probability distribution of depth location is spindle-shaped with poorly resolved depth location represented by broad maximum likelihood area. In such situation the determined depth location will contain large uncertainties and the final depth will be set as the average value rather than the maximum likelihood one.

Fig. 8

An example of the 2D probability distribution of relocated event M1.4 occurred on 11.07.2015. The earthquake was recorded by nine stations (denoted with green triangles, red triangles denote unavailable stations) and was localized in the center of the network. The picture in the left shows X–Y plane, two others show X–Z and Y–Z plane

Fig. 9

An example of the 2D probability distribution of relocated event M0.1 occurred on 16.10.2013. The earthquake was recorded by five stations (denoted with green triangles, red triangles denote unavailable stations) and was located in the center of the network. The picture in the left shows X–Y plane, two others show X–Z and Y–Z plane

In the second case, events occurred outside of the network, or at the edge of it, and was recorded by numerous stations. This situation belongs to the minority of SENTINELS network. An example of such case is plotted in Fig. 10. The event occurred on 10th January 2017, M_L 0.2, and was recorded by 10 stations. The biggest probability distribution in the X–Y plane is quite narrow, but the distribution in Z plane is elongated. In such situation, the accuracy of Z coordinate will contain higher uncertainties, therefore this event could be treated as an example of poor location. In Fig. 11 the results of the relocation of an earthquake recorded by only four stations are shown (the event occurred on 10th February 2014, with M_L = − 0.1). It is an extremely small event below the completeness level of the catalog and it is resulting in poor location results. The probability density function of the epicenter location is widely dispersed and fan-shaped, in a direction opposite to seismic stations. Two clusters of higher probability can be separated, they are stretched to the north and northwest, the accuracy of epicenter coordinates will have a bigger error than the previous ones. Such location is not reliable, however, because that event is below the completeness level of the data set it is excluded from the final data. The probability density function of depth is widely spread. The precision of the determined location is smaller when events occurred outside or near the edge of the network. If the magnitude of an earthquake is small it is more likely that it will be recorded by fewer stations and the location result will be less reliable. Classical approximation of location errors by the error ellipse (covariance matrix) brakes down in case of small events. We excluded from the final analysis events with a magnitude below completeness level and these, which had too big residual values of Rml and Rav, which resulted with a dataset covering 148 events from 25.08.2013 to 15.11.2017 with magnitude M_L ≥ 0. In the next step the values of three main TRMloc parameters: Cp characterizing origin time variance in seconds, Cme and Cmz were specified. They characterize a priori the variance for epicenter (X–Y) and depth (Z) coordinates, quantifying a priori expected location uncertainties. Influence of the two latter parameters in the final result is much smaller than Cp, in such a case, their values were set only once and remained unchanged during calculations, respectively 10000 m for Cme and 5000 m for Cmz. The relatively high value of these parameters provides a very weak a priori constraints on the final solutions.

Fig. 10

An example of the 2D probability distribution of relocated event M0.2 occurred on 10.01.2017. The earthquake was recorded by ten stations (denoted with green triangles, red triangles denote unavailable stations) and was localized outside of the network. The picture in the left shows X–Y plane, two others show X–Z and Y–Z plane

Fig. 11

An example of the 2D probability distribution of relocated event M-0.1 occurred on 10.02.2014. The earthquake was recorded by four stations and was localized outside of the network. The picture in the left shows X–Y plane, two others show X–Z and Y–Z plane

The likelihood variance Cp, determining location errors has the biggest influence on the final result. If Cp was set significantly larger than the variance of real errors, we end up with the likelihood function which is apparently consistent with the statistic of real errors. If Cp is taken too small the likelihood function will fail to describe the real errors, becoming too narrow with respect to the statistics of real errors (Debski 2015a). To find the value of these parameters, TRMloc was ran into a loop with Cp set in a range between 0.001 to 1 with asynchronous step given by the formula, ${{Cp\left( i \right)} \mathord{\left/ {\vphantom {{Cp\left( i \right)} {Cp\left( i \right) + 1}}} \right. \kern-0pt} {Cp\left( i \right) + 1}}$ equal to 1.2. The first value of that series are given by: 1, 0.8, 0.64, 0.512, …, 0.001. All together it gives back 31 datasets to analyse, holding 195 events each. Analysis of inversion parameters, evidence (Z) and entropy (H), provide dependence of location results on Cp, and estimation of the final value of Cp parameter. Evidence (Z) and entropy (H) are the simplest meta-characteristics of a posteriori distribution, easily calculated within a probabilistic approach. A place of flattening of function H = H(Cp) for larger Cp determine the maximum value of Cp below which, the a priori setting becomes irrelevant. Original linear decrease breaks at some point below which H(Cp) decrease faster than predicted Gaussian distribution, it is also applied for the Z = Z(Cp) curve. The breaking point is a point of estimated Cp value. To estimate the optimal value of Cp a derivate $H' = {{dH} \mathord{\left/ {\vphantom {{dH} {dlog\left( Z \right)}}} \right. \kern-0pt} {dlog\left( Z \right)}}$ was calculated. If the well pronounced maximum is visible, it will be taken as optimum Cp value. Functions of evidence Z = Z(Cp) and entropy H = H(Cp) were plotted in Fig. 12. A plot of evidence Z shows that there is no single point where Z function is flattening. First evidence Z values drop with the Cp equal to 0.08. Most of if drop when Cp is equal to about 0.01, but some of it did not drop till the end of the scale at point 0.001. A plot of entropy H show a high dispersion of value, it is difficult to determine the point of Cp breakdown, nevertheless, it could be estimated at the point of Cp equal to 0.1. Plots of calculated derivative H’ for each cases are shown in Fig. 13. It is impossible to determine the global maximum of these functions, that will be an optimal Cp value due to high dispersion. Half of derivative H’ value starts to growth after reaching Cp equal to 0.1, part of them constantly decrease, and a few of them show a local maximum. Analysis of individual H’ function shows that, if H’ has a local maximum it is in the Cp range between 0.01 and 0.1. Analysis of evidence Z, entropy H and derivative H’ shows that in the case of SENTINELS data it is very hard to find optimal Cp value, because plots of entropy H and derivative H’ were overlapping, optimal Cp value was estimated on the base of evidence Z plot, but events with constant H’ were skipped following approach described by Debski (2015a). However, it may be more save to consider every case separately, which is more convenient when the individual event is processed, but less handy when many events are processed at once, therefore we can assume the calculated Cp as an optimum value for this data set is 0.068. This value corresponds to Cp series taken into calculation.

Fig. 12

A dependence of evidence Z (left plot), and entropy H (right plot) for all calculated cases. Cp values are on a logarithmic scale, as log(Z) values

Fig. 13

A plot of derivative H′ = d(H)/d(log(Z)). Cp on a logarithmic scale

An overview of relocation results obtained with TRMLoc was plotted in Figs. 14, 15 and 16. The first one shows a map view of epicenter locations obtained with routine LocSat procedure using the IASPEI velocity model and TRMloc average solutions obtained for local 3D velocity model. LocSat location results were highly concentrated in a single cluster in the network center, some of the events were gathering in East–West line in the south part of the research area. Average TRMloc solutions of the main cluster are located more to the NE with average epicenter location shift around 2 km and 90% of the events location change was below 5 km. Only few events located at the network edges had bigger location change (Fig. 14, upper histogram insert). The depth of the relocated events was mainly deeper, average depth change was around 1.5 km and most of the events depth change was not bigger than 3 km (Fig. 14, lower histogram insert). For maximum likelihood solutions of TRMloc, the cluster is slightly moved to N (Fig. 15), but events were more concentrated within the main cluster than in case of average TRMloc solutions. Epicenter location and depth changes between primary LocSat location and TRMloc maximum likelihood were smaller than in case of an average TRMloc solution. Similarly to the LocSat relocation solutions obtained with local gradient velocity model the average solutions of TRMloc main cluster was moved further to NE (Figs. 6 and 14). However, LocSat relocated event solutions were deeper than the primary locations—average depth change was around 5 km (Fig. 6, lower histogram). The TRMloc maximum likelihood depth solutions were at a similar depth as primary LocSat depth locations. Main cluster events were more dispersed between 7 and 10 km depth in comparison to primary 8–10 km, while average TRMloc depth locations were deeper than both previous solutions with the main cluster located between 9 and 11 km depth. The difference in location was observed in both TRMloc average and maximum likelihood solutions consistently. Analysis of the uncertainty of relocation results for each direction shows that a mean error value for X-axis direction was about 2 km, in Y-axis about 1.5 km and about 2.7 in depth. The values seem to be big, but we have to keep in mind that events located in the center of the network were about 20% of all events. Most of the events have small magnitudes and were located only by a few stations, sometimes outside the network. When we look at the uncertainties for relatively strong, recorded by 9 stations and well-located event inside the network (see Fig. 8), the epicentre location uncertainty is less than 300 m, while depth location uncertainty is less than 2 km. On the other hand the smaller events located inside the network by fewer stations have epicenter location uncertainty around 2 km and depth location uncertainty is larger than 5 km (Fig. 9). However, it is better to have more stations available in terms of the epicenter and depth location accuracy, even when the event is located on near the edge of the network. It can be exemplified with the event of M_L 0.2 shown in Fig. 10, its epicenter location uncertainty is about 800 m, depth location uncertainty is around 4 km. The magnitude of the event obviously plays a role, because it is much easier to properly pick P wave arrival for bigger magnitude event. The average normalized RMS errors are smaller for both TRMloc solutions, than LocSat primary locations (see RMS histograms in Figs. 14 and 15).

Fig. 14

Relocation results of average solutions obtained with TRMloc on a map view (white dots denote relocated events, grey dots denote primary location). Histograms of normalized RMS errors for primary locations obtained with IASPEI velocity model (upper plot) and histogram of normalized RMS errors for average solutions of relocation (lower plot) with use of 3D velocity model. Main calculation parameters were set to Cme = 10000 m, Cmz = 5000 m, and Cp = 0.0687 s

Fig. 15

Relocation results of maximum likelihood solutions obtained with TRMloc on a map view (brown dots denote relocated events, grey dots denote primary location). Histograms of normalized RMS errors for primary locations obtained with IASPEI velocity model (upper plot) and histogram of normalized RMS errors for maximum likelihood TRMloc location solutions with depth change after relocation (lower plot) with use of 3D velocity model. Main calculation parameters were set to Cme = 10000 m, Cmz = 5000 m, and Cp = 0.0687 s

Fig. 16

Relocation results of maximum likelihood solutions (brown dots) and average (white dots) obtained with TRMloc on a map view. Histograms of normalized RMS errors for maximum likelihood solutions (upper plot) and average solutions (lower plot) obtained with TRMloc and 3D velocity model. Main calculation parameters were set to Cme = 10000 m, Cmz = 5000 m, and Cp = 0.0687 s

7 Focal Mechanism Reliability Validation

Focal mechanisms were estimated using the moment tensor inversion method. MT was calculated using inversion of the P wave amplitudes in the time domain (Wiejacz 1992; Kwiatek and Martinez-Garzon 2016). Uncertainties of the estimated moment tensors can be estimated through the normalized for the number of stations root-mean-square (RMS) error between theoretical and estimated amplitudes (Stierle et al. 2014a, b). There are 23 moment tensor solutions calculated for the events of 0.2 ≤ M_L ≤ 2.1, located in the center of the network and recorded by at least eight stations (IS-EPOS, 2017 and Leptokaropoulos et al. 2018). Therefore there was no significant change in location after the relocating procedure. All solutions were recalculated after the relocation with the use of the average 3D velocity model. There were no significant differences in the nodal plane and MT components in the final solutions. Most of the events are strike-slip or normal faults with the strike-slip component. The MT decomposition results (Fig. 17) confirmed the triggering origin of the seismicity in the area of Czorsztyn lake. Shearing component is the most dominant part of the mechanisms, but the non-DC parts of MT are observed—most of the events have at least 20% of non-shearing components (Lizurek 2017).

Fig. 17

Histogram of DC components of full MT solutions obtained for Czorsztyn lake related seismic events

The synthetic tests of the noise influence and take off angle uncertainties on the results of full MT solutions were carried out to investigate the reliability of the MT results. In all tests set up of 11 stations was used. Synthetic amplitudes and polarities were generated for a priori known fault orientations which were typical for the Czorsztyn lake area: normal (strike/dip/rake: 33°/58°/− 104°) and strike-slip (188°/85°/124°). The mechanism for the synthetic data was assumed as pure shearing (99.99% of the DC full MT component). The synthetic input data were amplitude and polarities of P-wave displacement following the proper focal sphere quadrants to the configuration with the amplitudes following the P and T axes orientations for the assumed fault. The RMS error of all MT solutions for such prepared synthetics was smaller than 0.01. Take off angle uncertainty was set up to 15° due to lateral anisotropy of the 3D velocity model and limitations of 1D velocity model used in routine MT inversion in FOCI. For all cases bootstrap amplitude resampling tests of noise and take off angle uncertainty influence were conducted (Kwiatek and Martinez-Garzon 2016). In amplitude resampling procedure, the random noise was introduced to input amplitude data. The amount of noise added to the original amplitude was specified as a random number drawn from Gaussian distribution with mean 0 and standard deviation of 0.1, 0.2 and 0.3. The noise to input amplitude data reached a factor of 20, 40, and 60% respectively.

Results of the noise influence on normal fault solutions of MT for SENTINELS network are shown in Fig. 18. The bigger the noise, the bigger the fraction of solutions of full MT which tend to deviate from the assumed geometry. However, over 90% of the solutions are similar to the assumed geometry. The noise increase is also followed by RMS error increase up to 0.7 value in case of the highest noise contamination. The number of solutions with high RMS error (> 0.4) is larger than 20 out of 100 solutions only in case of the highest noise contamination. The spurious non-DC components are also increasing with increased noise, but only up to 25% of the solutions in case of the highest noise contamination is characterized by the dominance of false non-DC components.

Fig. 18

Results of the bootstrap noise influence tests of normal fault solutions. From the left noise contamination up to 20, 40, and 60%

Similar results were obtained for strike-slip geometry and take off angles resampling (Figs. 19 and 20). In the case of this faulting geometry set up the influence of noise is bigger. More spurious non-DC components and more outlying nodal plane orientations are observed (mainly in 60% noise contamination). Nevertheless, DC components are dominant in the majority of cases—around 70% of the obtained solutions. The influence of take-off angles uncertainties is negligible in both assumed fault geometries. The spurious non-DC components are not larger than 20% in 99% of the cases and the RMS values are less than 0.2 with 0.1 being mean RMS value. It proves the SENTINELS network is reliable for the MT inversion purposes and confirms the reliability of full MT solution obtained for Czorsztyn lake related seismic events. It confirms also shearing as the main component of the mechanisms obtained in this area.

Fig. 19

Results of the bootstrap noise influence tests of strike-slip fault solutions. From the left noise contamination up to 20, 40, and 60%

Fig. 20

Results of the takeoff angle resampling test for the strike-slip mechanism. Clockwise from the left upper corner: Hudson plot of full MT decomposition, DC component histogram, RMS error, a beach ball with all solutions

8 Conclusions

SENTINELS network is located in a very complex geo-tectonic setting. Until 2015 there had been very few information about seismic velocity and tectonic structures were unconfirmed. The velocity model obtained during IS-EPOS project did not cover all area of SENTINELS network. The missing data were extrapolated from existing information and velocity model was smoothed to fulfill TRMloc relocation software demands.

TRMloc requires to be run several times to calculate stations residuals, analyse the 3D probabilistic distribution and to estimate final parameters. Station residuals help to show whether P waves were picked correctly or if there is systematic timing error on the seismic station. If earthquakes were recorded only by few stations or occurred outside seismic network location, results contain bigger errors, or the software has a problem to estimate one final value like in case shown in Fig. 11. Because SENTINELS network is a surface network, location errors in depth are greater than in XY plane, especially when a priori main cluster of solutions was estimated at the depth of 8 km, what is close to half network span. Average normalized RMS errors were smaller for both TRMloc solutions than for LocSat primary locations (Figs. 14, 15), while maximum likelihood solutions RMS errors were smaller than the average solution RMS errors. Average solutions obtained with TRMloc are deeper than the ones from routine LocSat estimations, while maximum likelihood solutions are at a similar depth as primary routine solutions. Most of the events with a relatively large change of the location after using TRMloc are the ones with small magnitudes close to completeness level and originally located outside of the network or close to its edges. Above mentioned differences in locations are probably caused by using a 3D velocity model by TRMloc and 1D IASPEI91 model by LocSat in primary solutions. However average TRMloc epicenter location solutions are similar as the LocSat solutions obtained with an average 1D model derived from the 3D local model, the main difference is in depth locations, which in case of the latter method are about 3 km deeper than the average TRMloc solutions. The 3D model is much more complex and provides more information about velocities. This can be seen in Fig. 4. The latter 1D gradient model is only an estimation of the complex CMP velocity distribution. The variance of velocity for the depth of 5 km is between 4600 to 5600 m/s, when the gradient model has estimated velocity equal to 5200 m/s. A similar situation happened when we compared results from LocSat when the top layer of the velocity model describing the first 20 km was a single value. Replacing this value by gradient model moved depth location of the events deeper.

Comparing final depth from TRMloc plotted in Fig. 21, with the seismic cross-section with tectonics, it shows that locations were correlated with tectonics only by the elongated shape of the main cluster, in both TRMloc solutions event clouds are following PKB_N and further north discontinuities. However, they are located more to N and deeper than these discontinuities in Fore Magura–Silesian nappes formation. The most of the events from the main cluster of maximum likelihood solutions are stretched in depth from 7 to 10 km, and are located to the north of PKB Northern fault, in Fore Magura–Silesian nappes formation, between Grybów Unit base and North European Platform base. Average locations of events from the main seismic cluster are deeper, between the border of the base in Fore Magura–Silesia nappes formation with some of them located inside North European Platform. Both TRMloc solutions are more scattered than the original LocSat locations, which show the 3D velocity model influences the location solutions. Lateral variability of 3D model results with more scattered seismic event locations, while 1D initial model LocSat location tends to cluster the events. Reliability of the full MT solutions obtained after relocation was tested with synthetic tests designed for two main faulting geometries: normal and strike-slip mechanisms. Tests took into account discrepancies of the velocity model and noise in the data. Results of the tests prove the SENTINELS network is reliable for the MT inversion purposes. The full MT solutions obtained for Czorsztyn lake related seismic events are stable and the RMS errors are relatively low even in case of high noise contamination. Focal mechanisms are mainly strike-slip or normal with some strike-slip component with the nodal planes striking mostly NNE–SSW or NW–SE. The shearing was the main component of the mechanisms obtained in this area. The strike-slip events solutions are in agreement with the main discontinuities orientations parallel to PKB northern (PKB_N) fault, while normal faulting mechanism orientation is hard to associate with the main discontinuities. The latter ones may be some secondary discontinuities not recognized in the seismic sounding results interpretation, but following the direction of the boundary between North European Platform and Fore Magura–Silesia nappes formation. Strike-slip events and the main discontinuities have similar orientations, but the events are not located exactly on the discontinuity planes (Fig. 21). This discrepancy might be explained in two ways. First, events are located on smaller discontinuities parallel to the main ones, which are too mature to be reactivated, while the smaller ones are triggered. A second explanation might be connected with location uncertainties which on average are around 1.5 km in epicenter location and 2.7 km in depth. Taking this into account we cannot rule out that the main cluster of seismicity is not connected with the main discontinuities plotted in Fig. 21.

Fig. 21

TRMloc relocated solutions (green dots denote maximum likelihood locations, red dots denote average locations) and primary LocSat locations (black dots) at the background of the seismic cross-section (location of the section is shown in Fig. 5) and interpreted discontinuities (based on Dec 2015). The white missing part of the profile denotes the lake location, which caused the gap in the coverage of receivers

The rise of seismic activity from 2013 to the beginning of 2014 was an exception from the normal seismic activity state observed in this area since its impoundment. However it wasn’t the first episode of increased seismic activity observed in Czorsztyn lake vicinity, the first rise of seismic activity was observed in November 2011 and the next in spring 2013. SENTINELS network started in August 2013, therefore the detailed studies and location of the events could be conducted for data recorded after the biggest activity was reported. None of these periods of increased seismic activity were associated with water level changes in the lake (Białoń et al. 2015). The small water volume in the reservoir, long gap between the lake impoundment and increase of seismicity, and no clear evidence on direct influence of the water level changes on seismicity in the PKB were shown. Therefore seismicity of Czorsztyn reservoir should be considered as delayed response seismicity type triggered on existing discontinuities following Simpson et al. (1988) classification of reservoir-induced seismicity.