Have the 1999 Izmit–Düzce earthquakes influenced the motion and seismicity of the Anatolian microplate?

SUMMARY

In the current plate tectonics paradigm, relative plate motions are assumed to remain unperturbed by temporal stress changes occurring during the seismic cycle, whereby stress slowly built up along tectonic plate boundaries is suddenly released by rapid fault slip during earthquakes. However, direct observations that could challenge such a tenet have not been identified so far. Here we show that the rigid motion of the whole Anatolian microplate, measured using space geodetic techniques, was altered by the stress released during the 1999 Izmit–Düzce earthquakes, which ruptured along the North Anatolian Fault. This kinematic change requires a torque change that is in agreement with the torque change imparted upon the Anatolian microplate by the Izmit–Düzce coseismic stress release. This inference holds across realistic ranges of data noise and controlling parameters, and is not hindered by active deformation in western Anatolia. These results suggest the existence of a whole-plate kinematic signal associated with the stress released by large earthquakes.

Plate motions, Seismic cycle, Satellite geodesy, Dynamics of lithosphere and mantle

1 INTRODUCTION

The theory of plate tectonics (Wilson 1965; McKenzie & Parker 1967; Morgan 1968; Le Pichon 1968), which describes the rigid motions of Earth’s |$\sim \, 100$|-km-thick outer shell (e.g. Gordon & Jurdy 1986; DeMets et al. 2010), revolutionized the geosciences by offering a framework to understand geological observations, including the global distribution of earthquakes, which occur predominantly along tectonic plate margins. In this framework, relative motions between adjacent tectonic plates, which today can be directly measured via space geodetic techniques (e.g. Dixon 1991; Stein 1993; Argus et al. 2011), result over periods of few decades in the slow buildup of stress along their crustal interfaces, followed by its sudden release through earthquakes. In fact, from the viewpoint of plate dynamics, the entire earthquake cycle can be thought as comprising two phases in which two torques are built upon a plate: a first phase in which a plate experiences a torque whose magnitude slowly grows as the stress builds up along a portion of its margins during the interseismic period. In the subsequent phase, corresponding to the earthquake itself, the plate experiences a sudden, additional stress over the finite size of the rupture area that is equal in magnitude but opposite in direction to the one developed during the first phase. Such a stress counters the one that has grown during the slow interseismic phase—hence the term stress drop or stress release over the rupture area. To first degree, the torque associated with the earthquake is similar in magnitude but opposite in direction to the one developed during the interseismic phase. This way, the net torque change upon the plate from the beginning to the end of the earthquake cycle remains negligible.

It is commonly accepted that earthquakes do not alter the rigid plate motions that generate the stress buildup, but instead induce post-seismic deformation in the form of afterslip and/or viscous stress diffusion in the near-field of the event (Pollitz et al. 2006). This may also be evident in the form of small (<100 km) block deformation along plate margins, as inferred in the vicinity of the 2011 M_W 9.0 Tohoku earthquake epicentre (Meade & Loveless 2017). Such a paradigm is ubiquitous, for example, in models of the seismic cycle (see Govers et al. 2018, for a review) or the slip deficit of locked crustal faults (e.g. Savage 1983). It stems from the tenet that torques associated with the release of stress during earthquakes (Allmann & Shearer 2009) remain small compared to the viscous resistance exerted by Earth’s asthenosphere at the plates base, which must be overcome in order to change plate motions over time (e.g. Iaffaldano 2014).

Recent analyses of long-wavelength glacial rebound data provide tighter constraints on the viscosity and thickness of the asthenosphere (Paulson & Richards 2009; Richards & Lenardic 2018), and thus on the resistance the latter exerts against plate motion changes. This, in the context of recently discovered microplates (rigidly-moving tectonic units whose typical length does not exceed |$\sim \, 1000~$| km, Wallace et al. 2005) calls for a revision of the tenet above, because forces arising from the rupture of locked plate margins during large earthquakes can be comparable to the viscous resistance exerted by the asthenosphere on the limited basal area of microplates or relatively small plates. Numerical synthetic tests of microplate viscoelastic dynamics over timescales comparable to the typical cycle of large earthquakes provide theoretical support to this (Martin de Blas & Iaffaldano 2019), which nonetheless remains to be tested by actual observations.

Here, we resort to publicly-available Global Positioning System (GPS) data to investigate whether the rigid motion of the Anatolian microplate (AT) was different before and after the 1999 Izmit–Düzce seismic events, and whether any difference is consistent with the torque change imparted upon the microplate by the coseismic stress released over the rupture area. We use GPS data collected at stations deployed within the core of AT to estimate its rigid-body instantaneous rotation, expressed in the form of an Euler vector, prior to and after the 1999 Izmit–Düzce events, and isolate a temporal change of the microplate motion. The GPS data utilized here consist of both continuous and campaign (survey mode) measurements. As the latter ones are more prone to be affected by local processes unrelated to rigid plate motions (i.e. data noise), we perform a suite of tests to rule out with significant confidence that the observed change in AT motion from before to after the 1999 Izmit–Düzce events owes to such noise. Next, we address whether the stress released over the rupture area during the 1999 Izmit–Düzce earthquakes generated the torque required to explain the observed AT rigid motion change. Lastly, we compare Gutenberg–Richter statistics of the seismicity along the AT margins before and after the 1999 Izmit–Düzce events, and correlate their temporal changes with the modified loading rates caused by its changed motion relative to the adjacent tectonic plates.

2 THE ANATOLIAN MICROPLATE TECTONIC SETTING

The Anatolian region (McKenzie 1972), a prominent example of microplate, has been the focus of numerous studies of seismicity, crustal strain accumulation and associated hazard (e.g. England et al. 2016; DeVries et al. 2017). Its tectonic structure features a very low seismicity core of continental lithosphere bounded by the North Anatolian Fault (NAF) to the north (Barka 1992), the East Anatolian Fault (EAF) to the east (Bulut et al. 2012), the Cyprus Arc to the south (Wdowinski et al. 2006) and the Western Anatolian Extensional Province (WAEP) to the west (Fig. 1). Previous studies have established, on the basis of survey-mode and continuous GPS measurements gathered since the 1990s, that AT moves rigidly at rates up to 30 mm yr⁻¹ counterclockwise relative to Eurasia (EU, Reilinger et al. 1997, 2006; McClusky et al. 2000; England et al. 2016), and that internal deformation is <2 mm yr⁻¹ (McClusky et al. 2000). Deformation in the WAEP (Aktug et al. 2009) occurs in the form of diffuse north–south stretching (Bird 2003) at rates of around 10 ± 5 mm yr⁻¹. This means that surface velocities in that region depart by up to 10 ± 5 mm yr⁻¹ from the velocity associated with the AT rigid plate motion. Such a pattern is also evident from analyses of recently collected InSAR images (Weiss et al. 2020).

Figure 1.

GPS-measured surface velocities relative to the ITRF08 reference frame. In red are velocities obtained from episodic measurements between 1994 and 1998. In blue are velocities obtained from continuous measurements between 2013 and 2015. Black ellipses indicate 95 per cent confidence velocity uncertainties. Empty squares show GPS station positions according to the colour legend above. Green dots are M_W > 5 earthquakes occurring between 1975 and the present day. Earthquake epicentres come from the International Seismological Centre (ISC) seismic catalogue—ISC 2020 Online Bulletin https://doi.org/10.31905/D808B830. Purple dots are the epicentres of the 1999 Izmit and Düzce earthquakes. Major active faults are shown in light orange (Emre et al. 2018). Continents are shown in grey. Contours of tectonic plates are shown with solid black lines. The Anatolian zone featuring N–S stretching is in black hatched. AR = Arabia plate, AT = Anatolian plate, EU = Eurasia plate and NU = Nubia plate.

Most of the seismicity in AT occurs along the NAF, an active right-lateral strike-slip fault that separates AT from EU, and that staged over 15 earthquakes of magnitude M_W ≥ 6.5 since the beginning of the last century (DeVries et al. 2017; Barka 1992; Stein et al. 1997). The most recent of these is the 2014 M_W = 6.9 Samothraki–Gökçeada event along the North Aegean Trough (Saltogianni et al. 2015). Significant seismicity also occurs along the EAF (McKenzie 1972; Bulut et al. 2012), which separates AT from the Arabian plate (AR), and which recently staged the 2020 M_W = 6.8 Elaziǧ seismic event (Doğru et al. 2020). The largest earthquake registered along the AT margins in the past 50 yr is the 1999 August 17 M_W = 7.5 strike-slip event, which ruptured a 120-km-long segment of the NAF (Reilinger et al. 2000; Tibi et al. 2001) near the city of Izmit. This event was followed shortly after, on 1999 November 12, by an M_W = 7.2 strike-slip earthquake that ruptured the adjacent 40 km of the NAF, east of the Izmit rupture segment and next to the city of Düzce (Ayhan et al. 2001; Utkucu et al. 2003). These two events are often referred to as the 1999 Izmit–Düzce earthquakes because of their proximity in time and space (Fig. 1).

3 TEMPORAL CHANGE OF ANATOLIAN RIGID MOTION

3.1 Available GPS data and inferred Euler vectors

Prior to the 1999 Izmit–Düzce earthquakes, tectonic velocities inside AT are constrained through publicly-available, survey-mode GPS solutions collected in three different campaigns during 1994, 1996 and 1998 (Reilinger et al. 1997; McClusky et al. 2000; Reilinger et al. 2006). Among these, we preselect stations that are located inside AT, sufficiently away from its tectonic margins (Bird 2003) and the WAEP (McClusky et al. 2000; Nocquet 2012), and that are thus not affected by interseismic elastic strain accumulation. We use the GAMIT/GLOBK software (Herring et al. 2018) to process data from these stations and obtain horizontal velocities as well as uncertainties (Fig. 1) relative to the International Terrestrial Reference Frame 2008 (ITRF08, Altamimi et al. 2012). Next, we perform a final selection of GPS sites (red squares in Fig. 1) based on (i) their velocity uncertainties and (ii) the residual velocities upon inversion for the AT Euler vector. Specifically, we request that, upon inferring the Euler vector, station residuals are smaller than twice the associated standard deviation (95 per cent confidence). Table 1 reports observation days for each GPS station, as well as surface velocities and associated standard deviations relative to ITRF08.

Table 1.

Availability of GPS stations before the 1999 Izmit–Düzce seismic event. Columns are: (1) sites acronyms, (2–3) longitude and latitude of site location (in decimal degrees), (4–6) days of the year for which measurements are available and (7–10) station east/north velocities (v_e/v_n) and standard deviations (σ) in  mm yr⁻¹.

GPS stations			Stations availability			GPS velocities
Site	Lon	Lat	1994	1996	1998	ve	vn	|$\sigma _{v_e}$|	|$\sigma _{v_n}$|
ANKA	32.726	39.871	279–280	267–268	–	5.9	8.6	2.3	2.5
ANKS	32.759	39.888	258–280	266–271	–	4.1	11.9	1.1	1.1
MELE	33.191	37.378	258–280	267–269	261–263	14.3	12.7	1.9	1.9
MERS	34.552	36.900	274–276	263–266	257–259	17.3	13.9	2.1	2.2
YOZG	34.813	39.801	270–272	259–261	–	8.4	12.8	1.9	2.0
ANKR	32.758	39.887	Continuous IGS station			4.1	11.9	1.1	1.1

The motion of AT after the 1999 Izmit–Düzce earthquakes is constrained by GPS data measured via continuously operating stations (Blewitt et al. 2018) that are publicly available via the Nevada Geodetic Laboratory repository. We select 14 stations (blue squares in Fig. 1) following the same criteria utilized for the motion of AT before the 1999 Izmit–Düzce earthquakes. While the associated position time-series are available for the period between 2009 and 2015, we deliberately truncate them to the most recent 3-yr available interval—longer than the 2.5-yr minimum standard commonly used to obtain tectonically representative velocities (Blewitt & Lavallée 2002; Herring et al. 2018)—that is, from 2013 to 2015. This choice is motivated by the fact that tomography studies (e.g. Priestley & McKenzie 2013) constrain the Maxwell time interval of the asthenosphere underneath central AT, where stations are located, to be around 7 yr (see the Supporting Information for details). Choosing stations (i) away from margins and (ii) whose position time-series begin around two Maxwell time intervals after 1999 (i.e. starting from 2013) allows us to minimize the imprint of earthquake-induced transient viscoelastic motions onto the surface tectonic velocities (Pollitz et al. 2006; Govers et al. 2018; Martin de Blas & Iaffaldano 2019), evident in post-seismic studies of the Izmit–Düzce earthquakes (Reilinger et al. 2000; Bürgmann et al. 2002b; Hearn et al. 2009) as well as in tectonic studies using data collected over periods across or immediately after these earthquakes (Aktuğ et al. 2013). For the selected stations, we obtain horizontal velocities and associated uncertainties relative to ITRF08 using the MIDAS software (Blewitt et al. 2016, Fig. 1). Table 2 reports the GPS measuring times for every station, as well as surface velocities and associated standard deviations.

Table 2.

Availability of continuous GPS stations after the 1999 Izmit–Düzce seismic event. Columns are: (1) sites acronyms, (2–3) longitude and latitude of site location (in decimal degrees), (4–5) starting/ending time of GPS measurements (in decimal years) and (6–9) station East/North velocities (v_e/v_n) and standard deviations (σ) in mm yr⁻¹.

GPS stations			Stations availability		GPS velocities
Site	Lon	Lat	Start	End	ve	vn	|$\sigma _{v_e}$|	|$\sigma _{v_n}$|
AKDG	35.872	39.660	2013.00	2015.83	6.0	17.5	0.5	0.6
ANKR	32.758	39.887	2013.00	2015.83	2.4	10.6	1.2	1.2
BOGZ	35.255	39.194	2013.00	2015.83	6.6	19.1	0.8	0.9
CIHA	32.922	38.650	2013.00	2015.83	6.5	11.3	0.5	0.6
FEEK	35.912	37.815	2013.00	2015.83	11.4	18.5	0.9	0.7
GEME	36.081	39.185	2013.02	2015.83	5.1	16.9	0.7	1.2
HALP	34.183	37.445	2013.00	2015.83	11.9	15.7	0.6	1.1
KAYS	35.524	38.708	2013.00	2015.83	9.5	17.2	0.6	0.6
KIRS	34.155	39.165	2013.01	2015.83	5.0	13.1	0.5	0.5
KKAL	33.518	39.843	2013.00	2015.83	4.1	12.1	0.6	0.7
KLUU	33.065	39.079	2013.00	2015.83	6.1	10.6	0.7	0.6
NEVS	34.703	38.617	2013.00	2015.83	8.1	16.7	0.5	0.4
NIGD	34.679	37.959	2013.07	2015.83	10.4	16.0	0.6	0.6
POZA	34.872	37.422	2013.00	2015.83	12.7	17.0	1.0	0.8

Next, we use the two sets of horizontal velocities and uncertainties for the 1994–1998 and 2013–2015 time intervals to infer Euler vectors and associated covariances that are representative of the AT/ITRF08 rigid motion before and after the 1999 Izmit–Düzce earthquakes. We do so utilizing the commonly used minimization of the sum of squared velocity misfits (e.g. Nocquet et al. 2001). Specifically, for the period before the 1999 Izmit–Düzce earthquakes we draw 10⁶ samples of each horizontal velocity from their mean values and standard deviations. This provides us with 10⁶ sets of horizontal velocities, from which we infer an ensemble of 10⁶ Euler vectors representative of the AT rigid motion before the 1999 Izmit–Düzce earthquakes. Through the same procedure, we infer an equally large ensemble of Euler vectors for the period after the earthquakes. Lastly, from both ensembles we calculate mean Euler vectors and covariances of their Cartesian components. Table 3 presents both Euler vectors relative to ITRF08, together with entries of their covariance matrices. Fig. 2 illustrates the inferred Euler vectors (poles and angular velocities) as well as the AT surface velocity fields associated with them. The inferred Euler vectors indicate that after the 1999 Izmit–Düzce earthquakes, AT features a different rigid motion than before 1999. The change between the 1994–1998 and 2013–2015 AT/ITRF08 Euler vectors, which expresses the temporal variation in AT rigid motion, concerns mostly the Euler pole, which wanders ∼90 km southward, and to a much lesser degree the angular velocity, which remains steady within its uncertainty. The fact that the Euler pole wanders beyond the 95 per cent confidence ellipses (which encompass 95 per cent of the Euler vector ensembles mentioned above) means that the change in AT rigid motion is warranted by the data at the 95 per cent confidence level. This plate motion change can also be visualized by looking at the surface velocity fields calculated at arbitrary locations for the periods before and after the earthquakes (see Fig. 2). Additionally, Fig. S6 in the Supporting Information shows the difference between these two velocity fields, thus illustrating a speed-up of the Anatolian plate. Lastly, Fig. S7 in the Supporting Information shows the difference between the GPS velocities for the period 2013–2015 and the predicted velocity field for the period 1994–1998 calculated at the locations of the GPS sites used for the period 2013–2015.

Figure 2.

Euler vectors for the AT rigid motions before (in red) and after (in blue) the 1999 Izmit–Düzce events. Euler vectors are inferred—through minimization of the sum of squared residuals—from ensembles of surface velocities drawn from the GPS-measured velocities and their uncertainties. Euler poles (calculated as ensemble averages) are shown in the map: dots with error ellipses (95 per cent confidence) are relative to ITRF08. Triangles without ellipses are relative to fixed EU, and are calculated by adding the rigid motion of EU relative to ITRF08 (Altamimi et al. 2012) to the AT/ITRF08 Euler vectors. Inset shows the ensemble distributions of angular velocities relative to ITRF08. Green contours highlight the regions where the most recurrent 95 per cent of poles would fall if one adds noise following a Gaussian distribution whose standard deviation equals to 1 mm yr⁻¹and assuming that noise is spatially correlated (we use a cross-site correlation coefficient of 0.5). Arrows within AT show predicted velocities relative to fixed EU at arbitrary locations. Empty squares, continents, plate margins and acronyms are as in Fig. 1.

Table 3.

Euler vectors before and after the Izmit–Düzce seismic event. We present the Euler pole location (in decimal degrees), the angular velocity (in |$^{\circ }\, \mathrm{Myr^{-1}}$|⁠) and the associated covariance matrix C elements (in |$10^{-8}\, \mathrm{rad\:^2\:Myr^{-2}}$|⁠).

Time period	Lon	Lat	ω	Cxx	Cxy	Cxz	Cyy	Cyz	Czz
[yr]	[° E]	[° N]	[° Myr−1]	[10−8 rad 2 Myr−2]
1994–1998	28.658	41.383	1.710	983	645	950	426	625	921
2013–2015	28.636	40.651	1.749	124	85	120	59	83	116

GPS station ANKR is the only one featuring a continuous position time-series record that bridges across the time periods before and after the Izmit–Düzce earthquakes analysed here. Fig. S8 in the Supporting Information shows its velocity between 1996 and 2016, calculated using 3-yr-long intervals. The initial and final velocity of ANKR (in red in Fig. S8, Supporting Information) indeed reflect the observed change of AT rigid motion illustrated by the 1994–1998 and 2013–2015 Euler vectors (Fig. 2): a change in direction of motion, with the velocity magnitude remaining steady. It is out of the scope of this study to test hypotheses on the possible origin of the temporal variations of ANKR’s velocity throughout the period from 1996 to 2016. However, if one were to assume that the 1996–2016 velocity record of ANKR is sufficiently representative of the temporal evolution of the AT rigid motion throughout this time period, then ANKR’s record may be utilized to explore how the response of the laterally heterogeneous, viscoelastic asthenosphere transitions from elastic in the short term to viscous in the longer (i.e. longer than its Maxwell interval) term—a problem for which an analytical solution exists only in the simpler case of laterally uniform model asthenosphere (for more details, see Martin de Blas & Iaffaldano 2019).

Upon inferring an Euler vector from surface velocities via the minimization of the sum of squared velocity misfits, it is standard practice to verify a posteriori that in fact the observed input velocities represent a rigid motion on the sphere (i.e. one that is adequately described by Euler vectors). This is achieved by verifying that the differences between input observed velocities and those predicted by the output Euler vector (commonly referred to as residual velocities) are acceptably small relative to the uncertainties associated with the input data. That is, one would sample distributions of the input velocities using such uncertainties, obtain distributions of Euler vector and residuals, and find that the ranges where most samples of the residuals distributions (e.g. 95 per cent) as well as their mean values fall encompass zero. This allows one to conclude that residuals are zero within the associated uncertainties, and that the most–recurrent part of the inferred Euler-vector distribution coincides with the true value of the Euler vector. A more challenging scenario is one where the input velocity measurements are a combination of a rigid plate motion plus other components not associated with rigid plate motions, for instance arising from hydrological or atmospheric loadings (e.g. Blewitt & Lavallée 2002), that affect each available data point (in the following we refer to any velocity component that is not a realization of rigid motions as random noise). The challenge in this case comes from the fact that one generally does not know a priori whether any random noise component is present in the data, nor what its magnitude might be. On the one hand, if the unknown noise is larger than the known uncertainties on the input velocity measurements, then the residuals would not be zero within the associated uncertainties, and thus one would readily conclude on the non-rigidity of the input velocity field due to the presence of noise. If, on the other hand, the unknown noise is less than the uncertainties on input velocity measurements, the residuals would nonetheless be zero within the associated uncertainties, and the inferred Euler vector would still provide an adequate description of the input velocities to the extent allowed by the uncertainties of the latter ones. However, not knowing how random noise compares to the input uncertainties leaves one wondering how much has noise impacted on the accuracy of the Euler-vector inference. In other words, larger-than-noise uncertainties on input velocities might in principle allow noise to go undetected in the analyses of the residuals, while still driving the Euler-vector inference somewhat away from the value one would otherwise infer from non-noisy velocity measurements.

From the AT/ITRF08 Euler vectors (⁠|$\vec{\omega }$|⁠), we predict station–site velocities (⁠|$\vec{v}_{P} = \vec{\omega }\times \vec{r}$|⁠, |$\vec{r}$| being the station-site position vector) and then calculate residual velocities (⁠|$\vec{v}_{R}$|⁠) as the difference between the GPS observed velocity (⁠|$\vec{v}_{O}$|⁠) and its prediction from the Euler vectors at each station site—that is, |$\vec{v}_{R} = \vec{v}_{O} - \vec{v}_{P}$|⁠. Table 4 and Fig. 3 report residual velocities for the period before the 1999 Izmit–Düzce earthquakes, while Table 5 and Fig. 4 report residual velocities for the period after the seismic events. In both cases station residuals are smaller than or at most equal to their uncertainties. This indicates that the velocity fields measured by geodetic stations over both time intervals are in fact adequately described through the use of Euler vectors—that is, they are velocity fields representing rigid motions over the Earth’s surface. However, the fact remains that the station–site velocities we utilize to constrain the 1994–1998 AT/ITRF08 Euler vector are obtained from campaign measurements as opposed to continuous ones, and therefore feature standard deviations that are in some cases up to four times larger than those of the velocities utilized for the 2013–2015 Euler vector, which are instead inferred from continuous position time-series. Such a difference is reflected in the larger covariance and residual uncertainties associated with the 1994–1998 Euler-vector inference, compared to the 2013–2015 one. On this basis, as well as on the basis of the considerations above, it would not be unreasonable to hypothesize that the larger uncertainties on the campaign measurements utilized for the 1994–1998 interval might potentially allow any random noise present in the data to affect the accuracy of the 1994–1998 Euler–vector inference (i.e. the distance between such an inference and the true value of the 1994–1998 Euler vector) without being detected. In other words, one may suspect that noise in the 1994–1998 campaign measurements might be such that the actual change in AT rigid motion is in reality less than the one inferred here, or even that the AT rigid motion has in reality not changed at all from before to after the 1999 earthquakes.

Figure 3.

Residual velocities reported in Table 4. Ellipses are 95 per cent confidence. Continents, plate margins and plate acronyms are as in Fig. 1.

Figure 4.

Residual velocities reported in Table 5. Ellipses are 95 per cent confidence. Continents, plate margins and plate acronyms are as in Fig. 1.

Table 4.

Residual velocities at stations used to constrain the AT Euler vector before the 1999 Izmit–Düzce event. Columns are: (1) sites acronyms, (2–3) longitude and latitude of site location (in decimal degrees), (4–7) station east/north residual velocities (R_e/R_n) and standard deviations (σ) in  mm yr⁻¹ and (8) correlation coefficient (ρ_en) between R_e and R_n.

Site	Re	Rn	|$\sigma _{R_e}$|	|$\sigma _{R_n}$|	ρen
ANKA	0.1	−1.6	1.8	2.1	−0.1
ANKS	−1.8	1.7	1.2	1.2	−0.01
ANKR	−1.8	1.7	1.2	1.2	−0.01
MELE	0.1	1.4	1.6	1.7	0.05
MERS	1.3	−0.8	1.4	1.8	−0.2
YOZG	2.0	−2.5	1.7	1.7	0.1

Table 5.

Same as Table 4, for stations used to constrain the AT Euler vector after the 1999 Izmit–Düzce event.

Site	Re	Rn	|$\sigma _{R_e}$|	|$\sigma _{R_n}$|	ρen
AKDG	1.3	−1.1	0.5	0.6	−0.02
ANKR	−1.0	−0.0	1.1	1.0	−0.2
BOGZ	0.4	2.1	0.8	0.8	0.02
CIHA	−1.2	0.2	0.5	0.6	−0.003
FEEK	0.4	−0.1	0.8	0.6	−0.06
GEME	−1.3	−2.2	0.7	1.0	0.04
HALP	−0.0	1.4	0.6	1.0	0.03
KAYS	1.6	−0.5	0.6	0.6	−0.002
KIRS	−1.1	−1.1	0.5	0.5	0.008
KKAL	0.4	−0.5	0.6	0.7	−0.01
KLUU	−0.1	−0.8	0.7	0.6	0.001
NEVS	0.1	1.1	0.5	0.5	0.001
NIGD	0.1	0.5	0.6	0.6	−0.003
POZA	0.6	0.9	0.9	0.8	−0.03

3.2 Testing the potential impact of data noise

In order to address this, we perform two semi-synthetic tests on the available GPS data. Their semi-synthetic nature stems from the fact that we utilize noise values generated in a computer environment (so that we can store them and perform statistical analyses a posteriori), but we apply them to the AT rigid motions inferred from the actual GPS measurements. In the first test, we presume that the Euler vectors inferred for both time intervals are accurate within the associated covariances, and explore how large would noise need to be in order for the two inferences to start overlapping, in which case one could no longer exclude with absolute certainty that the Euler vector has remained steady through time. In the second, more stringent suite of tests, we start from the assumption that the difference between the less constrained 1994–1998 Euler vector and the 2013–2015 one is only apparent (i.e. not real) and owes to random noise in the 1994–1998 campaign measurements being present but undetected due to the larger measurement uncertainties. The combination of both tests aims at demonstrating that, despite the larger standard deviations of the 1994–1998 station-site velocities, the Euler vector difference illustrated in Fig. 2 is highly unlikely to owe to data noise, and shall thus be interpreted as a real change in the AT rigid motion.

We start by adding random noise values to the station-site velocities before and after the 1999 Izmit–Düzce seismic event, and then recalculate the AT Euler vector. The value of noise added to each station velocity is randomly selected from a Gaussian distribution whose standard deviation is 1 mm yr⁻¹. This means that noise in both the east and north components of the site velocities can be up to ∼3 mm yr⁻¹ (three standard deviations comprising 99.8 per cent of the Gaussian distribution). In addition, we take noise that is spatially correlated across sites (Wdowinski et al. 1997), and assume that the cross-site correlation coefficient is equal to 0.5, which has been previously determined to be appropriate for the inter-stations distances at hand (e.g. Wang et al. 2012, 2018). We repeat such test 10⁶ times in order to generate an ensemble of noise-related realizations of the Euler vectors before and after the 1999 Izmit–Düzce earthquakes. Next, we infer the contours/intervals where 95 per cent of the ensembles of Euler poles and angular velocities fall. Contours for the Euler poles do not overlap (see Fig. 2) even for a high level of noise of up to ∼3 mm yr⁻¹ in the GPS network and a cross-site correlation coefficient of 0.5. Furthermore, the direction of Euler-pole change inferred from geodetic data from 1994–1998 to 2013–2015 is highly skewed relative to the main axes of the noise-related ellipses (Fig. 2), making the observed Euler-pole change unlikely to be the result of data noise.

For the second suite of tests, we proceed by assuming that the AT Euler vector has not changed from before to after the 1999 earthquakes, but has remained constant and equal to the one we infer for the 2013–2015 time interval, which is based on a larger number of continuous station-site velocities. Initially, we perform an F-ratio test for two Euler-vector solutions: the one obtained using velocities between 2013 and 2015, and the one obtained using all velocities from the two time periods—consistent with the assumed steadiness of the AT motion. This test provides the residual probability that the null hypothesis (i.e. that the AT Euler vector remained steady through time, and thus that the two sets of observed velocities can be described by the same Euler vector) can be accepted. We obtain an F-ratio of 1.75 that, given the degrees of freedom associated with the two data sets, corresponds to a residual probability of 11.6 per cent. Next, we proceed by accepting the null hypothesis despite its low probability, and calculate surface velocities at each of the positions of the 1994–1998 station–sites (6 positions, see Table 1) using the 2013–2015 Euler vector. We take these as the velocities that campaign observations would yield in the absence of any noise. We then add to each station-site velocity an ensemble of 10⁶ values of noise that have been randomly sampled from a Gaussian distribution whose standard deviation is a free parameter in range from 0.5 to 10 mm yr⁻¹ (referred to as level of noise). As in the previous test, we assume noise to be correlated across sites with correlation coefficients equal to 0.5 and 0.9. We generate noise ensembles independently for each station-site velocity. The level of noise increases at steps of 0.25 mm yr⁻¹. For each level, we generate ensembles of 10⁶ noisy station-site velocities, and calculate 10⁶ realizations of noisy Euler vector. Next, we quantitatively assess the extent to which the 10⁶ samples of noisy Euler vectors fall near the 1994–1998 Euler pole inferred from the observed campaign velocities (in red in Fig. 2). Specifically, we count what fraction of such ensemble (i.e. how many of the 10⁶ samples) falls within a contour/interval where the most-recurrent 20 per cent and 30 per cent of 1994–1998 Euler pole/angular velocity inferred from real data lie (referred to in the following as target ranges). Given the size of the ensembles (10⁶), such fraction represents an estimate of the probability that noise in the data utilized for the 1994–1998 time interval is responsible for the Euler-vector inference being different from that obtained using the 2013–2015 time-interval data. Fig. 5 reports such estimates. The fact that, for a given target range, the probability initially increases with the assumed level of random noise to then decrease should not be surprising: as the noise level increases, the ensemble of noisy Euler vectors lies initially close to the 2013–2015 Euler vector, to then migrate towards the 1994–1998 one, to finally move progressively away from it. The probability is around 4 per cent in the most optimistic case (i.e. for the largest values of target range and correlation coefficient), and is maximum for noise levels in range from 1.5 to 2 mm yr⁻¹.

Figure 5.

Percentage of the ensemble of noisy Euler vectors that falls within the 20 per cent and 30 per cent target ranges plotted against the assumed level of noise, for assumed cross-site correlation coefficients of 0.5 and 0.9 (see the main text for details on noise tests).

These results indicate that if one starts from the presumption that the AT Euler vector has not changed from before to after the 1999 earthquakes, and that the different Euler-vector solution obtained from the 1994–1998 and 2013–2015 sets of data owe to noise in the former set, then one is bound to accept the following: (i) there is ∼10 per cent probability that the two data sets can be adequately described by one single Euler vector. (ii) In the unlikely event that this was actually the case, the probability that noise is responsible for the apparently different Euler-vector solution obtained from the 1994–1998 data set is at most ∼4 per cent. In order to accept the hypothesis of Euler vector steadiness given the available data, both events should occur. However, the probability of both events occurring is the product of the probabilities estimated above—that is, ∼0.4 per cent. On this basis, we conclude that having used campaign data for the 1994–1998 time interval affects the precision of the inference (i.e. the covariance of the inferred Euler vector, also reflected in the red confidence region/range in Fig. 2) rather than its accuracy (i.e. the proximity of the inferred Euler vector to the true one). Therefore, we regard the Euler-vector change illustrated in Fig. 2 as representative of a real change in the AT rigid motion.

4 AT MICROPLATE TORQUE VARIATIONS

In order to determine whether the observed rigid motion change of AT results from the coseismic stress drop occurred over the Izmit–Düzce rupture area, we estimate the torque variation experienced by AT as a result of the 1999 Izmit–Düzce rupture, and compare it with estimates of the torque variation required upon AT in order to change its Euler vector to the degree evidenced by the geodetic data.

4.1 Torque variation required upon AT

For the latter torque-variation calculation, we build on previous studies (Iaffaldano & Bunge 2015; Martin de Blas & Iaffaldano 2019) that derived an analytical equation linking the torque variation |$\Delta \vec{M}$| experienced by a tectonic plate of basal area S to the resulting change of Euler vector |$\Delta \vec{\omega }$|⁠, which represents a rigid motion change in an absolute reference frame (e.g. fixed with the relatively more sluggish deeper mantle). Specifically, by writing two initial equations for the torque balance of a tectonic plate at two different points in time, and then taking the difference between them (so that only the terms that have actually changed through time survive in the resulting equation) one obtains

$$\begin{eqnarray*} \Delta \vec{M} = \int _{S} \frac{\mu _A}{H_A} \cdot \vec{r} \times [\vec{v}_p(\vec{r},t_2) - \vec{v}_p(\vec{r},t_1) ] \cdot dS \end{eqnarray*}$$

(1)

where |$\vec{v}_p(\vec{r},t)$| is the generic surface velocity at position |$\vec{r}$| (within the plate) and time t (see Iaffaldano & Bunge 2015; Martin de Blas & Iaffaldano 2019, for more details and equation development). If tectonic plates behave perfectly rigidly, then |$\vec{v}_p(\vec{r},t) = \vec{\omega }(t) \times \vec{r}$|⁠, |$\vec{\omega }(t)$| being the plate Euler vector at time t. In the more realistic case of a tectonic plate that features regions of internal deformation (as, for instance, the N–S stretching in the WAEP), the surface velocity is the summation of a rigid component and a deformation component (⁠|$\vec{v}_d$|⁠) — that is, |$\vec{v}_p(\vec{r},t) = \vec{\omega }(t) \times \vec{r} + \vec{v}_d(\vec{r},t)$|⁠. On the basis of earlier studies (see Gordon 1998, for a comprehensive review), it is reasonable to expect that |$\vec{v}_d = 0$| over the largest part of a plate, while |$\vec{v}_d \ne 0$| only in the deforming part of a plate area S_d < <S. The torque-variation equation thus becomes

$$\begin{eqnarray*} \Delta \vec{M} = \int _{S} \frac{\mu _A}{H_A} \cdot \vec{r} \times [(\vec{\omega }(t_2) \times \vec{r}) - (\vec{\omega }(t_1) \times \vec{r}) ] \cdot dS + \int _{S_d} \frac{\mu _A}{H_A} \cdot \vec{r} \times [\vec{v}_d(\vec{r},t_2) - \vec{v}_d(\vec{r},t_1) ] \cdot dS_d \end{eqnarray*}$$

(2)

or

$$\begin{eqnarray*} \Delta \vec{M} = \int _{S} \frac{\mu _A}{H_A} \cdot \vec{r} \times [\Delta \vec{\omega }\times \vec{r}] \cdot dS + \int _{S_d} \frac{\mu _A}{H_A} \cdot \vec{r} \times \Delta \vec{v}_d(\vec{r}) \cdot dS_d \end{eqnarray*}$$

(3)

where |$\Delta \vec{\omega }= \vec{\omega }(t_2) - \vec{\omega }(t_1)$| and |$\Delta \vec{v}_d(\vec{r},t_2,t_1) = \vec{v}_d(\vec{r},t_2) - \vec{v}_d(\vec{r},t_1)$|⁠. The latter integral on the right-hand side of eq. (3) represents the contribution to the torque-variation estimate arising from temporal changes in the pattern of internal deformation. Its expression carries the assumption that |$\Delta \vec{v}_d$| is the same at the surface and at the lithosphere/asthenosphere boundary. In other words, it carries the assumption that deformation evident at the surface propagates vertically down to the asthenosphere. This is arguably an unlikely scenario, which in cases such as deformation in the WAEP stands at odds with the observed depths of intraplate earthquakes. One may account for this by introducing a positive multiplying coefficient c ≤ 1. Instead, the former integral on the right-hand side of eq. (3) represents the contribution to the torque-variation estimate arising from plate rigidity. Since the Euler-vector change |$\Delta \vec{\omega }$| is independent of |$\vec{r}$|⁠, the whole integral can be written as the result of a linear map

$$\begin{eqnarray*} \int _{S} \frac{\mu _A}{H_A} \cdot \vec{r} \times [\Delta \vec{\omega }\times \vec{r}] \cdot dS = \underbrace{ {\begin{pmatrix}\int _{S} \frac{\mu _A}{H_A}(y^2+z^2)dS & -\int _{S} \frac{\mu _A}{H_A}xydS & -\int _{S} \frac{\mu _A}{H_A}xzdS \\ -\int _{S} \frac{\mu _A}{H_A}xydS & \int _{S} \frac{\mu _A}{H_A}(x^2+z^2)dS & -\int _{S} \frac{\mu _A}{H_A}yzdS \\ -\int _{S} \frac{\mu _A}{H_A}xzdS & -\int _{S} \frac{\mu _A}{H_A}yzdS & \int _{S} \frac{\mu _A}{H_A}(x^2+y^2)dS \end{pmatrix}}}_{\mathbf {P}} {{\begin{pmatrix}\Delta \omega _x \\ \Delta \omega _y \\ \Delta \omega _z \end{pmatrix}}} \end{eqnarray*}$$

(4)

where the operator |$\mathbf {P}$| links the vectorial space of rigid Euler-vector changes to that of torque variations acting upon a tectonic plate. Therefore,

$$\begin{eqnarray*} \Delta \vec{M} = \mathbf {P} \Delta \vec{\omega }+ \int _{S_d} \frac{\mu _A}{H_A} \cdot \vec{r} \times c\Delta \vec{v}_d(\vec{r}) \cdot dS_d \end{eqnarray*}$$

(5)

Eq. (5) above holds true when the intervals of time over which |$\Delta \vec{M}$| and |$\Delta \vec{\omega }$| are accomplished are (i) short enough (in a geological sense) that the plate has not changed shape, but at the same time (ii) long enough that the elastic behaviour of the viscoelastic asthenosphere can be neglected. Indeed such conditions apply to the case of AT for Euler vectors calculated using geodetic data collected in the periods 1994–1998 and 2013–2015, since these are more than two Maxwell time-intervals apart from each other. Furthermore, eq. (5) indicates that plate non-rigidity would impact on the torque-variation estimate to a non-negligible degree only if, across the interval of time at hand, the pattern of deformation changed (i.e. |$\Delta \vec{v}_d \ne 0$|⁠) to a degree that compensates for the fact that c ≤ 1 and S_d < <S. Table 6 reports the entries of the operator |$\mathbf {P}$| specific of AT when |$\mu _M = 1.5 \times 10^{21}\, \mathrm{Pa \cdot s}$| and μ_A is set variably to 5 × 10¹⁹ and |$1 \times 10^{20}\, \mathrm{Pa \cdot s}$|⁠. In using eq. (5), we implement here the findings of previous studies of long-wavelength glacial rebound data (Paulson & Richards 2009; Richards & Lenardic 2018), which constrain the cube of the asthenosphere thickness (H_A) to be proportional to the viscosity contrast between the asthenosphere (μ_A) and the lower part of the upper mantle (μ_M). That is, |$H_A^3 = a \cdot (\mu _A/\mu _M)$|⁠, with |$a = 4.67 \times 10^5~\, \mathrm{km^3}$|⁠. Setting the global average of μ_A to |$5 \times 10^{19}\, \mathrm{Pa \cdot s}$| and that of μ_M to |$1.5 \cdot 10^{21}\, \mathrm{Pa \cdot s}$| yields H_A ∼ 150 km. Instead, H_A ∼ 190 km when the global average of μ_A is set to |$1 \times 10^{20}\, \mathrm{Pa \cdot s}$|⁠. Furthermore, we allow for asthenosphere lateral variations of viscosity compatible with the lateral temperature variations mapped by the tomography model PM_v2_2012 by Priestley & McKenzie (2013, see the Supporting Information for more details and images).

Table 6.

Entries of the linear operator |$\mathbf {P}$| in eq.(4). Units are |$10^{39}\, \mathrm{Pa \cdot s \cdot m^3}$|⁠.

μA	|$\mathbf {P}_{11}$|	|$\mathbf {P}_{12}$|	|$\mathbf {P}_{13}$|	|$\mathbf {P}_{22}$|	|$\mathbf {P}_{23}$|	|$\mathbf {P}_{33}$|
|$5 \times 10^{19}\, \mathrm{Pa \cdot s}$|	1.56	−0.77	−1.13	2.27	−0.72	1.71
|$1 \times 10^{20}\, \mathrm{Pa \cdot s}$|	3.37	−1.68	−2.43	4.89	−1.57	3.72

Next, we draw 10⁶ samples of the AT Euler vectors that describe the rigid motion respectively before and after the 1999 Izmit–Düzce events. We use these to build an ensemble of geodetically measured AT Euler-vector changes that we then map, through eq. (5), into an ensemble of torque variations required to generate the temporal change of AT rigid motion inferred from geodetic data. In order to estimate the impact of deformation in the WAEP onto the required torque variation, we perform two additional calculations where we assume that the pattern of N–S stretching observed in the WAEP increases/decreases by up to 7 mm yr⁻¹ from 1994–1998 to 2013–2015. The western Anatolia deformation pattern appears to have remained unchanged when comparing geodetic estimates of crustal deformation performed before (Reilinger et al. 1997) and after (Nocquet 2012) the end of the last century. These estimates feature an uncertainty of ±5 mm yr⁻¹. Therefore, calculations that assume an increase/decrease of deformation rates up to 7 mm yr⁻¹ arguably encompass what the change in deformation pattern might have been in reality.

4.2 Torque variations resulting from the 1999 Izmit–Düzce rupture

Before expressing in analytical form the torque variation imparted by an earthquake upon a tectonic plate, it is appropriate to make the following clarifications: (i) earthquakes are currently viewed as processes where slip occurs on a fault of finite area—hence the term finite-fault solution used when imaging the amount of slip from seismic and/or geodetic observations. This is different from the early view of earthquakes as processes having an equivalent force system represented by a double couple (i.e. no-net-torque) in an intact (i.e. unfaulted) medium, which served the purpose of explaining the compression/decompression observed radiation pattern while casting the pattern source as a simple point source (e.g. Julian et al. 1998). With a view of faults having a finite area, earthquakes do impart torques to tectonic plates, because slip occurs in a given direction on a finite-size fault of given orientation. (ii) One shall expect that the earthquake stress drop triggers an elastic-body response as well as a rigid-body response of the adjacent plates sharing the rupturing fault, which is then transferred to the viscoelastic asthenosphere. However, observations constrain the amount of slip associated with an earthquake, not the stress drop. This means that when using an elastic model to fit observations, typically made in the vicinity of a ruptured fault, one actually infers the portion of stress drop that triggered the elastic response, not necessarily the entire stress drop value, nor the portion that would trigger a rigid-body response. This study is concerned with testing the existence of the latter one. Therefore, in order to circumvent the issue above, we only make the assumption that the relationship between stress drop over the rupture area and amount of coseismic slip is linear (Madariaga 1977), but keep the stress drop value as a free parameter in our formulation.

We calculate torque variations upon AT caused by the 1999 Izmit–Düzce rupture (⁠|$\Delta \vec{M}_e$|⁠) by integrating the coseismic stress drop occurred during the Izmit–Düzce earthquakes (Reilinger et al. 2000; Bürgmann et al. 2002a), which lies within the NAF strike-slip rupture planes and is parallel to the local AT/EU relative motion, over the finite rupture area. Coseismic stress is released quasi-instantaneously compared to the time interval over which geodetic velocities are calculated, and is opposite to the stress that has been building up slowly before the 1999 Izmit–Düzce events, when the sections of the NAF that ruptured in 1999 were locked. We assume that the stress drop over the rupture area scales linearly with the amount of coseismic slip that has been estimated through independent geodetic data collected in the vicinity of the rupture (Reilinger et al. 2000), and that are not used here to determine the AT Euler vectors. In quantitative terms, we implement the following equation

$$\begin{eqnarray*} \Delta \vec{M}_e = \int _{\Sigma _r} \vec{r} \times \bigg [ \Delta \sigma _{\mathrm{ max}} \cdot \frac{s(\vec{r})}{s_{\mathrm{ max}}} \cdot \frac{\vec{v}_{AT/EU}(\vec{r})}{| \vec{v}_{AT/EU}(\vec{r}) |} \bigg ] \cdot d\Sigma _r \end{eqnarray*}$$

(6)

where Δσ_max is the maximum stress drop over the rupture area Σ_r, |$s(\vec{r})$| is the slip at position |$\vec{r}$| over Σ_r, s_max is the maximum slip, while |$\vec{v}_{AT/EU}(\vec{r})$| is the motion of AT relative to EU at position |$\vec{r}$|⁠. We keep Δσ_max, occurring where slip is maximum, as a free parameter that we vary from 2 to 10 MPa, in line with global statistics of strike-slip earthquakes (Allmann & Shearer 2009). We vary the maximum depth of rupture (i.e. the depth extent of Σ_r) from 10 to 15 km (Reilinger et al. 2000). We use lateral extent and strike of rupture at the surface (i.e. intersection of Σ_r and Earth’s surface) put forth in previous studies (Reilinger et al. 2000; Bürgmann et al. 2002a), and assume that the strike may vary by up to ±5°. We randomly select values of the controlling parameters within the ranges/sets above, and integrate the infinitesimal torque associated with the infinitesimal stress-drop vector over the rupture area surface. This allows us to generate an ensemble of 10⁶ samples of torque variations caused upon AT by the 1999 Izmit–Düzce rupture. Through the same approach, we account for the effect of the smaller 2014 M_W = 6.9 Samothraki–Gökçeada event, which occurred during the period of continuous GPS recordings used here to calculate the AT Euler vector. This event featured a stress drop 10 times smaller than that of the 1999 Izmit–Düzce events (Saltogianni et al. 2015; Konca et al. 2018).

4.3 Comparison of torque variations

We express the torque variations computed above in terms of their magnitude and axis, the latter displayed as a pole located where the torque-variation axis intersects the Earth’s surface (Fig. 6). We find that the torque variations required to explain the AT rigid motion change are in good agreement with those generated upon AT by the 1999 Izmit–Düzce earthquakes. In fact, 99 per cent of the ensemble poles associated with the torque variation caused by the 1999 Izmit–Düzce rupture fall near the most likely pole associated with the torque variation required upon AT, and anyways within the region where the most-likely 68 per cent of the latter ensemble falls, which corresponds to 1.1 per cent of the Earth surface. This means that, had the AT Euler vector change been the result of a random process, then the probability that such a process would result in a torque-variation pole falling within the 68 per cent confidence contour of the torque-variation pole required to explain the geodetically inferred kinematic change would be as small as 1.1 per cent. Required torque-variation magnitudes (inset in Fig. 6) are in the order of 10²² N·m, less than the typical torques driving/resisting the motions of larger tectonic plates (Copley et al. 2010; Iaffaldano 2014), but in line with the fact that microplates offer a limited basal area to asthenospheric viscous resistance. Their distribution trades off with the average viscosity of the asthenosphere. Assuming the latter to be |$5\times 10^{19}\, \mathrm{Pa \cdot s}$| (Fig. 6) results in a good agreement with a scenario where the maximum stress drop over the rupture area is in range from 2 to 6 MPa. Instead, assuming the average asthenosphere viscosity to be |$1\times 10^{20}\, \mathrm{Pa \cdot s}$| (Fig. S9, Supporting Information) results in a good agreement with a scenario where the maximum stress drop over the rupture area is in range from 6 to 10 MPa. Thus, the distribution of required torque-variation magnitude is in good agreement with the magnitude of the torque variation provided by the 1999 Izmit–Düzce earthquakes within a realistic range of average asthenosphere viscosity values. This inference also holds when we account for the impact that deformation in the WAEP might have on the estimate of the torque variation required upon AT.

Figure 6.

Comparison of the distributions of torque-variation ensembles (i) required in order to explain the geodetically observed Euler-vector change illustrated in Fig. 2 (in green and light/dark grey) and (ii) modelled from the rupture characteristics of the 1999 Izmit–Düzce events (in red/blue). Distributions of torque-variation magnitudes are in the inset. Torque-variation poles are in the map: dots are the ensemble averages, contours/hatched regions show the area where the most recurrent 68 or 95 per cent of the ensembles fall (see the figure legend for details). Green parentheses show the probability that a point randomly selected on Earth’s surface would fall within the hatched area. Required torque variations are calculated using an average asthenosphere viscosity of |$5\times 10^{19}\, \mathrm{Pa \cdot s}$|⁠. The average viscosity of the lower part of the upper mantle is set to |$1.5 \times 10^{21}\, \mathrm{Pa \cdot s}$|⁠. Continents, plate margins and acronyms are as in Fig. 1.

5 DISCUSSION AND CONCLUSIONS

The AT rigid motion change evidenced by the analysis of geodetic velocities would impact, in turn, on the loading rates (i.e. the local fault-plane-parallel velocities) along the crustal boundaries separating AT from the adjacent EU, AR and Nubia (NU) plates. We calculate AT margins loading rates as follows: we subtract the geodetically measured Euler vectors for the motions of EU, AR and NU relative to ITRF08 (Altamimi et al. 2012) from the ensembles of Euler vectors of the AT/ITRF08 motion before/after the 1999 Izmit–Düzce earthquakes (Fig. 2). This allows us obtaining ensembles of Euler vectors for the AT/EU, AT/AR and AT/NU rigid motions before/after the 1999 Izmit–Düzce earthquakes (Fig. 2 and Figs S10 and S11, Supporting Information). We use these to calculate ensembles of the surface velocities (before/after the 1999 events) along the NAF, EAF and Cyprus Arc, and project these velocities onto the associated margin planes to obtain ensembles of loading rates before/after the 1999 Izmit–Düzce earthquakes. The regime (i.e. strike-slip, normal, or reverse) and orientation (i.e. dipping toward/away from the AT interior) of the NAF, EAF and Cyprus Arc are taken from the tectonic model PB2002 (Bird 2003). The dip angle of strike-slip margins is randomly chosen across the ensemble in range from 70° to 90°. The dip angle of normal margins is randomly chosen in range from 40° to 70°. The dip angle of reverse margins is randomly chosen in range from 20° to 40°. Lastly, we subtract the loading rate calculated for the period before the 1999 Izmit–Düzce earthquakes from that calculated for the period after, in order to obtain an ensemble of temporal variations of loading rates. Fig. 7(a) shows the mean (main panel) and the values of the lower/upper extremes (insets) of the ranges within which the most recurrent 68 per cent of the ensemble falls. Within the 68 per cent confidence associated with the AT Euler vectors, we find that AT plate boundary loading rates increased by up to 2.5 ± 1 mm yr⁻¹ along the EAF and the NAF, and decreased by around 1.5 ± 1 mm yr⁻¹ along the Cyprus Arc.

Temporal changes of loading rates predicted from Euler vectors should not be expected to perfectly match estimates of loading rates inferred from geodetic observations taken in the vicinity of the AT margins, as local margin deformation impacts on the latter, but is not captured in the former ones. Nonetheless, it is reasonable to expect that, to first degree, these estimates align with each other. Recent studies (Ergintav et al. 2014; Ozarpaci et al. 2021) that utilized geodetic data collected near the 1999 Izmit–Düzce rupture segment since 2014 found some evidence that the local loading rate increased by around 2 mm yr⁻¹ relative to estimates preceding the 1999 events. The inference of the present study is based on independent data, and is indeed compatible with these findings.

To first order, these changes appear to correlate with temporal changes in the number of moderate/large (i.e. M_W ≥ 5) earthquakes (Ekstrom et al. 2012) that have occurred along these margins from before to after the 1999 Izmit–Düzce earthquakes (Fig. 7). To illustrate this, we infer Gutenberg–Richter statistics along the tectonic margins of AT by resorting to the International Seismological Centre (ISC) catalogue of earthquakes—ISC 2020 Online Bulletin https://doi.org/10.31905/D808B830. Specifically, we take all events with an available estimate of M_W larger than or equal to 3 that have occurred between 1970 and 2020. From these, we remove those events that could be regarded as aftershocks of the large earthquakes (M_W ≥ 6) in the set. Specifically, we define an aftershock as any earthquake (i) within three rupture lengths from the epicentre, (ii) whose magnitude is at least 0.5 smaller than the main shock and (iii) occurring at the latest three months after the main shock. Furthermore, we only select the earthquakes that occurred within 50 km on either side of the NAF, the EAF and the Cyprus Arc. We then re-group these events within 10-yr-long continuous time intervals between 1970 and the first half of 1999, and between 2000 and 2020. These time frames divide the available catalogue into periods before/after the 1999 Izmit–Düzce earthquakes. A linear regression of the log₁₀ counts over these time intervals is representative of the data only for M_W ≥ 5 (Fig. 8). This indicates that the completeness threshold of the catalogue for the specific region and time intervals analysed is M_W = 5. For each tectonic segment and within each 10-yr-long time interval, we count the number of events whose magnitude is equal to or larger than values between 3 and 7. Figs 7(b)–(d) report, for each tectonic segment, the log₁₀ of average (solid line) and minimum-to-maximum range (hatched/light-coloured areas) of the 10-yr-long counts before/after (green/black) the 1999 Izmit–Düzce earthquakes. The comparison of loading rates and seismicity along the AT margins features some degree of uncertainty propagating from the range of values of the parameters involved in the calculation (i.e. Euler vectors, dip angles, seismic catalogue temporal coverage). Nonetheless, it suggests that seismicity tends to increase where loading rates also do—an inference that holds even when one limits the analyses to earthquakes occurred within 25 km on either side of the NAF, the EAF and the Cyprus Arc (Fig. S12, Supporting Information).

Figure 7.

(a) Temporal variation of margins loading rates from before to after the 1999 Izmit–Düzce earthquakes, calculated from the difference of the projection of along-margin surface velocities (relative to neighbouring plates) onto the margin planes. Insets show the limits of the most recurrent 68 per cent of loading-rate temporal variations calculated from the ensemble of Euler vectors. Plate margins and acronyms are as in Fig. 1. (b)–(d) Gutenberg–Richter statistics of the seismicity along the Cyprus Arc (segment AB), the EAF (segment BC) and the NAF (segment CD) from before to after the 1999 Izmit–Düzce events.

Figure 8.

Completeness tests of the portion of the ISC catalogue used here, which includes earthquakes whose M_W ≥ 3, between 1970 January and 2020 April, between 21.70°E and 43.34°E, between 32.49°N and 44.52°N.

Lastly, we acknowledge that processes such as surface mass changes or post-seismic viscous relaxation might impact the surface velocities within AT. While testing alternative dynamic scenarios is outside the scope of this study, we note that they involve viscoelastic responses that would unfold over timescales in the order of the Maxwell time interval. For most of AT, this is estimated to be ≤7 yrs when the average asthenosphere viscosity is assumed to be |$5 \times 10^{19}\, \mathrm{Pa \cdot s}$|⁠, and ≤15 yr when the average viscosity is assumed to be |$1 \times 10^{20}\, \mathrm{Pa \cdot s}$|⁠. The fact that the geodetic observations utilized here were performed at least 15 yr (i.e. 1–2 Maxwell time intervals) apart from each other (from 1994–1998 to 2013–2015) indicates that, if processes other than those considered here were at work, they are of second-order importance, as their effects would have decayed significantly across the time period covered by the geodetic observations used in this study. We therefore conclude that the coseismic stress release associated with the Izmit–Düzce earthquakes modified the kinematics of the Anatolian microplate by a small, yet observable amount. This microplate kinematic change, in turn, corresponds to a torque change that is in remarkable agreement with the torque change imparted upon Anatolia by the Izmit–Düzce coseismic stress release. This result is in line with early theoretical findings (Martin de Blas & Iaffaldano 2019) and suggests that microplate kinematics can be affected in a geodetically measurable way by large earthquakes along their boundaries. At the same time, we acknowledge that the validity of the modelling assumptions made here, particularly those made in estimating the torque variations imparted by large earthquakes to tectonic plates, shall be corroborated with further studies focused on other tectonic settings.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Figure S1. Depth-averaged asthenosphere viscosity μ_A underneath AT, for the case where the global average of μ_A is set to 5 × 10¹⁹ Pa · s. Lateral variations of μ_A are inferred from lateral temperature variations in the seismic tomography model of Priestley and McKenzie (2013). Plate margins and acronyms are as in Fig. 1 in the main text.

Figure S2. As Supplementary Fig. 1, for the case where the global average of μ_A is set to 1 × 10²⁰ Pa · s.

Figure S3. Depth-averaged Young’s modulus of the asthenosphere underneath AT, inferred from lateral temperature variations in the seismic tomography model of Priestley and McKenzie (2013). Plate margins and acronyms are as in Fig. S1.

Figure S4. Depth-averaged asthenosphere Maxwell time interval underneath AT, for the case where the global average of μ_A is set to 5  × 10¹⁹ Pa · s. Plate margins and acronyms are as in Fig. S1.

Figure S5. As Fig. S4, for the case where the global average of μ_A is set to 1 × 10²⁰ Pa · s.

Figure S6. Difference between velocities calculated at arbitrary locations for the periods 2013–2015 and 1994–1998. Continents, plate margins and plate acronyms are as in Fig. 1 in the main text.

Figure S7. Difference between the GPS velocities for the period 2013–2015 and the predicted velocity field for the period 1994–1998 calculated at the locations of the GPS sites used for the period 2013–2015. Ellipses are 95 per cent confidence. Continents, plate margins and plate acronyms are as in Fig. 1 in the main text.

Figure S8. Evolution of the (a) surface velocity magnitude and (b) azimuth of the GPS continuous station ANKR. Open circles represent the velocity magnitude/azimuth across time, while vertical lines show the associated 68 per cent confidence ranges. Grey dashed lines show the time intervals for which these velocities have been obtained. In red are velocities corresponding to the periods used to calculate the Euler vectors before/after the Izmit–Düzce earthquakes.

Figure S9. As Fig. 6 in the main text, for the global average of μ_A set to 1 × 10²⁰ Pa · s.

Figure S10. As Fig. 2 in the main text, but using AR as reference-frame plate.

Figure S11. As Fig. 2 in the main text, but using NU as reference-frame plate.

Figure S12. As Fig. 7 in the main text, but considering earthquakes within 25 km on either side of the AT margins.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

All data and methods used in this study come from referenced sources and are described herein. The authors are grateful to the editor (Bert Vermeersen), Jean–Mathieu Nocquet and Marco Bohnhoff for constructive reviews that substantially improved the manuscript. JMB acknowledges support from a fellowship (LCF/BQ/EU17/11590058) from “la Caixa“ banking foundation (ID 100010434), as well as support from the Department of Geosciences and Natural Resource Management at the University of Copenhagen. GI acknowledges support from a Villum Fonden research grant (00023121). EC acknowledges support from the Institut Universitaire de France.

DATA AVAILABILITY

GPS data used here are openly available from the UNAVCO (http://www.unavco.org), Nevada Geodetic Laboratory (http://geodesy.unr.edu/) (Blewitt et al. 2018) and IGS (www.igs.org) archives. Survey-mode GPS data between 1994 and 1998 were processed using the GAMIT/GLOBK software (Herring et al. 2018), publicly available at MIT (http://geoweb.mit.edu/gg/). Surface velocities for the period 2013–2015 were obtained using the MIDAS software (Blewitt et al. 2016), which is publicly available at the Nevada Geodetic Laboratory (http://geodesy.unr.edu/). Some of the figures were generated using Generic Mapping Tools (Wessel et al. 2013).