Earthquake Characteristics and Structural Properties of the Southern Tyrrhenian Basin from Full Seismic Wave Simulations

Modelling the response of seismic wavefields to sharp lateral variations in crustal discontinuities is essential for seismic tomography application and path effects correction in earthquake source characterization. This is particularly relevant when wavefields cross back-arc oceanic basins, i.e. mixed continental-oceanic settings. High-frequency (>0.05 Hz) seismic waves resonate and get absorbed across these settings due to a shallow Moho, crustal heterogeneities, and energy leakage. Here, we provide the first high-frequency wave-equation model of full seismograms propagating through realistic 3D back-arc basins. Inversion by parameters trial based on correlation analyses identifies P-, S- and coda-wave as attributes able to estimate jointly 3D Moho variations, sediment thickness, and earthquake source characteristics using data from a single regional earthquake. We use as data waveforms produced by the Accumoli earthquake (Central Italy, 2016), propagating across the Southern Tyrrhenian basin and recorded across Southern Italy. The best model comprises a deep Moho (∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}18 km) in the middle of the basin and a crustal pinch with the continental crust in Sicily. The deep Moho corresponds to the Issel Bridge, a portion of continental crust trapped between the Vavilov and Marsili volcanic centres. The Accumoli earthquake is optimally described at a depth of 7.3 km using a boxcar with rise time of 6 s. Our results show that the early S-wave coda comprises trapped and reverberating phases sensitive to crustal interfaces. Forward modelling these waves is computationally expensive; however, adding these attributes to tomographic procedures allows modelling both source and structural parameters across oceanic basins.


Article Highlights
• Modelling the high-frequency seismic response to sharp crustal variations in mixed continental-oceanic settings • High-frequency forward model of the seismic wavefield propagating through a realistic 3D back-arc basin (South Tyrrhenian) • Improved estimate of Moho depth, sediment thickness and seismic source features using a single regional earthquake

Introduction
Back-arc basins characterized by mixed continental-oceanic crust and high-absorption magmatic systems linked to subduction (Sartori 2003) are challenging environments for fullwaveform and adjoint tomography techniques (Blom et al. 2020;Magnoni et al. 2022). These settings distort seismic wavefields significantly, dispersing high-frequency information. Consequently, standard inversions of phases and amplitudes measured on the recorded wavefield produce inconsistent results with standard seismological approaches (Cormier and Sanborn 2019). Only waveforms at frequencies lower than 3 Hz are consistently recorded across oceanic basins at the regional scale. However, waves recorded down to 0.05 Hz show resonant scattering properties due to strong interactions with large-scale velocity discontinuities of sizes comparable with seismic wavelengths (Cormier and Sanborn 2019). Strong scattering from the Moho (Sanborn et al. 2017;Sanborn and Cormier 2018) and loss of energy due to leakage into the mantle (Wegler 2005;Margerin 2017) trade-off with the standard description of wave attenuation as a combination of scattering and absorption (Morozov and Safarshahi 2020;Nardoni et al. 2021), making the application of any standard waveform-based tomography unfeasible. These effects are a primary reason why high-frequency full-waveform inversions across basins are problematic, and models are generally inverted using few coherent attributes. Seismic forward modelling allows discriminating the stochastic signatures produced by crustal highly scattering and absorption media from coherent reverberations due to shallow Moho. These reverberations likely have sensitivity to Moho depths, generally measured using teleseisms ) or the ambient noise field (Manu-Marfo et al. 2019). Modelling and eventually including them in full-waveform inversions (Fichtner et al. 2009) could provide new tools to image oceanic basins.
Wave-equation modelling is the primary forward modelling strategy in seismology. Recent advances in parallel computing and numerical methods enable large-scale 3D numerical simulations of the seismic wavefield by taking advantage of fast computing codes adaptable to the studied area (Komatitsch et al. 2010;Maeda et al. 2017;Magnoni et al. 2022). These simulations tackle fundamental mode surface waves and the topography effects directly (Takemura et al. 2015), accounting for the complete physics of 3D wave propagation (Chen et al. 2007;Fichtner et al. 2009;Romanowicz et al. 2020). Advances in modelling also enable implementing the interaction between Rayleigh and body waves during scattering by including the contribution of the fundamental mode surface wave (Xu et al. 2022). However, no wave-equation model has ever tackled strong high-frequency reverberations from shallow Moho variations and statistical variations in crustal heterogeneity across basins. Both these phenomena generate stochastic wavefields (coda waves) that are primary attributes to image back-arc basins (Sanborn and Cormier 2018). Coda waves are thus not part of full-waveform inversion strategies, limited to minimizing misfits of coherent amplitudes and phases.
Radiative-transfer theory is the standard for forward modelling such as high-frequency wavefields, sometimes combined with wave-equation modelling across continental and volcanic settings (Przybilla et al. 2006;Obermann et al. 2013). Radiative transfer theory has effectively reproduced seismic amplitudes across 3D basins, recognizing those with the highest sensitivity to basin structures: P-wave, S-wave, and early S-coda amplitudes (Sanborn et al. 2017;Sanborn and Cormier 2018;Cormier and Sanborn 2019). However, this technique offers no full phase-amplitude modelling suitable to fit high-frequency attributes and, eventually, full seismogram data. In this study, we use constraints from coda-attenuation imaging and radiative-transport-based forward modelling (Nardoni et al. 2021) to perform such wave-equation forward modelling of high-frequency waves crossing a 3D basin.
The OpenSWPC framework (Maeda et al. 2017) can include constraints from radiative transfer models (Nardoni et al. 2021) in a wave-equation finite-difference code. Developed by the heterogeneous Earth community (Sato et al. 2012) and tailored to work in high-attenuation and high-contrast lithospheric media, OpenSWPC includes statistical fluctuations in velocity and heterogeneous topography, typical of a mixed oceanic and continental crust. With its efficient parallelization, it represents the ideal tool to solve the primary issue associated with high-frequency wave-equation modelling across basins: computational time. Computational time remains the primary challenge to modelling highfrequency seismic wavefields across 3D realistic basin structures.
The Italian peninsula and the back-arc Tyrrhenian basin are ideal for testing and discriminating the potential of full-waveform modelling in mixed continental-oceanic settings. This basin is the result of multi-stage extension following the Ionian slab rollback (Malinverno and Ryan 1986;Doglioni et al. 1999;Faccenna et al. 2001;Lo Bue et al. 2021). Its irregular extension has supposedly entrapped the Issel Bridge (a portion of thinned continental crust) within the Vavilov and Marsili volcanic centres (Sartori 2003;Cocchi et al. 2009). Extended seismic networks ( Fig. 1) have already enabled the use of state-of-the-art full-waveform codes to image the lithosphere across the Mediterranean (Blom et al. 2020). Receiver functions (Piana Agostinetti and Amato 2009;Di Stefano et al. 2011) and noise-based surface-wave imaging (Manu-Marfo et al. 2019) are the primary source of information for Moho depths across this basin.
The station coverage in the region remains too sparse to discriminate Moho variations across the Issel Bridge and at the interface between the basin and the continental crust with these approaches (Manu-Marfo et al. 2019;Nardoni et al. 2021;Magrini et al. 2022). In such mixed continental-oceanic settings, full-waveform inversion studies have imaged the crust and the upper mantle at large epicentral distances (>500-600 km) using lowfrequency (<0.05 Hz) body and surface waves (Fichtner and Villaseñor 2015;Simutė et al. 2016). None of these applications comprises interface modelling, as the low-frequency smooth gradients necessary to full-waveform inversions only provide information on seismic velocities. Considering regional distances (<600 km) across the Tyrrhenian basin allows investigating higher-frequency waves, where the crustal heterogeneity and, thus, the dependency on Moho and sediment structures become relevant.
Seismic waveforms propagating across the basin are sensitive to bathymetry, whose effects trade-off with those from the Moho and earthquake source. The properties of the crust are still dominant up to 1 Hz and at regional distances, as assessed by Furumura and Kennett (1997). Alvarado et al. (2007) performed forward waveform modelling to sample the crust, finding that synthetic regional seismograms (epicentral distances from 100 to 350 km) are sensitive to crustal parameters, especially the crustal thickness. However, in the frequency band 0-2 Hz, Fu et al. (2002) investigated the energy attenuation of regional phases across a randomly rough topography via numerical modelling, showing 1 3 that topographic scattering attenuates regional waves. As the irregular topography and its stochastic fluctuations disperse reverberating phases (Takemura et al. 2015), Zhang et al. (2022) propose a method to suppress the trade-offs caused by the topographic scattering in imaging the interfaces with the receiver functions. After including topography in the modelling, reverberating phases (e.g. Lg-wave) within the S-coda (Sens-Schönfelder et al. 2009) become an additional tool for estimating Moho depths across basins using regional earthquakes. Additionally, the high-velocity contrast at the Moho and the interference of other phases, such as lower crustal (Pg) and upper mantle (Pn) phases, with the critically reflected Moho rays (PmP) change coherent phases and coda amplitude and arrivals. Their potential for estimating Moho depths has never been investigated in the past.
In this work, we demonstrate that high-frequency inversions of structural and source parameters across mixed continental-oceanic settings are feasible using: (1) regional observations from a single earthquake and (2) both coherent and stochastic seismic attributes. We establish the sensitivity of direct and coda waveforms recorded at frequencies up to 0.3 Hz to the complexity of basins, highlighting the physics of the existing trade-offs. First, the sensitivity of phase arrivals and amplitudes to both propagation medium and source has been assessed using wave-equation modelling on a limited number of models available in the literature. Then a best model has been defined through correlation analyses (Alvarado et al. 2007;Lee et al. 2014) (Sects. 2.2, 2.3 and 3.2) using seismograms from arrays surrounding the basin (Fig. 1). The final result is a model of the basin structure and source characteristics that best correlates with the existing data. Once applied to specific modelled The black dot is the modelled source, and the red triangles indicate the location of the receivers in the simulation. The dashed green line represents the northsouth profile taken as a reference to characterize the seismic source and velocity profiles and visualize the correlation between the data and the synthetic data corresponding to different Moho depths attributes, this approach allows estimating earthquake source characteristics, sediment thickness and Moho depths using a single regional earthquake.

Methods
We use Open-source Seismic Wave Propagation Code (OpenSWPC, Maeda et al. 2017) to model seismic waves in a 3D isotropic viscoelastic medium at the regional scale. The numerical simulation based on a parallel architecture solves the equations of motion described by the velocity and stress components in the 3D Cartesian coordinate system with viscoelastic constitutive equations, using a finite difference method. The code is equipped with a frequency-independent attenuation model based on the generalized Zener body and an efficient perfectly matched layer for absorbing boundary conditions. OpenSWPC is the only open-access wave-equation modelling code that allows including statistical fluctuations of seismic velocities within layered regional structures. Also, the code allows users to conduct seismic wave propagation simulations that explicitly target media with a mixed oceanic and continental crust: a crucial advantage over other codes for these settings.

Model Parameters Setup
The model covers the area of the Tyrrhenian basin and the Italian peninsula (Fig. 1). The volume used for the 3D simulation is 600 km by 750 km in the horizontal plane and 130 km in depth. It is discretized with a grid spacing of 0.5 km in the horizontal and vertical directions. Due to the simulation grid and the cluster memory, the maximum frequency of the synthetic wavefield is 0.33 Hz.
The initial layered model is built using the data available for the topography ( (Table. 1). We used , a and a Von Kármán distribution function to implement stochastic random heterogeneities in the crustal structure. The heterogeneities are kept fixed, while the velocity of each layer has been varied to fit phases arrivals and amplitudes (Table. 2). We locate the seismic receivers using the INGV network coordinates (Fig.1). To investigate the sensitivity of the full-wavefield to velocity and source variations (see Sects. 2.2 and 2.3), we consider the seismic receivers located in Sicily (Fig.1).

Source Implementation and Starting Models
We implement a seismic source corresponding to the earthquake that occurred in Accumoli (central Italy, 2016) and characterized by a constant magnitude M w of 6.1 (Tinti et al. 2016). The seismic sources are described in terms of moment tensor, source time function and depth using far-field and long-period data (Hanks and Wyss 1972;Ferreira and Woodhouse 2006). Besides, at regional scales, both body and surface waves have been used to provide source parameters for intermediate earthquakes (M∼ 6) using a point source Table 1 Starting models. representation (Wallace et al. 1981). In the moment tensor representation strike, dip, and rake describe the focal mechanism. The earthquake source is assumed point-like, as the wavelengths are longer than the fault dimensions and much shorter than the source-receiver distance (Aki and Richards 2002). Several studies estimated the source parameters assuming different moment rate functions (MRFs, Tinti et al. 2005;Cirella et al. 2018;Supino et al. 2019) to describe the fault displacement duration (rise time) and velocity. Tinti et al. (2016) modelled the fault slip velocity by imposing a boxcar source-time function. Nevertheless, the MRFs duration and shape affect the synthetic seismogram envelope (Takemura et al. 2020;Krischer et al. 2017). We test the source kinematics assuming a boxcar function and then a triangle function and explore different source parameters (i.e. hypocentral depth, focal mechanism and rise time, see Supplementary Materials), keeping fixed the magnitude M w at 6.1 (Tinti et al. 2016). The OpenSWPC forward simulations allow the evaluation of the source properties from the fit with real data recorded across a sourcereceiver path (Takemura et al. 2020). This fit is performed both qualitatively, by looking at recordings and synthetics, and quantitatively, by single waveform and 2D spatial correlation analysis (Alvarado et al. 2007).
We define different source models (Table 1) by varying the focal depth, the focal mechanism, the moment rate function, and the rise time separately in each model. Tinti et al. (2016) and Cirella et al. (2018) obtained an average rise time smaller than 2 s for each sub-faults of the modelled fault of the Accumuli earthquake. As the majority of the seismic moment is released between 2 and 6 s, we set the boxcar function with rise time T R = 4 s (Tinti et al. 2016). We vary the hypocentral depth (relative to the sea level), strike, dip, and rake separately, based on previous studies (INGV Istituto Nazionale di Geofisica e Vulcanologia, GCMT Global Centroid Moment Tensor Catalog, USGS Earthquake Hazards Program, Tinti et al. 2016) (Source A and B, Table 1). The changes in earthquake depths are primarily observed on S-wave phases, with minimal changes in P-and coda arrivals (Fig. S1). We compute the correlation coefficients between the synthetic and real seismic envelopes at each station ( Fig. 1) and map the percentage differences between the investigated source parameters. We select a depth of 7.3 km as it shows the highest correlation values for the time interval 0-200 s comprising the P arrival and the crustal S waves packet (Fig. S1d). Among the three proposed radiation patterns, the USGS estimation appears most inconsistent with the synthetic waveforms and is excluded from further analyses (Fig. S2). Again, we select the INGV estimate, which shows the best fit to the S-wave peak ground velocity.
The wavefield primarily depends on the source time function shape. Synthetic seismograms were thus obtained using both the boxcar (Fig. 2 Table 1) and the triangle ( Fig. S3-Source D) functions for different rise times. Source duration effects produce a decrease in seismogram amplitudes and smoothing with time as the duration of the source increases (Matsuzawa et al. 2009). Only direct waves are usually considered in the estimation of source characteristics. Therefore, we set a boxcar function with an initial duration of T R = 4 s to model the source (Fig. 2) as it gives higher correlation values.

Starting Velocity and Sediments Models
After setting the starting source parameters fitting at best compressional and regional phases amplitude, we investigate how the larger-scale variability in crustal structure affects the seismic wavefield (Furumura and Kennett 1997). We explore the effects of structural models differing from the starting model (compare Tables 1 and 2) performing an inversion by parameters trial focused on velocity variations and interfaces. These structural models were proposed by Nardoni et al. (2021) based on radiative-transfer simulations. The source-station geometry (north-south) allows the discrimination of radial (north-south) from transverse (west-east) components easily.

Results
Based on the comparison between different velocity models, Fig. 3a shows the best agreements between the synthetic and recorded phases (see Fig. S4 for the other models). As expected, the P-phase (head wave) on the north-south and vertical components shows sensitivity exclusively to the P-wave velocity ( V P ) variations for the crust and mantle. The vertically polarized S wave (SV) is more difficult to detect on radial and vertical components as reverberations between the free surface and the Moho increase P-coda. This likely explains the difference in phase and amplitude on north-south and  (Table 1). a: Comparison between the station recording (coordinates: 37.89 • , 13.30 • ; triangle contoured by the black edge) and the synthetic waveforms that are obtained using different rise times for the boxcar source time function. The depth is set to 7.3 km. b and c: Enlargements of the two insets in panel (a), which are indicated by the dashed red and blue lines, respectively. They show the comparison between envelopes in the two time windows corresponding to compressional waves arrival (dashed red line) and from S-wave arrival (dashed blue line), respectively. d: Percentage differences of the correlation coefficients computed for a rise time of 4 s with respect to 8 s in the time window 0-200 s vertical components between data and synthetics around 150 s (Fig. 3a). Instead, we can satisfactorily reproduce the arrival of the east-west component, which best describes the horizontally polarized S wave (SH). The primary controller of the east-west waveforms is thus the S-wave velocity ( V S ) in the crust and upper mantle. Velocity model 1.3 (Fig. 3), corresponding to higher V P and lower V S velocity in the crust, yields the best fit between the recordings and the synthetic seismograms for both P-and S-waves.

Sensitivity of the Coda to Sediments and Adjustment of Rise-Time Variations
Model 1.3 does not reproduce coda amplitudes, which have shown greater sensitivity to sediments (Nardoni et al. 2021) and are also affected by changes in source rise times (Fig. 2). We thus model the effects of sediment thickness variations on the full waveform ( Fig. 3b and c). A thick sedimentary layer primarily affects the S-coda wave-packet propagating and reverberating inside the crust by delaying the maximum amplitude peak (Borleanu et al. 2017). Thinner sediments instead dampen coda reverberations (Fig. 3b), with the east-west component most sensitive to small sediment thickness variations. For a thicker and low-velocity sedimentary cover (Molinari and Morelli 2011), reverberating and slower waves dominate the later part of the seismogram (Fig. 3c) delaying the maximum amplitude and increasing the coda duration (Kaneko et al. 2019). Considering the results with evidence from literature (Moskalenko 1992;Pepe et al. 2005;Pontevivo and Panza 2006), thinner sediments (velocity Model 1.3- Table 2) appear to characterize the continental areas. At the same time, there is no sedimentary cover detectable under the sea at these frequencies. Two-km-deep continental sediments do not produce a good fit to coda amplitudes (Fig. 3c). We thus decreased their quality factor and thickness (from 2 to 1 km) across the continental crust according to the results of Nardoni et al. (2021) (Fig. S5).  These variations do not change the P-and S-wavefields as much as the rise time of the source time function (compare Figs. 2-S6) but drastically improve coda-wave fitting. So far, we selected the best rise time in the starting source model by looking at direct phases without considering the coda as a relevant attribute. In Fig. S6, we used velocity model 1.3, thinner sediments and lower quality factor while changing the moment rate function duration. Setting the best velocity model for the phase arrivals, we show that the coda is very sensitive to the duration of the moment rate function. In the final model (Table 2) used to estimate Moho variations, the rise time is thus set to 6 s for Model 1.3 (Table 2). This source model allows a good fit for direct phase arrivals and coda amplitudes.

Sensitivity to the Moho
As the lateral heterogeneity of the Moho is the other relevant factor affecting wave propagation, different crustal models are also implemented. The results are compared by varying Fig. 3 a: Comparison of the three seismogram components (the receiver location is shown in Fig 2) with the synthetic waveforms obtained implementing velocity Model 1.3. Layers velocity is shown in Table 2. The crust thickness is set according to Manu-Marfo et al. (2019), the sediments are 2-km thick below the continental regions, and there are no sediments below the sea (as explained in section 2.1). b: Comparison between synthetic seismograms obtained for 2 km of sediments (red) and thicker sedimentary cover (4-5 km-blue) set according to Molinari and Morelli (2011) (Model 2), assuming velocity Model 1.3. c: Comparison between synthetic seismograms obtained for sediment thickness as in Model 2, where the sediment layer is characterized by low velocities, V P = 2.9 km/s and V S = 1.5 km/s. The black waveforms are the receiver recordings the thickness of the crust. We implement different Moho models, including those estimated by Manu-Marfo et al. (2019) and Molinari and Morelli (2011) (EPcrust) (Fig. 4). Across the oceanic basin, the models revealed a thinner crust across the Tyrrhenian sea, giving values of about 10 km of Moho depth. However, a region separating the two oceanic subbasins, the Marsili basin in the east and the Vavilov basin in the west, has been geologically identified as a portion of continental crust (Issel bridge, Sartori 2003). Assuming that the Vavilov and Marsili Plains are separated by a thicker continental crustal, we implement two other crustal models by deepening the depth of the Moho of the EPcrust model across the southern Tyrrhenian Sea. Figure 4 shows the differences between the four Moho models, pointing out the sharp crustal thickness variation across the basin. Figures S7 and 5 compare the recordings at the station in Sicily with the synthetic seismograms obtained for the thicker and thinner crustal models and the model including an 18 km deep Moho in the centre of the basin (MohoEP18). The comparison is performed in the P-and S-coda time windows, between the origin time and a lapse time of up to 200 s. The corresponding correlation coefficients (Fig. 5) discriminate the effects of the crustal An increase in the oceanic crust thickness leads to delayed arrivals on the three components (Fig. 5). As the propagation is mainly north-south (Fig. 1), the east-west component isolates the SH waves and is reproduced effectively (Fig. 5a).
Wave polarization helps discriminate the effects of a horizontal interface such as the Moho for the peculiar geometry of the chosen north-south propagation. To perform a similar correlation analysis at the INGV network, we rotate all WE and SN seismic components, obtaining radial and transverse components for each source-station path (Fig. 5). For each station in Fig. 1, the synthetic seismogram is then compared to data (Figs. 5 and S7). We estimate the correlation coefficients in the time window 0-200 s after filtering the seismograms in the frequency band 0.05-0.33 Hz and compare the results with the other two Moho models for the basin (Molinari and Morelli 2011;Manu-Marfo et al. 2019). The correlation is mapped in space with a regionalization approach, averaging the coefficients corresponding to different ray paths and crossing each cell (Fig. 6).

Fig. 5 a:
Comparison between the recordings and the best Moho model, MohoEP18 (Fig. 4), for each seismogram component recorded at the station in Sicily. b and c: Comparison between the synthetic seismograms obtained for a crustal model with/without the pinch and the recordings in two time windows corresponding to compressional arrivals and crustal phases, respectively. The data are filtered in the frequency band 0.05-0.33 Hz. The numbers indicate the correlation coefficients between the data and the three synthetic seismograms corresponding to the pinch model (blu waveform), the flat 10 km and 30 km Moho (violet and red waveforms) 1 3 Fig. 6 Regionalization of the positive correlation coefficients ( > 0 ) between synthetic waveforms and recordings at the stations located across Italy (Fig 1)

Comparison Between Flat and Realistic Interfaces: Topography and Correlation Analysis
Because the free surface plays a relevant role in the propagation of regional phases as Lg and Rg, we created three two-layered crustal models by varying the flat-Moho depth and modelled the wavefield in the case of a flat free-surface and the topography (ETOPO1-Amante and Eakins 2009). Figure S8 shows the results of the wavefield simulations for the radial, transverse and vertical components recorded at a seismic receiver located in Sicily ( 37.89 • , 13.30 • ) across the NS profile (Fig. 1). The topography has marginal effects on the transverse component and all direct phase arrivals and amplitudes; however, it substantially changes the early coda of radial and vertical components, damping and dispersing Moho reverberations and increasing late coda amplitudes. Topography is thus included in the final models. The qualitative and quantitative comparison between different realistic Moho models for all the available seismic receivers has been performed to identify the quantities that can constrain the depth of the discontinuity through the inversion. Then, we performed additional simulations using flat Moho models and estimated the correlation coefficients between the synthetic seismograms (those obtained for the flat interfaces and the MohoEP18 cases) and the recordings. The three components are divided into two time windows (Fig. 5). The first window comprises the compressional arrival (i.e. the P head wave), while the second includes the S-wave and the following crustal, reverberating and converted waves. Figure 5 also shows the correlation coefficients obtained by comparing the data with synthetic models obtained embedding a crustal pinch between the basin and Sicily (Fig. 5a) and two flat Moho depths at 10 km and 30 km (Fig. 5b and c). These tests are initially computed only for the station in Sicily, directly south of the earthquake source. The correlation coefficients computed in the two time windows clarify how a pinch is necessary to fit low-frequency data for this source-station path (Fig. 5a).
We performed the final correlation analysis at all the array stations using a single time window of 200 s that comprises P-, S-, and early coda windows. We estimate the correlation coefficient between the synthetic and the recorded waveforms for each of the crustal models available (Fig. 6). By regionalizing the correlation coefficients in space, we constrain the validity of the crustal models for the basin. The MohoEP18 model, presenting a crustal pinch, is the one that best explains data when considering all three components in the Issel Bridge, the region between Vavilov and Marsili.

Discussion
The mixed continental-oceanic crust in basins produces wave reverberations due to sediments, crustal extension and topography in the early coda wavefield (Fig. 7). Source parameters estimations are thus generally based on near-source data (Tinti et al. 2016) as modelling these phases is challenging. Our study demonstrates that once full-wavefield simulations are available, regional-scale recordings can discriminate between different models of source characteristics (Fig.2 and Supplementary) after including crustal velocities and thicknesses derived from previous amplitude-dependent modelling (Nardoni et al. 2021). The full wavefield simulations and use of later coda waves enable the assessment that the optimal source function for the Accumoli earthquake is a boxcar with a rise time of 6 s, a depth of 7.3 km, and the radiation pattern previously defined by INGV (Table 2).  (Table 1). The insights in each panel show the details of the propagation effects due to the Moho interface and the topography. a: The P head wave is produced, and the Moho slope causes energy leakage into the mantle. b: After 66 s, the P-wavefront travels unperturbed through the flat Moho. Coda waves diffuse under central Italy while reverberating P-waves (see insight) hinder the detection of S-wave phases on vertical and radial components. c: Reverberating P-and S-waves propagate through the flat Moho, while the pinch north of Sicily focuses shear waves in a waveguide. The lack of strong heterogeneity in the mantle produces recognizable P-and S-waves and a relatively small coda All regional seismic phases show high sensitivity to crustal structures (Furumura and Kennett 1997). The Moho depth changes the arrivals and amplitudes of complex crustal phases, especially Lg and P/SV conversions and surface waves (Sens-Schönfelder et al. 2009). Figures 5 and S7 show data fits across the full seismogram recorded at a single receiver located in Sicily: this source-station path crosses the region characterized by the largest difference between the three available Moho models. The first Moho (Moho10 km-Manu-Marfo et al. 2019) is very shallow (approximately 10-12 km thick) below the southern Tyrrhenian basin and offshore the Calabria coast, reflecting the effect of the eastward retreat of the Ionian slab, which resulted in the opening of the Southern Tyrrhenian basin (Malinverno and Ryan 1986;Faccenna et al. 2001;Doglioni et al. 1999;Lo Bue et al. 2021). The EPcrust model, similarly to Di Stefano et al. (2011), has Moho depths varying from an average of 12-20 km in Sicily and the volcanic area of the Eolian Islands. By investigating a single phase primarily dependent on the vertical component of motion, Moho10 only offers excellent vertical correlation. EPcrust is derived from the integration of larger-scale models (Molinari and Morelli 2011) and produces average correlations across the entire basin, marginally better than Moho10 when using the transverse component. Only MohoEP18 (Fig. 5) offers satisfactory correlations with data across the three components and especially between Vavilov and Marsili (Fig. 5).
MohoEP18 has been established using the full wavefield: its best fit to data is a consequence of considering a wider time window and information discarded by analytical techniques. A thin crustal pinch is necessary to fit data across the Marsili plain. The results imply the need for a deep oceanic Moho across the basin, specifically between Vavilov and Marsili, to produce the highest and statistically valuable correlations for radial, transverse, and vertical components. The best fit across the Issel Bridge (Sartori 2003;Cocchi et al. 2009) is evidence of the need for thinned continental crust between the volcanic centres. Detecting the Bridge confirms that thinned continental crust is still present between the two volcanic centres associated with two different extension episodes: Vavilov, during Tortonian to Pliocene WE migration and Marsili, formed by Pleistocene NW-SE migration (Sartori 2003;Cocchi et al. 2009). We can still not recover the eastern structure of the basin satisfactorily due to more complex continental structures (i.e. faults, thicker sediments- Owens et al. 1984) that are not implemented in the modelling. Figure 3a shows that reverberations in thicker sediments dominate the transverse component in the early coda time window, producing the primary difference in the very late coda for the vertical and radial components. Transverse components appear most sensitive to the low-velocity sedimentary layers, making them ideal markers of mid-crustal boundaries induced by volcanism across the Aeolian Islands, which is not implemented in the model. Schneider et al. (2013) demonstrate that the transverse component is dominated by dipping structures, such as slabs, or anisotropy, and not by horizontal discontinuities. Instead, horizontal interfaces significantly influence the north-south and vertical plane, i.e. the P and SV polarizations having a component in the vertical direction. Figure S8 demonstrates the influence of water embedded within realistic bathymetry on the radial and vertical components. The comparison in Fig. S8 suggests that the sea column is primarily responsible for the attenuation of these components. Low P-wave velocity in the sea relative to the surrounding bedrock causes wave energy trapped within the basin and decreased amplitude at the receiver, as stressed by Maeda et al. (2013). Water acts primarily on the Rayleigh modes (i.e. radial and vertical components- Petukhin et al. 2010), while the effects on Love waves are weaker and due to conversions of Rayleigh into Love waves (e.g. due to transversely inclined shallow interface of oceanic sediments). Instead, the water does not affect the propagation of direct P/S or converted waves. Implementing topography and water-filled bathymetry is thus a necessary step to obtain feasible Moho depths when using full-wavefield information (Zhang et al. 2022).
The results demonstrate that the thinning of the crust and the pinch slope are responsible for reverberations, energy leaked into the mantle and converted waves (Furumura and Kennett 1997). The evolution of the wavefield in Fig. 7 clarifies the strong interaction of P and S waves with the topography and Moho. As shown in Fig. 7a, the interfaces generate P head waves, P-SV conversion, and energy loss into the mantle. Once the P waves propagate and reverberate across the crustal waveguide, there are still several conversions between P and S waves due to the topography and presence of a solid-liquid interface (Fig. 7b). After 130 s (Fig. 7c), the coherent S-wave packet still propagates under Sicily, while the wavefield has reached diffusion across the continental Moho of Central Italy. Figure 7c confirms that the time window 130-200 s can be used to estimate the location and shape of the crustal pinch, where shear waves are guided, producing waveforms with unique characteristics. We can thus best map the Moho depth across mixed continental-oceanic settings from regional observations using two attributes: the first-wave arrival, sensitive to the continental Moho depth (Fig. 7a), and the crustal coda waveguide associated with an oceanic Moho (Fig. 7b).
Our work definitively clarifies how the properties of a basin affect seismic attributes that are suitable to invert for the Moho depth, using a single distant earthquake. Across mixed continental-oceanic settings, reverberating waves inhibit standard full-waveform inversion that uses misfits based on coherent phases and amplitude. The forward modelling and inversion by parameters trial we propose in this paper are necessary to set the correct attributes to use in a full-waveform inversion for the basin structure and earthquake source simultaneously. The correlation analysis, a standard approach in seismology (Alvarado et al. 2007;Lee et al. 2014), is key for understanding how to couple the coherent and the stochastic information in full-waveform inversions. Simulating full seismic wavefields up to 0.3 Hz across realistic basin structures is a computationally heavy task, to our knowledge, never attempted before. Previous full-waveform modelling across basins was limited to lower frequencies (Fichtner and Villaseñor 2015;Simutė et al. 2016). Even by disregarding variations in source characteristics, which are efficiently constrained using near-source data (Tinti et al. 2016), and concentrating only on Moho depth, the minimization of appropriate cost functions requires a large number of simulations and, at present, is well-beyond our computational capacity.

Conclusions
We perform high-frequency (up to 0.3 Hz) wave-equation modelling to explore the effects of the 3D structural variations on the propagation of the full wavefield across the Southern Tyrrhenian Sea. Forward modelling of P-, S-and early coda waves, produced by the Accumoli earthquake in Central Italy, depends on topography and sediment, crustal thicknesses, velocities, and velocity fluctuations. The models are compared qualitatively and inverted by parameter trial with phases and amplitudes observed at the regional INGV network. The sensitivity of the seismic wavefield to source properties and crustal discontinuities at depth (i.e. crustal and sediment thicknesses and velocities) in the frequency band 0.05 − 0.33 Hz provides refined propagation corrections to estimate earthquake source properties. The optimal earthquake source model is a boxcar source time function with a rise time of 6 seconds, depth of 7.3 km, and normal faulting with a strike, dip and rake of 155 • , 49 • and − 87 • , respectively.
A thin sedimentary layer is necessary to reproduce the data avoiding reverberations from trapped waves. The simulations highlight the existence of a crustal waveguide in the early S-coda due to a laterally varying Moho interface. Its detection trades off with topography, which can disperse and scatter the waves later in the seismogram. The inversion by trial uses a correlation analysis. The best Moho model results in a thinned continental crust between the Vavilov and Marsili volcanic centres, the Issel Bridge, a product of the changes in the back-arc extension caused by the Ionian slab. We still cannot recover the full wavefield propagating through the Italian continental crust as we neglect broader structures such as extended faults or sediments. However, the simulations show that the central portion of regional seismograms is sensitive to crustal thinning. This sensitivity proves the need for crustal pinching between the basin and the Sicilian crust. The modelled regional phases and amplitudes are an ideal forward model for future inversions of Moho depths, especially across mixed continental-oceanic regions, using recordings from a single high-magnitude earthquake with dense arrays. An advantage of this technique is that it can be applied when high-magnitude earthquakes recorded at wider arrays on the continent nucleate across oceanic arcs or in areas with no near-source seismic deployments. The reproduced features are essential for assessing past back-arc extension following subduction in the Southern Tyrrhenian sea: similar modelling can be implemented through any interpreted geodynamic model after transformation to seismic velocities, allowing the assessment of present-day data.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.