Geodynamic Modeling in Central America Based on GNSS Time Series Analysis—Special Case: The Nicoya Earthquake (Costa Rica, 2012) ^†

Abstract

GNSS systems allow precise resolution of the geodetic positioning problem through advanced techniques of GNSS observation processing (PPP or relative positioning). Current instrumentation and communications capabilities allow obtaining geocentric and topocentric geodetic high frequency time series, whose analysis provides knowledge of the tectonic or volcanic geodynamic activity of a region. In this work, a GNSS time series study is carried out through the use and adaptation of R packets to determine their behavior, obtaining displacement velocities, noise levels, precursors in the time series, anomalous episodes and their temporal forecast. Statistical and analytical methods are studied; for example, ARMA, ARIMA models, least-squares methods, wavelet functions, Kalman techniques and CATS analysis. To obtain a geodynamic model of the Central American region, the horizontal and vertical velocities obtained by applying the above methods are taken, choosing the velocity with the least margin of error. Significant GNSS time series are obtained in geodynamically active regions (tectonic and/or volcanic).

1. Introduction

In the area of Central America, one of the most interesting geodynamic zones in the world is observed, with the convergence of five tectonic plates: the South American, North American, Nazca, Cocos, and Caribbean plates. The evolution of these plates is well known, except for the origin of the Caribbean plate, which is still the subject of debate [1,2,3].

There are two main models for the origin of the Caribbean plate. The first is one that contemplates the origin of the plate “in situ” between the North and South American plates [2]. The second and most commonly accepted model is that the Caribbean plate originated in the Pacific Ocean during the Late Cretaceous, composed primarily of a large igneous province due to hot spot magmatism and normal-thickness oceanic crust that migrated eastward to its present position between the American plates [1,4,5].

The Caribbean plate is an enclosed oceanic basin composed primarily of the Caribbean Greater Igneous Province and large areas of ocean crust of normal thickness within the Venezuelan and Colombian basins [4]. It is bounded to the north and south by two continental plates, the North and South American plates, which are characterized by large strike-slip fault systems, although convergence also occurs on a smaller scale [6]. It is bounded to the east and west by two main subduction zones, the Lesser Antilles and Central America, respectively, where the oceanic lithosphere of the Atlantic Ocean to the east and the Pacific Ocean to the west subducts beneath the Caribbean plate [7]. In addition, this plate subducts under northwestern South America, reaching a depth of 600 km [4].

At the western edge of the Caribbean plate lies the Nicoya Peninsula, where the Cocos Plate subducts northeastward beneath the Caribbean Plate along the Middle America Trench at about 8 cm/yr, with a range of obliquity of 10° counterclockwise from the normal direction of the trench [8,9]. The movement of the plates along the Middle American Trench is partitioned. In the normal direction of the trench, the subduction velocity is 74–84 mm/yr [10], while in the forward arc the movement is 8–14 mm/yr [8,11,12]. The subducting Cocos plate forms on both the rapidly extending East Pacific Ridge (EPR), with relatively uniform seafloor topography, and the slowly extending Cocos-Nazca Ridge (CNR), with a relatively rugged seabed [13,14].

The Nicoya Peninsula is elongated in a northwest–southeast direction and extends within 60 km of the Trench [15], and is located over the shallow portion of the subduction interface that generates earthquakes, called megathrust [16]. There are only a few subduction zones in the world where there is land access of as close as 50–60 km to a deep trench [17] and where the megathrust seismogenic zone is covered by land rather than ocean [16].

In this area, large megathrust earthquakes with magnitudes greater than 7 have historically occurred with a well-defined seismic cycle of about 50 years, 1853, 1900 and 1950 ($M_w=7.7$), 2012 ($M_w=7.6$) [14,16,18], in addition to the record of other smaller nearby events of magnitude ($M_w=7$) in 1978 and 1990 [14]. Slow slip events (SSEs) are common below the Nicoya Peninsula [10,19]. These events are basically largely aseismic slip that occurs at the plate boundary for weeks to months [20,21] releasing a fraction of the accumulated stress aseismically or weakly seismically, often accompanied by low-frequency earthquakes and non-volcanic seismic tremors [22].

On 5 September 2012, in the North Pacific region of the Nicoya peninsula, a $M_w=7.6$ earthquake was recorded SW of Sámara, which had a large number of aftershocks in the following months (42 in September to the SSW of the peninsula from Nicoya; 24 in October to the WSW of the Nicoya Peninsula; 10 in November on the coast of the Nicoya Peninsula). Most of these earthquakes have their hypocenters at depths between 15 and 20 km. This earthquake on 5 September and its aftershocks were caused by the subduction process of the Cocos plate under the Caribbean plate, a process that has generated other historic earthquakes in Guanacaste such as the 1950 Nicoya earthquake ($M_s=7.7$). There was damage to homes and buildings. The solution of the focal mechanism carried out by the USGS for the Samara earthquake shows a pure inverse type mechanism that confirms its relationship with the subduction process of the Cocos plate [23].

Ref. [17] relocated the earthquake of 5 September 2012 using data from the local seismic network. He located the hypocenter at $9.76$° N, $85.56$° W at 10 km offshore, 13 km deep at megathrust, with seismic moment being $3.5\times {10}^{20}$ Nm, giving $M_w=7.6$. The joint finite fault inversion of GPS data, seismometers, accelerometers, teleseismic P-waves, and GPS static offsets revealed that the coseismic rupture propagated downward from the hypocenter with a rupture velocity of $3.0$ km/s and a total source duration of 21 s.

2. Data Collection

GNSS technologies and permanent station networks have created a very relevant terrestrial reference framework and tool for the study of deformations of the Earth’s crust due to tectonic forces. These technologies are of great interest for geodynamics and deformation studies. Although strain is a more objective indicator than displacement because no frame of reference is required [24], GNSS techniques make it possible to quantify with guarantee the displacements of the stations that occur during earthquakes and relate them to other unaffected areas; as a consequence, the horizontal and vertical movements can be measured in faults and tectonically active regions. The GPS system has proven to be a very effective tool to carry out deformation studies due to its high precision and accuracy [25].

Nevada Geodetic Laboratory, (http://geodesy.unr.edu/, accessed on 2 March 2023) provides all available raw GPS geodetic data from over 17,000 stations around the world; these stations form the MAGNET network. Managing this large amount of data has led to the development of new processing strategies, automated systems, algorithms and robust estimation techniques. From these data, time series of Cartesian coordinates (every 24 h), topocentric coordinates (every 24 h and every 5 min), tropospheric delay and predictions are available. The NGL laboratory provides a multitude of graphs and their corresponding displacement velocities, and indicates with dashed lines the instant of time of the earthquakes that have occurred and the equipment change events near the station.

Plots of the unfiltered and filtered data, the results of fixing the series against an African plate, are provided. The NGL Lab also provides detrending series plots and predictions, filtered against the African and Eurasian Plate [26]. Every week the daily position coordinates of about 10,000 stations are updated. Every day, position coordinates are updated every 5 min for more than 5000 stations. Every hour, we update the 5-min position coordinates for about 2000 stations.

The NGL laboratory routinely updates station velocities, which can be used to image deformation rates of the Earth’s surface for a variety of interdisciplinary applications. These velocities are robustly estimated using the asymmetry-adjusted mean interannual difference (MIDAS software), a median-based GPS station velocity estimator that is insensitive to outliers, seasonality, step functions (abrupt changes) arising from earthquakes or equipment changes, and variability of statistical data [27]. Still, at NGL, for cases where an earthquake of magnitude greater than 6.9 has occurred, close enough to the station, we solve an exponential function defined by $A(1-exp(-(t-t_0)/\tau ))H(t-t_0)$ where $t_0$ is the earthquake time, $\tau$ is the relaxation time, A is the decay amplitude and H is the Heaviside step function. In these cases, we retrend after removing the exponential terms to obtain a self-consistent model for the time series.

Figure 1 shows the vector model developed from the velocities extracted from NGL, where the arrows mark the magnitude and direction in which each permanent GNSS station is moving per year.

3. Methodology

In this paper, the behavior of the time series of the stations near the earthquake that occurred in Nicoya on 5 September 2012 of magnitude 7.6 will be analyzed (see Figure 2).

The topocentric time series provided by the NGL laboratory was extracted. In this work, time series will be distinguished into three phases: preseismic, coseismic and postseismic. The preseismic phase is formed by the data prior to the moment of the earthquake; the coseismic phase comprises the data from the occurrence of the earthquake until the earth’s crust recovers; and the postseismic phase is characterized by being the data that goes from when the trend becomes linear again (recovers its preseismic behavior) until the end of the available data. Various analytical and statistical techniques have been applied to these data [28]. To obtain the velocities corresponding to the preseismic and postseismic phases, the CATS adjustment will be applied, which provides a model formed by the sine and cosine functions that adapt to the values of the time series; through the given model the displacement velocities are obtained from the series. The linear fit will be used to obtain the velocities of the coseismic series.

4. Analysis of the Series Affected by the Earthquake

The series affected by earthquakes have a common behavior between them (see Figure 3):

Following the methodology, the displacement speeds of each of the defined phases are obtained, thus giving three vector models of horizontal displacement. In this way it will be possible to see the magnitude and direction of the displacement of the GNSS-GPS permanent stations during the different phases. In addition, the stress–strain models will be obtained, thus seeing the zones of maximum geodesic deformation for each of the phases. For this, the “Q-Str2-Models” plugins available in QGIS will be implemented [29].

5. Conclusions

All the time series close to the Nicoya earthquake on 5 September 2012 present a behavior similar to that shown in Figure 3, in which three phases are distinguished: preseismic, coseismic and postseismic. The travel speeds were calculated for each station shown in Figure 2, in each of the phases. In this way, the graphs of Figure 4, Figure 5 and Figure 6 are obtained. In Figure 4a, the vector displacement model can be seen in the preseismic phase, seeing that the stations near the coast present a similar behavior. Figure 4b shows the maximum geodesic deformation that occurred in the preseismic phase, observing a greater deformation in the area where the SAJU, PNE2 and GRZA stations are located. In Figure 5a we have the vector model for the case of the coseismic phase; it can be seen how the direction and magnitude of the horizontal displacement changes. Figure 5b shows how the geodesic deformation increased compared to the values obtained in Figure 4b. However, said deformation is caused in the same zone by also adding the zones of the HATI, IND1, QSEC and BON2 stations. In Figure 6a you can see the moment in which the Earth’s crust again has a linear trend; even so, it does not present a behavior as homogeneous as that visible in Figure 4a. Figure 6b shows the maximum geodesic deformation of the region; it can be seen how the PNE2, SAJU, GRZA, IND1 and QSEC stations continue to form part of the areas with the greatest deformation; in addition, this deformation can be seen in a downtown station such as CDM1.