Multiscale Post-Seismic Deformation Based on cGNSS Time Series Following the 2015 Lefkas (W. Greece) Mw6.5 Earthquake

Abstract

In the present work, a multiscale post-seismic relaxation mechanism, based on the existence of a distribution in relaxation time, is presented. Assuming an Arrhenius dependence of the relaxation time with uniform distributed activation energy in a mesoscopic scale, a generic logarithmic-type relaxation in a macroscopic scale results. The model was applied in the case of the strong 2015 Lefkas Mw6.5 (W. Greece) earthquake, where continuous GNSS (cGNSS) time series were recorded in a station located in the near vicinity of the epicentral area. The application of the present approach to the Lefkas event fits the observed displacements implied by a distribution of relaxation times in the range τ_min ≈ 3.5 days to τ_max ≈ 350 days.

1. Introduction

Earthquake rupture creates static stress changes within the lithosphere that are large enough to produce transient deformation that can be observed on the Earth’s surface. Satellite geodesy using GNSS techniques for observing Earth’s surface displacements have advanced in such a level that the post-seismic deformation can be monitored with high accuracy [1].

To study post-seismic relaxation, various mechanisms have been proposed, such as those of poroelasticity of the rock medium [2,3], shallow after-slip on the main faulting plane or nearby faults [4,5], and viscoelastic relaxation [6]. The after-slip mechanisms are localized and follow logarithmic decay [4,7], while the viscoelastic one is supported by a bulk deformation process and exhibits exponential decay [8,9,10].

In [11], it was demonstrated that, during the first 6 months following the 1999 Chi-Chi earthquake (M = 7.6, Taiwan), both the post-seismic relaxation and the cumulative number of aftershocks followed the same temporal dependence, which reduced to a logarithmic one after a long time from the main event. In addition, logarithmic deformation could be explained by an exhaustion process, which is based on the assumption that the number of sources of strain-producing events evolve according to a thermally activated process, in which the activation energy is stress dependent [12].

The recent strong earthquakes present examples where the macroscopic post-seismic deformation is interpreted in view of the above relaxation mechanisms, which are based on the assumption of a single relaxation time estimated for each post seismic relaxation observed [13,14,15,16,17,18].

Experimental results and theoretical models suggest that post seismic relaxation is thermally activated and enhanced by the presence of fluids [19]. Deformation is controlled by a variety of temperature-dependent healing processes that are related to the crack size [20,21]. Since different crack sizes have different activation energy, the associated relaxation processes operate at different timescales [22], indicating an intrinsic organization [23] with a distribution of relaxation times, instead of the use of a single one. The latter implies that post-seismic relaxation could be viewed as the superposition of several different relaxation times, each one activated in different mesoscopic domains of the macroscopic relaxed seismic volume. Thus, the macroscopic relaxation observed results from the superposition of different relaxation times on the microscopic–mesoscopic scale, which obviously operates on different spatiotemporal scales. The recovery process that returns the deformed volume to equilibrium after a strong earthquake could be characterized as a slow dynamics process [24], because an instantaneous perturbation causes a transient response, which moves the system to a state of equilibrium [25,26,27,28,29,30]. It is, thus, possible that relaxation, which varies with the logarithm of time, can be interpreted as the result of the superposition of relaxation mechanisms with a spectrum of relaxation times.

The purpose of this work is to give a description of the post seismic deformation relaxation. It is based on the thermodynamic Arrhenius model for the distribution of the relaxation times [31] and describes the post-seismic recovery as a slow thermodynamic process taking place in a mesoscopic scale. The model predicts a physically accepted finite value just after the co-seismic perturbation, a logarithmic time dependence for the time period between the minimum and maximum relaxation time, and an exponential decay while approaching an equilibrium state a long time after the main event. Based on that, the explanation that relaxation occurs as a logarithm of time does not depend on the details of Earth’s physical and chemical properties, but is based on a generic, or even a universal, multiscale property of the relaxation mechanism, as resulted from the Arrhenius expression, which connects relaxation time with the activation energy of the relaxation process.

To test the proposed approach for the post-seismic relaxation, positioning data from the continuous GNSS station PONT [32] and www.geodesy.gein.noa.gr (accessed on 30 December 2020) located in the southern part of Lefkas Island were used. This cGNSS station is located in the close vicinity of the earthquake epicenter of the 17 November 2015 Mw6.5, shallow event (~10 km depth), that occurred on the island, causing severe damages and human losses.

2. A Thermodynamic Model of Post Seismic Relaxation

A post-seismic transient deformation stage, where strain varies with time at a decreasing rate, occurs following the instantaneous deformation during an earthquake. Transient deformation is followed by the return to steady-state deformation, where the strain rate is a constant after some time following the main event.

Following [33], we define a relaxation function R(t), where t is the time after the main seismic event, as a perturbation of the displacement given by the generic expression:

$$u\left(t\right)=u_{\infty }+S R\left(t\right)$$

(1)

where $u_{\infty }$ is the equilibrium value of displacement after a long time and S is a scaling factor. Notice that $R\left(t\to \infty \right)\approx 0,$ while $u_{\infty }=u\left(t\to \infty \right)$ is the remaining distortion at the end of the relaxation process. The initial value $u_o$ of the post-seismic deformation immediately after the co-seismic perturbation δu_c is given from Equation (1) as $u_o=u\left(t=0^{+}\right)=u_{\infty }+S R_O$, where $R_0=R\left(t=0\right)$. The scaling factor S of the relaxation constrains the rheological properties of the medium [34,35] and depends on the rheology of the crust and upper mantle, as well as on the local conditions around the site (i.e., crustal heterogeneities). In addition, S depends on the stress released in the medium and reflects the response of the medium (i.e., of the magnitude of the earthquake event) after the co-seismic displacement due to earthquake.

In most of the cases a simple relaxation mechanism with a single relaxation time τ that corresponds to a relaxation function exp (−t/τ) (i.e., not varying logarithmically with time) is the simplest rheological model introduced to explain post-seismic deformation, known as the linear Maxwell model [6]. In this case:

$$R_o=1,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{u}_o=u_{\infty }+S{ }_{}\mathrm{and} u\left(t\right)=u_{\infty }+S e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{$\tau $}\right.}$$

(2)

while the normalized displacement is $\mathrm{D}\left(t\right)=\frac{u\left(t\right)}{u_o}=\frac{1+\frac{S}{u_{\infty }}{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{$\tau $}\right.} }{1+\frac{S}{u_{\infty }}}$. In Figure 1, we plot the displacement as a function of normalized time $\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{$\tau $}\right.$. In the very early stage, expanding the exponential in Equation (2) leads us to $u\left(t\right)=u_{\infty }+S\left(1-\frac{t}{\tau }\right)\mathrm{or}\mathrm{D}\left(t\right)=1-\frac{\frac{S}{u_{\infty }} }{1+\frac{S}{u_{\infty }}} \frac{t}{\tau }$, a linear expression of normalized displacement with time, which for $\frac{S}{u_{\infty }}\gg 1,$ leads to D = 1 − (t/τ) with the slope equal to 1/τ, and independent of S.

We proceeded to study a post-seismic relaxation process with a spectrum of relaxation times within the range defined by the minimum and maximum values ${\tau }_{min}$, ${\tau }_{max}$ of relaxation time, respectively. If f(τ) is the distribution function of the relaxation time, then:

$$R\left(t\right)={{\displaystyle \int }}_{{\tau }_{min}}^{{\tau }_{max}}f\left(\tau \right)e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{$\tau $}\right.}d\tau$$

(3)

The expression (3) is based on the assumption that the total macroscopic relaxation observed is a result of the superposition of relaxation processes operating with different relaxation times τ. This approach could be viewed as a generalization of the Burgers body rheological model [8,36,37], expanding the two discrete relaxation times to a continuum spectrum in the range values ${\tau }_{min}$, ${\tau }_{max}.$

Assuming a relaxation process with activation energy E, then following Arrhenius’ law [31], the relaxation time is given by:

$$\tau ={\tau }_0{e}^{\raisebox{1ex}{$E$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}$$

(4)

with ${\tau }_0$ a constant proexponential factor, T the temperature, and k_B the Boltzmann constant. If the activation energy is distributed following a distribution function N(E), then the distribution function f(τ) for the relaxation times satisfies the expression f(τ) =N(E)(dE/dτ). Using Expression (4), dτ/dE= (${\tau }_0$/k_BT)exp (E/k_BT) = τ/k_BT, hence:

$$f\left(\tau \right)=\frac{k_{B}TN\left(E\right)}{\tau }$$

(5)

In the case of a uniform distribution of the activation energy, bounded by a minimum E_min and maximum E_max activation energy, N(E) has a constant value N_o, which, according to Equation (5) for a temperature T, leads to a distribution of relaxation times $f\left(\tau \right)=(k_{B}T{{\rm N}}_0)/\tau ,$ which is proportional to 1/τ. It is evident that the minimum ${\tau }_{min}$ and maximum ${\tau }_{max}$ relaxation times are related to the corresponding minimum and maximum activation energies by the expressions ${\tau }_{min}={\tau }_0{e}^{\raisebox{1ex}{$E_{min}$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}$ and ${\tau }_{max}={\tau }_0{e}^{\raisebox{1ex}{$E_{max}$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}$ and, thus, $\Delta {\rm E}=E_{max}-E_{min}=kT \ln (\frac{{\tau }_{max}}{{\tau }_{min}}$).

The macroscopically observed relaxation is given by the superposition of relaxation processes, expressed as:

$$R\left(t\right)=C{{\displaystyle \int }}_{{\tau }_{min}}^{{\tau }_{max}}\frac{1}{\tau }{e}^{-t/\tau }d\tau$$

(6)

where C = $T N_o$ is a constant for a given T and N_o, which is omitted in the calculation merging it with the scaling parameter S. We note that N_o is a constant value for a uniform distribution of energy N(E) between its minimum E_min and maximum E_max value. The weight factor 1/τ in Equation (6) is a result of the application of Arrhenius’ law due to mesoscale relaxation processes operating on different timescales. In Figure 2, we present examples of the relaxation function with time on a logarithmic scale along the horizontal axis with a constant ${\tau }_{min}$ and various values of ${\tau }_{max}$. In a logarithmic scale, the curve in Figure 2 decreases almost linearly, up to a time t ≈ τ_max, and at the end, it flattens out. Thus, the relaxation function R(t) is proportional to ln (t) for t << τ_max. Based on that, the maximum relaxation time τ_max could be estimated as the time where the relaxation curve flattens out.

The integral in Equation (6) is the exponential integral E₁(x) [38], which can not be solved analytically but behaves like a negative exponential for very large values of the argument and decays on a logarithmic mode in the intermediate range of the relaxation curve R(t) while deviating from a logarithmic time dependence for t < τ_min. Furthermore, as results from Equation (6) when t = 0, R(t) can be calculated analytically:

$$R_0=R\left(t=0\right)=\ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)$$

(7)

supporting the physical condition to have a finite relaxation function for t = 0. Since $\Delta {\rm E}=E_{max}-E_{min}=kT\ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)=kT{R}_0$, the finite value of the relaxation function for t = 0 is required in order for the energy range $\Delta {\rm E}$ to be a finite quantity. Taking the derivative of R(t) with respect to t, we form the following equation:

$$\frac{dR}{dt}=\frac{1}{t}\left(e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{min}$}\right.}-e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{max}$}\right.}\right)$$

(8)

which, for τ_min$\approx$ t≪τ_max, expands exponentially and leads to ${\left[\frac{dR\left(t\right)}{dt}\right]}_{t=0}=-\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right)$ and

$$R\left(t\right)\approx \ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)-\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right)t=R_0-\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right)t$$

(9)

In the observed cases where τ_max ≫ τ_min, the slope of R(t) vs. t at t = 0 is well approximated as ${\left[\frac{dR\left(t\right)}{dt}\right]}_{t=0}=-\frac{1}{{\tau }_{min}}$, while the ratio of the slope to R(0) follows the expression: $\frac{1}{R\left(t=0\right)}{\left[\frac{dR\left(t\right)}{dt}\right]}_{t=0}=\frac{1}{{\tau }_{min}} \ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)=\frac{R_o}{{\tau }_{min}}$.

Introducing expression (9) in (1), the very early stage displacement immediately following the co-seismic one is given by the expression:

$$u\left(t\right)=u_{\infty }+S{R}_0-S\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right)t=u_o-F_{1}t$$

(10)

where the slope of Equation (10) is $F_1=S\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right)\approx S/{\tau }_{min}$ for ${\tau }_{min}\ll {\tau }_{max}$ and $u_o=u_{\infty }+S{R}_0$.

In the case when τ_min ≪ τ_max and the time t is such that τ_min ≪ t ≪ τ_max, we could approximate exponentials such as $e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{min}$}\right.}\approx 0 \mathrm{and} e^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{max}$}\right.}\approx 1$ and Equation (8) results to:

dR/dt ≈ −1/t

which, after integration, gives:

R(t) ≈ R₁ − ln(t) for τ_min≪τ≪τ_max

(11)

with R₁ an integration constant. Equation (12) reveals that the relaxation function R(t) exhibits a logarithmic behavior in an intermediate time interval bounded by the relaxation times ${\tau }_{min}$ and ${\tau }_{max}$ with a slope −1 in a graph of R(t) versus ln (t). It becomes evident that when multiplying R(t) with a scale factor S in Equation (1), there is a modification of the slope of the postseismic relaxation versus ln (t) related to the particular characteristics of earthquake event and the seismotectonic settings and the rheology of the deformed zone as well. Substituting Equation (11) into (1), displacement obtains logarithmic behavior:

$$u\left(t\right)=u_{\infty }+S{R}_1-S\ln t$$

(12)

for τ_min ≪ τ ≪ τ_max with a slope F₂ = S.

Defining τ* as the point where the transition from the linear to logarithmic behavior starts, there is a continuation of displacement (as given by Equations (10) and (12), respectively), that permits the estimation of R₁ as:

$$R_1=\ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)+\ln {\tau }^{ *}-\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right){\tau }^{ *}$$

(13)

In order to estimate the order of magnitude of R₁, we calculate its value for two reasonable values of τ*, as that of ${\tau }^{*}={\tau }_{min}$ and ${\tau }^{*}=2{\tau }_{min}$. For ${\tau }^{*}={\tau }_{min},$(i.e., the transition in the logarithmic range starts in the very early stage of deformation), we have $R_1=\ln {\tau }_{max}-1+\frac{{\tau }_{min}}{{\tau }_{max}},$ which for $\frac{{\tau }_{min}}{{\tau }_{max}} \ll 1,$ leads to:

$$R_1=\ln {\tau }_{max}-1$$

giving $R=-\ln \left(\frac{t}{{\tau }_{max}}\right)-1$. Since R(t) is positively defined, the above expression is valid for t < τ_max/e ≈ 0.37${\tau }_{max}$. For a more reasonable estimation, we set ${\tau }^{*}=2{\tau }_{min}$, leading to $R_1=\ln {\tau }_{max}+\ln 2-2(1-\frac{{\tau }_{min}}{{\tau }_{max}}),$ which for $\frac{{\tau }_{min}}{{\tau }_{max}}\ll 1,$ leads to:

$$R_1=\ln {\tau }_{max}+\ln 2-2,$$

with $R=-\ln \left(\frac{t}{{\tau }_{max}}\right)+\ln 2-2.$ From the positive definition of R(t), we conclude that the above expression is valid for $t<\frac{{\tau }_{max}}{2-\ln 2} \approx 0.76{\tau }_{max}$. A comparison between Equations (10) and (13) indicates that the slope of the logarithmic part $F_2=S$ is scaled with that of early stage deformation behavior $F_1=S\left(\frac{1}{{\tau }_{min}}-\frac{1}{{\tau }_{max}} \right).$ For τ_min ≪ τ_max, this scaling can provide an estimation of τ_min, since $\frac{F_2}{F_1}={\tau }_{min}$.

Analyzing the relaxation function (3) for the case τ_min ≪ τ_max ≪ t, the long time behavior of R(t) can be estimated, replacing the upper limit τ_max with infinity so that [33]

$$R\left(t\right)\approx \frac{{\tau }_{max}}{t}{e}^{-\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{max}$}\right.}$$

(14)

indicating that, for t ≫ τ_min, the relaxation function is controlled only by τ_max and not τ_min. In that case, the relaxation time τ_max is predominant, and R(t) leads to $R\left(t\right)=e^{\raisebox{1ex}{$(-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{max}$}\right.)},$ implying that the term τ_max/t is associated to the relaxation processes with relaxation time τ_min < τ < τ_max. We note that, since the time dependence of 1/t is slow compared to $e^{\raisebox{1ex}{$(-t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{max}$}\right.)},$ the relaxation function could be well approximated as an exponential a long time after the main event. Furthermore, the exponential relaxation, as presented in Equation (2), is a special case of the general expression of the relaxation function (3) and results when τ_min and τ_max have similar values, i.e., when τ_min = τ₀ − Δ and τ_max= τ₀ + Δ, which leads to $R\left(t\right)\approx \frac{2\Delta }{{\tau }_o}{e}^{-\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{${\tau }_0$}\right.}$ presenting an exponential decay over time.

3. Seismotectonic Setting, Data Selection, and Analysis

Intense seismic activity in the Central Ionian Islands of Cephalonia and Lefkas during the last few years, along with several large events (M > 6) over the last century [39,40], clearly highlight the key role of this area in the kinematic process in the Eastern Mediterranean. The area is dominated by a major regional seismogenic zone known as the Cephalonia Transform Fault (CTF) (Figure 3), as well as several other faulting features, resulting in generation of the highest seismic energy release in Europe [41,42]. The broader area is characterized by a complex kinematic and geodynamic regime. The area links the oceanic subduction boundary in the south to the collision between the Apulian microplate and Eurasia in the north [43,44,45]. This complicated tectonic regime creates a stress field that generates intense seismic activity and ground deformation on both a local and regional scale.

Based on seismological data and geodetic measurements, the CTF zone is described as a dextral strike-slip faulting type system. This system is composed of two separate segments that run along the western shores of the Cephalonia and Lefkas Islands. The southern segment, offshore of Cephalonia, extends for approximately 90 km on a SW–NE direction and ends on the northern part of the island. The northern segment, offshore of western Lefkas Island, strikes a NE–SW direction for ~50 km, ending close to north-western Cephalonia. This segment dips to the ESE and exhibits dextral strike-slip motion combined with a small thrust component [45].

Over the last twenty years, several strong earthquakes have devastated the islands of Lefkas and Cephalonia. In August 2003, a strong event (M_w6.3) occurred off-shore of the north-western coast of Lefkas Island [46,47,48,49,50]. In late January–early February 2014, two large-magnitude earthquakes (M_w5.9 and M_w6.1) were located in the western part of Cephalonia [51,52]. However, these events were not linked directly with the CTF zone.

This study focuses on a strong seismic event that took place in near proximity to the PONT continuous GNSS station. On 17 November 2015, a strong earthquake (M_w6.4) occurred in the southern part of Lefkas Island [53,54] just north of Cephalonia. The main event was followed by a significant post-seismic activity that extended southwards towards Cephalonia. The hypocenter of this earthquake was located at a depth of ~10 km, and the fault plane solutions demonstrated a NNE–SSW striking dextral strike-slip seismic fault with a reverse component that dipped east at a high angle of about 70° ± 5°. The after-shock seismic activity extended along three distinctive clusters north and south of the main event [54]. The northern cluster seemed to extends towards the epicentral area of the 2003 strong earthquake, the central one was located in the vicinity of the main event, while the southern one extended towards Cephalonia Island and about 15 km southwest of the epicenter. The spatial coverage of the aftershocks indicated less dense distribution close to the epicentral area, while away from the main shock, there was higher concentration of post-seismic activity. The latter may indicate that, during the main event, higher energy was released in the main ruptured fault, which comprised part of the major CTF zone, while the post-seismic activity was triggered in nearby secondary and adjacent local faults, away from the main rupture plain [54]. The post seismic events extended in a NE–SW direction for approximately 50 km, extending north towards the 2003 earthquake and south towards the northern Cephalonia. The southern seismic cluster occurred in the same area that was also activated during all the previous intense seismic periods in 2003 and in 2014 [48,49].

3.1. GNSS Data Analysis

The epicentral area of the 2015 Lefkas earthquake, located to the southern part of the island, is remarkably close to the local continuous GNSS station PONT [32]. This station has been operating in the area since 2007. The GNSS data are processed daily, using 30 s RINEX files, in order to calculate the daily coordinates of the station (Figure 4). From the time series of PONT station coordinates, the station velocity was estimated for the period 2009–2015. Prior to the 2015 strong event, the station exhibited the anticipated regional motion of the area, and the velocity was calculated as V_East = (20.45 ± 0.03) mm/yr, V_North = (6.81 ± 0.02) mm/yr, and V_Up = (−0.71 ± 0.08) mm/yr with respect to IGb08 reference frame (http://igscb.jpl.nasa.gov/network/refframe.html) (accessed on 30 December 2020). The earthquake of 17 November 2015 caused a co-seismic displacement of the PONT station of 360 mm to south, 201 mm to the west, and a subsidence of 73 mm (Figure 4).

The files from PONT station were processed using the Bernese Software version 5.2 [55]. The software allows the estimation of epoch wise receiver coordinates in precise point positioning mode (PPP), as well as in the double difference mode. For the present research, coordinates were calculated after the earthquake of 2015 on a static mode. To get the best results, the PPP technique was used as a first step to get a priori coordinate values, which were, consequently, introduced in the more precise double difference method. In the processing procedure, GNSS data from 25 IGS and EUREF cGNSS stations from Europe, Greece, Cyprus, and Turkey were used to define the local reference frame. Ambiguities were solved using several resolution strategies depending on the baseline length between the stations: the wide-lane and narrow-lane for the medium baseline (<10–200 km) and the quasi-ionosphere-free (QIF) strategy for long baselines (<1000–2000 km). For the troposphere modelling, the dry and wet Vienna mapping function (VMF) of Bohm et al. (2006), the dry and wet Niell model [56], as well as the Saastamoinen model [57] were the models used to analyze the GNSS microwave data in different stages of the processing. Effects of the solid earth tides were taken into account and corrections of the ocean tide loading were applied based on the FESS2004 model. In analogy to site displacements caused by ocean tide loading, the loading effect of the atmospheric tides were taken into account in the processing of the GNSS data. The latter was done by introducing a specific file in the processing containing atmospheric tidal loading coefficients from a global grid-based model [58]. Moreover, the absolute antenna phase center corrections were used in the processing and to better improve the final solutions files containing the precise orbital parameters; the Earth orientation parameters and the satellite clock corrections were taken from the Center for Orbit Determination in Europe (CODE) (DOI:10.7892/boris.75876.4). On a weekly basis, the calculated coordinates were evaluated for the repeatability error, and in case of large deviations from the weekly solution, the values were excluded. The latter resulted in an unweighted RMS error, in mm, for each component that varied from 0.5 to 1.01 in the east component, 0.8 to 2.1 in the north one, and 1.8 to 4.6 for the up component. However, that was not applied for the first ten days after the main event. The daily coordinates were estimated with respect to IGb08 reference frame.

3.2. Post Seismic Time Series Analysis

The post seismic displacement evolution is made up of two contributions: the displacement u(t) due to the relaxation one on the fault zone and the long term V_pst t, where V_pst is the long term loading velocity a long time after the earthquake when the system has reached a new equilibrium state. In our case, V_pstEast = (17.48 ± 0.07) mm/yr, V_pstNorth = (3.72 ± 0.07) mm/yr, and V_pstUp = (−0.44 ± 0.18) mm/yr, which was removed from the observed postseismic data. To apply our approach, we analyzed the post seismic relaxation observed after the 17 November 2015 Mw6.5, Lefkas island earthquake. In Figure 5, the North–South displacement, which was the event’s predominant post seismic displacement, presented as a function of time for a period 1000 days after the event. For energy reasons, the minimum relaxation time τ_min cannot be equal to zero in order to have a converged relaxation function. From Figure 5, the transition from the linear to logarithmic part is observed, permitting the estimation of τ_min ≈ 3–5 days. Since the longest relaxation time τ_max is the time where the relaxation stops, from Figure 5, we estimate τ_max ≈ 350 days. The aforementioned estimations of τ_min and τ_max result in a value R_o ≈ 4.2 to 4.7. Using the observed values of $u\left(t\right),$ we estimate that an average value of displacement in the range from 300 to 600 days is about −26 mm, a value that could be used to approximate $u_{\infty }.$ In Figure 5, the observed displacement, along with the theoretical one, calculated using Equation (1), with S = 4.8 mm, $u_{\infty }=-26 \mathrm{mm}$, τ_min ≈ 3.5 days, and τ_max ≈ 350 days (i.e., R_o ≈ 4.6), is presented as a very good agreement. Furthermore, the time evolution of displacement suggests a slope F₂ ≈ −4.21 for the logarithmic part and an approximate value F₁ ≈ −1.09 for the short linear early part of deformation. Comparing the above estimation with Equations (10) and (13), we conclude that, since τ_min << τ_max, the scaling of the slopes leads to τ_min ≈ 3.8 days, in agreement with our estimation from Figure 5.

A function to present the slow dynamics of postseismic relaxation could be defined in terms of the asymptotic long-term value $u_{\infty }$ as [59]:

$$\Omega \left(t\right)=\frac{u\left(t\right)-u_{\infty }}{u_o-u_{\infty }}$$

(15)

Substituting Equation (10) in (15) leads us to:

$$\Omega \left(t\right)=\frac{R\left(t\right)}{R_o}$$

(16)

In Figure 6, the experimentally estimated Ω(t) function for the first 400 days of the post-seismic relaxation after the Lefkas main event is given, along with the calculated one using Equation (16) for τ_min ≈ 3.5 days and τ_max ≈ 350 days, respectively. A comparison of the observed and calculated values is given in Figure 7, along with a dichotomous line.

Furthermore, we note that ${\int }_0^{\infty }\Omega (t)dt=\frac{{\tau }_{max}-{\tau }_{min}}{R_o},$ which, after using the experimental values of Ω, leads to ${\tau }_{max} \approx 360$ days in quite good agreement with the value used in Figure 5 and Figure 6. We note that the law values of $\frac{R\left(t\right)}{R_o}$ in Figure 7 correspond to the logarithmic part of relaxation, while the high values are related with the part just after the rupture where the uncertainty in the estimated displacement is higher.

4. Discussion

Most of the models used to describe post-seismic relaxation are based on the assumption of a single relaxation time estimated for each relaxation mechanism proposed [4,11,14,60,61]. A co-seismic deformation leaves the system in a metastable state taking place in the broad earthquake zone. In the present work, the post-seismic relaxation is viewed as a superposition of several different relaxation times, each one activated in different mesoscopic domains of the macroscopic relaxed seismic volume. The macroscopic relaxation observed is a result of the superposition of different relaxation times on the microscopic–mesoscopic scale. This slow thermodynamic recovery could be explained by a physical model based on the Arrhenius Equation [62].

The lithosphere as a dynamical system under a post-seismic relaxation is a complex systems with strong interactions between its mesoscopic “subvolumes” [63] and, thus, could be viewed as a macroscopic system driven by hierarchically constrained dynamics [28] involved in the macroscopically observed relaxation process as a sequential series of correlated mesoscopic processes.

The view of hierarchical constrained dynamics could be a physical pathway to describe many of the complex systems that present logarithmic relaxation [64,65,66]. In the case of post-seismic phase, the applied tectonic load drives the system to fracture, and its relaxation creates a reorganization of the stress field in some selected mesoscopic regions of the system. The deformed subvolumes nearest to the fractured region are rearranged and the system reorganized. These new stress changes will affect the deformed subvolumes outside the area nearest to the fracture, permitting their relaxation, and so on. The above picture suggests that the relaxation in the fractured fault system is slow because it consists of a large number of correlated and strongly interacting reorganization steps with an increasing relaxation time. It seems to be appropriate to describe macroscopical post-seismic relaxation in terms of a lithosphere that presents a hierarchical pattern, with the faster relaxation to be associated with the subvolume close to earthquake focus regions and, due to strong correlations and long-range interactions, the slower relaxation time to be related with the “outer” deformed, but in the broad vicinity of epicenter, regions. The latter is in agreement with that proposed for complex relaxed systems with cluster interactions, where the relaxation time τ scales with the relaxation volume as τ~V^α, where α is an exponent related macroscopically with the fractal geometry of the subvolume distribution and the dynamical nature of relaxation [67,68,69]. This type of hierarchically constrained dynamics is able to derive the logarithmic relaxation frequently observed and associate it with thermally activated processes. We note that this behavior can result either from the parallel relaxation of subvolumes or by a sequential relaxation series of strongly correlated complex events [28].

Our results are in agreement with that presented in [70], where a Burridge-Knopoff (BK) model is used to describe the time dependence of the mean stress, assuming that a uniform statistical distribution of stress values is acting on each block (i.e., deformation in subvolumes) in the relaxing lithospheric system.

In [70], the post-seismic normalized stress relaxation of the strongly interacted elements of the BK model presents a stress relaxation that could be viewed as hierarchically constrained since blocks increase the stress on their nearby elements when they slip, and as a consequence, the probability that their neighbors will slip increases as well. They conclude that the normalized mean stress <σ> during the post-seismic relaxation in a BK model with a uniform distribution of stresses [71] is:

$$<\sigma >={{\displaystyle \int }}_0^1\sigma e^{-\nu t\exp (-\beta {\epsilon }_m\left(1-{\sigma }^2\right))}d\sigma$$

(17)

where ε_m is the maximum size of energy barrier, ν is an attempt frequency to jump over the energy barrier, and β = 1/k_BT. Approximation of (17) leads to the following behavior in short, intermediate, and long time periods, respectively:

$$\underset{t\to 0}{\lim }<\sigma >=\frac{1}{2}-\frac{1}{2\beta {\epsilon }_m}\left(1-\exp \left(-\beta {\epsilon }_m\right)\right)\nu t$$

(18)

$${<\sigma >}_{int}\approx \frac{1}{2}-\frac{1}{2\beta {\epsilon }_m}(\gamma +\ln \left(\nu t\right)$$

(19)

where γ is the Euler–Mascheroni constant (γ ≈ 0.577).

$$\underset{t\to \infty }{\lim }<\sigma >=\frac{\exp \left(\beta {\epsilon }_m\right)}{2\beta {\epsilon }_m}\frac{\exp (-\nu t\exp \left(-\beta {\epsilon }_m\right)}{\nu t}$$

(20)

which, for very long time periods, t is proportional to exp(−t/τ) with τ = exp(β${\epsilon }_m$)/ν.

A comparison of Equations (18)–(20) with Equation (10), (12), and (14), respectively, implies that $S ~\frac{1}{2\beta {\epsilon }_m}$, supporting the hypothesis that the scaling factor depends on temperature and the energy barrier that must be overcome before the relaxed subvolume moves to a new stress stage. In addition, ${\tau }_{max} ~\frac{1}{\nu }\exp \left(\beta {\epsilon }_m\right) \mathrm{and} {\tau }_{min} ~\frac{1}{\nu }$ with $\ln \left(\frac{{\tau }_{max}}{{\tau }_{min}}\right)~\beta {\epsilon }_m$ being an expression similar to that defined in the present approach using the Arrhenius expression, if ${\epsilon }_m=E_{max}-E_{min}$ is the maximum size of energy barrier.

The relaxation function, as presented in Equation (6), explains the logarithmic relaxation observed. In addition, it results in finite values, both just after the perturbation (i.e., t = 0) and for a long time after the earthquake event (t → ∞). The relaxation function R(t) can be explained by Arrhenius’ law, presenting a universal behavior of the relaxation process. Following Equation (6) and the ideas of hierarchically constrained dynamics, the faster relaxation mechanisms, related to small relaxed volumes (small τ), offer a more significant contribution to the relaxation process than the slower relaxation mechanisms (large τ related with the broader relaxed area), as indicated by the presence of term 1/τ in the expression of R(t), which permits greater weight to fast relaxation processes. We note that the superposition of relaxation mechanisms, as expressed in Equation (6), implies that, for early time periods, the fast relaxation mechanism associated with small relaxed volumes that have a relaxation time close to τ_min dominates, while for late time periods, the slow relaxation mechanism, related to the relaxation of volumes farther from the epicenter and present relaxation time close to τ_max, contributes most. It is reasonable to assume that the maximum relaxation time depends on the perturbation that triggers the relaxation and, thus, on the deformed volume related with the main event, along with the geotectonic conditions in the deformed region, such as stress, temperature, and fluids.

In addition, measurements of τ_min could give information on the fastest relaxation mechanism, which is scaled with the smallest deformed subvolume, even if it is sometimes not clearly observed due to the immediate transition of the deformation process into the logarithmic range, resulting in an unreliable estimation of τ_min.

Applying the generic hierarchically constrained dynamics view to the North–South displacement of the post-seismic relaxation of Lefkas 17 November 2015, Mw6.5 earthquake, we conclude that this approach describes with good agreement the displacement as estimated using cGNSS recordings during the post-seismic relaxation period, leading to the estimation of τ_min ≈ 3.5 days and τ_max ≈ 350 days, respectively.

This is a rather new approach that is used to describe post-seismic relaxation based on the assumption of a spectrum of relaxation times that are superposed, each of them activated in different mesoscopic domains of the macroscopic relaxed seismic volume. This slow thermodynamic recovery (i.e., a logarithmic time dependence) could be explained by a physical model based on the Arrhenius equation. With such an approach, we are led to finite values for t = 0, without assuming any empirical constitutional equation.

5. Concluding Remarks

In the present work the post-seismic relaxation is viewed as a superposition of several different relaxation times, each one activated in different mesoscopic domains of the macroscopic relaxed seismic volume. The macroscopic relaxation that is observed is a result of the superposition of different relaxation times on the microscopic–mesoscopic scale. This approach is a new one since most of the models used to describe post-seismic relaxation are based on the assumption of a single relaxation time estimated for each relaxation mechanism proposed. The slow thermodynamic recovery of the post-seismic relaxation could be explained by a physical model based on the Arrhenius equation and an appropriate distribution of relaxation times.

In summary, we can state that, in the present work:

• A multiscale post-seismic relaxation mechanism based on the existence of a distribution in relaxation time is presented.
• Assuming an Arrhenius dependence of the relaxation time with the uniform distributed activation energy in a mesoscopic scale, we conclude there is a generic logarithmic type relaxation on a macroscopic scale.
• The hierarchically constrained dynamics model could be used to understand the evolution of post seismic relaxation.
• The model was applied in the case of 2015 Lefkas (Greece) Mw6.5 earthquake, where cGNSS time series were recorded in a station located in the vicinity of the epicentral area. The application of the present approach to the Lefkas event fits the observed displacements, implying a distribution of relaxation times in the range τ_min ≈ 3.5 days to τ_max ≈ 350 days.

These results offer a stimulating finding for further studies that will investigate post-seismic deformation timeseries for other strong earthquakes, which may provide new insights into the underlying dynamics of earthquake relaxation processes, along with the possible parameters that influence this behaviour (i.e., fault type, local conditions, regional stress field, etc.).