Influence of Upper Mantle Anisotropy on Isotropic P‐Wave Tomography Images Obtained in the Eastern Mediterranean Region

Seismic body‐wave tomography studies typically assume an isotropic upper mantle, possibly mapping anisotropy into artificial isotropic velocity anomalies in the resulting images. The Eastern Mediterranean with its oceanic, continental, and extinct subduction systems, as well as dense station coverage, provides an ideal setting to explore this issue. To examine the influence of seismic anisotropy, our study deals with both synthetic and real data inversions in which realistic seismic anisotropy models derived from 3D mantle convection simulations and shear wave splitting measurements are taken as a priori constraints. Spatial large‐scale velocity perturbations are mostly consistent between models derived with and without considering anisotropy. Small differences in the magnitude (up to 2%) and shape of velocity perturbations occur, and some structures are less diffuse when including anisotropy. Additionally, good backazimuthal coverage of teleseismic events and a larger data set improve the resolution of our model with respect to previous tomography studies and allow us to better interpret first‐order isotropic velocity anomalies. Key features, such as the half‐arc subducting oceanic plate in the southern Aegean and a wide and deep tear in the slab beneath southwestern Turkey, are clearly visible in all models. Our final tomography images also provide evidence for a shallow horizontal tear in the northern Hellenides and a vertical tear between two parts of the Cyprian slab. In eastern Anatolia, slab‐related high‐velocity anomalies are absent due to the continental collision and break‐off.


Introduction
Seismic tomography is a well-established method that is frequently used to investigate upper mantle subduction tectonics and kinematics by constraining variations in seismic velocity. Over the last two decades, there have been many P and/or S wave tomography studies conducted to understand the complex tectonic setting beneath the Eastern Mediterranean region (e.g., Biryol et al., 2011;Blom et al., 2020;Çubuk-Sabuncu et al., 2017;Piromallo & Morelli, 2003;Portner et al., 2018) (Figure 1). Although these studies employed different data sets, the resolved images of the lithosphere and asthenosphere have implied very similar features in the mantle such as subducting and fragmented slabs as well as hot upwelling material. However, it is noted that some anomalies, for instance, those underneath the slab or small velocity perturbations, are often left out from interpretations, which mostly focus on the large-scale characteristics of seismic velocity anomalies.
In the Eastern Mediterranean region, the Hellenic (Aegean) subduction zone, which accommodates the convergent plate motion between the Nubian and Eurasian plates, has been active since the Eocene (e.g., Jolivet et al., 2015) and has produced present-day seismic anisotropy in various directions and strengths. Confal et al. (2018) numerically simulated the strain-induced lattice preferred orientations (LPO) of A-type olivine in the upper mantle for a hypothetical subduction zone resembling the Hellenic subduction system. The modeling results are consistent with observations on seismic anisotropy, i.e., station averaged fast polarization directions (FPDs) primarily based on shear wave splitting (SWS) measurements, which indicate the presence of strong anisotropy evidenced by splitting time delays (TDs) of up to 2 s (e.g., Confal et al., 2016;Evangelidis et al., 2011;Paul et al., 2014). These measured time delays between fast and slow shear waves are proportional to the thickness of the anisotropic layer and the strength of the anisotropy. For teleseismic isotropic tomography studies with incoming waves that have steep incidence angles, the propagation direction would be subparallel to the anisotropic slow direction in layers with horizontally aligned fast and intermediate directions. Consequently, this produces delayed arrival times that can be mapped by isotropic tomography into a low-velocity anomaly as shown by Bezada et al. (2016). Conversely, mantle structures with vertically aligned fabrics would produce high-velocity perturbations ( Figure 2). Seismic velocity perturbations of up to ±3-4% in various tomographic studies (e.g., P and S waves and surface waves) are, in general, attributed to the lateral variations in temperature and partial melt, or the composition of crustal and mantle rocks. Velocity anomalies caused by anisotropy could interfere with real isotropic velocity anomalies and lead to incorrect interpretations of the thermophysical parameters. Sobolev et al. (1999) used synthetic data to show that in isotropic P-wave tomography, significant artifacts can occur for some specific anisotropy structures, especially in regions with a dipping olivine a-axis. Lloyd and van der Lee (2008) examined the influence of anisotropy on both S and Rayleigh wave derived isotropic tomographic images obtained for North America. They were able to quantify the anisotropic bias in images and concluded that the magnitude of the bias decreased with increasing depth location of the anisotropic material. Eakin et al. (2010) Taymaz et al. (2007) and Jolivet et al. (2013) and references therein. Black-pink arrows represent extensional regimes (about 15 mm/yr in the northern Aegean and Gulf of Corinth). The triangles represent stations and their respective number of waveforms used in this tomography study. Inset in the lower right shows averaged GPS vectors (black arrows in mm/yr) of the plates with respect to a stable Eurasian plate after Reilinger et al. (2006) and Le Pichon and Kreemer (2010). Anisotropy parameters from direct S wave and SKS wave studies in the Eastern Mediterranean were retrieved from the splitting database of Wüstefeld and Bokelmann (2007) and Paul et al. (2014). anisotropic time corrections to P-wave traveltime residuals that were calculated based on SKS splitting parameters for the western United States. Their results highlighted the importance of the magnitude of anisotropy rather than its orientation, under the assumption that the axis of symmetry is horizontal. Furthermore, Eken et al. (2012) show that large-scale anisotropy related to fabrics of the continental mantle lithosphere (typically characterized by dipping symmetry axes) contaminated tomographic images in some parts of their models beneath the Fennoscandian shield and concluded that this effect should not be ignored. Both O'Driscoll et al. (2011) and Eken et al. (2012) reported that magnitude and anisotropic orientation could influence P-wave traveltime residuals and noted that tomography methods need to be improved, especially for subduction regimes. These studies also emphasize the importance of azimuthal coverage of earthquake distribution. Calculating azimuthal terms using observed SKS splitting parameters may provide a quick method as a first-order approximation. However, this approach could have disadvantages, primarily because SKS phases arrive subvertically and are polarized in nearly horizontal directions, which makes them sensitive to horizontal anisotropy only (O'Driscoll et al., 2011). Directly inverting for anisotropy and isotropic velocity (e.g., Eberhart-Phillips & Reyners, 2009;Ma et al., 2019;Wang & Zhao, 2013;Wei et al., 2019;You & Zhao, 2012) would be a solution to overcome some of these limitations. However, data availability issues mainly related to the lack of dense permanent station networks and poor azimuthal coverage make the inversion even more underdetermined, and a loss in resolution is usually the consequence when solving for additional unknowns. More recently, Bezada et al. (2016) tested the effect of anisotropy on a synthetic inversion experiment. A strain-induced LPO of upper mantle aggregates in a 3D subduction model was established prior to the tomographic inversion. They observed that artificial velocity anomalies produced by anisotropy could be up to hundreds of kilometers wide and that including estimates of the anisotropy field as a priori constraints could be useful in reducing these artifacts. Estimates of the anisotropy field could come, for instance, from numerical models that simulate the geodynamic evolution of the region and consequent mantle flow and LPO formation.
In this study, we invert a large number of teleseismic P-wave traveltime residuals from several permanent and temporary seismic networks that have operated in the Aegean and Anatolia. The larger and more comprehensive data set gives our model improved resolution compared to other available tomographic models (e.g., Biryol et al., 2011;Piromallo & Morelli, 2003;Portner et al., 2018), which results in novel insights into the tectonic features, of Greece and the Aegean, improving the knowledge of the velocity structures present in the upper mantle of the Eastern Mediterranean and Anatolia. In addition to the new isotropic model, to better characterize the 3D P-wave velocity anomalies, we apply a correction to our observed P-wave traveltime residuals using anisotropy parameters inferred from numerical models as well as SWS measurements. Discrepancies between the corrected and uncorrected tomographic models reveal the role of seismic anisotropy in changing the velocity perturbations in the upper mantle beneath this active tectonic region.

Tectonic Setting of the Region
The Eastern Mediterranean has been affected by ancient and current subduction systems since the Cretaceous (Görür, 1988), impacting the crust and generating destructive earthquakes in the region (Taymaz et al., 2004). This study investigates the more recent subduction process starting in the Oligocene (Jolivet, 2001), from about 30 Ma (million years ago) until the present. Currently, the Nubian and Arabian plates are actively moving toward the Eurasian plate (Reilinger et al., 2006), while the Anatolian plateau is moving to the west and the Aegean is characterized by an extensional regime due to subduction rollback (e.g., Le Pichon & Kreemer, 2010). Continental lithosphere is subducting at the northwestern end of the Hellenic subduction system (Evangelidis, 2017;Pearce et al., 2012), while southeast of the Kefalonia Transform Fault (KTF) the Ionian oceanic lithosphere is subducting in a half-arc. The subduction rates at the northern Hellenic trench (continental part) are much slower than at the southern trench (oceanic part) (e.g., McClusky et al., 2000), which offsets the southern trench by 70-85 km at the KTF (Pearce et al., 2012). In contrast, below Cyprus, the subducting slab is no longer active and is considered to be in transition to continental collision similar to eastern Anatolia (e.g., Feld et al., 2017), where collision is forming mountain belts in southeastern Anatolia and the Caucasus (Tan & Taymaz, 2006). The Hellenic and Cyprus subduction systems have experienced tearing in several locations. A geological and tectonic reconstruction study by Jolivet et al. (2013) dated the tearing in southwestern Turkey to ∼15 Ma, with subsequent slab break-off in eastern Anatolia occurring around 10 Ma and the most recent tear in the Gulf of Corinth starting around 5 Ma.
Body and surface wave tomography studies (Biryol et al., 2011;Çubuk-Sabuncu et al., 2017;Govers & Fichtner, 2016;Salaün et al., 2012;Portner et al., 2018;Taymaz, 1996;Wei et al., 2019) indicate low-velocity anomalies and separation of the Aegean and Cyprian slab in southwestern Turkey, at depths between 50 and at least 300 km. They interpret this as a north-south fragmentation of the Hellenic slab, which started in the Eocene-Miocene as a result of the different retreat rates of the trenches. The rollback in the Aegean and the steepening of the Cyprus slab started subsequently. Even thermal anomalies in western Anatolia, extending from the upper mantle into the crust, might be generated by the rollback and slab tearing (Roche et al., 2019). Biryol et al. (2011) andPortner et al. (2018) found evidence for a smaller tear in the Cyprian trench, defining an eastern and a western Cyprian slab and linked the tearing to volcanism in the Central Anatolian Volcanic Province. Biryol et al. (2011) found a vertical tear reaching ∼200 km depth, while Portner et al. (2018) interpreted a horizontal tear propagating from east to west in the Cyprian slab. To the south of Cyprus, the Eratosthenes Seamount collided with the island and at present appears to block the subduction process (Schattner, 2010, and references therein). Gürer et al. (2018) dates another extensional and rollback phase related to subduction in central Anatolia to about 80-43 Ma, implying that some slab fragmentations may be older than previously thought. The latest strong extensional phase (past 10-15 million years) in the Aegean (Wortel & Spakman, 2000) as well as mantle flow through the tear and around the slab (Jolivet et al., 2018) could produce strong anisotropy, which may affect imaging the velocity structures of the upper mantle. Eastern Anatolia and Arabia started colliding about 30-35 Ma (Jolivet & Faccenna, 2000); after the closure of the Bitlis suture (16 Ma, Govers & Fichtner, 2016), the slab's dipping angle steepened. The slab most likely broke off at 10-15 Ma (Jolivet et al., 2013), making asthenospheric inflow possible. Domal uplift of the Eastern Anatolian Plateau (EAP) accompanied by volcanism (e.g., Keskin, 2003) began after the slab break-off.

Anisotropy in the Region
Several SWS studies have investigated lateral variations in the direction and strength of anisotropy in various tectonic provinces of the study area, for example, eastern Turkey (e.g., Paul et al., 2014;Sandvol et al., 2003), north-central Turkey (e.g., Biryol et al., 2010;Paul et al., 2014), western Turkey (Paul et al., 2014, the Aegean Sea and Greece (e.g., Confal et al., 2016;Evangelidis et al., 2011;Olive et al., 2014;Paul et al., 2014), and along the Cyprus trench (Yolsal-Çevikbilen, 2014). Anisotropic behavior in the lithosphere has been well constrained by the anisotropic inversion of Pn traveltime residuals in the whole region (Mutlu & Karabulut, 2011) and through the analysis of the directional dependence of teleseismic receiver functions collected at the central NAFZ (Licciardi et al., 2018). Paul et al. (2014) explained the pervasive pattern of the NE-SW fast-axis orientations from the northern Aegean Sea to eastern Anatolia as the result of instantaneous density-driven mantle flow in the asthenosphere with additional local effects, such as slab rollback in the Aegean Sea and a slab window beneath southwestern Anatolia. Later, Confal et al. (2018) were able to simulate such regional coherency in splitting measurements using a 3D petrological-thermo-mechanical model, where a significant S-N asthenospheric mantle flow resulted from various effects including, primarily, the Nubian-Eurasian plate convergence with slab rollback in the Aegean Sea, a tear in the African slab, and detachment occurring within the Arabian plate (break-off) ( Figure 1). Observed splitting time delays of up to 2 s suggest that this primarily sublithospheric anisotropy might significantly alter isotropic velocity models at these depths. In the Aegean back-arc region, FPDs are predominantly perpendicular to the trench and TDs are large (1.5 ± 0.4 s Confal et al., 2016;Paul et al., 2014). In Anatolia, FPDs are NE-SW oriented, while in southwestern Anatolia and the Peloponnese, the pattern becomes more complex because of suspected slab tearing (Paul et al., 2014) and local changes due to return flow (Olive et al., 2014), respectively.
Numerical 3D synthetic anisotropy calculations of the Eastern Mediterranean and Anatolia by Confal et al. (2018) show an overall consistency with existing shear wave splitting measurements in the Eastern Mediterranean region. Near the subducting slab, however, they also suggest strong vertically directed mantle flow; thus, vertical anisotropy is likely to be present. In section 4, we explore the effect of such an anisotropy field on the tomographic imaging.
Recently, Wei et al. (2019) conducted an anisotropic inversion of P-wave traveltime data based on global catalogs from beneath the Eastern Mediterranean and Middle East. They later synthetically calculated path-integrated SKS splitting parameters based on the vertically stratified FPDs of P-waves resolved from the inversion. A comparison between their synthetic SKS splitting parameters and those observed from previous SWS studies revealed similar patterns; for example, NE-SW-oriented FPDs in most of Anatolia, a circular pattern around the tear in southwestern Anatolia, and trench-parallel FPDs in the fore-arc and subslab regions. It is only in the central-south Aegean that the synthetically estimated TDs are very small, a result that contradicts to the seismic observations reported in Paul et al. (2014) and Confal et al. (2016), and the numerical models of Confal et al. (2018). Such inconsistencies may stem either from the damping and smoothing parameters used in the regularization of the inverse problem, or the insufficient amount of P-waves or SKS-phase data recorded in those regions.

Data
In this study, we invert relative P-wave traveltime residuals to obtain a model of 3D P-wave velocity perturbations in the upper mantle. Measurements were made on records from 686 broadband and short-period seismic stations located between 20-48°EW and 33-43°NS in a region covering the Aegean, Anatolia, Greece, and Georgia ( Figure 1). The seismic stations belong to several seismographic networks with digital seismic waveform recordings extracted for the years 2005-2010 and 2013-2015. Most stations utilized in this study are operated by the Regional Earthquake-Tsunami Monitoring Center (KOERI-RETMC). We also used data recorded by several regional seismic stations operated by the Greek National Observatory of Athens (NOA-IG). Furthermore, three-component data from 72 broadband seismic stations in central Anatolia, which operated between 2013 and 2015 as part of the Continental Dynamics-Central Anatolian Tectonics project (CD- CAT Sandvol, 2013), were also used.
We selected 1,135 teleseismic events with magnitudes Mw >5.0 and epicentral distances between 30°and 90°. Following a visual inspection process, we consider the P-wave arrival times for 935 earthquakes for further analysis, since some of the observed data had unacceptably low signal-to-noise ratio. To obtain relative arrival times and uncertainty estimates, we used multichannel cross-correlation (VanDecar & Crosson, 1990) with three center frequencies, 1 Hz (22.9% of data), 0.5 Hz (32.6%), and 0.3 Hz (44.6%). This enabled us to obtain up to three sets of relative delays for each event. After cross-correlating the data, 557 good events producing 107,283 frequency-dependent delays were used for the inversion ( Figure S1 in the supporting information). The backazimuthal distribution is somewhat uneven (with more events between 5°and 110°t han between 111°and 4°) due to epicentral distance limitations, yet events of nearly all backazimuths are represented in this study ( Figure S1).

Method
We use the hybrid ray-tracing method of Bezada et al. (2013) for isotropic delay calculations. This method combines velocity-dependent iterative ray tracing and approximate finite-frequency Born kernels (Bezada et al., 2013;Schmandt & Humphreys, 2010). Traveltimes outside the model volume are calculated using a 1D model, whereas the times inside the box are calculated using a 3D velocity model applying the graph-theory-based method of Toomey et al. (1994). As a background model, for the initial ray paths and the calculation of sensitivity kernels, we employ the AK135 1D velocity model of Kennett et al. (1995). For the inversion, we use a study area approximately 2,000 km (E-W) by 500 km (N-S), with a 30 to 50 km vertical (increasing with depth) and 42 to 56 km horizontal (increasing toward the edges of the model) node spacing, resulting in 69×39×24 (x y z) nodes (64,584 nodes in total, Figure 3c). Our model extends from 750 km depth to the Earth's surface. To overcome the limited resolution and contamination of the model resulting from unimaged anomalies in the lithosphere above 50 km, the surface wave velocity model of Delph et al. (2015) is used as a constraint for about 80% of our study area. The S wave velocities are converted to P-wave velocities with an average V P /V S ratio of 1.8, adopted from a recent receiver function study by Schiffer et al. (2019). The compiled data are implemented into our starting model to represent the crustal structure prior to the inversion. High damping values prevent model velocities in the crustal block from significantly changing in subsequent iterations (Bezada et al., 2013). We prefer this approach over static corrections derived from crustal thickness because it accounts for the effect of crustal velocity variations on ray tracing and includes the effect of lateral variations in sedimentary basin thickness. The crustal structure in our final model closely resembles the model of Delph et al. (2015) but with slightly poorer resolution. In regions where we do not have crustal structure information, no prior constraints are placed on the inversion. In the present work, we will discuss only the structures and anomalies deeper than 60 km, since the focus of this study is on upper mantle structures and anisotropy. In order to regularize the inversion, constraints on the model norm and model roughness are imposed on the model. After the first iteration, the delays are recalculated in two additional iterations, sufficient for the model to stabilize. In addition to P-wave velocity model parameters, we inverted for station and event parameters.
Damping and smoothing parameters for the inversion were carefully selected and applied in such a way that the norm and roughness of the total model (rather than the updates to the model) are minimized in each iteration. We calculated the variance reduction and L2 norm for damping values of 1 to 11 and smoothing values of 1 to 14 resulting in 154 different models. We selected the best parameters using an L-curve for the first iteration (see Figure S4 and Table 1). Additionally, we looked at the misfit with norm and mean curves of the relative residuals ( Figure S3). We calculated hit quality maps for each depth layer ( Figure  S5), to evaluate the quality of the sampling of the model space. Our hit quality metric is based on the number of rays in a node and takes the backazimuthal distribution into account (Schmandt & Humphreys, 2010). The hit quality is good (>0.4) in nearly the whole study area in depths ranging from 100 to 300 km. It is noted that underneath two offshore regions in the Aegean that have small station coverage is the hit quality lower at about 0.3-0.4. Depths greater than 300 km and regions toward the edges of the model lose resolution. The resolution of the checkerboard test is very good in the upper 300 km ( Figures S6 and S7). Toward the edges of the model, and with depth, some smearing and amplitude loss of perturbations is visible ( Figure S6). In the eastern part of the model, data resolution seems to be problematic, resulting in smearing and difficulties in resolving small structures ( Figure S7).

P-Wave Tomography Corrected for Anisotropy
To include anisotropy in our velocity perturbation calculations, we used an approach that has been described and successfully tested on a synthetic data set in Bezada et al. (2016). Delay calculations are processed similarly to fully isotropic cases (see section 3), except including prescribed a priori anisotropy. The assumed anisotropy field, inferred from various sources (e.g., numerical model or shear wave splitting measurements), is interpolated into the ray-tracing grid and included in the forward traveltime calculations. As a result, each delay has a unique correction that depends on its corresponding ray path. The 3D ray-tracing method is limited to anisotropy with a hexagonal symmetry for simplicity (Toomey et al., 1994); therefore, the elastic tensors from the numerical model are approximated to the best-fitting hexagonal symmetry tensors with D-REX (Browaeys & Chevrot, 2004;Kaminski et al., 2004). In the first iteration, the prescribed anisotropy field is included in the reference velocity model and the input delays are recalculated. In principle, this removes the effect of the assumed anisotropy from the delays. For every following iteration, delays are recalculated by considering the updated isotropic structure from the previous iteration and the anisotropy field. We define a priori anisotropy fields for two different cases ( Figure 3): 1. Numerical model: Anisotropy parameters based on the numerical model of the region by Confal et al. (2018) are introduced to the starting model prior to the tomographic inversion. The strength of seismic anisotropy and its direction are calculated (D- Rex Faccenda & Capitanio, 2013;Kaminski et al., 2004) using the mantle flow and temperature yielded from the numerical reconstruction of the tectonic evolution of the study area over the past 22 million years. When including the 3D anisotropy of the numerical model in the inversion, some geometric adjustments were made to fit the real geometry of the region. The anisotropy taken directly from the numerical model (Confal et al., 2018) resulted in calculated SKS splits that were larger than the observed splits by a factor of ∼2, most notably in the Aegean region; thus, the anisotropy fractions were halved. The modeled half-arc region of the Aegean was flipped to complete the actual full arc of the Hellenic trench. To fit the geometry of the slab, the model was sheared with depth and rotated toward the east by 5°. The EW direction was scaled to 65% of its original size and the NS to 35%, respectively ( Figure 3a). 2. Observational constraints (SKS-wave measurements): We assume that the shear wave splitting is the result of a single layer of anisotropy with a laterally variable but vertically homogeneous fast polarization direction and anisotropy strength (Figure 3b). At each point, the fast direction of propagation is the same as the average FPD observed directly above, and the percentage of anisotropy is such that integrated over the thickness of the layer (placed at 90 to 350 km depth), it results in the observed TDs. The SWS parameters of Paul et al. (2014) were used since they cover the entire study area.

P-Wave Tomography With Isotropic Starting Model
Our preferred model uses damping (4) and smoothing (8) values determined through an L-curve analysis and achieves a variance reduction of 74.68% and L2 norm of 3.33 in the first iteration of the isotropic model; the same damping and smoothing values are used for the inversions corrected for anisotropy.
In western Greece and the southern Aegean (between ∼20°and 28°longitude), a high-velocity semicircularshaped structure that follows the Hellenic subduction trench is visible (Marker 1, Figures 4a, 5a, 7a1, and S9 at 60-690 km). In the west, north of the KTF, inside the high-velocity structure, an embedded low-velocity structure appears to exist parallel to the trench (NW-SE directed) present above 160-km depth (Marker 2, Figure 4a). In the southern Aegean, the high-velocity structure is approximately 200-300 km thick. Below 90km depth, beneath Crete, it dips north with a 60°angle until it flattens (10-20°dip angle) after the ∼410 km discontinuity (Marker 1, Figures 7a1 and S10). In the east, the high-velocity structure is not connected to the Cyprian slab at shallower depths (Marker 3, Figure 4a). First-order low-velocity structures (about 50-100 km wide and 100-200 km long) at 60-to 450 km depths occur in the model (Marker 3/4, Figures 4a, 5a, and 7b1). They appear to reside inside or in close proximity to a large deep gap in the high-velocity structure. The high-velocity anomaly underneath Cyprus and western Anatolia dips more steeply and appears to be fragmented. A western Cyprian slab (WCS) fragment is thin at shallow depths and becomes a thick bulge at 190 to 230 km depth (Marker 5, Figure 4a). An eastern Cyprian slab (ECS, Marker 6, Figure 4a) starts south of Cyprus but connects with the WCS at about 250-km depth (Marker 7 in Figure 4a). Toward the east, beneath the Bitlis Zagros suture zone, some high-velocity fragments 100-200 km wide at 200 to 400 km depth can be seen (Marker 8, Figures 5a and 7d1). Beneath the central Aegean, the uppermost mantle is mostly free of intense anomalies. There is a high-velocity anomaly at depth shallower than ∼200 km beneath the region of Istanbul and northern Anatolia (Marker 9, Figures 4a and 7c1). Farther south, in central and eastern Anatolia, three large and deep low-velocity structures can be found down to a depth of 300 km, separated in the deeper parts but connected in the upper 50 km of the model (Marker 10, Figures 4a and 7c1/d1).

P-Wave Tomography With an Anisotropic Starting Model
When performing the inversion with a priori constraints on anisotropy, most first-order velocity perturbations (i.e., subducting and fragmented slab, tear and break-off in the slab, and significant low-velocity zones in eastern Anatolia) are similar to those resolved after isotropic inversions. Nevertheless, some important features have significant differences in the shape and magnitude of the velocity variations (Figures 4-7, marked with a green circle and an uppercase letter).
A comparison of synthetic tests (see Figure S8) suggests that even minor features can be resolved from all of our inversions and that the effect of anisotropy on the delay times surpasses the impact of presumable noise, especially with a standard deviation of 0.1 s, which can be regarded as a reasonable range for data error (e.g., Timkó et al., 2019). When introducing noise with a standard deviation of 0.2 s, the range of velocity perturbation differences, compared to the purely isotropic case, becomes similar to the case obtained when anisotropy was added. Our synthetic tests imply that anisotropy-induced artifacts are most prominent and distinct with respect to the random noise effect mainly around the subduction zone, where mantle flow is likely to be more vigorous.

Main Differences of the Tomography Including Anisotropy From a Numerical Model
When comparing the isotropic models, derived with and without using a priori constraints on anisotropy from the numerical model of Confal et al. (2018) (Figure 6), the maximum discrepancies for velocity perturbations are −5.8% (the isotropic model is slower) and 5.2% (the isotropic model is faster). More than 52% of all nodes are faster in the isotropic inversion, while nearly 48% are slower. The average of the differences is the same with 0.16% for high-velocity and low-velocity. Only 3.23% of the model is more than 1% different (1.57% isotropic faster, 1.66% isotropic slower). Most first-order structures are in the same locations but may differ in their geometry.
In anisotropically corrected images, the low-velocity structure below western Greece seems to be stronger and larger with respect to the one resolved in the isotropic model. To a depth of 230 km, the high-velocity anomaly in the southern Aegean becomes thinner or splits in two parts (Marker A, Figure 4b). The subslab low-velocity anomaly is slightly stronger when we introduce numerically modeled anisotropy  Figures 5b and 7a2). In northeastern Turkey, low-and high-velocity perturbations exhibit prominent variations (±2-3%) without a clear pattern (Marker C, Figures 4b, 6a, and 7d2). The low-velocity anomalies between the Hellenic and Cyprus high-velocity perturbations differ in their shape (Marker D, Figure 4b, 5b, and 6a). South of Cyprus, there is an additional small low-velocity zone to the east of the high-velocity zone (Marker E, Figure 4b), which is dismissively small in the isotropic model. Low-velocity anomalies underneath the Central and Eastern Anatolian Plateau (CAP and EAP) appear to be 1-3% more intense between 230 and 270 km depth in the inversion that includes anisotropy from the numerical model (Marker F, Figure 6a).

Main Differences of the Tomography Including Anisotropy From SKS Measurements
In addition to using a starting model with numerically derived estimates of anisotropy, we also invert our observed traveltime residuals with a second anisotropic starting model, which is based on the available SKS splitting measurements for the entire region reported in Paul et al. (2014). A comparison between uncorrected isotropic tomography images and those corrected for anisotropy using SKS splitting measurements ( Figure 6) displays values ranging from a minimum of −5.44% (isotropic slower) to a maximum of 4.54% (isotropic faster). The fraction of slower nodes in the isotropic version is 47% (difference average ¼ −0.18%) while fast nodes represent 53% (difference average ¼ 0.17%). Of the nodes, 3.59% are more than 1% different (1.83% isotropic faster, 1.76% isotropic slower). The low-velocity structure underneath western Greece appears to be even stronger and larger than in the tomography that includes anisotropy from the numerical model (Marker G, Figures 4c and 6b). Starting at a depth of 125 km, the high-velocity anomaly in the southern Aegean appears to be thicker, and there is no gap visible in the slab as seen in the model corrected for numerical anisotropy (Marker H, Figures 4c and 6b). In the northern Aegean and the rollback area, the isotropic tomography is about 1-2% slower at depths down to 230 km (Marker I, Figures 4c, 6b,

Journal of Geophysical Research: Solid Earth
and 7a3). A low-velocity anomaly present south of Istanbul appears to be stronger at 90 km depth (Marker J, Figure 4c), and in the upper 200 km, another low-velocity anomaly in western Anatolia is weaker than in the model corrected for numerically calculated anisotropy (Marker K, Figure 4c). High-and low-velocity perturbations in northeastern Anatolia differ by ±2-3% from the isotropic model as in the model with numerical anisotropy, but in different locations (Marker C, Figures 4c, 6b, and 7d3).

Discussion
In this study, we showed that anisotropy can affect velocity perturbations in some areas, with local differences of up to 2% compared to an isotropic tomography approach. While most discrepancies between the models seem small if our primary target is to detect and prove the existence of a slab, others might be relevant for the interpretations of mantle heterogeneities. In cases where the exact geometry or presence of relatively small-scale anomalies, like tears and fragments of the slab, are sought and if estimating the temperature of anomalies is the goal, then small variations might be meaningful. In the following subsections, we will discuss the stable and prominent features of this P-wave tomography that are visible in all of the models and examine the differences between isotropic models and models corrected for anisotropy.

Hellenic Subduction System and Aegean
In isotropic tomography models (Biryol et al., 2011;Portner et al., 2018, and this study), the dipping high-velocity anomaly, which is commonly explained by a subducting crust and lithosphere, appears to be thicker than those expected by other methods (8-to 10-km oceanic crust, 70-to 80-km Nubian lithosphere Kind et al., 2015;Pearce et al., 2012). The smearing of high-velocity anomalies is highly possible and most likely stems from an insufficient number of crossing rays in a given cell, which could make the slab appear thicker (Bezada et al., 2016). On the other hand, the checkerboard test ( Figure S6) shows no smearing for the upper 300 km. In our tomograhic images obtained after anisotropy-correction using mantle flow modeling results, the slab appears to be more defined and some high-velocity perturbations north of the slab in the southern Aegean are reduced (Marker A, Figure 4, 90 to 160 km). Based on numerical models (e.g., Confal et al., 2018;Faccenda, 2014) we know that vertically/subvertically aligned anisotropy exists below and inside the slab, and in the fore-arc and back-arc regions above the mantle wedge. However, considering the effect of anisotropy, the geometry does not change drastically. Therefore, smearing and difficulties to resolve the half-arc geometry of the slab in tomography models might still be a problem, as discussed in Portner et al. (2018). Bezada et al. (2016) described a trench-parallel low-velocity artifact below a hypothetical slab structure and a decrease in the strength of the anomaly when anisotropy is included. Piromallo and Morelli (2003) ascribed the observed low-velocity anomaly present below the Aegean slab to a possible mantle upwelling. We observe this low-velocity anomaly in all three models (similar to Biryol et al., 2011;Portner et al., 2018) as it appears even more intense in the models corrected for anisotropy (e.g., Marker B, Figures 5b and 7a2). The literature implies the low-velocity anomaly may be due to hot uprising material or subparallel mantle flow. In a recent tomography study of Wei et al. (2019), this low-velocity anomaly is less intense and only seen close to the Cyprian trench. Very strong, trench-parallel FPDs in their study, identified as mantle flow, and a larger study area could explain the abundance of this low-velocity anomaly compared to observations in our study. VanderBeek and Faccenda (2020) found low-velocity artifacts, especially beneath the slab, to be very persistent, when performing realistic anisotropic inversions. Subslab anisotropy might be difficult to correct for, but a better representation and inclusion of stations from the south of Crete could reduce low-velocity artifacts.
In the north-western segment of the Hellenic trench, high-velocity anomalies are less strong or totally absent at shallow depth (<160 km), contrary to a P-wave tomography model that was originally developed by Amaru (2007) Özbakır et al. (2020) interpreted the tear as a semihorizontal proto-STEP (Subduction-Transform-Edge-Propagator) tear until a depth of about 130 km in their full-waveform tomography. Their findings fit very well with our results, in which we see a low-velocity structure parallel to the trench in 60 to 150 km depth. The slab fragmentation seems to have proceeded from the NW to SE ending at the KTF, which marks the border of the continental and oceanic subduction regimes. Contrary to the interpretations of some studies (Evangelidis, 2017;Guillaume et al., 2013;Jolivet et al., 2013), this study shows that this disruption in the western Hellenic slab does not seem to be a vertical tear but rather a horizontal (Hansen et al., 2019) or semihorizontal (Özbakır et al., 2020) shallow break-off. Nevertheless, we support the theory that this recent break assisted the fast retreat of the Hellenic slab in the southern Aegean.
The subduction angle of the Hellenic slab in the southern Aegean is approximately at 60°, similar to Portner et al. (2018) but steeper than what Biryol et al. (2011) suggested. Even though the resolution decreases with depth, it is sufficient to show that the dip angle decreases to about 10-20°; it appears to become stagnant ( Figures 5 and 8c), but others, for example, Portner et al. (2018), have shown the Hellenic slab to extend further into the mantle.

Slab Tear of Southwestern Anatolia and the Cyprus Subduction System
The tear beneath western Turkey is very well developed from shallow depths to nearly 600 km, similar to earlier studies (Biryol et al., 2011;Portner et al., 2018). Low-velocities within the tear zone to the north and south (Marker 4, Figures 4 and 7) can be interpreted as volatile rich mantle inflow from below the slab. The volume and the exact locations of the anomalies differ in our three models depending on the anisotropy model used for the corrections. Paul et al. (2014) attributed the NW-SE oriented FPDs inferred from SKS splitting measurements in south-western Anatolia to the effect of toroidal mantle flow. In fact, numerical models in Confal et al. (2018) have since shown mantle flow with an anticlockwise toroidal pattern moving through a tear toward the west, a result of the retreat of the Hellenic trench. Toward the east, the tear is connected to the smaller Cyprian tear (Marker 11, Figures 4a and 8b/c). As reported in Biryol et al. (2011) and Wei et al. (2019), the Cyprian slab is resolved as being fragmented in this study, but not by a purely horizontal tear as recently suggested by Portner et al. (2018). In our model, the shallow Cyprian tear is probably connected to the larger tear beneath western Turkey between 60 and 100 km depth that separates the western and eastern Cyprian slab (WCS and ECS) down to a depth of about 250 km. Our results correlate well with the locations of detected seismicity in the region ( Figure S11) that indicate no shallow events south of the Antalya basin but the presence of deep ones in the WCS. A lack of earthquakes correlates with the Paphos Transform Fault (PTF). The ECS is thin (similar to Piromallo & Morelli, 2003, at 150 km depth) but continuous, starting south of Cyprus and connecting with the WCS at a depth of about 250 km. An absence of deeper earthquakes suggests that the ECS is not moving except for in the upper 70 km ( Figure  S11). It is likely that the subduction process was disrupted about 3 Ma, when the Eratosthenes Seamount collided with the Cyprus slab (Schattner, 2010, and references therein). Portner et al. (2018) suggested buoyancy-driven tearing due to the denser oceanic plate of the WCS and the continental ESC. East of Cyprus, strike slip faults mark the end of the subduction due to continent-continent collision and a transition to slab break-off (Schattner, 2010, and references therein). Westward-directed mantle flow through the various tears resolved in this study, could facilitate the westward escape movement (Schildgen et al., 2014) of the Anatolian microplate.

North, Central, and Eastern Anatolia
Beneath the Istanbul and Pontides block, north of the NAFZ, a prominent shallow (depth < 150 km) high-velocity structure (Marker 9/IST, Figure 4a, 7c1, and 8), similar to those seen in the models of Biryol et al. (2011) and Portner et al. (2018), is clearly visible and could be related to Neotethyan sutures (e.g., Biryol et al., 2011;Okay & Tüysüz, 1999;Salaün et al., 2012). Similar to Portner et al. (2018), we see a relatively sharp border from high averaged velocity perturbations in the west to low averaged velocity perturbations in the east in central Anatolia ( Figure S2)

Journal of Geophysical Research: Solid Earth
and this study) possibly facilitated the upwelling of hot asthenospheric material and active volcanism. The recent uplift in this region could be a consequence of the rebound of the subducting lithosphere after the shallow break-off (Delph et al., 2017;Schildgen et al., 2014), equally asthenospheric inflow beneath central Anatolia could also be responsible for the uplift (e.g., Cosentino et al., 2012;Govers & Fichtner, 2016). Velocity anomalies in our study suggest that the Cyprian slab has not totally broken-off beneath central Anatolia, which is in contrast to the latest tomography findings of Portner et al. (2018). More precisely, our models show that it is connected with the gently dipping deeper parts of the Hellenic slab ( Figure 8b). Therefore, we suggest recent uplift is not responsible for the rebound but is mostly the result of asthenospheric melt inflowing through the tear between the Cyprian slabs. The trench retreat slowed down (Schildgen et al., 2014); therefore, extension is absent in this region and thinning of the crust has not occurred.
The slab detachment beneath the Bitlis Zagros suture zone seems to be at an advanced stage, with no large continuous slab and close to no earthquake activity below 20 km ( Figure S11). In central Anatolia, the relic slab lies at about 400-km depth and is connected to the Hellenic slab. Biryol et al. (2011) suggested a stagnated slab at the mantle transition zone, while Portner et al. (2018) did not find any evidence of a slab present beneath central and eastern Anatolia. In our tomography, it appears that the Cyprian slab (which reaches 400 to 500 km depth) never sunk as deep as the Hellenic slab (until 500 to 660 km depth), although they are connected. The slab may be stuck, and this could be preventing it from sinking deeper due to the collision with the Eratosthenes Seamount (e.g., Biryol et al., 2011;Schattner, 2010). Although resolution decreases with depth and toward the east, Bitlis slab fragments beneath eastern Anatolia (roughly cubic high-velocity anomalies with side lengths of 50-100 km) can be distinguished (Figure 8c).
Given the insufficient data coverage in eastern Anatolia, our model resolution of this region might be relatively poor and may generate artificial anomalies. While the normalized hit quality values of 0.4-0.5 indicate reasonably good data resolution ( Figure S5), the checkerboard test shows smearing in this region, implying that the results may not be sufficiently robust for interpretation.

Uncorrected Isotropic Model Versus Models Corrected for Anisotropy
Previous studies have shown that whether anisotropy is neglected in tomographic inversions or some effort is made to account for its effect, the main features recovered are generally stable, with second-order differences between the models often present (e.g., Bezada et al., 2016;Eken et al., 2012;Lloyd & van der Lee, 2008;O'Driscoll et al., 2011;Sobolev et al., 1999). Even when anisotropy is explicitly solved for in the traveltime inversion, a similar scale of differences is found between isotropic and anisotropic models (e.g., Ishise & Oda, 2005;Koulakov et al., 2009;Tian & Zhao, 2013). This is generally consistent with what we observe in our study, although the differences we detect are perhaps smaller than expected.
Several factors can affect the relative performance of seismic tomography with and without consideration of anisotropy, including the azimuthal coverage of the events and the characteristics of the anisotropic structure beneath the study area (e.g., Bezada et al., 2016;Lloyd & van der Lee, 2008;Sobolev et al., 1999). Although the backazimuthal distribution of the events used in this study is uneven, our synthetic inversions ( Figures S6 and S7) exhibit a promising model resolution performance, without much smearing down to 300-400 km, except for eastern Anatolia, where the station density is relatively sparse. Earlier, Bezada et al. (2016) concluded azimuthal coverage would not have a severe impact on anisotropy-related artifacts for teleseismic data. For these reasons, it seems unlikely that azimuthal coverage is a dominant factor in the outcome of our inversions. Regarding the spatial character of anisotropy variations, Grésillaud and Cara (1996) revealed that laterally homogeneous anisotropy showing long-wavelength variations over a large study area would not bias isotropic inversions drastically. However, our study region has a long history of several subduction events and is presently subject to extensional and compressional tectonics as well as an active subducting slab. This tectonic complexity is likely to produce strong lateral variations in anisotropy with relatively short wavelengths, which could reasonably be expected to have a significant impact on teleseismic P-wave traveltime inversions. At the same time, subduction (where the lateral variations in anisotropy have the shortest wavelengths) occupies a relatively small fraction of our entire model space, which may limit the influence of anisotropy on the results. Another important factor is the dip of the fast axes of anisotropy. Sobolev et al. (1999) tested the effect of seismic anisotropy caused by various cases developed under different tectonic conditions. They predicted the highest deviations to result from isotropic models for the hypothetical scenarios with dipping axes of symmetry (e.g., a 3-5% change in amplitude) and with sublithospheric mantle flow (a 2-3% change in amplitude). VanderBeek and Faccenda (2020) in a more recent study have shown that the anisotropic bias in isotropic inversion of teleseismic P-wave delays would not be removed efficiently in cases where azimuthal anisotropy (horizontal fast axes) dominates. The anisotropic inversion of P-wave traveltime residuals can be carried out when the fully anisotropic tensor is approximated to a hexagonal symmetry, where the axis of symmetry is often assumed to be horizontal (e.g., Wei et al., 2019). A recently developed method  allows for arbitrarily oriented anisotropy in 3D. By applying it to the northern Fennoscandian shield, it was found that the best-fitting fast axes were not always horizontal . The numerical model (Confal et al., 2018) employed in estimating the anisotropic contribution to traveltime delays in this work includes dipping anisotropy above and beneath the active subducting slab, regions where we observed the most significant discrepancies between the isotropic and anisotropy-corrected models.
After correcting for the effects of anisotropy, the biggest discrepancies in our models (up to ±2%) compared to the isotropic model lie, for instance, in Greece, in the downgoing slab and in central Anatolia. Velocity anomalies in western Greece are up to 2% lower (e.g., Marker G, Figure 6a By applying anisotropy corrections derived from the numerical model of Confal et al. (2018), we were able to only marginally improve the total variance reduction of the tomography, which can be regarded as an overall measure for the goodness of data fit, by 1.6% (e.g., an increase from 74.32% to 75.9%, Table 1). Similarly, Eken et al. (2012) reported that applying numerically calculated anisotropic correction terms to the observed data for the Fennoscandian shield would produce an improvement of a few percent in the variance reduction. , despite explicitly inverting for anisotropy with an arbitrarily oriented axis of symmetry, obtained only a small improvement of up to 7% for the variance reduction compared to isotropic inversions. Results from completely synthetic tests designed to focus on subduction zones, presented in Bezada et al. (2016), suggest a much better increase in variance reductions should be possible. Counterintuitively, our synthetic analyses show that isotropic inversions in which anisotropy has ( Figure  S8b) and has not been accounted for ( Figure S8e) yield very similar variance reductions, meaning that even anisotropic traveltimes can be mapped to isotropic velocity variations very successfully. We attribute this mismatch with Bezada et al. (2016) to be primarily due to the fact that our study area is much larger and more complex than the isolated synthetic slab explored in the earlier study. Our model space covers an active subducting African Plate along the Hellenic and Cyprus trenches and the entirety of Anatolia which experiences present-day extensional and collisional tectonics as well as sublithospheric mantle flow evidenced by the azimuthal anisotropy. While lateral changes in the strength and dip of anisotropy have short wavelength changes near trenches, the extend of azimuthal anisotropy is fairly homogeneous over all of Anatolia, which may help explain why there is relatively little bias in the isotropic inversions.
Another important factor to consider is our treatment of the shallower velocity and anisotropy structure. The numerical model does not include possible "frozen-in" dipping anisotropy that may be present in the continental lithosphere given the complex deformation history of the area. Regions of high-velocity perturbations in northern Anatolia (Marker 9/IST, Figures 4a, 7c1, and 8) can be related to relatively thick and old lithosphere. Similar anomalies in the study region have also been reported in the early tomography studies of Biryol et al. (2011) and Salaün et al. (2012) and interpreted as the presence of a remnant slab underneath a highly deformed basement of the Pontides and Istanbul Zones, which are characterized by the accretion of continental terrains during the Cimmerian orogeny (Bozkurt, 2001). If any "frozen-in" fabric with a dipping axis of symmetry exists in this region, then ignoring it within our correction procedure will lead to biased inversions (Sobolev et al., 1999) and thus contribute to lowering the total variance reductions with and without corrections. In contrast, Paul et al. (2014) suggested that such complicated anisotropic structures could not be considered for the Aegean-Anatolia, which has a relatively thin mantle lithosphere, 10.1029/2019JB018559

Journal of Geophysical Research: Solid Earth
given the regionally coherent pattern in the lateral variation of SKS splitting parameter observations. Previously, Plomerová et al. (2011) and Eken et al. (2010) showed that the actual orientation of complicated anisotropic structures beneath the Fennoscandian shield would require a joint inversion of multiple data sets, that is, P-wave traveltime residuals and SKS splitting parameters. Thus, more sophisticated modeling studies are further needed to clarify whether a subregional dipping anisotropy is present in this particular region. Crustal anisotropy that can be produced by the layering of sediments, oriented cracks of variable length and width, or the foliation of rock complexes (Babuska & Cara, 1991) is neglected in our model since previous studies in the region based on local shear wave splitting measurements (e.g., Eken et al., 2013;Hurd & Bohnhoff, 2012;Peng & Ben-Zion, 2005) or receiver function analyses/models (e.g., Licciardi et al., 2018;Vinnik et al., 2016) have indicated its relatively insignificant contribution compared to the predominant mantle anisotropy. If strong and laterally variable crustal anisotropy is present in the study area, it would be an additional source of variance that neither of our models would be able to explain. Additionally, although we make efforts to account for isotropic crustal structure, the model we use for this purpose is imperfect. Any potential bias that is introduced by these a priori constraints on the crustal part of the starting model will be identical in both uncorrected and anisotropy-corrected models and equally contribute to the postinversion variance in both cases. We note that none of these crustal effects were factors considered in the synthetic study of Bezada et al. (2016) which may also help explain why our anisotropic corrections do not achieve the improvement in variance reduction that is suggested by that study.
Our findings suggest that while upper mantle anisotropy in tectonically active regions such as subduction zones could bias P-wave inversions, the effect is relatively small and comes into play only when interpreting second-order features and the amplitude of specific anomalies. Although these may seem like details in the context of the broader model, they may be important components of the overall interpretation of results, especially if they relate to slab tears and other spatially small features. On the other hand, the fact that changes in the model due to anisotropy are not very significant adds further credibility to the interpretation of the features resolved by this study and previous studies based exclusively on isotropic tomography (e.g., Biryol et al., 2011;Portner et al., 2018).

Conclusion
This study contributes to the knowledge of the active tectonic evolution of the Eastern Mediterranean region with a very high-resolution 3D velocity model of the upper mantle as well as to the discussion on the effect of anisotropy on isotropic tomography models. Compared to recent isotropic P-wave tomography models, our model benefits from improved station coverage in the western Aegean, and we present new findings, especially about western Greece, the Aegean, and Cyprus. Our tomography results suggest the presence of a horizontal tear in western Greece as part of the northern Hellenic slab segment. A very deep and pronounced vertical tear in western Anatolia and a subhorizontal tear between the western and eastern part of the Cyprian slabs can be observed in our results. The Cyprian tear reaches about 200 km depths between the thin eastern and deeper western Cyprian slab. The dip of the slab in the Aegean has flattened from 410 to 660 km depth, while in central Anatolia, the slab appears to not reach deeper than 500 km.
Including anisotropy, especially with a significant component of plunging axes of symmetry, in anisotropy-corrected tomography is very challenging, since anisotropy parameters are derived from numerical models that are only approximations to true mantle structure. The configuration of these geodynamic numerical models is based on prior first-order interpretations of tomography models themselves. In general, corrected and uncorrected tomographic images indicate similar features that are very robust and can be therefore interpreted with a high degree of certainty, while some smaller features are less stable. Even though the tomography results change weakly by including anisotropy, the changes we see highlight the fact that caution should be taken when using these images to interpret the physical state of the upper mantle or linking them to active tectonics such as volcanism (O'Driscoll et al., 2011) or active subducting slabs (Sobolev et al., 1999). The discrepancies we observe between isotropic and anisotropy-corrected tomography models are similar to those between the results of different tomography studies previously conducted in the region and depend on model parametrization. Our findings show that the influence of anisotropy in the study region should be taken into account, especially when interpreting small-scale features but that first-order features observed in isotropic models are robust.

Journal of Geophysical Research: Solid Earth Data Availability Statement
The IRIS Data Management Center (https://www.iris.edu/hq/) was used to access seismic waveforms and the United States Geological Survey (USGS) to obtain earthquake hypocentral details respectively. This work includes data from various permanent and temporary regional networks ( Digital waveform recordings were retrieved from EIDA, an initiative within ORFEUS (hhtp://www.orfeus-eu.org/data/eida/networks/). Data from prior SKS splitting measurements were taken from Wüstefeld and Bokelmann (2007). The figures are made using the Generic Mapping Tools (Wessel & Smith, 1998), GeoMapApp (http://www.geomapapp. org), Matlab, and Inkscape softwares. The tomography results can be downloaded on OSF under a project with the same title as this paper (https://osf.io/bu46n/).