Tomographic Image Interpretation and Central-Western Mediterranean-Like Upper Mantle Dynamics From Coupled Seismological and Geodynamic Modeling Approach

The Central-Western Mediterranean (CWM) is one of the most complex tectonic setting on Earth. Episodes of slab rollback, break-off and tearing, the opening of back-arc extensional basins (i.e., Liguro-Provencal, Alborean, Algerian and Tyrrhenian basins), the presence of large mountain ranges, active volcanoes and violent earthquakes have made the Mediterranean an ideal environment to study a wide range of geodynamic processes and an important target for seismological studies (e.g, seismic tomography). Here we build a geodynamic model which, although it does not reproduce its exact tectonic structure (e.g., due to the limits of the numerical method, approximations in the initial setup, etc), presents multiple and geometrically complex subduction systems analogous to those found in the CWM. The tectonic evolution of this model is estimated with petrological-thermo-mechanical 3D simulations, then, we dynamically compute the upper mantle fabrics and seismic anisotropy as a function of the strain history and local P-T conditions. After comparing the model with SKS splitting observations in order to quantify the discrepancies with the true Central-Western Mediterranean, we use the elastic tensors predicted for the modeled configuration to perform 3D P-wave anisotropic tomography by inverting synthetic P-wave delay times. Using the geodynamic model as reference, we evaluate the capabilities of a recently developed seismic tomography technique to recover the isotropic anomalies and anisotropy patterns related to a complex subduction environment in different conditions, such as poor data coverage and bad data quality. We observe that, although P-wave tomography still remains a powerful tool to investigate the upper mantle, the reliability of the retrieved structures strongly depends on data quality and data density. Furthermore, the recovered anisotropic patterns are consistent with those of the target model, but in general an underestimation of the anisotropy magnitude in the upper mantle is observed. In the light of future developments, our study suggests that by combining micro- and macro-scale geodynamic simulations and seismological modeling of seismic anisotropy it will be possible to reproduce, at least to a first order, the tectonic evolution of real study regions (e.g., the Mediterranean) thus providing fundamental constraints on the processes that have contributed in shaping their current geological scenario.


INTRODUCTION
Since the early 1990s numerous seismological studies have been carried out to image the Earth's upper mantle and seismic tomography proved to be a fundamental tool for constraining the past and present-day mantle dynamics and structure (Liu and Gu, 2012;Rawlinson et al., 2014;Van der Meer et al., 2018;Romanowicz, 2021). Tomographic methods (e.g. P-, S-and surface-wave tomography) yield wave velocity models that are commonly used to infer distributions in physical and chemical properties affecting seismic-wave propagation such as density, temperature, melt fraction and volatile content.
At the same time, petrophysical analysis of exhumed rock samples and micromechanical laboratory experiments (Ribe, 1989;Savage, 1999;Blackman and Kendall, 2002;Kaminski et al., 2004;Karato et al., 2008;Karato et al., 2008;Long and Becker, 2010;Faccenda, 2014;Skemer and Hansen, 2016) have shown that the development of mineral and compositional fabrics mainly associated with rock deformation can create significant directional variations in seismic velocities known as seismic anisotropy. Although the presence of seismic anisotropy in Earth's upper mantle is well-established, scientists have often assumed the Earth's interior as seismically isotropic. This approximation certainly simplifies the computational approach but at the same time it can introduce notable imaging artefacts and, consequently, errors in the interpretation of the tomographic results (Kendall, 1994;Blackman et al., 1996;Blackman and Kendall, 1997;Sobolev et al., 1999;Lloyd and Van Der Lee, 2008;Menke, 2015;Bezada et al., 2016;VanderBeek and Faccenda, 2021).
Recently, VanderBeek and Faccenda (2021) and Wang and Zhao (2021), have independently developed a methodology to invert for P-wave isotropic (mean velocity) and anisotropic (magnitude of hexagonal anisotropy, azimuth and dip of the symmetry axis) parameters. When tested on a relatively simple, 3D geodynamic model of subduction, VanderBeek and Faccenda (2021) found that the new inversion technique produces a much more accurate reconstruction of the upper mantle isotropic and anisotropic structures. In contrast, ignoring for seismic anisotropy (isotropic approximation) or allowing for only azimuthal variations in seismic velocity (i.e., no dipping fabrics) generates strong imaging artifacts. From these tests it follows that taking into account seismic anisotropy can provide new insights into the 3D upper mantle structure and dynamics. Despite these encouraging results, it remains unclear whether isotropic and anisotropic structures of the Earth's mantle can be simultaneously recovered by P-wave anisotropic inversions in real and more complex tectonic settings.
Along with seismic imaging techniques, over the last decades numerical geodynamic modelling became an essential approach for understanding the long-term and deep evolution of a wide range of geological processes, which otherwise would remain unconstrained due to the lack of geological data (Gerya, 2019). Owing to the development of increasingly high performance computers and more advanced numerical techniques, it is nowadays possible to simulate the multiscale tectonic evolution of 3D complex settings for 10s or 100s of millions of years (van Zelst et al., 2021). However, despite being a powerful tool, numerical modeling is also affected by several limitations that could potentially bias the final output, such as uncertainties in the employed initial model geometry, physical parameters (mainly viscosity), chemical compositions, and limited computational power.
In order to test the limitations of the tomographic and numerical modeling methods, a promising approach is combining micro-and macro-scale geodynamic modeling simulations of mantle flow to predict mantle isotropic and anisotropic structures and then perform seismological synthetics (Faccenda andCapitanio, 2012, 2013;Hu et al., 2017;Confal et al., 2018;Zhou et al., 2018;Lo Bue et al., 2021). We decided to apply this combined methodology to the Central-Western Mediterranean (CWM) region (Supplementary Figure S1). In the last 20-30 million years, this area has experienced a complex tectonic activity characterized by backarc extension related to slab retreat in the Liguro-Provençal, Alborean, Algerian and Tyrrhenian basins and episodes of slab break-off, lateral tearing and interactions between slabs Platt and Vissers, 1989;Lonergan and White, 1997;Carminati et al., 1998;Wortel and Spakman, 2000;Rosenbaum et al., 2002;Faccenna et al., 2004;Mauffret et al., 2004;Spakman and Wortel, 2004;Jolivet et al., 2006;Faccenna et al., 2007;Jolivet et al., 2008;Vignaroli et al., 2008;Jolivet et al., 2009;Carminati et al., 2012;Faccenna et al., 2014;van Hinsbergen et al., 2014;Király et al., 2018;Van Hinsbergen et al., 2020) and a wealth of geological and geophysical data are available. Numerous tomographic models and geodynamic studies focusing on the CWM upper mantle are available, which can be used here to test the reliability of our approach.
We first extend the modeling methodology of Lo Bue et al. (2021) to create a geodynamic model that resemble observed slabs morphology and anisotropic mantle fabrics of the CWM. The geodynamic model is then exploited as synthetic case study to test the capabilities and limitations of P-waves isotropic and anisotropic inversions in recovering complex geological scenarios using the methodology of VanderBeek and Faccenda (2021).
In this work, we attempt to answer some fundamental questions. How well does P-wave anisotropic tomography recover the modeled isotropic and anisotropic structures? How reliable are the inferred anisotropic patterns with respect to the upper mantle fabrics? Which are the main artefacts introduced in the tomographic image when neglecting seismic anisotropy? To which extent vertical smearing bias the inverted structures when only using teleseismic P-wave travel times?

Geodynamic Numerical Modelling
We construct a 3D petrological-thermo-mechanical numerical model of the Central-Western Mediterranean convergent margin using I3MG (Gerya, 2019), which is based on the finite difference method (FDM) combined with a marker-in-cell (MIC) technique. The mass, momentum and energy conservation equations are solved on a staggered Eulerian grid while the physical properties are interpolated to the Lagrangian markers for advection. The Earth's mantle is treated as a highly viscous incompressible medium. Visco-plastic deformation is simulated by combining a Drucker-Prager yielding criterion with dislocation, diffusion and Peierls creep mechanisms.
In this paper we refer to our geodynamic model as Model CWM (Central-Western Mediterranean Model). This model is an updated version of the Reference Model CM of Lo Bue et al. (2021). Here, the computational domain has been enlarged and has dimensions of 3700 x 700 x 2200 km (373 x 101 x 229 nodes) along the x − y − z coordinates, with y being the vertical direction. As in Model CM, subduction modeling is self-consistent, driven only by internal buoyancy forces. Velocity boundary conditions are free slip everywhere. We impose a constant incoming heat flux of 2 mW/m 2 at the bottom boundary, while the top boundary is characterized by a constant temperature of 273 K. The side boundaries are insulating. The models account for frictional and adiabatic heating, and for thermal and dynamic effects of phase changes (except that the medium is assumed to be incompressible).
We used the MATLAB toolbox geomIO (Bauville and Baumann, 2019) to create the 3D initial temperature and compositional fields. The tectonic plate geometry has been designed according to the paleogeographic and tectonic reconstructions at~30 Ma proposed by Lucente and Speranza (2001); Lucente et al. (2006);Faccenna et al. (2014);van Hinsbergen et al. (2014); Romagny et al. (2020) although some simplifications were required due to limitations imposed by numerical modeling.
In the initial setup ( Figure 1), a subducting oceanic plate, that represents the Ionian Ocean, is surrounded by lateral continental blocks corresponding to the Adria, Africa, Iberia and European plates. The position of the plates in the Oligocene-Miocene was adapted from a reconstruction of van Hinsbergen et al. (2014). It is worth noting that, not having applied a convergence rate between the plates (self-consistent subduction), their relative position slightly differs from the present-day one. However, an initial geometry defined in the Oligocene-Miocene should not have a strong impact on mantle flow directions and splitting parameters as the slow Africa-Europe plates convergence has not caused a drastic change in plates position over this time span.
In Model CWM, we considered a more realistic paleo-tectonic configuration of the region, which is characterized by the presence of multiple subducting slabs rather than a single one as in Model CM.
Subduction in the Ionian plate occurs along two trenches as in Romagny et al. (2020). The longest one stretches from the Alps to the southeast of the Baleares and is associated with a slab dipping 40°NW and extending down to 300 km in the upper mantle. A second one is placed in the Alboran domain, where a slab with the same dipping angle but an opposite vergence extends down to 350 km in the upper mantle (Romagny et al., 2020). Throughout the manuscript we use "Ionian slab and Ionian trench" to indicate the former subduction zone and "Alboran slab and Alboran trench" when referring to the latter. It is worth noting that to trigger slab roll-back self-consistently, the Ionian trench has been positioned further south as in Lo Bue et al. (2021) and the initial depth of the Alboran and Ionian slabs has been increased compared to tectonic reconstructions. This could cause a difference in rates of slabs retreat when compared to those reported in the literature and be representative of a more recent stage of the Central-Western Mediterranean history rather than the 30 Ma assumed here.
Two large collisional suture zones, are present in the continental Adria and European plates as in Faccenna et al. (2014). To the north, we find the Alpine trench with its characteristic arcuate shape and, to the east of the model, the Dinaric-Hellenic trench that extends from Eastern Alps to the southernmost tip of the Hellenic peninsula . In both trenches the slab dips almost vertically into the upper mantle to a depth of about 350 km to simulate locked collision zones. In this area slabs extended down to 350 km depth to model flow barriers due to the presence of subducted slab.
In the area of the model corresponding to the present-day Apenninic chain, we use the same initial configuration as in Lo Bue et al. (2021), characterized by lithospheric heterogeneities which are fundamental for the development of key tectonic features such as a prolonged eastward retreat of the Ionian plate and the formation of a slab window below the modelled central Apennines. The Adria plate consists of a thin continental lithosphere in the Umbria-Marche region and of a stiffer continental promontory in its central portion corresponding to the Abruzzo-Laziale platform (Calcagnile and Panza, 1980;Geiss, 1987;Lucente and Speranza, 2001;Panza et al., 2003;Lucente et al., 2006;Miller and Piana Agostinetti, 2012;Maino et al., 2013;Lo Bue et al., 2021). The African plate structure, nearby Sicily and the Sicilian Channel area, is characterized by a slightly thinner margin (Arab et al., 2020;Lo Bue et al., 2021).
The initial lithosphere thermal structure was modeled using the half-space cooling equation (Turcotte and Schubert, 2014), while the underlying asthenosphere consists of a 0.5 K/km constant adiabatic temperature gradient. The thermal age of the Ionian oceanic plate is 80 Myr, while that of the two slabs is 70 Myr to simulate partial heating upon subduction. The age of the continental plates (Africa, Africa eastern margin, Iberia, Adria and Adria promontory) is 150 Myr while an age of 90 Myr was imposed for the thinned portion of Adria continental lithosphere. To activate a self-consistent subduction, the Ionian plate north of the two trenches is composed of a young lithospheric portion (1 Myr -the young age is justified by assuming a well-developed continental rifting system North of the Balearics and Corsica-Sardinia block). Furthermore, rheologically weak zones (constant viscosity of 10 18 Pa s and constant density of 3200 kg/m 3 ) have been inserted 1) on the slabs top surface to lubricate the initial contact between the overriding and the subducting plates, and 2) around southwest Iberia and northwest Africa to facilitate the Alboran trench retreat (e.g., Chertova et al. (2014)). The plates thermal structures and flow law parameters have been tuned to allow a self-consistently subduction and simultaneously to reproduce the main tectonic events as close as possible to the geological reconstructions. This may cause a too weak rheology and faster rates of mantle convection once self-sustained subduction has started due to the non-linear viscous behaviour of the mantle.
The density is computed using the thermodynamic databases generated with PERPLE_X (Connolly, 2005) and tested by Mishin et al. (2008) for a pyrolytic mantle composition. The continental crust density is calculated as being that of the mantle minus 400 kg/m 3 , except for the Adria thin margin where we subtract 200 kg/m 3 to model a less buoyant continental lithosphere. Instead, for the crust of the Adria promontory we use a constant value of 2700 kg/m 3 . More details about the physical parameters used in the geodynamic model can be found in Lo Bue et al. (2021).

Predicting Mantle Anisotropy and SKS Splitting
The development of seismic anisotropy in the upper mantle is calculated using a modified version of D-Rex (Kaminski et al., 2004), which incorporates the deformation mechanisms inducing LPO (plastic deformation, dynamic recrystallization and grainboundary sliding) and accounts for deformation history and nonsteady-state evolution of geodynamic systems (Faccenda and Capitanio, 2013;Faccenda, 2014).
A large number of Lagrangian particles representing mineral aggregates are regularly distributed throughout the numerical domain (25 km reciprocal distance along the three directions, for a total of 364,672 aggregates). Each particle consists of 1024 randomly oriented crystals, which results in an initially isotropic upper mantle. We use a harzburgitic upper mantle composition (70% olivine and 30% orthopyroxene modal abundance) and a more fertile pyrolitic mantle composition in the transition zone (60% spinel and 40% majoritic garnet) (Faccenda, 2014). The Eulerian velocity field obtained by the macro-flow simulation is then used to passively advect the particles and LPO is generated at each time step through the re-orientation of such particles in response to the gradients in the velocity field. Since SKS splitting parameters are mostly sensitive to the upper mantle (Sieminski et al., 2008), we only model the anisotropy from the Moho to the 410 km discontinuity. We use the same dimensionless crystallographic parameters as in (Rappisi and Faccenda, 2019) with the nucleation rate λ = 5, the grain-boundary-mobility M = 1 and the threshold volume fraction χ = 0.9, which generate weaker fabrics and seismic anisotropy more consistent with the observations.
Synthetic SKS splitting parameters are computed using the software package FSTRACK (Becker, 2006). Through the stiffness matrix the code recovers the elastic tensors for each aggregate and then, below each station and down to 400 km, it builds a vertical stack of horizontal layers (minimum thickness of 25 km) where the elastic tensor of each layer is radially averaged within a distance of 50 km. Next, assuming an incident plane wave (5°for typical SKS arrivals) into the mantle over a range of frequencies from 0 to 25 Hz, using the inverse Fourier transform, it computes a pulse seismogram that will be further filtered to construct SKS waves (i.e. from 0.1 to 0.3 Hz). Finally, by applying the cross-correlation method of Menke and Levin (2003) and averaging all the fast azimuths and delay times at each station measured by rotating the vertical stack of elastic tensors by 5°intervals around the y-axis, the SKS splitting parameters are determined. The software for computing mantle aggregates fabrics and SKS splitting can be found in the open source software package ECOMAN (https://newtonproject.geoscienze. unipd.it/ecoman/).

3D P-Wave Anisotropic Tomography
We use the anisotropic seismic imaging method by VanderBeek and , that solves for perturbations to P-wave slowness and three additional parameters that define the anisotropic magnitude, azimuth, and dip in a hexagonally symmetric medium. The tomographic algorithm does not require an anisotropic starting model which could potentially distort the results if not close enough to the true solution as in the case of the anisotropic imaging method of Munzarová et al. (2018). Additionally, changes in elevation and surface velocity are explicitly addressed in teleseismic imaging using 3D ray tracing through a user-defined 3D velocity model that incorporates elevation (Toomey et al., 1994).
Ray theoretical travel-times are estimated 1) with the shortestpath algorithm (Moser, 1991) through Model CWM described in the previous section using the mantle aggregates full elastic tensor at~21 Myr, and 2) with the tau-p method (Crotwell et al., 1999) outside the study area using a 1D radial Earth velocity model. The geodynamic model was centered on 42°N 12.5°E to match the main seismic structures with the real positions observed in current tomographic images.
Partial derivatives of the travel-times with respect to the model parameters are computed along the discretized ray paths. The LSQR method (Paige and Saunders, 1982) is used to solve the resulting linear system of equations relating changes in model parameters to changes in travel-times. To regularize the ill-posed inverse problem, damping and smoothing constraints are used. The choice of the regularization parameters that limit the norm of the model perturbational vector and enforce the Laplacian spatial smoothness of the model perturbations, thus controlling the length and the roughness of the solution vector relative to the length of the data residual vector, i.e. λ d and λ s respectively, is discussed in Section 2.3.1.

Starting Model, Discretization and Regularization
We use a regular grid with uniform 10 km node spacing for the forward calculation of travel-times. The initial mantle velocity model is the isotropic 1D AK135 model (Kennett et al., 1995). We applied an Earth flattening transform (Müller, 1971) to account for Earth's curvature in our Cartesian model domain.
Perturbations to the three anisotropic parameters and the mean P-wave slowness (i.e. inverse of velocity, u = 1/v) are solved on a coarser regular grid with 40 km node spacing and then, at each iteration, linearly mapped to the finer model used for travel-times calculation. Model CWM was considered down to 700 km depth, however, to limit the number of inversion parameters, anisotropic perturbations were restricted to the upper 400 km where there is the best ray crossing and mineral physics predicts mantle anisotropy to be most significant (Karato et al., 2008).
To resemble realistic conditions a first inversion was performed using delay times calculated through our Model CWM with the same distribution of sources (Supplementary Figure S2A) and receivers (Figure 2A), and the same regularization parameters as in the anisotropic tomography model of the Central Mediterranean by Rappisi et al. (2022) (Test 1). For this first test normally distributed errors with a standard deviation of 450 ms was added to the seismic data.
Next, several sets of inversions were run imposing a 1-sigma error of 125 ms applied to synthetic data and a smoothing-todamping ratio (λ s /λ d ) of 100 and damping values (λ d ) of 1,2, . . . 9,10 with different synthetic datasets. For these sets of tests we used: 1) same sources (Supplementary Figure S2A) and station array as in Test 1 (Test2); 2) an ideal on land receivers distribution ( Figure 2B) (Test 3); 3) an ideal marine and on land receivers distribution ( Figure 2C) (Test 4). In the last two cases (Test three and 4) the receivers are equally spaced 75 km apart and the teleseismic events are placed at a distance from the center of the domain from a minimum of 35 up to a maximum of 110°, every 10°of azimuth (Supplementary Figure S2B), guaranteeing a perfectly homogeneous azimuthal events distribution, thus removing any bias associated with preferential sampling of certain back azimuths.
In addition we performed 4) purely isotropic inversions in order to evaluate the effect of neglecting seismic anisotropy on the tomographic image (Test5), and 5) an inversion where the Model CWM is considered to be isotropic below 200 km to address vertical smearing of anisotropic structures (Test 6).
We constructed L-curves (Aster et al., 2018) plotting the squared model norm (|dm| 2 ) as a function of the squared norm of the delay time residuals normalized by the estimated data uncertainty (χ 2 ) for different values of damping factor (λ du ) (Supplementary Figure S3). Ideal solutions are considered those near the corner of the L-curve where an increase in model norm does not result in an appreciable decrease in data residuals. For each test, convergence is usually reached before or at iteration 3. All the inversions performed are summarized in Table 1.

Reliability of the Tomographic Results
To explore possible trade-offs between isotropic and anisotropic parameters, a synthetic inversion was aimed at reconstructing the isotropic component of model CWM. To test if velocity anomalies present in our preferred isotropic model could yield erroneous anisotropy, delays predicted through this model-not considering the anisotropic components-were inverted for both isotropic and anisotropic parameters. The result is showed in Figure 3. Isotropic anomalies were faithfully recovered with minimal anisotropic perturbations throughout the entire study area (generally < 1%) with the exception of higher-magnitude anisotropic perturbations ( < 2%) in the Southern Tyrrhenian sea, Ionian sea and Sicilian Channel. Permuted data test similar to Spakman (1991); Bijwaard et al. (1998); Rawlinson and Spakman (2016) was performed in order to assess the model amplitude error as the anomaly amplitudes are interpreted in terms of geodynamic features and errors could potentially bring to wrong interpretations. The inverted dataset is the data vector of the last iteration of Test 2, randomly permuted. The "permuted dataset" can be considered noise that has the same average, standard deviation and distribution of the delay times "not permuted" used for Test 2. By permuting the data vector order, there is no more correlation between the delay times and the raypaths. The starting model is the tomography obtained at the last iteration of Test 2. This choice allows to guarantee that the ray geometry of the permuted data test is comparable to that in Test 2. The result is showed in Figure 4. We do not observe regions with systematic anomaly patterns, on the contrary, random anomalies are recovered with delay time residual χ 2 ≈ 8.8 and the RMS amplitude variations ≈ 0.99%. This value is much bigger than the data residual observed for Test 2 (i.e. χ 2 ≈ 0.95; Supplementary Figure S3A) and for Test 3, 4 and 5 as well  Isotropic restoration synthetic test. Anisotropic inversion of purely isotropic synthetic data calculated through our model CWM, i.e., non taking into account the anisotropic patterns. While no anisotropic structures have been considered when performing the forward problem, the inversion does introduce some anisotropic perturbations. Anisotropy is represented by ellipse symbols where the major axis of the ellipse parallels the fast-direction and the minor axis scales linearly with the symmetry axis dip into the view plane such that fabrics parallel and normal to the crosssections plot as lines and circles, respectively. Anisotropic perturbations were restricted to the upper 400 km. See legend. Areas of poor data coverage are masked in grey.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 884100 6 (i.e. 0.85 < χ 2 < 3.6; Supplementary Figure S2B-D). The random distribution of the retrieved anomalies and the low data fit obtained for the permuted data test compared with the higher data fit of the correlated (i.e. the data vector not permuted) data test suggest that the results of our tests (i.e. Test 1-6) are reliable.

Geodynamic Evolution of the Central-Western Mediterranean (Model CWM)
In this section the geodynamic evolution of Model CWM is addressed (Supplementary Movie S1). The following discussion focuses on a mere description of our geodynamic model evolution. Analogies and differences between Model CWM and the real tectonic evolution of the Central-Western Mediterranean region will be discussed in Section 4.1. The Alpine and Dinaric slabs have been included in the Model CWM only to evaluate their influence in the mantle flow below the Adria plate surrounding regions, but they also represent important targets to be recovered by our tomographic inversions. Their geodynamic evolution is characterized by further slab verticalization and final break-off.
The slabs negative buoyancy drives the evolution of the two active oceanic subductions. The oceanic plates progressively sink down to the mantle transition zone and after bending start to rollback accompanied by a stretching of the overriding lithosphere. The Ionian slab migrates south-eastward while the Alboran slab south-westward. The rollback of the two slabs evolves with episodes of lateral tearing, segmentation and break-off when trenches impact with a continental margin.
The tectonic evolution of the Ionian slab is similar to that of Model CM (Lo Bue et al., 2021). In a few million years (~4 Myr), the western part of the Ionian trench collides with the African plate inducing slab tearing along the passive margin and subduction of continental crust fragments. The tear propagates along the African margin, favouring the eastward slab rollback ( Figure 5A). Subsequently, the northeastern edge of the trench reaches the thin northwestern Adria margin progressively causing subduction of the Adria continental crust, slab lateral tearing along the oceanic-continental lithosphere transition, and the formation of a curved trench due to the variations of buoyancy along it. Meanwhile, the Alboran slab rapidly rolls back westward, accommodated by lithosphere tearing along the African and Iberian margins   2021), here, this phenomenon also occurs on the side of the African continent. This leads to a final geometry of the Ionian slab characterized by the presence of two wide windows, one below the area corresponding to the current Central Apennines and one beneath the north-eastern African margin ( Figure 5D).
After~20 Myr ( Figure 5D-F), slab remnants are found in model areas corresponding to the present-day Northern Apennines, Southern Tyrrhenian sea, Alboran sea and Kabylides. At 20 Myr, the Ionian slab (portions beneath the Northern Apennines and south of the Tyrrhenian Sea -Supplementary Figure S4A) and the Alboran slab extend continuously from the surface down to the mantle transition zone.
At~25 Myr, the Northern Apenninic and Kabylides slabs hang down to~150 km depth (Supplementary Figure S4B). The first one extends further deeper from~180 km down to about~660 km depth while a horizontal segment of the Kabylides slab is still joined to the Ionian slab from~180 km down to about~300. The remaining portion of the Ionian slab (south of the Tyrrhenian Sea) instead extends continuously from the surface down to the mantle transition zone. The Alboran slab is already detached. At~30 Myr, the Calabrian slab appears still anchored to the surface (Supplementary Figure S4C) showing clear evidences of breakoff. The model evolves with the complete detachment of all the slabs ( Figure 5F).

Upper Mantle Flow, LPO, and Synthetic Seismic Anisotropy
The subduction and rollback of the Ionian and Alboran slabs in the Model CWM induces a complex flow in the surrounding mantle characterized by the presence of poloidal and toroidal components (Supplementary Movie S1). The initial sinking of the two slabs (i.e. Ionian and Alboran) generates a dominant poloidal flow component and mantle upwelling in the mantle wedge (i.e., arrows pointing downward or upward in correspondence of slabs and basins, respectively -Supplementary Movie S1). Subsequently, toroidal cells are also generated by slab rollback that forces the mantle to flow circularly around the edges of the two slabs and through the slab windows that are formed at later stages. The upper mantle fabrics patterns are reflected in those of the synthetic SKS splitting measurements shown in Figure 7. In the back-arc regions, the fast azimuths orient parallel to the trajectory of the Ionian and Alboran trenches migration. The delay times in these regions are very high (δt = 2-3.2 s), reflecting fabrics that are consistent within the entire upper mantle, and are reduced in the areas near the two trenches (δt = 1-1.5 s) due to the superposition of mantle domains with contrasting fabric patterns. In the fore arc regions, the teleseismic fast shear wave components align trenchparallel, while around the slabs edges, they form a circular pattern highlighting the underlying return flow (δt = 1-1.5 s).
Following the workflow described in Section 2.3.1, we first inverted a set of time delays computed through Model CWM using the distribution of sources and receivers as in Rappisi et al. (2022) (see Supplementary Figure S2A and Figure 2A; Test 1). We added random errors with a standard deviation of 450 ms to the data (i.e. a value corresponding to the amount of error usually encountered in real case studies). This solution reproduces realistic study conditions to test the ability of our method in recovering the main isotropic and anisotropic structures of the target ( Figure 6E-H).
The marine areas of the Tyrrhenian, Adriatic, Ionian Sea and Strait of Sicily are poorly sampled, resulting in a general loss of fast and slow anomaly amplitude. Nevertheless, the main isotropic structures (i.e. the Alpine, Northern Apennines, Calabrian and Dinaric-Hellenic slabs) are well recovered. The Northern Apenninic and Calabrian slabs are imaged as a single weak fast anomaly stretching along the N-S direction, while in the geodynamic model a~150 km wide window is present at shallow depth, i. e~100-200 km beneath central Italy ( Figure 5D, Supplementary Figure S4A and Figure 6A,B).
Test 1 exhibits a~-2% low velocity artefact in correspondence of the northern Tyrrhenian basin and Corsica-Sardinia block at 200 km depth ( Figure 6F), indicating some vertical smearing of the true low velocity structure confined in the upper 100 km of the domain. Anisotropy patterns are well recovered where seismic ray coverage is relatively abundant, e.g., the near-horizontal circular pattern of P-wave fast azimuths around the Western Alps in Southern France. Trench perpendicular steeply dipping fabrics are imaged above the Calabrian slab in the Tyrrhenian Sea, while E-W oriented fabrics are found in the Northern African margin.
The recovered isotropic and anisotropic structures from the model resulting from Test 2 ( Figure 6I-L) indicate that a better quality dataset increases the probability of better retrieving the magnitude and sharpness of the true anomalies. However, at the same time new and increased in magnitude tomography artefacts are observed. For example, the low velocity artefact at 200 km depth ( Figure 6J) in the Tyrrhenian sea, east of the Sardinia-Corsica block, is in Test 2 much stronger than it was in Test 1 ( Figure 6F) with an increase in magnitude of~1%. And also, a new 100 km wide low velocity artefact (~-2%) appears at 300 km depth south of Sicily ( Figure 6K). High velocity artefacts already observed in Test 1, such as the one in Spain and west of the Alps, are in Test 2 much bigger (i.e. joined in a single broader anomaly) and slightly stronger in amplitude, covering the entire southern portion of France at 100 and 200 km depth ( Figures 6I,J). Although the anisotropy patterns do not differ from the ones of Test 1, a reduction in their magnitude is observed above the Calabrian slab (i.e. in the Tyrrhenian Sea dipping fabrics). Probably due to the trade-off between isotropic and anisotropic components, it is worth noting that this reduction is associated with an increase in the magnitude of the isotropic fast anomaly. In Test 4, with an ideal distribution of sources (Supplementary Figure S2B) and marine and land receivers ( Figure 2C), the fast anomalies are better retrieved in terms of amplitude and spatial distribution (i.e. size and geographic position). However, many artefacts still persist, for example, at 100 km depth ( Figure 6M) the Calabrian fast anomaly exhibits a weak magnitude (i.e.,~1% vs~3% in the true model). More importantly, at 200 km depth ( Figure 6M) the slow velocity artefact located est of the Sardinia-Corsica block in the previous tests now affects the entire Liguro-Provencal and Tyrrhenian basins. Although the weak magnitude of the Calabrian fast anomaly, the gap between Northern Apennines and Calabrian slab is better retrieved at 100 and 200 km depth ( Figure 6M,N) with respect to the previous Test 1 and Test 2. The Alboran fast anomaly is recovered as well and placed in the correct geographic position ( Figure 6P). The high velocity artefact imaged in Test 2 beneath the southern France at shallow depth (100-200 km; Figure 6I,J), now disappears and the recovered model better resembles the target. The ideal station coverage also helps in retrieving anisotropic patterns and magnitude at every depth slice, including the well sampled marine areas where the NW-SE fabrics are now recovered in the Tyrrhenian Sea. Similar results, i.e. about sensitivity of teleseismic P-wave tomography under different conditions, have been previously described (e.g., Spakman and Nolet (1988); Lévěque et al. (1993); Rawlinson and Spakman (2016)).
Considering that placing marine receivers is a costly procedure, we also performed a set of inversions with an ideal distribution of on-land receivers only ( Figure 2B). The result is shown in Figure 8A-D. The main effect of having reduced the FIGURE 8 | Depth slices at 100 km, 200 km, 300 km and 400 km depth for the tomographic results from Test 3 (A-D), Test 5 (E-H) and Test 6 (I-L). Isotropic anomalies are plotted with respect to the starting model. Anisotropy is plotted using ellipses as described in Figure 6. Areas of poor data coverage are masked in grey.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 884100 number of receivers is the underestimation of anisotropy in poorly sampled areas. For example, the trench-parallel patterns bordering the eastern side of the Apenninic fast anomaly in Test 3 ( Figure 8A-D) appears weaker than it is in the true model and in Test 4 ( Figure 6A-D, M-P). Figure 8 also shows the results of Test 5 and 6 performed, respectively, in isotropic approximation ( Figure 8E-H; i.e. ignoring seismic anisotropy) or with a model that is anisotropic only in the top 200 km ( Figure 8I-L). In both cases, we observe that the isotropic solution contains a number of fast anomaly features broadly consistent with the true model ( Figure 6A-D). However, several slow velocity artefacts are imaged around and above the main slabs (i.e. Alpine and Calabrian slabs) especially when not considering seismic anisotropy ( Figure 8E-H). Figure 8I-L indicate that seismic anisotropy is retrieved also at depths below 200 km where the true model is instead isotropic. This suggests that when only using teleseismic P-waves anisotropic structures are vertically smeared similarly to isotropic anomalies.

How Well Does Model CWM Fit Seismological Observations?
In this work we have extended the modeling methodology of Lo Bue et al. (2021) to build a structurally complex geodynamic model which was then exploited to test the capabilities of anisotropic P-wave seismic tomography to recover a Mediterranean-like subduction environment. With respect to Model CM of Lo Bue et al. (2021), the new geodynamic model CWM has been updated by using a different paleotectonic configuration characterized by the presence of additional subducting plates (i.e., Alboran, Alpine and Dinaric subduction zones).
Similarly to previous studies (Luth et al., 2013;Jagoutz et al., 2015;Holt et al., 2017Holt et al., , 2018Király et al., 2018;Peral et al., 2020), here we notice that the presence of multiple subducting slabs influences the overall force balance, the geometry and kinematics of the subduction systems, as well as the mantle flow patterns. The inclusion of additional subduction zones has partly improved the prediction of the mantle dynamics leading to a better correspondence between the modeled and observed surface and deep isotropic structures, and seismic anisotropy patterns when compared to Model CM of Lo Bue et al. (2021). A quantitative comparison between predicted and observed SKS splitting measurements (Figure 7 Figure S5). In detail, the general pattern of the synthetic 3D anisotropy calculations matches the observed data with a relatively lower misfit angle in the model areas corresponding to the Iberian peninsula and the Apennines chain if compared to the Model CM (Lo Bue et al., 2021). In the Dinaric Alps area the synthetic fast splitting directions are trench parallel. Here, the average angular misfit between predicted and observed fast azimuths remains high, but with slightly reduced values of about 70, −,80°degrees, compared to 90°d egrees of Lo Bue et al. (2021). Lastly, we acknowledge that split times are not particularly well-fit considering that the residual mean and standard deviation are 667 and 592 ms, respectively (Supplementary Figure S5). However, the magnitude of the average misfit of the whole model is mainly due to the broad mismatch found in the south of the Iberian Peninsula.
Uncertainties in the initial model geometry and in the modeled mantle rheology are likely responsible for the major discrepancies. Other sources of mismatch could be related to the modeling of the mantle textures, and to the presence of fossil fabrics within the oceanic and continental subducted lithosphere that have not been included here. Furthermore, the employed free slip boundary conditions prevent lateral mantle flow across the bottom and vertical boundaries. The large discrepancy in the area of the southern Iberian Peninsula may be partly due to the fact that the model does not account for the Cenozoic Eurasia-Africa convergence, and the relative position of the African plate has remained fixed since~30 Ma differing slightly from the presentday one. As such, the Alboran arc is positioned further south than its present-day position (under Morocco).
Although Model CWM is based on paleogeographic and tectonic reconstructions of the region in the Oligocene-Miocene (Lucente and Speranza, 2001;Lucente et al., 2006;Faccenna et al., 2014;van Hinsbergen et al., 2014;Romagny et al., 2020), geometrical assumptions, that could potentially bias the final output, were required due to limitations imposed by numerical modeling. First, the initial geometry and thermal ages of the subducting slabs were partly simplified and this could strongly influence the comparison with seismological data. This could be the case of the Alps and Dinarides where a simplified initial portion of the subducted lithosphere was imposed to model flow barriers. We note that the detachment of Alpine slabs prior to collision is still debated in some areas, such as the Western Alps where Kästle et al. (2020) favor the interpretation of a recent European slab break-off, consistent with observations of strong exhumation and sedimentation that started around 2-7 Ma ago and is still ongoing (Escher and Beaumont, 1997;Kuhlemann, 2007;Fox et al., 2016;Nocquet et al., 2016). On the contrary, Zhao et al. (2016) document the lateral continuity of the European slab from the Western Alps to the Central Alps, and the downdip slab continuity beneath the Central Alps, ruling out the hypothesis of slab break-off to explain Cenozoic Alpine magmatism. The teleseismic P-wave tomography of Rappisi et al. (2022), referred as model ani-NEWTON21, exhibits a continuous slab beneath the Alps, divided into an Eastern, Central, and Western segment characterised by changes in dip. Similarly, the extent of the Dinaric slab at~30 Ma is largely debated. Post-collisional uplift and contemporaneous emplacement of igneous rocks (33-22 Ma) in the internal Dinarides may suggest either 1) "that the Oligocene-Miocene orogen-wide uplift was driven by post-break-off delamination of the Adriatic lithospheric mantle" Balling et al. (2021), or 2) verticalization of the Adria slab driven by slab pull and consequent upper plate extension, which is exactly what is modeled in our Model CWM. In conclusion, the available geophysical and geological data do not allow to discriminate between a model of post-collisional slab break-off and one of post-collisional slab verticalization (as modelled for the Alps and Dinarides in our Model CWM), as both would imply upper plate extension, uplift and magma emplacement (Faccenda et al., 2008(Faccenda et al., , 2009. Secondly, the tectonic reconstruction of Romagny et al. (2020) shows that~30 Ma the Mesozoic Tethyan lithosphere was consumed in two different trenches located from the Alps to the southeast of the Baleares and in the Alboran domain, respectively. Two incipient slabs (~150-200 km) were already subducted in the upper mantle. However, to trigger a "spontaneous" subduction system, the Ionian trench has been initially positioned further south with the slab extending to a depth of 300 km, while the Alboran one further west with a 350 km deep slab, in order to model a more developed subduction and increase the slab negative buoyancy. This could cause a difference in rates of Ionian and Alboran slabs retreat at the model early stage when compared to those reported in the literature. However, we note that in the reconstructions by Faccenna et al. (2014) and Romagny et al. (2020) the initial plate geometry at~23 Ma does not differ substantially from that at~35 Ma, and from our initial setup. This is likely related to the slow dynamics typical of incipient subduction systems.
Despite the modeling limitations, Model CWM reproduces several episodes of slab lateral tearing and break-off that have been proposed according to geological and seismological data.
On the contrary, the modeled Alboran slab, in addition to being in a wrong position (i.e. further south than its current position), possesses a morphology which is not entirely realistic. This is probably due to the imposed initially 350 km long slab and to the slab tearing occurring as soon as the trench interacts with the continental margins, thus preventing any arcuate shape of the margin. However, we note that the geometry and length of the Alboran slab in model CWM at 0 (initial conditions; Figure 1) and~20 Ma ( Figure 5D) are similar to those obtained by (Chertova et al., 2014)

How Well Does Tomography Recover the Target Model?
We performed seismological forward and inverse simulations by testing different types of data coverage and quality. To help the comparison between the different tests and evaluating their capabilities in recovering model CWM, Figure 9 shows the difference between the true model and the solution of Test 1, Test 2, Test three and Test 4, both in terms of isotropic and anisotropic structures (i.e. dlnVp = dlnVp true model -dlnVp tomography model). We observe that the average difference in retrieved isotropic velocity is in general low (~[−1%,+1%]) and gradually decreases moving from Test 1 to Test 4. For example, Figure 9A-D shows maximum values of dlnVp of~3% for Test 1, i.e. in the Apenninic slab at 400 km depth ( Figure 9D), that gradually decrease to~1% for Test 4. For Test 4 (i.e. test with perfect data coverage), higher values are observed in the western side of the model, with peaks of~1.5% in the Liguro-Provencal basin at 200 km depth ( Figure 9N). True and recovered anisotropy patterns are plotted in black and red, respectively, showing high degree of matching both in terms of azimuth and dip. With few exceptions of sparse differences in dip angles, no particular areas of discrepancy are identified.
From our results it emerged that even with a non-ideal sourcestation coverage the recovery of isotropic structures and anisotropic patterns is quite good, although anisotropy magnitude is overall underestimated (especially in poorly sampled areas). This suggests that the amount of mantle anisotropy could be higher than that retrieved by tomographic models with commonly uneven source-receiver distributions. The consequences of the inhomogeneous distribution of seismicity and stations on ray coverage and on retrieved tomographic images is known in isotropic tomographic models (e.g., Bozdağ et al. (2016); Boschi and Dziewonski (1999); Masters et al. (1996); Antolik et al. (2003); Dalton and Ekström (2006); Ruan et al. (2019)). Here we show that similar problems are also found in anisotropic seismic tomographic models (e.g., causing underestimation of anisotropy magnitude). In addition, it emerged that tomographic images calculated from data with a scarce seismic coverage are potentially affected by the presence of anomalies placed in a wrong geographic position. This is the case of the Alboran slab that in Figure 6E-L appears shifted toward the east. This kind of artefacts could bring errors in the tomographic model interpretation when fast anomalies are present close to the boundaries of the sampled area. Increasing data quality (i.e., decreasing data error; Test 2) helps in better retrieving the magnitude of the isotropic and anisotropic structures, but at the same time leads to an increase in the magnitude and size of the artefacts in poorly sampled areas ( Figure 6E-L). In addition, Figure 6M-P shows that ideal data coverage allows for a more accurate retrieval of anomaly magnitudes without increasing artefact amplitudes. However, we note that the higher number of receivers (e.g., in the Tyrrhenian and Liguro-Provencal basins) at 200 km depth amplifies the smearing of the upper low-velocity layer with respect to Test three where, on the contrary, the limited number of stations (i.e. limited rays) reduces this effect.
In the inversions where seismic anisotropy is ignored (Test 5), we observe that several slow anomalies appear in the tomographic sections ( Figure 8E-H). This is especially evident in the area north of the Alps ( Figure 8A) and below the Calabrian slab (i.e. in the Ionian Sea, Figure 8B-D). Considering that these anomalies are not present in the true model and indeed completely disappear in the anisotropic inversions ( Figure 6), it follows that they are seismic artefacts due to the isotropic approximation.
Lastly, the test carried out on the model isotropic only from 200 km depth down (Test 6), showed that both the isotropic and anisotropic features are subjected to vertical smearing (8g,h). This should be taken into account when interpreting teleseismic P-wave anisotropic tomography.
In order to further characterize model differences between true and tomographic models, we have computed the average misfit values for fraction of anisotropy (df), azimuth and dip angles (i.e. dψ and dγ) with increasing depth ( Figure 10A-C) with respect to the true values. We observe that the average solution gradually improves, better resembling the true model, with decreasing data error and improving data coverage. The higher values of misfit are in fact observed for the model obtained from the inversion performed with the bigger data error (i.e. 450 ms, Test 1) and the worst receiver distribution. This is true for both df and dψ, while for dγ is valid below~70 km depth ( Figure 10C). For all models the average azimuthal misfit is highest in the upper 50 km (due to the poor ray coverage at these depths by teleseismic P-waves), below which it rapidly decreases and remains roughly constant with depth except for a slight increase toward the bottom of the anisotropic domain. In contrast, the dip angle average misfit gradually increases with depth. The misfit curves for Test 4 and Test 6 show similar shapes but with shifted absolute values in the upper 150 km. This indicates that the presence of deeper anisotropy (Test 4) associated with poor vertical resolution deteriorates the quality of the retrieved shallower structures.

CONCLUSION
We applied the modeling methodology of Lo Bue et al. (2021) to simulate the geodynamic evolution over~20-30 Myr of a model that presents similar characteristics to those currently observed in the Central-Western Mediterranean region (e.g., detached or stagnating slabs, slab windows, etc). To quantify similarities and discrepancies between the obtained geodynamic model and the current tectonic setting, the model results were verified by comparing seismological synthetics (isotropic P-wave anomalies, P-wave anisotropy and SKS splitting) and major tectonic features (i.e., slab and trench geometry) with observations. This comparison confirms that, with respect to the previous study of Lo Bue et al. (2021), using a more complex initial geometry (i.e. including the Alboran, Alpine and Dinaric-Hellenic slabs) allows us to perform a step forward toward the better recovering of the mantle flow, overall evolution and current tectonic beneath this region. However, we note that model CWM is still far from reproducing the exact evolution and present-day tectonic setting of the area and further studies need to be performed in this direction. For example, next-level numerical studies should attempt to improve the model geometry by considering the Earth's sphericity and the Africa-Eurasia plates convergence.
Despite the several limitation of the numerical methods (e.g., Cartesian coordinates system, no plates convergence, no fossil LPO, fabrics within the lithosphere, free slip boundaries, no mantle in/outflow, etc.) and the assumptions necessary to start and drive the simulation self-consistently (e.g., initial slab depths, mantle rheology parameters, etc.), we observe that at~20 Myr model CWM exhibits interesting geological features resembling those found in the Central-Western Mediterranean (e.g., Calabrian slab continuous from the surface down to the base of the upper mantle, the presence of two wide windows in the Ionian slab, etc). For this reason, we used the modeled elastic properties at this stage (i.e. the elastic tensors at~20 Myr), to perform 3D P-wave anisotropic tomography using the approach proposed by VanderBeek and Faccenda (2021). Using the geodynamic model as reference model, we evaluated the capabilities of seismic tomography to recover a complex subduction environment in different conditions, such as poor  station coverage and bad data quality. From the seismological inversions and the comparison between purely isotropic and anisotropic solutions it emerges that 1) it is fundamental to invert for anisotropy to improve the reliability of the tomographic result and 2) even a non-ideal source-station coverage allows to recover isotropic structures and anisotropic patterns from teleseismic P-wave tomography. Anisotropy magnitude, although consistent with those of the synthetic target model, is overall underestimated in the upper mantle especially in poorly sampled areas. In light of this, it is recommended to increase the number of marine and land stations and improve the accuracy of teleseismic arrival time measurements. However, it should be noted that perfect coverage of receivers does not guarantee an ideal tomographic solution.
For example, Test 4, despite being performed with receivers distributed over the entire study area, presents various imaging artefacts. Future steps aiming at recreating a "perfect coverage" should be characterized by a good ray sampling, thus to include seismic rays that cover different directions in order to guarantee an excellent resolution (e.g., not only teleseismic events). Furthermore, although the synthetic inversions confirm that the developed methodology for P-wave anisotropic tomography is capable of retrieving with a good approximation the modeled upper mantle structures, the employed geodynamic simulations do not account for compositional variations, presence of fluids/melt and lithospheric fossil fabrics that can affect the seismic properties of natural tectonic settings. The presence of these further complexities remains to be tested, and it will be considered in future studies. The synthetic tomography results demonstrated that using a combination of geodynamic and seismological numerical modeling techniques could represent a powerful tool to investigate mantle dynamics. Although the modeling limitations, we obtained a 3D complex mantle structure that partly resemble some main characteristics of the actual presentday mantle in the Central-Western Mediterranean. This opens new perspectives towards the future possibility of creating models of the geodynamic evolution that can be constrained by the structure and mantle anisotropy obtained from P-travel time anisotropic tomography. To better constrain the initial tectonic configuration, an interesting future development would be to formulate a fluid dynamic inverse problem to reproduce unknown mantle flow back in time from seismic tomographic observations of the mantle and reconstructions of past plate motions using variational data assimilation (Bunge et al., 2003). Adjoint modeling is in fact a great opportunity to produce realistic mantle retrodictions models. However, the development of testing of this methodology is still far from being applicable to complex 3D tectonic settings such as the Mediterranean. This technique has been successfully applied to reproduce the recent dynamics in the South America and North America subduction zones (Hu et al., 2017;Zhou et al., 2018), However, we believe that exhumation back in time of the slabs stagnating in the mantle transition zone is non-trivial when the plate convergence rate is quite small (basically, by inverting gravity there is not easy way to exhume back at the surface these slabs), which is one of the reason why we choose to model forward in time the Central-Western Mediterranean dynamics. Reuber and Simons (2020), although showing the potentials (and limitations) of this technique on quite simple 2D and static (not dynamic) model configurations, concluded that the method "needs to be thoroughly tried and tested on real-world examples". Adding mantle fabrics to improve the mantle flow retrodictions is an ongoing research activity in our group, that we hope to include in our models in the near future.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors conceived the study. RLB performed the geodynamic and seismological (i.e. strain-induced fabric estimation and SKS splitting calculations) numerical modelling. FR performed the seismological forward and P-wave tomography. BPV and MF supervised the findings of this work. All the authors contributed equally to the discussion of the results and to the conclusions of this study.

FUNDING
This study is supported by the ERC StG 758199 NEWTON.

ACKNOWLEDGMENTS
T. Gerya provided the I3MG code used for the subduction modeling. The modified version of the D-Rex code used for the fabric modeling, and the routines used to calculate SKS splitting parameters, P-wave and Rayleigh wave anisotropy can be found inside the ECOMAN software package (https:// newtonproject.geoscienze.unipd.it/ecoman/). The MATLAB toolbox geomIO used to define the geometry of the model initial setup can be found at https://geomio.bitbucket.io/. Paraview was used for graphic visualization of the model output (https://www.paraview.org/). Tomographic maps were created using Generic Mapping Tools (Wessel et al., 2019) with colormaps developed by Crameri (2018a,b). The manuscript was significantly improved thanks to the constructive feedback and comments from the editor, Wim Spakman and Claudia Piromallo.